Spatio-Temporal Modeling of SST for the Assessment of Climate Risk Over Aquaculture in the Coast of the Valencian Region

Received: Revised: Accepted: Published: Copyright: © 2025 by the authors. Submitted to Journal Not Specified for possible open access publication under the terms and conditions of the Creative Commons Attribution (CC BY) license. Article Spatio-temporal modeling of SST for the assessment of climate risk over aquaculture in the coast of the Valencian Region Laura Aixalà1, Irene Lopez-Mengual1Javier Atalah2, Juan Aparicio1, David Ballester3, David Conesa4, Aitor Forcada3, Jonatan Gonzalez-Monsalvo1, Antonio López-Quílez4, Pablo Sanchez-Jerez4and Xavier Barber1* 1Centre of Operations Research, Miguel Hernández University of Elche, Spain 2Environment Aquaculture Interactions, Cawthron Institute, New Zeland 2Department of Marine Sciences and Applied Biology, University of Alicante, Spain 2Statistics and Operations Research Department, University of Valencia, Spain *Correspondence: [email protected] Abstract 1 Climate change poses significant risks to Mediterranean aquaculture, with sea surface tem2 perature (SST) identified as a critical stressor affecting cultivated species. This study aims 3 to assess climate-related risks for coastal aquaculture in the Valencian Community (Spain) 4 by analyzing SST spatiotemporal variability and predicting future trends. A multi-method 5 approach was employed, combining ARIMA models for 10-year predictions at eight coastal 6 locations, Bayesian hierarchical models (BHM) fitted via INLA for spatiotemporal analysis 7 of maximum SST and temperature range (2000–2024), and Generalized Additive Models 8 (GAM) to evaluate relationships with climate indices (NAO, AMO, ENSO). Results revealed 9 a consistent warming trend since the 1990s, with ARIMA predictions indicating maximum 10 SST values of 27.2 ± 0.1 °C in September over the next decade. The spatiotemporal model 11 showed effective spatial correlation ranges of 246 km for maximum SST and 207 km for 12 SST range. Anomalous warming years (2003, 2006, 2018, 2023–2024) coincided with doc13 umented marine heatwave events. The GAM explained 98.2% of deviance, with AMO 14 showing significant influence ( p< 0.001) while ENSO was not statistically significant. 15 Notably, the area north of San Antonio Cape exhibited lower warming trends, suggesting 16 potential climate refuge characteristics. Southern locations (Altea, Campello) currently 17 experience the highest temperatures, but projections indicate Valencia and Sagunto will 18 become the warmest areas. These findings provide essential information for marine spa19 tial planning and recommend a precautionary approach when considering aquaculture 20 relocation towards northern coastal areas. 21 Keywords: Time series; ecosystems; nonlinear methods; short-term and long-term dynam22 ics, SST; marine spatial planning; climate change. 23 1. Introduction 24 Achieving a smarter and sustainable aquaculture and fishing industry, equipped 25 with technological advancements tailored to climate change and its monitoring, stands 26 as a predominant objective in marine science research and development. The inadequate 27 planning could become economically unsustainable due to climate stressors and the actual 28 climate variances have impacted aquaculture production negatively around the world ([ 1 ]). 29 There has been a notable increase in the amplitude of the daily sea surface temperature 30 Version December 16, 2025 submitted to Journal Not Specified https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 2 of 24 (SST) seasonal cycle over many ocean basins since 1950 ([ 2 ]). Also, many studies have tried 31 to analyze these trends, highlighting that some areas are more vulnerable than others. Our 32 area of interest, located in the Mediterranean Sea, is one of these sensitive regions ([ 3 ]). This 33 region is recognized as the largest climate change hotspot, experiencing warming rates up 34 to 20% higher than the global average ([4]). 35 In order to establish a better spatial planning for the sustaintable acuaculture, the 36 Available Zones for Aquaculture (AZAs) defined by the FAO ([ 5 ]) summarize the main 37 oceanographic features necessary for developing sustainable aquaculture along the Mediter38 ranean coast. However, under a climate change scenario, these conditions will evolve in 39 the coming decades. The magnitude and frequency of these changes will expose cultivated 40 species to varying levels of stress, which may directly affect their health and survival 41 ([ 6 ],[ 7 ]). Relocating fish farms to areas less sensitive to climate change is one option for 42 mitigating the effects of temperature increase ([ 8 ],[ 9 ]), and efficient Marine Spatial Planning 43 must be performed. It is crucial to consider the factors that have the greatest impact on 44 the survival of cultured species. Studies about growing projections in a climate change 45 scenario in the Mediterranean Sea ([ 10 ]) alert about future areas that will be more suitable 46 for aquaculture, and it is expected that the western Mediterranean becomes the suitable 47 one. In order to archive that, we need to pay attention to studies like [?], which summarize 48 about all the stressors in the cultivated species in a warming climate change scenario. 49 Some studies confirm a persistent warming trend in daily sea surface temperature 50 (SST) data derived from satellites throughout the Mediterranean region. This increase 51 in temperature is observed on various temporal scales, ranging from daily to monthly, 52 seasonal, and decadal assessments, affecting not only surface waters but also extending 53 to deeper and intermediate areas ([ 3 ], [ 11 ]). While there is some controversy surrounding 54 this issue, [ 8 ] emphasized the potential effects of increased surface temperature on coastal 55 aquaculture, acknowledging both negative and positive impacts. Several studies have 56 demonstrated that SST variability is significantly correlated with fish growth,(Islam et al., 57 2020; Haberle et al., 2024). Islam et al. (2021) conducted a comprehensive study on the 58 effects of climate change on aquaculture. Additionally, high to moderate temperatures may 59 have beneficial effects on growth models, but they can also potentially cause muscular 60 dystrophies and neurological impairments, among others (Madeira et al., 2020). For this 61 reason is difficult to set up a braking point for the SST in the aquaculture context. The season 62 of growth and climate risk may differ from one year to another, sometimes overlapping 63 among them. This is supported by the study of [ 11 ] which shows how the range of seasonal 64 patterns is changing due to the effects of climate change in the Mediterranean Sea. 65 Within this context, the present study aims to focus on the risks associated with coastal 66 warming for aquaculture. Therefore, it is essential to understand how temperatures have 67 evolved during key periods in the past to better predict their future trends. To achieve this, 68 the focus must be on potentially cultivable areas, specifically in seashore areas within 50 69 km from the coast to facilitate the establishment of logistical bases. The primary objective 70 is to generate 10-year future predictions of SST and conduct a temporal analysis from 1990 71 to the present, including interannual variability. The locations studied for this objective are 72 several points where aquaculture activity is currently taking place. Secondly, we aim to 73 delve into the spatio-temporal variability of the maximum SST recorded in the area, as well 74 as the temperature range (difference between maximum and minimum values). In order to 75 improve and expand the avaliable information from the data, a complex prediction model 76 will increase the resolution of the area and the coverage of the data near to the weak points 77 as the coastal margins.This spatio temporal model will cover the ummarizing, this study 78 aims to provide a localized and detailed assessment of climate risks associated with high 79 temperatures in a region where aquaculture plays a significant role. By addressing both 80 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 3 of 24 historical variability and future projections, the research seeks to equip stakeholders with 81 the insights needed to mitigate risks and adapt to changing environmental conditions. 82 2. Material and Methods 83 2.1. Study area 84 This study focuses on the East Mediterranean coasts of Spain, a Region named Va85 lencian Community with located between the Catalonia and Murcia region (Figure 1). In 86 terms of aquaculture production in Spain, this region was the second largest producer 87 of cultivated fish in 2022, with a production of 13,975 t. In 2023, it became the leading 88 producer, achieving a 52% increase, which corresponds to a total production of 21,227 t.([? 89 ]). Initially, the investigation examined eight individual time series from different locations 90 within the Valencian Community, near the municipalities of Altea, Burriana, Campello, 91 Denia, Guardamar, Sagunto, Valencia, and Vinaroz (Figure 1). The selected location, except 92 Denia, were chosen because they are key areas for current fish aquaculture development. 93 Denia was included to understand the behaviour of the region north of San Antonio Cape. 94 This area’s orography and the current system’s characteristics set up a barrier, which can 95 change the environmental factors from north to south ([ 12 ]). Conducting separate analyses 96 for these locations allowed us to obtain comparable results regarding the serial behaviour 97 inherent in each place. 98 Figure 1. Area of study located in the western Mediterranean (left) and locations details (right). 2.2. Data sources 99 Ecosystem phenomena databases are commonly evaluated through the system100 atic collection of repeated measurements of environmental variables [ 13 ]. ). This 101 study employs the Reprocessed Mediterranean Sea Surface Temperature (SST) dataset 102 SST_MED_SST_L4_REP_OBSERVATIONS_010_021 , which provides a stable, long-term SST 103 time series developed specifically for climate studies. The product consists of daily, night104 time, optimally interpolated (L4) satellite-based SST estimates that are largely free from 105 diurnal variability. It spans from January 1, 1981, to approximately one month before 106 present and is provided at a 0.05° resolution gridcovering the entire Mediterranean Sea 107 and a portion of the adjacent North Atlantic ([ 14 ];[ 15 ]). The dataset is constructed from a 108 consistent reprocessing of the ESA Climate Change Initiative (CCI) Climate Data Record 109 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 4 of 24 (CDR) v3.0, which covers up to 2021, and is extended with an Interim Climate Data Record 110 (ICDR) for data beyond 2022. 111 Climate indices were obtained from various sources. The North Atlantic Oscillation 112 (NAO) was retrieved from the National Oceanic and Atmospheric Administration (NOAA) 113 Physical Sciences Laboratory (PSL) as a daily index, which was subsequently averaged 114 to generate a monthly index. These indices are derived from 500 mb geopotential height 115 patterns, calculated using the NCEP/NCAR Reanalysis 1 dataset. The index is obtained 116 by subtracting the area-averaged values over 55–70°N and 70–10°W from those over 117 35–45°N within the same longitudinal range. A spectral truncation at total wavenumber 118 10 is applied to the height fields to isolate large-scale teleconnection patterns, using the 119 1981–2010 period as the reference climatology ([ 16 ]). The monthly Atlantic Multidecadal 120 Oscillation (AMO) index was also obtained from NOAA PSL. As the version based on 121 the Kaplan SST dataset is currently unavailable, the index based on NOAA ERSSTv5 was 122 used instead. Combined SST data (HadISST1 and NOAA OI SST) were used to compute 123 anomalies in the North Atlantic (0°–60°N, 80°W–0°) relative to the 1901–1970 climatology. 124 The global signal (60°S–60°N) was removed by subtracting annual means, and a 13-point 125 running mean filter was applied to isolate decadal-scale variability. The resulting series is a 126 corrected AMO index, representative of internal multidecadal variability in the Atlantic 127 ([ 17 ]). El Niño–Southern Oscillation (ENSO) has several indicators. For this case, the 128 Niño Index (ONI) calculated from the equatorial Pacific difference of surface temperature 129 has been chosen . Warm and cold phases are defined as a minimum of five consecutive 130 3-month running averages of SST anomalies. This dataset, obtained from NOAA’s Climate 131 Modeling branch, was originally provided as quarterly moving averages ([ 18 ]). To integrate 132 it into the analysis, each value was assigned to the central month of the corresponding 133 trimester, creating a monthly time series. The dates of all datasets were harmonized using 134 the lubridate library ([ 19 ]) in r , and the time series were merged into a unified database 135 using date-based joins, ensuring temporal consistency. 136 2.3. Spaitotemporal statistical methodology 137 2.3.1. ARIMA 138 Sea surface temperature prediction methods can be broadly categorized into three 139 groups: numerical methods, data-driven methods, and a combination of both ([ 20 ]). Data140 driven methods center their focus on thorough analysis of observed patterns and rela141 tionships, enabling them to efficiently estimate and predict accurate SST values. These 142 methods learn from the available data to infer the response variable and can operate with143 out extensive domain knowledge of ocean and atmospheric processes, making them less 144 sophisticated than numerical methods. However, this simplicity allows them to excel 145 in high-resolution SST predictions at smaller scales ([ 20 ]. Techniques such as traditional 146 statistical methods and machine learning and artificial learning approaches are commonly 147 employed to estimate and predict SST ([20],[21], [22],[23]) 148 The temporal series of each location (Figure 1) has been extracted and analyzed to 149 obtain the monthly maximum SST values from January 1990 until March 2022. After that, 150 these values have been modelled by ARIMA models ([ 24 ]) and 10-year predictions have 151 been calculated. 152 In particular, the model implemented is ARIMA(1, 1, 1)(0, 1, 1)12:153 (1−ϕ1L)(1−L)(1−L12)yt=c+ (1+θ1L)(1+θ12L12)εt, (1) 154 where L is the lag operator, ϕi are the autoregressive model parameters, θi are the 155 parameters of the moving average part, and εt are the error terms. εt are generally assumed 156 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 5 of 24 to be independent, identically distributed variables sampled from a normal distribution 157 with zero mean. 158 Leting ∇kyt=yt−yt−k be the k differentiation of the observation k . Then the model 159 ARIMA(1, 1, 1)(0, 1, 1)12 is equivalent to: 160 ∇yt=∇yt−12 +ϕ1(∇yt−1− ∇yt−13) + εt+θ1εt−1+θ12εt−12 +θ1θ12εt−13. (2) 161 The ARIMA model used does not specifically predict the maximum temperature 162 for each year. Instead, is designed to model and predict monthly values of sea surface 163 temperature (SST), particularly the monthly maximum SST values. This model includes 164 components for seasonal differencing as well as non-seasonal differencing. Summarizing, 165 it focuses on capturing recurrent monthly and seasonal patterns rather than providing a 166 single annual prediction. The model for each location has been implemented through the 167 function arima() from the stats R package ([25]). 168 As a path to better analyze the SST gradient over time at a global scale, more sophis169 ticated modelling techniques were employed. Specifically, data from this century was 170 used, as it revealed the highest recorded values in recent years. To capture comprehensive 171 spatio-temporal patterns spanning the entire study area, a cohesive approach was adopted, 172 treating the entire dataset as a single entity. The data analyzed for this approach is located 173 along the coast, as visually depicted in Figure 1(-0.98°to 1.38°E and 37.85°to 40.52°N), 174 with a designated coastal buffer highlighted in blue, extending 50 km offshore. 175 2.3.2. Bayesian hierarchical model 176 For the analysis, a Bayesian hierarchical model (BHM) ([ 26 ]) was employed, well177 suited for accounting for the complex correlations within the geostatistical data and their 178 temporal evolution. This holistic perspective provided valuable insights into areas at risk 179 during different periods. The BHM was used to analyze, between June and October, the 180 maximum SST values and the annual range between maximum and minimum SST values 181 for each year. Those period range has been selected conducting a preliminary study of the 182 extreme temperatures over the late spring, summer and early autumn. This representation 183 of the 95th percentile of SST over the study area provides a preliminary view of how 184 extreme values have shifted over the decades. Although June and October are not among 185 the peak months for these extremes, they have shown a progressive increase over time, 186 reaching extreme temperatures above 25°C. For this reason, these two months have been 187 included in the study period of the BHM. 188 Let Yit be the values of the annual SST maximum or range in summer-autumn season 189 measured in the locations i=1,...,I where t=2000,...,2021. We assume: 190 Yit ∼Gammaµit,ϕ2 e, log(µit)=β0+w(si,t)+α(t), (3) 191 where βo is the intercept and ϕ2 e is the precision parameter of the measurement error. 192 ϕ2 e is assumed to follow a PC prior. w(si , t) is a spatially correlated at location si random 193 effect that changes in time t with first order random walk. Specifically, each w(si , t) follows 194 a zero-mean Gaussian distribution temporally independent but spatially dependent at each 195 time with Matérn covariance function. α(t) is the temporal random effect that assumes 196 that the data is independent and identically distributed. To avoid parameter identification 197 problems the maximum and range models have been fitted by using different percentages 198 of the whole dataset (from 10% of the data to 80%). 199 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 6 of 24 The models have been fitted in a Bayesian framework through the Integrated Nested 200 Laplace Approximation implemented in the r package r-INLA ([?]). As INLA can not fit 201 continuous Gaussian Fields (GFs) it approximates them using a Stochastic Partial Differen202 tial Equation approach (SPDE). Hence, a separable space-time model is defined as a SPDE 203 model for the spatial domain and a random effect for time dimension. 204 To assess temporal trends in sea surface temperature (SST) over the study period, the 205 Sen’s slope estimator ([ 27 ]) was used, a non-parametric measure of trend that is robust 206 to outliers. Monthly SST data were reshaped into a spatial-temporal matrix, with spatial 207 coordinates (x, y) and columns corresponding to each month/year. Sen’s slope (°C°C 208 per year) was calculated for each spatial cell using the sens.slope() function from the r209 package trend ([ 28 ]). The slope estimate considers all observations in the time series of each 210 cell without assuming a normal distribution. Additionally, the Mann-Kendall test([ 29 ];[ 30 ]) 211 was applied to assess the statistical significance of the trends, as previous works as [ 31 ]. 212 This non-parametric test detects the presence of monotonic trends in time series without 213 requiring normality. Cells with significant positive or negative trends at =0.05=0.05 are 214 reported. 215 2.3.3. Generalized Additive Models 216 A Generalized Additive Model (GAM) was used to capture potentially nonlinear 217 relationships between sea surface temperature (SST) and various climate forcings (AMO, 218 ENSO, NAO). This approach allows for the modeling of smooth and flexible effects in con219 tinuous variables, as well as the inclusion of seasonal and temporal components ([ 32 ];[ 33 ]). 220 The implementation was carried out using thethe mgcv package ([ 34 ]) in r , with REML 221 estimation employed to avoid overfitting ([ 35 ]). Temporal lags were also considered for 222 certain indices (e.g., a onemonth lag for NAO) to account for potential delayed effects on 223 SST ([ 36 ]). Observations with missing values (NA) were removed to ensure model integrity. 224 The model follows the equation: 225 SSTt=β0+factor(Yeart) + f1(Month) + f2(AMOt) + f3(ENSOt) + f4(NAOt−1) + εt(4) 226 where fi represent smooth functions (splines or tensor product interactions) estimated 227 with penalization. To capture annual seasonality, a cyclic smooth term over the month of 228 the year was included ( f1(Month) ). Alternative models were also tested, incorporating 229 year either as a smooth effect ( f2(Year) ), as a categorical factor ( factor(Year) ), or as part of 230 interaction terms with the climate variables (f3(ClimateVar ×Year)). 231 3. Results 232 3.1. ARIMA: a comparison between locations 233 The data shows a consistent temperature increase since the 1990s, with a noticeable 234 change in the rate of increase after the year 2000 (see Figure 9). This observation highlights 235 the impact of climate change on the world’s oceans. Notably, after the year 2000, while 236 temperatures continued to rise, the rate of increase appeared to slow down compared to the 237 preceding decade. Regarding the inter-annual variation, the summer months are expected 238 to have abnormally high temperatures in specific years (2002, 2017) with a larger range 239 of temperatures meanwhile, the winter seasons have less variation in temperatures over 240 the years. Another relevant observation from the seasonal analysis is that the temperature 241 during the spring and autumn months has increased from 3 to 4 degrees over the years. 242 Not only that, but the changes are stratified, they are not as random as the peaks recorded 243 in summer and winter, but they are more systematic and progressive. 244 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 7 of 24 Figure 2. Seasonal time series of SST at each location, showing observed data and ARIMA-based 10-year prediction band (pink). The shaded area indicates the forecast uncertainty. The plot highlights seasonal variation, with summer months showing higher inter-annual variability, while winter months remain relatively stable. The predictions’ red band shows how the upcoming years using the ARIMA model 245 reveal a clear upward trend in temperatures across all locations for the next 10 years (Figure 246 9). This trend shows how the autumn season is going to change over the years, and in the 247 next 10 years, the SSTs will reach a maximum of 27,2± 0,1 ºC in September, 24,7± 0,08 ºC in 248 October and 21,3± 0,08 ºC in November. 249 To gain insights from the data, we performed a decomposition of the series and 250 predictions, showcasing the trends for each location in Figure 10. The results reveal a 251 consistent upward trend in SST over time. Additionally, an interesting correlation emerges 252 when examining temperature variations with latitude. The points observed in the southern 253 part of San Antonio Cape (Campello, Altea, Guardamar) will reach higher temperatures 254 than those in the northern part (Valencia, Sagunto and Burriana). Another exception can 255 be noted in the northern locality, Vinaroz, which trend is not as high as the others. Also, 256 the warmest locations initially observed are the municipalities of Altea and Denia, but the 257 data suggests that Sagunto and Valencia will surpass current predictions and become the 258 subsequent warmest regions. 259 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 8 of 24 Figure 3. Trend decomposition of SST time series and 10-year ARIMA predictions for each location. Points represent monthly SST observations, lines the trend component, and the shaded area the prediction interval. Seasonal and latitudinal variations are captured in the trend component The models parameters are slightly different for each location (Table 6). Hence, follow260 ing the equation 2, the temperature differentiations ∇yt−1 (temperature difference between 261 month t and its previous month t− 1) and ∇yt−13 (temperature difference between month 262 t and its previous month the year before t− 13) have less relevance in the northen locations 263 estimations as ϕ1 is minor. Moreover, the parameter θ1 which represents the magnitude 264 of the effect of the error the month before has practically the same effect in all locations 265 (with a slightly minor effect in Campello and Guardamar. On the other hand, the error 266 terms εt−12 (twelve months before) and εt−13 (thirteen months before) have pretty much 267 the same effect in all locations with the minimum in Altea and the maximum in Valencia as 268 the parameter θ12 has the smaller and maximum values respectively. Both parameters θ1269 and θ12 have negative effect at temperature. 270 ϕ1θ1θ12 Vinaroz 0.300 -1.000 -0.922 Burriana 0.326 -1.000 -0.916 Sagunto 0.333 -1.000 -0.929 Valencia 0.350 -1.000 -0.930 Denia 0.387 -1.000 -0.920 Altea 0.404 -1.000 -0.885 Campello 0.416 -0.997 -0.925 Guardamar 0.386 -0.994 -0.900 Table 1. ARIMA( 1, 1, 1 )( 0, 1, 1 )12 parameter values. ϕ1 represents the magnitude of the autoregressive effect from the previous month, θ1 the impact of the previous month’s error term, and θ12 the effect of the error 12 months prior. Northern locations show lower ϕ1 values, indicating smaller month-to-month autocorrelation, while θ1 and θ12 are similar across sites, with small differences in magnitude. 3.2. Spatio-temporal modelling 271 In general, the spatio-temporal evolution of SST in the Mediterranean is studied 272 using reanalysed satellite and in situ data. These models offer high reliability, as they 273 incorporate multiple levels of validation. However, they also present certain limitations, 274 such as limited spatial resolution and reduced accuracy near the coastline. The advantages 275 of the model used in this study include its ability to directly model spatial correlation, 276 capture information about temporal evolution, and account for potential interactions with 277 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 9 of 24 spatial variability. In addition to providing a comprehensive estimate of the uncertainty, it 278 also allows an accurate quantification of the reliability of each prediction. 279 Our model 3adapts one part of the original data to build the model and the rest for 280 validation. Based on the summary statistics for each value of the data percentage, it has 281 been decided to proceed with 40% of the data for the model that estimates the monthly 282 maximum SST and 35% of the data for the model that estimates the monthly range SST . 283 A summary of the estimated values for each model parameters: β0 , ϕe , Range for i and 284 sd for i are displayed in Table 7(maximum SST) and Table 8(range SST). The logarithm 285 of the annual maximum SST between June and October is on average 5.706 meaning that 286 the mean yearly maximum SST is 27.03 ºC estimated with an acceptable residual standard 287 deviation ( ϕe ). The effective range of the Gaussian field, the distance in which the spatial 288 effect stabilizes for higher distances, has a mean of 246 km with a small sd with a mean of 289 16.75 Km. 290 On the other hand, the logarithm of the maximum annual difference in SST is on aver291 age 1.293 meaning that the mean yearly difference is 3.64 ºC estimated with an acceptable 292 residual standard deviation ( ϕe ). The effective range of the Gaussian field, which denotes 293 the distance at which the spatial effect stabilizes for greater distances, exhibits an average 294 of 207 Km with a narrow standard deviation of 7.69 Km. 295 mean sd q0.025 q0.5 q0.975 β05.706 0.0024 5.701 5.706 5.710 ϕe0.002 0.000021 0.002 0.002 0.002 Range for i246.365 16.756 212.799 246.700 278.490 sd for i0.0046 0.0003 0.0041 0.0046 0.0051 Table 2. Posterior summaries of the maximum SST spatio-temporal model (40% of the data). β0 is the intercept on the log scale, ϕe is the residual standard deviation, “Range for i” indicates the effective spatial range of the Gaussian field (km), and “sd for i” the standard deviation of the spatial effect. The results show a moderate spatial correlation (246 km) and low residual uncertainty. mean sd q0.025 q0.5 q0.975 β01.293 0.111 1.074 1.293 1.512 ϕe0.005 0.000085 0.005 0.005 0.004 Range for i207.767 7.695 194.379 207.162 224.596 sd for i0.298 0.010 0.280 0.297 0.321 Table 3. Posterior summaries of the SST range spatio-temporal model (35% of the data). The parameters are analogous to Table \ref{tab:st_param_max}. The effective spatial range is slightly smaller (207 km), indicating stronger localized spatial correlation for SST variability. The results of the INLA models leave the general evolution of the maximum trends in 296 the summer-autumn season and the spatial effect of each value in the map. The effective 297 range of maximum SST results (Table 7) indicates an influence area of 246 Km. This 298 influence is a lightly higher value compared with the influence area 207 Km of the SST 299 range (Table 8). The greater spatial correlation in the SST range and maximum trends could 300 be associated with regional phenomena that impact uniformly over large areas of the ocean. 301 Figure 11 show the prediction maps of the maximum annual between June and October 302 for the years since 1995 to 2024. The general trend shows the maximum temperatures 303 throughout the time series ranges between 27 and 28°C. A slight upward trend can be 304 observed, with annual increases between 0.02 and 0.08°C (see Figure 13 part a). Other 305 studies as [ 37 ] have analysed both the seasonal component and the sum of the trend and 306 seasonal components in the Western Mediterranean. Significant trends were found for 307 the study area during summer, while trends were not significant in autumn. Overall, a 308 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 16 of 24 spatial variability. In addition to providing a comprehensive estimate of the uncertainty, it 405 also allows an accurate quantification of the reliability of each prediction. 406 Our model 3adapts one part of the original data to build the model and the rest for 407 validation. Based on the summary statistics for each value of the data percentage, it has 408 been decided to proceed with 40% of the data for the model that estimates the monthly 409 maximum SST and 35% of the data for the model that estimates the monthly range SST . 410 A summary of the estimated values for each model parameters: β0 , ϕe , Range for i and 411 sd for i are displayed in Table 7(maximum SST) and Table 8(range SST). The logarithm 412 of the annual maximum SST between June and October is on average 5.706 meaning that 413 the mean yearly maximum SST is 27.03 ºC estimated with an acceptable residual standard 414 deviation ( ϕe ). The effective range of the Gaussian field, the distance in which the spatial 415 effect stabilizes for higher distances, has a mean of 246 km with a small sd with a mean of 416 16.75 Km. 417 On the other hand, the logarithm of the maximum annual difference in SST is on aver418 age 1.293 meaning that the mean yearly difference is 3.64 ºC estimated with an acceptable 419 residual standard deviation ( ϕe ). The effective range of the Gaussian field, which denotes 420 the distance at which the spatial effect stabilizes for greater distances, exhibits an average 421 of 207 Km with a narrow standard deviation of 7.69 Km. 422 mean sd q0.025 q0.5 q0.975 β05.706 0.0024 5.701 5.706 5.710 ϕe0.002 0.000021 0.002 0.002 0.002 Range for i246.365 16.756 212.799 246.700 278.490 sd for i0.0046 0.0003 0.0041 0.0046 0.0051 Table 7. Posterior summaries of the maximum SST spatio-temporal model (40% of the data). β0 is the intercept on the log scale, ϕe is the residual standard deviation, “Range for i” indicates the effective spatial range of the Gaussian field (km), and “sd for i” the standard deviation of the spatial effect. The results show a moderate spatial correlation (246 km) and low residual uncertainty. mean sd q0.025 q0.5 q0.975 β01.293 0.111 1.074 1.293 1.512 ϕe0.005 0.000085 0.005 0.005 0.004 Range for i207.767 7.695 194.379 207.162 224.596 sd for i0.298 0.010 0.280 0.297 0.321 Table 8. Posterior summaries of the SST range spatio-temporal model (35% of the data). The parameters are analogous to Table \ref{tab:st_param_max}. The effective spatial range is slightly smaller (207 km), indicating stronger localized spatial correlation for SST variability. The results of the INLA models leave the general evolution of the maximum trends in 423 the summer-autumn season and the spatial effect of each value in the map. The effective 424 range of maximum SST results (Table 7) indicates an influence area of 246 Km. This 425 influence is a lightly higher value compared with the influence area 207 Km of the SST 426 range (Table 8). The greater spatial correlation in the SST range and maximum trends could 427 be associated with regional phenomena that impact uniformly over large areas of the ocean. 428 Figure 11 show the prediction maps of the maximum annual between June and October 429 for the years since 1995 to 2024. The general trend shows the maximum temperatures 430 throughout the time series ranges between 27 and 28°C. A slight upward trend can be 431 observed, with annual increases between 0.02 and 0.08°C (see Figure 13 part a). Other 432 studies as [ 37 ] have analysed both the seasonal component and the sum of the trend and 433 seasonal components in the Western Mediterranean. Significant trends were found for 434 the study area during summer, while trends were not significant in autumn. Overall, a 435 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 17 of 24 Figure 11. Predicted maximum monthly SST (June–October) across all locations from 1995 to 2024, obtained from the spatio-temporal INLA model. Colors represent SST in °C. The model captures both spatial and temporal variation, with the shaded areas (or color gradient) reflecting the uncertainty of the predictions. A slight upward trend is observed in maximum SST, particularly during 2015–2024. The effective range of the spatial Gaussian field is 246 km, indicating the distance over which spatial correlation stabilizes. general summer warming trend of 0.016±0.002°C per year was reported. However, the 436 increase in maximum temperatures became more evident during the most recent decade 437 (2015–2024), rising from around 28°C to nearly 30°C. It’s quite evident that the most 438 transcendental increase has been observed between the years 2023 and 2024 Those two 439 years have registered the records of temperature anomalies in the past decades, directly 440 related with one El Niño event ([ 38 ]). It can also be observed that in certain anomalous 441 years, such as 2003, 2006, and 2018, the mean maximum temperature increases by 1 to 2°C 442 relative to adjacent years. The maximum SST values predicted by this model are higher 443 than those reported in other studies, such as[ 11 ] , where, for example, in 2003 and 2018 444 the maximum temperature reached in the western Mediterranean was 27.5ºC. In contrast, 445 our model records maximum values of up to 29ºC for the same years. Both the anomalous 446 years and the period of intensified warming have also registered a higher frequency of 447 marine heatwaves (MHWs) and extreme marine summers (EMS) (See Table9)448 The predicted maximum annual difference (Figure 12) depends directly from the 449 minimum and maximum predicted SST for the hole period. A complementary analysis 450 with the previous SST maximum predictions (Figure 11) help us to understand how the 451 range has evolved over the years. Along the entire data frame, the mean temperature 452 ranges differed from the 5 to 7 degrees. Exceptionally, the range increases up to the 8 453 degrees of difference in the previously mentioned anomalous years (2003, 2006, 2018). 454 These increases are directly related to the increase in the predicted temperature and linked 455 with different heating events (see Table 9) . However, in some other years, the widening of 456 the range is not due to an increase in the maximum, but rather to a rise in the minimum 457 temperatures. For example, the year 2013 shows a predicted temperature (see Figure 11) 458 lower than that of adjacent years, while the temperature range spans from 8 to 1ºC. Similar 459 patterns were observed in 1998 and 2008. The Sen’s slope analysis in Figure 13 (b) gives 460 us a general evolution of the ranges changes along the area. The northern part of the 461 buffer appears to be more affected by the increase in range, while the southern part has 462 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 18 of 24 Figure 12. Predicted range of monthly SST (difference between maximum and minimum SST) for June–October, 1995–2024. Colors indicate the magnitude of SST variability (°C), with higher ranges corresponding to greater intra-annual variation. The spatial model accounts for temporal and spatial correlation, with an effective range of 207 km. Years of interest Impact type Reference 2018, 2003 MHWs with medicanes [39] 2022, 2016 MHWs in the West Med [40] 2003, 2020 EMS in the West Med [41] 2008-2017 Higher SST trend in the West Med [42] 2003, 2018, 2015 MHWs in the West Med [43] 1998, 2003, 2012, 2015, 2022 Sumer and Autumn MHWs in the Med [44] Table 9. Bibliographic recompilation of different studies about extreme sea surface temperatures and their effects in the Mediterranean Sea not shown any increase—in fact, the range has even decreased in some areas. However, no 463 statistically significant trends were detected in any of the examined time series according 464 to the Mann-Kendall tes 465 4.3. Generalized additive model 466 The GAM, assuming a Gaussian distribution with identity link, explained 98.2% of 467 the total deviance, with an adjusted R² of 0.979, indicating an excellent fit to the monthly 468 sea surface temperature (SST) data. 469 Among the smooth terms disposed on the 10, month showed a highly significant effect 470 (p< 0.001), effectively capturing intra-annual seasonality. The AMO also had a strong and 471 significant influence on SST (p< 0.001). In contrast, ENSO was not statistically significant 472 (p= 0.44), and NAO with a one-month lag showed a marginal effect (p= 0.093) (14). 473 Compared to the baseline year (1995), SSTs were significantly cooler in 2010 (–0.94 °C, 474 p= 0.002), 2013 (–0.61 °C, p= 0.031), and 2024 (–0.91 °C, p= 0.022). The effect in 2005 was 475 marginally significant (–0.56 °C, p= 0.051). This cooling effect is also reflected on the spatio476 temporal modelling (11 and 12), and gives some proxy about the geographic areas where the 477 cooling effects are evident and if there is a pattern on a specific area. For example, in 2010 478 and 2013, the northest area of the San Antonio Cape seems to be the area where its observed 479 this cooling. Not only in this years, but also in the 1997,2002,2007,2014,2015,2016 and 2020 480 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 19 of 24 Figure 13. Sen’s slope estimates for the long-term trends in SST maxima (panel a) and SST ranges (panel b) across the study area. The slopes quantify annual change (°C/year), highlighting areas with significant warming or increased variability. This analysis complements the spatio-temporal model predictions, confirming the seasonal and inter-annual patterns of warming and variability. Figure 14. Estimated smooth effects of cyclic month, AMO, ENSO index, and lagged NAO on monthly SST, obtained from the GAM. Each panel shows the partial effect of the corresponding predictor, with the solid line representing the estimated smooth function and the shaded area its 95% confidence interval (calculated using the Bayesian covariance matrix of the smooth). The y-axis indicates the deviation in SST (°C) relative to the model intercept. A smooth is considered statistically significant if its confidence band does not overlap zero and/or if the associated approximate significance test (summary output) yields p<0.05. https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 20 of 24 Term EDF Ref.df F p-valor s(month)7.688 8.000 1527.703 p < 0.001 s(amo_value)1.000 1.001 36.031 p < 0.001 s(enso_value)1.000 1.001 589 0.443 s(nao_lag1)1.849 2.388 2.264 0.930 Table 10. Summary of smooth effects from the GAM of monthly SST.The estimated degrees of freedom (edf) indicate the complexity of the smooth function(EDF ≈ 1 corresponds to a linear effect, higher values indicate non-linear responses).F-tests are approximate significance tests for each smooth term. SST shows a strong seasonal cycle ( s(month) ) and a significant positive association with AMO. ENSO and lagged NAO did not contribute significantly. Figure 15. Forest plot of estimated annual effects (±95% confidence intervals) of year on monthly sea surface temperature (SST), relative to the baseline year 1995. Points represent model coefficients from the GAM, with vertical lines indicating approximate 95% confidence intervals calculated as estimate ± 1.96 × SE. Years shown in red are significantly lower than the baseline ( p< 0.05.), based on parametric tests in the GAM summary. Non-significant effects are shown in black. The dashed horizontal line marks zero, i.e., no difference from the reference year. events if the mean effect cooling effect is not significant but this specific area remain as a 481 climate refuge. Along the same lines of identifying temperature-dampening areas, the Sen’s 482 slope 13 results on the map indicate that the area north of the cape—particularly along 483 the coastline—shows a less pronounced SST warming trend compared to the surrounding 484 region. 485 5. Conclusion 486 The importance of spatially-aware models, such as INLA, lies in their ability to account 487 for spatial correlations when estimating SST maxima. This is particularly relevant in the 488 Mediterranean, where both marine heatwaves (MHWs) and the maximum temperatures 489 reached in various areas are influenced by underlying climate and oceanographic factors. 490 Accurately capturing this spatial influence is therefore essential for predicting which regions 491 are more likely to experience cumulative effects during MHW events, for example. When 492 referring to the biological limits of a species, we are not only considering thresholds that 493 cause immediate mortality, but also those that may lead to metabolic stress. Such limits may 494 involve high rates of environmental change or seasonal extremes which, when exceeded 495 for prolonged periods, can result in physiological disorders and metabolic failure. 496 In the future, oceanographic dynamics, climate indexes and different climate change 497 scenarios will play a key role in determining which SST trends are going to dominate 498 in the Mediterranean sea. It’s because of that, that spatio temporal regional studies are 499 necessary in order to generate more accurate predictor models, with relevant and variable 500 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 21 of 24 scenarios adapted to the global climatic system. The most relevant information obtained 501 from these spatio-temporal models, as applied to our study area, is as follows: (A) Those 502 regions that are currently experiencing problems with high temperatures (like Campello 503 and Altea) may not be as troublesome in the future as those that expect less vertical water 504 movement (like Valencia and Sagunto). (B) The spatial variability of the SST along the 505 years leave us a footage of variable warming trends, which with a better understanding of 506 seasonal hydrodynamic can define climate refuges in the north of the San Antonio Cape. 507 (C) The climate indices provide valuable insights into the variability and trends of SST in 508 the Mediterranean Sea; however, their use does not currently allow for accurate predictive 509 modelling under extreme conditions. 510 Therefore, a precautionary approach is recommended when considering the migration 511 of aquaculture production zones towards the northern Valencian coast, as this region is 512 projected to experience more accelerated warming patterns over the next 10 years compared 513 to the southern areas. 514 Author Contributions: Conceptualization, L.A. and I.L.-M.; methodology, L.A., I.L.-M., A.L.-Q. 515 and X.B.; software, L.A. and I.L.-M.; validation, L.A., I.L.-M. and X.B.; formal analysis, L.A. and 516 I.L.-M.; investigation, L.A., I.L.-M., A.L.-Q., X.B., A.F. and P.S.-J.; resources, X.B.; data curation, L.A. 517 and I.L.-M.; writing—original draft preparation, L.A. and I.L.-M.; writing—review and editing, J.A. 518 (Javier Atalah), J.A. (Juan Aparicio), D.B., D.C., J.G.-M., A.L.-Q., P.S.-J. and X.B.; visualization, L.A. 519 and I.L.-M.; supervision, X.B.; project administration, X.B.; funding acquisition, X.B. All authors have 520 read and agreed to the published version of the manuscript. 521 Funding: This study forms part of the ThinkInAzul programme (https://thinkinazul.es/) and 522 was supported by MCIN with funding from the European Union NextGenerationEU (PRTR-C17.I1) 523 and Generalitat Valenciana (THINKINAZUL/2021/044-TOWARDS and THINKINAZUL/2021/021524 MODESTA). This paper is part of the project PID2022-136455NB-I00, funded by Ministerio de 525 Ciencia, Innovación y Universidades of Spain (MCIN/AEI/10.13039/501100011033/FEDER, UE) 526 and the European Regional Development Fund. This research was also partially supported by the 527 CIPROM/2024/34 grant, funded by the Conselleria de Educación, Cultura, Universidades y Empleo, 528 Generalitat Valenciana. 529 Institutional Review Board Statement: Not applicable. 530 Informed Consent Statement: Not applicable. 531 Data Availability Statement: The environmental data used in this study are available from the 532 Copernicus Marine Service (CMEMS). Mortality data are subject to confidentiality agreements. 533 Acknowledgments: We thank the General Directorate of Agriculture, Livestock and Fisheries of the 534 Generalitat Valenciana for providing the mortality datasets. 535 Conflicts of Interest: The authors declare no conflicts of interest. 536 Abbreviations 537 The following abbreviations are used in this manuscript: 538 539 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 22 of 24 SST Sea Surface Temperature ARIMA Autoregressive Integrated Moving Average BHM Bayesian Hierarchical Model INLA Integrated Nested Laplace Approximation GAM Generalized Additive Model NAO North Atlantic Oscillation AMO Atlantic Multidecadal Oscillation ENSO El Niño–Southern Oscillation ONI Oceanic Niño Index MHW Marine Heatwave EMS Extreme Marine Summer SPDE Stochastic Partial Differential Equation AZA Available Zones for Aquaculture REML Restricted Maximum Likelihood FAO Food and Agriculture Organization NOAA National Oceanic and Atmospheric Administration CMEMS Copernicus Marine Environment Monitoring Service 540 References 541 1. FAO. El estado mundial de la pesca y la acuicultura (SOFIA); Number 2024 in El estado mundial de la pesca y la acuicultura (SOFIA), 542 FAO: Roma, Italia, 2024; p. 278. 543 2. Jo, A.; Lee, J.Y.; Sharma, S.; Lee, S.S. Season-Dependent Atmosphere-Ocean Coupled Processes Driving SST Seasonality Changes 544 in a Warmer Climate. Geophysical Research Letters 2024,51.https://doi.org/10.1029/2023GL106953.545 3. López García, M.J. Recent warming in the Balearic Sea and Spanish Mediterranean coast. Towards an earlier and longer summer. 546 Atmósfera 2015,28, 149–160. https://doi.org/https://doi.org/10.20937/ATM.2015.28.03.01.547 4. Cherif, S.; Doblas-Miranda, E.; Lionello, P.; Borrego, C.; Giorgi, F.; Iglesias, A.; Jebari, S.; Mahmoudi, E.; Moriondo, M.; Pringault, 548 O.; et al. First Mediterranean Assessment Report - Chapter 2: Drivers of Change, 2022. https://doi.org/10.5281/zenodo.7100601. 549 5. Macias, J.; Avila Zaragozá, P.; Karakassis, I.; Sanchez-Jerez, P.; Massa, F.; Fezzardi, D.; Yücel Gier, G.; Franiˇcevi´c, V.; Borg, J.; 550 Chapela Pérez, R.; et al. Allocated zones for aquaculture: a guide for the establishment of coastal zones dedicated to aquaculture in the 551 Mediterranean and the Black Sea; Number 97 in General Fisheries Commission for the Mediterranean. Studies and Reviews, FAO: 552 Rome, 2019; p. 90. 553 6. Rosa, R.; Marques, A.; Nunes, M.L. Impact of climate change in Mediterranean aquaculture. Reviews in Aquaculture 2012, 554 4, 163–177. https://doi.org/10.1111/j.1753-5131.2012.01071.x.555 7. Islam, M.J.; Kunzmann, A.; Slater, M.J. Responses of aquaculture fish to climate change-induced extreme temperatures: A review. 556 Journal of the World Aquaculture Society 2021,53, 314–366. https://doi.org/10.1111/jwas.12853.557 8. López-Mengual, I.; Sanchez-Jerez, P.; Ballester-Berman, J.D. Offshore aquaculture as climate change adaptation in coastal 558 areas: sea surface temperature trends in the Western Mediterranean Sea. Aquaculture Environment Interactions 2021.https: 559 //doi.org/10.3354/aei00420.560 9. Stavrakidis-Zachou, O.; Lika, K.; Anastasiadis, P.; Papandroulakis, N. Projecting climate change impacts on Mediterranean finfish 561 production: a case study in Greece. Climatic Change 2021,165.https://doi.org/10.1007/s10584-021-03096-y.562 10. Haberle, I.; Hackenberger, D.K.; Djerdj, T.; Bavˇcevi´c, L.; Geˇcek, S.; Hackenberger, B.K.; Marn, N.; Klanjšˇcek, J.; Purgar, M.; Ili´c, 563 J.P.; et al. Effects of climate change on gilthead seabream aquaculture in the Mediterranean. Aquaculture 2024,578, 740052. 564 https://doi.org/10.1016/j.aquaculture.2023.740052.565 11. Pastor, F.; Valiente, J.; Palau, J.L. Sea Surface Temperature in the Mediterranean: Trends and Spatial Patterns (1982–2016). Pure 566 and Applied Geophysics 2020.https://doi.org/10.1038/npg.els.0003276.567 12. Pinardi, N.; Masetti, E. Variability of the large scale general circulation of the Mediterranean Sea from observations and modelling: 568 a review. Palaeogeography, Palaeoclimatology, Palaeoecology 2000,158, 153–173. 569 13. Lange, H. Time-series Analysis in Ecology. eLS. Wiley 2006.https://doi.org/10.1038/npg.els.0003276.570 14. Pisano, A.; Buongiorno Nardelli, B.; Tronconi, C.; Santoleri, R. The new Mediterranean optimally interpolated pathfinder AVHRR 571 SST Dataset (1982–2012). Remote Sensing of Environment 2016,176, 107–116. DOI PRODUCT: https://doi.org/10.48670/moi-00173, 572 https://doi.org/10.1016/j.rse.2016.01.019.573 15. Embury, O.; Merchant, C.J.; Good, S.A.; Rayner, N.A.; Høyer, J.L.; Atkinson, C.; Block, T.; Alerskans, E.; Pearson, K.J.; Worsfold, 574 M.; et al. Satellite-based time-series of sea-surface temperature since 1980 for climate applications. Scientific Data 2024,11. 575 https://doi.org/10.1038/s41597-024-03147-w.576 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 23 of 24 16. NOAA/ESRL Physical Sciences Laboratory. Daily North Atlantic Oscillation (NAO) Index Time Series. https://psl.noaa.gov/ 577 data/timeseries/daily/NAO/, 2025. Accessed: 2025-04-08. 578 17. NOAA Physical Sciences Laboratory. Atlantic Multidecadal Oscillation (AMO) Index. https://psl.noaa.gov/data/timeseries/ 579 AMO/, 2024. Accessed: 2025-04-08. 580 18. NOAA National Centers for Environmental Information. El Niño / Southern Oscillation (ENSO): Sea Surface Temperature 581 Monitoring. https://www.ncei.noaa.gov/access/monitoring/enso/sst, 2025. Accessed: 2025-04-08. 582 19. Grolemund, G.; Wickham, H. lubridate: Make Dealing with Dates a Little Easier, 2011. R package version 1.9.3. 583 20. Xiao, C.; Chen, N.; Hu, C.; Wang, K.; Xu, Z.; Cai, Y.; Xu, L.; Chen, Z.; Gong, J. A spatiotemporal deep learning model 584 for sea surface temperature field prediction using time-series satellite data. Environmental Modelling Software 2019.https: 585 //doi.org/10.1016/j.envsoft.2019.104502.586 21. Xiao, C.; Chen, N.; Hu, C.; Wang, K.; Gong, J.; Chen, Z. Short and mid-term sea surface temperature prediction using time-series 587 satellite data and LSTM-AdaBoost combination approach. Remote Sensing of Environment 2019.https://doi.org/10.1016/j.rse.20 588 19.111358.589 22. Dunstan, P.K.; Foster, S.D.; King, E.; Risbey, J.; O’Kane1, T.J.; Monselesan, D.; Hobday, A.J.; Hartog1, J.R.; Thompson, P.A. Global 590 patterns of change and variation in sea surface temperature and chlorophylla. Scientific Reports 2018.https://doi.org/10.1038/s4 591 1598-018-33057-y.592 23. Kisi, O.; Choubin, B.; Deo, R.C.; Yaseen, Z.M. Incorporating synoptic-scale climate signals for streamflow modelling over the 593 Mediterranean region using machine learning models. Hydrological Sciences Journal 2019.https://doi.org/10.1038/s41598-018-3 594 3057-y.595 24. Box, G.E.; Jenkins, G.M. Time Series Analysis: Forecasting and Control; Holden-Day: San Francisco, 1970. 596 25. Robitzsch, A.; Luedtke, O. Functions for the STARTS Model, 2022. R package version 1.3-8, https://doi.org/10.32614/CRAN. 597 package.STARTS.598 26. Good, I.J. Some history of the hierarchical Bayesian methodology. Trabajos de Estadistica Y de Investigacion Operativa 1980, 599 31, 489–519. https://doi.org/10.1007/bf02888365.600 27. Sen, P.K. Estimates of the regression coefficient based on Kendall’s tau. Journal of the American Statistical Association 1968, 601 63, 1379–1389. 602 28. Brugnara, S. trend: Non-Parametric Trend Tests and Sen’s Slope Estimator, 2024. R package version 1.2.1. 603 29. Mann, H.B. Nonparametric tests against trend. Econometrica 1945,13, 245–259. 604 30. Kendall, M.G. Rank Correlation Methods, 4th ed.; Charles Griffin, London, 1975. 605 31. Hirsch, R.M.; Slack, J.R.; Smith, R.A. Techniques of trend analysis for monthly water quality data. Water Resources Research 1982, 606 18, 107–121. 607 32. Hastie, T.J.; Tibshirani, R.J. Generalized Additive Models; Chapman & Hall: New York, 1990. 608 33. Wood, S.N. Generalized Additive Models: An Introduction with R, 2nd ed.; Chapman & Hall/CRC: Boca Raton, FL, 2017. https: 609 //doi.org/10.1201/9781315370279.610 34. Wood, S.N. mgcv: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation, 2025. R package version 1.9-3. 611 35. Wood, S.N. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear 612 models. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2011,73, 3–36. https://doi.org/10.1111/j.1467-9868. 613 2010.00749.x.614 36. Tian, B.; Fan, K. A Skillful Prediction Model for Winter NAO Based on Atlantic Sea Surface Temperature and Eurasian Snow 615 Cover. Weather and Forecasting 2015,30, 197–205. https://doi.org/10.1175/waf-d-14-00100.1.616 37. Pisano, A.; Marullo, S.; Artale, V.; Falcini, F.; Yang, C.; Leonelli, F.E.; Santoleri, R.; Buongiorno Nardelli, B. New Evidence of 617 Mediterranean Climate Change and Variability from Sea Surface Temperature Observations. Remote Sensing 2020,12.618 38. Copernicus Climate Change Service. Global Climate Highlights 2024, 2025. Accedido: 2025-07-02. 619 39. Androulidakis, Y.; Pytharoulis, I. Variability of marine heatwaves and atmospheric cyclones in the Mediterranean Sea during the 620 last four decades. Environmental Research Letters 2025,20, 034031. https://doi.org/10.1088/1748-9326/adb505.621 40. Atalah, J.; Ibañez, S.; Aixalà, L.; Barber, X.; Sánchez-Jerez, P. Marine heatwaves in the western Mediterranean: Considerations for 622 coastal aquaculture adaptation. Aquaculture 2024,588, 740917. https://doi.org/10.1016/j.aquaculture.2024.740917.623 41. Denaxa, D.; Korres, G.; Flaounas, E.; Hatzaki, M. Investigating extreme marine summers in the Mediterranean Sea. Ocean Science 624 2024,20, 433–461. https://doi.org/10.5194/os-20-433-2024.625 42. García-Monteiro, S.; Sobrino, J.; Julien, Y.; Sòria, G.; Skokovic, D. Surface Temperature trends in the Mediterranean Sea from 626 MODIS data during years 2003–2019. Regional Studies in Marine Science 2022,49, 102086. https://doi.org/10.1016/j.rsma.2021.102 627 086.628 43. Simon, A.; Plecha, S.M.; Russo, A.; Teles-Machado, A.; Donat, M.G.; Auger, P.A.; Trigo, R.M. Hot and cold marine extreme events 629 in the Mediterranean over the period 1982-2021. Frontiers in Marine Science 2022,9.https://doi.org/10.3389/fmars.2022.892201. 630 https://doi.org/10.3390/1010000 Version December 16, 2025 submitted to Journal Not Specified 24 of 24 44. Martínez, J.; Leonelli, F.E.; García-Ladona, E.; Garrabou, J.; Kersting, D.K.; Bensoussan, N.; Pisano, A. Evolution of marine 631 heatwaves in warming seas: the Mediterranean Sea case study. Frontiers in Marine Science 2023,10.https://doi.org/10.3389/ 632 fmars.2023.1193164.633 Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual 634 author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to 635 people or property resulting from any ideas, methods, instructions or products referred to in the content. 636 https://doi.org/10.3390/1010000