Temperature Trends and Seasonality in Neritic and Transitional Waters of the Southern Bay of Biscay from 1998 to 2023: Relationships with Local Factors and Climate-Ocean Teleconnections

1. Introduction

Temperature is a key driver of biological processes from molecular to ecosystem levels and varies over daily to multidecadal time scales. Much effort is being devoted to the development of theories and the prediction of biological responses to climate change across different levels of biological organization, and on different spatial and temporal scales [1,2,3]. A major focus of attention is the analysis of spatial and temporal temperature patterns in terrestrial and marine ecosystems in the context of the expected climate warming throughout the 21st century [4]. Globally, sea surface temperature (SST) is increasing [5], although there are some exceptions [6] and marked interregional differences in warming rates [7], which differ not only between ocean basins [8] but also at regional scales [9,10,11]. Coastal warming is highly heterogeneous, with temperature trends differing across continental shelves, and from coastal to open-ocean locations [12,13,14]. Beyond these systems, studies in Australia [15] and South Florida [16] have shown that estuarine waters are warming faster than adjacent shelf areas and at a greater rate than predicted by global ocean and atmospheric models, suggesting that these models may not accurately predict estuarine temperature trends [15]. A recent global assessment of estuarine surface water warming rates has shown that less than half of the studied estuaries have exhibited warming trends over the recent 30-yr period (1985-2022). Additionally, estuarine water warming trends have shown significant differences across hemispheres, continents, and between small and large estuaries [17].

Moreover, climate change affects not only mean SST trends but also variability, extreme values and seasonality, although these aspects remain less explored [6]. Models predict regional variability in SST warming in large marine systems such as the Baltic Sea, with clear latitude dependent seasonal warming differences [18,19]. Both predicted and observed water temperature changes also show spatial and seasonal differences in estuaries [20,21].

However, aquatic systems have a tridimensional structure and functioning, which cannot be shaped only from changes observed at surface. Studies on long-term water temperature variability with depth are more limited due to the scarcity of data at subsurface and bottom layers, but the existing evidence suggests that warming trends differ across depths, due to the influence of different forcing factors at each layer [22]. Similarly, the intensity, duration and spatial extension of heatwaves vary with water depth [23,24]. These differences are expected to be larger in coastal areas and estuaries, which are naturally prone to experiencing greater temperature fluctuations at smaller spatial (regional and local) and temporal scales than offshore environments, and therefore they are more likely to undergo increased ecological stress [25].

Estuaries are ecologically and economically valuable systems [26] but due to their high warming response to the current climate crisis [15,16], the biodiversity and uniqueness of estuarine biota, as well as the commercial and recreational services of estuaries are threatened. The rise in temperature is a common driver of water quality impairment in estuaries [27]. It can enhance the severity and extent of hypoxia [28] as well as the level of pathogens and infectious diseases [29], and increase the frequency, severity, and duration of harmful algal blooms [30]. Extreme thermal events can also affect fisheries, due to catastrophic fish mortality [31]. However, estuaries vary widely in size, morphology and hydrodynamics, and therefore encompass a broad range of vulnerabilities to climate warming [15]. Beyond atmospheric drivers, estuarine water temperature is affected by external factors such as upstream river temperature and discharge, marine water intrusions and coastal upwelling, as well as by internal processes such as mixing, dispersion, and water residence time e.g., [32,33,34,35]. This complexity in temperature patterns across spatial and temporal scales leads to different vulnerabilities to warming in the context of climate change, with implications for estuarine functions now and in the future [36]. Therefore, disentangling the spatial and temporal patterns and rates of temperature change in different types of estuaries and identifying their main driving factors is essential to improve our understanding of future environmental and biotic changes. Such insights are vital for informing management actions in these dynamic systems.

The aim of this study was to analyse water temperature variability at interannual and seasonal scales, as well as along axial and vertical spatial gradients, in two estuaries of the Basque coast, the estuary of Bilbao and the estuary of Urdaibai (Bay of Biscay), during the last 26-yr period (1998-2023). The specific objectives were the following

a) To identify distinct water mass clusters based on temperature variability similarities.

b) To describe long-term water temperature trends and assess rates of change.

c) To establish seasonal patterns, phenological changes, and seasonal and monthly differences in interannual temperature variations.

d) To determine the relationship of air temperature, river flow and upwelling events and the main climate teleconnections with water temperature variability at the different spatial and temporal scales studied.

e) To assess the role of estuarine morphology in shaping thermal conditions in the estuaries under study.

2. Materials and Methods

2.1. Study Area

The nearby estuaries of Bilbao (43° 23′ N, 03° 07′ W) and Urdaibai (43° 22’ N, 02° 43′ W) are located in the Basque coast (inner Bay of Biscay), within the eastern North Atlantic temperate region (Figure 1a). They share a temperate oceanic climate characterized by mild winters, warm summers, and frequent rainfall [37]. However, they differ largely in morphology, topography, and watershed characteristics [38]. The estuary of Bilbao is a highly man-modified system currently turned into an artificial channel 25 m wide and 0.5 m deep in the upper fluvial reaches, 270 m wide in the widest middle reaches and 10 m deep in the deepest lower reaches. This channel meets the Abra harbour, a semiconfined funnel-shaped embayment that increases in depth from around 12 m in the inner Abra to 32 m at the coastline (Figure 1b). Extensive land reclamation has greatly reduced intertidal areas, and dredging has deepened the estuary for navigation [39]. In contrast, the estuary of Urdaibai has remained much more natural, with an average depth of ~3 m at high tide and a greater proportion of intertidal and supratidal zones. However, a straight artificial channel was built alongside the natural meandering one in the upper estuary to facilitate navigation (Figure 1c).

Both estuaries experience a meso-macrotidal, semidiurnal tidal regime, but the estuary of Urdaibai has a much higher tidal prism volume/total volume ratio [40]. The estuary of Bilbao receives greater freshwater inflows, leading to a stratified system with a two-layer circulation [41] and a well-defined estuarine plume drifting offshore [42], whereas the estuary of Urdaibai is mostly a well-mixed system. Under usual river discharge conditions, both estuaries are sea-dominated systems with high salinity (euhaline > 30) water masses reaching the upper estuary below the halocline in the estuary of Bilbao and filling the outer half of the estuary of Urdaibai at high tide. However, the torrential nature of Basque coastal rivers [43] means that extreme flood events can for brief periods of time flush freshwater to the estuary mouths.

In the adjacent marine region, downwelling events prevail in autumn and winter due to southwesterly winds, while frequent episodes of weak upwelling occur in spring and summer under easterly and northeasterly winds [44].

2.2. Data Acquisition

Water temperature data were obtained from an ongoing hydrography and plankton monitoring program in the estuaries of Bilbao and Urdaibai. A Lagrangian sampling approach was used to assess changes in water masses of specific salinities instead of samplings at fixed locations. From 1998 to 2023, monthly measurements were taken at sites with salinities of 35, 34, 33, and 30 in the estuary of Bilbao and 35, 33, 30, and 26 in the estuary of Urdaibai, all of them below the halocline. These salinity sites were selected based on the axial distribution of water masses, ranging from the main salinity habitat (~35 salinity) for neritic plankton to brackish plankton habitats at varying salinities in each estuary. At each sampling site, water temperature was measured at 0.5 m intervals from just below the surface (where the sensor became submerged; hereafter referred to as surface (0.0 m)) to the bottom using different multiparameter meters (WTW and YSI) during the study period. To minimize tidal currents interference and short-term physical stressor variability, sampling was conducted, whenever possible, in the morning, close to the high tide slack, and during the second neap tide of the month.

To account for local influences, air temperature data for the marine area adjacent to the estuaries were obtained from the Copernicus Climate Data Store (https://cds.climate.copernicus.eu/datasets). Additionally, air temperature over land was sourced from three hydrometeorological stations (Figure 1a): Sondika (https://www.aemet.es), located between the basins of the two estuaries; Abusu (https://www.euskalmet.euskadi.eus/hasiera/), near the head of the estuary of Bilbao; and Muxika (https://www.euskalmet.euskadi.eus/hasiera/), near the head of the estuary of Urdaibai. River flow data were also obtained from the latter two gauging stations (https://www.uragentzia.euskadi.eus/), representing the combined flow of the Ibaizabal and Nerbioi rivers in the estuary of Bilbao and the Oka river in the estuary of Urdaibai. Upwelling index data were accessed from the Spanish Institute of Oceanography (http://www.indicedeafloramiento.ieo.es/interactivo.html).

Finally, teleconnection indices (NAO: North Atlantic Oscillation, EA: East Atlantic pattern, AMO: Atlantic Multidecadal Oscillation, and AO: Arctic Oscillation) were obtained from the NOAA Climate Prediction Center (https://www.cpc.ncep.noaa.gov/).

2.3. Data Selection

23 water temperature time series (combinations of sites and depths) were used for the present study. The selected temperature measurements included surface (0.0 m) and 2.5 m depths at all salinity sites in both estuaries, and 5.0 m and 7.5 m depths at all salinity sites in the estuary of Bilbao except for 30 salinity at 7.5 m, where this depth was generally unavailable. Due to the shallowness of the estuary of Urdaibai, data below 2.5 m were not consistently available. Similarly, in the estuary of Bilbao, measurements at 5.0 m and 7.5 m were occasionally missing due to shallow water depths or cable length limitation (8 m) of the multiparameter meter used until 2003. To fill these gaps at some depths, missing values were primarily replaced with values from the depth nearest to the missing depth. Likewise, for other occasional missing data, two gap-filling methods were applied: 1) When a predefined salinity was not found, or temperature data were missing due to instrument failure (overall <4.5%), missing values of temperature were estimated from polynomial regression models of temperature variation as a function of salinity. To this purpose data from other sampling sites of the same date were used; 2) When a sampling was cancelled, missing values were estimated using the average of the preceding and following sampling dates (< 2.5%).

Finally, the criterion to establish the months included in each season was that the three coldest and three warmest consecutive months in the highest salinity time series are winter and summer, respectively. As a result, January, February and March were considered as winter, April, May and June as spring, July, August and September as summer and October, November and December as winter.

2.4. Data Processing

First, the interannual trend component of the 23 water temperature time series, the afore mentioned local hydrometeorological factors and the teleconnection indices were extracted using the multiple seasonal decomposition function (mstl) from the forecast v.8.21 package in R [45].

According to similarities in interannual trends, the 23 water temperature time series were grouped into six distinct water mass clusters using Ward's hierarchical clustering method, based on Spearman’s rank correlation (r = 0.8), with clusters defined at a 0.3 height threshold. For each sampling date, mean temperatures were calculated for each cluster of time series. To identify breakpoints in the variations of these mean temperatures throughout the study period, a Cumulative Sum (CuSum) analysis was conducted for each water mass cluster.

Additionally, the monthly mean for the entire study period were calculated for each water temperature time series to establish the standard seasonal pattern.

Linear regression analysis was performed to estimate temperature change rates per decade and assess their statistical significance in the trend component of each water temperature time series. The same analysis was applied to temperature data for each water mass cluster split by seasons.

Finally, Spearman rank correlation analyses were conducted to explore relationships between water temperature and local environmental factors and North Atlantic climate teleconnections.

All statistical analyses were performed by means of R software [46].

3. Results

3.1. Water Temperature Time-Series Clustering

Hierarchical clustering analysis of the 23 water temperature time series identified six water mass clusters based on similarities in interannual trends (Figure 2) herein termed as:

1) BTD (Bilbao-Transitional-Deep). It includes time series from the intermediate and inner (transitional) waters of the estuary of Bilbao, specifically from below the halocline (deep). This water mass cluster comprises data from 2.5, 5.0, and 7.5 m depths at the 34 salinity site, 5.0 and 7.5 m at the 33 salinity site, and 5.0 m at the 30 salinity site.

2) UN (Urdaibai-Neritic). It includes time series from the neritic waters of the 35 salinity site (both surface and 2.5 m depths) of the estuary of Urdaibai.

3) BND (Bilbao-Neritic-Deep). It includes time series from the outer zone (neritic) of the estuary of Bilbao, specifically from below the halocline (deep). This group comprises data from 2.5, 5.0 and 7.5 m at the 35 salinity site.

4) UT (Urdaibai-Transitional). It includes time series from the transitional waters of the estuary of Urdaibai, comprising data from surface and 2.5 m depths at the 33, 30, and 26 salinity sites.

5) BTS (Bilbao-Transitional-Surface). It includes time series from the intermediate and inner (transitional) surface waters of the estuary of Bilbao. The group comprises data from 0.0 m depth at the 34, 33, and 30 salinity sites.

6) BTG (Bilbao-Transitional-Gradient). It is a mixed group from the estuary of Bilbao, which includes the time series from the 35 salinity site at surface and the 33 and 30 salinity sites at 2.5 m depth.

Figure 2. Similarity dendrogram of water temperatures from the 23 time series, constructed using Ward's hierarchical clustering method based on Spearman's rho correlations. Water masses (see text for definition of BTD, UN, BND, UT, BTS and BTG waters) are defined by a 0.3 height threshold (dotted line), with water clusters outlined by boxes.

Figure 2. Similarity dendrogram of water temperatures from the 23 time series, constructed using Ward's hierarchical clustering method based on Spearman's rho correlations. Water masses (see text for definition of BTD, UN, BND, UT, BTS and BTG waters) are defined by a 0.3 height threshold (dotted line), with water clusters outlined by boxes.

As shown in Figure 2, BND and UN were grouped into a single cluster at the set threshold and therefore they included data from both estuaries. However, this was an exception, as the rest of groups were related to a single estuary. Therefore, it was decided to treat BND and UN as separate clusters. At a lower similarity level, BTD merged with the neritic assemblage (BND and UN), while the other transitional waters (UT, BTS, and BTG) formed a separate group.

3.2. Interannual Trends in Water Temperature Across Water Masses

Interannual trends of water temperature at each depth, salinity site and estuary showed an overall significant increase in most cases during the study period (Figure 3 and Table 1).

Exceptions were observed in some surface waters of the estuary of Bilbao, where a clear downward trend (-0.3 °C decade^-1 at the 30 salinity site) or no trend (at the 33 and 35 salinity sites) was found. The strongest increases were observed at the 26 salinity site in the UT water masses and at the 33 salinity site in the BTD ones, with warming rates of 0.25 °C and 0.24 °C decade^-1, respectively. Additionally, UN waters showed a slightly higher increase (0.18-0.19 °C decade^-1) than BND ones (0.10-0.15 °C decade^-1).

However, the increase was not gradual. Water temperature trend showed a decrease across all time series from 2003-2006 to 2013-2015, followed by an increase in the last years of the study period. In addition, the warmest and coldest years varied between clusters. Annual mean values for each water mass (cluster) showed that 2023 was the warmest year in neritic waters (UN and BND) and 2022 in transitional waters (UT, BTG, BTS). In the first half of the study period temperature peaked in 2006 in neritic waters and 2003 in transitional waters. The coldest years were 2015 in neritic waters and 2013 in transitional waters. Temperature in BTD waters (nearest ones to the neritic waters of the estuary of Bilbao) were an exception because they aligned better with the neritic pattern (Figure 3).

The CuSum analysis of water temperature data calculated for the six water mass clusters (Figure 4) revealed five distinct periods: 1) 1998-2002, characterized by values below the mean in the main water masses (UN, BND, BTD, and UT); 2) 2003-2006, with values above the mean, especially in the earliest years in UT, BTG, and BTS, and in later years in BND, UN, and BTD; 3) 2007-2012, with values below the mean but closer to it in UT, BTG, and BTS waters; 4) 2013-2019, characterized by the lowest values below the mean across all waters; and 5) 2020-2023, with values significantly above the mean.

The linear interannual trends and the periods of increase and decrease identified by the CuSum analysis primarily reflected the interannual patterns of spring (April-June) water temperature, which showed stronger variations and higher rates of increase than other seasons temperatures (Figure 5 and Table 2). Indeed, spring showed the strongest warming across all water masses, with UT showing the largest spring increase (0.49 °C decade^-1) and BTS waters the smallest one (0.17 °C decade^-1). Autumn also showed considerable warming (0.28 to 0.22 °C decade^-1 in UN, UT and BTD clusters).

In contrast, summer and winter displayed greater variability, with summer showing both warming (e.g., in BTD waters 0.23 °C decade^-1) and cooling (e.g., minimum rise in UN waters: -0.10 °C decade^-1), while winter had the most inconsistent patterns, with predominantly negative trends. However, none of these relationships were statistically significant. Similarly, the interannual rates of change of water temperature at the monthly scale (Figure 6) were positive for all spring months, with May showing the strongest rates across most water masses (0.88 °C decade^-1 in UT, and between 0.5 and 0.6 °C decade^-1 in BND, BTD, and UN). In autumn months, highest rates were observed in November or December (e.g. 0.63 °C decade^-1 in December in UT). Negative change rates were primarily observed in summer and winter months, especially in August for UN (-0.37 °C decade^-1) and in March for BTS (-0.60 °C decade^-1).

3.3. Seasonal Patterns in Water Temperature Across Water Masses

Overall, as depth and salinity increased, winter temperatures were higher, summer temperatures lower, and the annual temperature range narrowed (Figure 7). An exception was observed in the estuary of Bilbao, where summer maxima were higher at the surface of the 33 salinity site compared to the 30 salinity site. Also, the estuary of Urdaibai showed higher temperatures at the same depths and salinity sites than the estuary of Bilbao. The ranges of average water temperatures of each water mass cluster (Figure 7) were the following: from 12.2-12.4 to 20.6-21.1 °C in BND waters, from 12.0-12.2 to 21.7-21.8 °C in UN ones, from 11.9-12.4 to 21.0-21.9 °C in BTD ones, from 10.1-11.6 to 22.5-23.4 °C in UT ones, and from 8.4-10.9 to 21.5-22.6 °C for the joint BTG plus BTS waters.

The timing of the annual water temperature minima and maxima, as well as the coldest and warmest periods varied between the different water masses. Coldest periods occurred from January to March, with annual minima in February for the BND, UN, and BTD waters, while for BTS and UT ones the coldest periods were from December to February, with annual minima in January. Annual maxima were skewed toward July in UN, UT and BTS waters (Figure 7). Spring water temperatures were always higher than autumn temperatures, with the smallest differences observed in deep water masses of the estuary of Bilbao (BND and BTD), and the highest in inner surface ones (BTS) (Figure 7). In the estuary of Urdaibai, such seasonal differences were also found from UN to UT waters with an increase in the difference between spring and autumn temperatures and decreases in the temperature difference between spring and summer, and between autumn and winter.

In all water masses temperature seasonality (Figure 8) was similar in the five time periods identified by the CuSum analysis. The most significant variation over time of the annual cycle was the delay of the coldest winter values towards March and the advance of the warmest summer values towards July, which led to a shortening of the seasonal warming period (although this process did not occur in a totally gradual manner). The delay of the winter minima was more evident in the neritic waters of both systems and the transitional waters of the estuary of Urdaibai, where a faster and greater warming from April to May was also observed in the last period (2020-2023). In contrast, the first period (1998-2002) showed the lowest May temperatures and the coldest springs in all water masses, except in BTS ones. The 2003-2006 period was characterized by the coldest winters (at least in January and February or February and March) in all water masses and warmest summers (except in BTD waters), and therefore, the highest temperature ranges. The 2007-2012 time span showed the coldest autumns in all water masses. The 2013-2019 period was characterized by water temperatures below the time series mean for all the seasons in all water masses, and by the lowest range of temperature variation between the warmest and coldest periods. In contrast, the 2020-2023 period showed water temperatures above the time series mean for all the seasons, with warmest springs and autumns in all water masses, warmest winters in UN, BND and BTD, and warmest summers in BTD.

3.4. Relationship of Water Temperature with Local Environmental Factors and Climate Teleconnections

Interannual trends of local environmental factors and teleconnection indices showed a significant increase over the study period (Figure 9 and Table 3). Regarding air temperatures, the rate of change varied considerably across meteorological station sources (Table 3).

At the monthly scale (Figure 10), warming rates were quite similar between over-land and over-ocean air temperature data. The assessment of the month specific rate of warming/cooling showed an overall increasing trend of warming rates from winter (with decreasing rates in January) to autumn, with highest warming rates in November-December. Similarly, at the monthly scale, river flow trends showed an almost coincident seasonal pattern between basins; with increasing rates in late autumn and early winter (from November to February) and clear decreasing rates in early spring (March-April) and early autumn (October). The upwelling index showed a somewhat opposite pattern to that of river flow, with a marked increase in April and October but a decline in January.

Regarding the relationship between environmental factors and water temperature (Table 4), air temperature showed the highest positive correlations across all cases, with over ocean air temperature data (ATc) being the most consistent. However, over-land-air temperatures showed stronger correlations during winter and/or summer in some water masses, the air temperature measured at the head of the estuaries (ATa/m) being the most strongly correlated factor with water temperature in all transitional water masses (BTD, UT, BTG and BTS). River flow showed significant negative correlations, except during the summer in the estuary of Urdaibai (UN and UT) and below the halocline in the estuary of Bilbao (BND and BTD), where no correlations were observed. The upwelling index was positively correlated with water temperature in spring in all water masses. In contrast, overall teleconnection indices showed fewer and less significant correlations. AMO index showed significant positive correlations with temperature in all water masses and generally in most seasons except summer. Negative correlations were found for NAO index, but only in autumn. EA pattern only correlated with water temperature in summer in deep waters of the estuary of Bilbao (BND and BTD), and in winter in non-deep transitional waters (UT, BTG and BTS). Positive correlations with AO were only observed for BTD and BTG in winter.

4. Discussion

4.1. Differences in Thermal Environments

Water masses differed between and within estuaries in their interannual variability, seasonal pattern and range of temperature. Neritic waters of both estuaries showed a very similar pattern, likely influenced by the connectivity of waters along the adjacent coastal area. However, the thermal features of transitional waters responded differently in each estuary depending on the water circulation and mixing patterns, which are shaped by the morphology of the estuaries and the river discharges. The higher river discharge and deeper water column generate vertical stratification in the estuary of Bilbao, and this causes differences in thermal features between surface and deeper transitional waters. The deep transitional waters in the estuary of Bilbao are formed by the intrusion of high-salinity marine water that flows landwards [41]. The limited air temperature and river flow influences allow these waters to retain a neritic thermal behaviour, as shown in the present work. The surface transitional waters, however, exhibited a different thermal environment due to the persistence and extent of the stratification that allowed a greater influence of air temperatures. In contrast, under conditions of low river discharge and shallowness, mixing dominates in the estuary of Urdaibai and a similar thermal behaviour over the entire column was observed, influenced also by air temperatures. The distinction between the neritic thermal environment and the thermal environments of transitional waters in both estuaries agrees with the differentiation of water masses due to the dominance of neritic and brackish zooplankton species, respectively [38,47]. The thermal differences we found followed a spatial pattern commonly observed in estuaries, where upper reaches are consistently warmer than the lower estuary [33]. However, results also revealed that the main waters of the estuary of Bilbao were, on average, colder and more thermally stable than those in the estuary of Urdaibai. This implies differences in the thermal habitats of brackish communities in each estuary, which could partly explain the observed between-estuary differences in the dominance and phenology of the main zooplankton brackish species [48].

4.2. Overall Trends and Change Rates

Water temperature increased during the 1998-2023 period in all the water masses, except for the lowest-salinity above-halocline waters of the inner estuary of Bilbao, where it decreased. Our data further suggests that warming trends may reverse if poorly mixed freshwater masses expand at surface towards lower reaches if river flows increase. In estuaries, freshwater discharge generally affects upper layers more than bottom layers, and freshwater runoff can influence the magnitude and spatial distribution of heating or cooling rates, especially in frontal areas [49]. This agrees with observations in the present study, where the correlation between water temperature and river flow was higher in transitional gradient waters than in bottom or inner surface transitional waters during the rainiest seasons in the stratified estuary of Bilbao. As observed in our study, sites within estuaries that deviate from the prevailing warming trend have been reported from many other estuaries [50,51]. Both of our systems showed the highest warming in transitional waters, with rates of 0.25-0.24 °C decade^-1 (considering data from all months in the series). The average warming rate from 1998 to 2023 in the main transitional waters of these two estuaries (0.21 °C decade^-1) was lower than the average increase of 0.37 °C decade^-1 reported for all modelled English estuaries between 1990 and 2022 [51]. However, our series lacks the colder years registered in the early 90’s (as shown by SST data from the nearby Aquarium of San Sebastian [52], which may bias the comparison towards lower warming rates. Nonetheless, interregional differences between our south-eastern Bay of Biscay estuaries and English estuaries are smaller than the intraregional differences (0.70-0.05 °C decade^-1) reported for the latter [51]. The higher warming rates in transitional water masses compared to the average 0.15 °C decade^-1 in neritic waters of our two systems is in agreement with the idea that estuaries are more sensitive to global warming than open sea waters [15,16]. Our results, however, also evidenced that warming rates in neritic waters decreased substantially from surface layers in the estuary of Urdaibai (0.19 °C decade^-1) to subsurface layers in the estuary of Bilbao (0.13 °C decade^-1). This is very likely due to the combined effect of depth and the west to east gradient of temperature increase that develops in the inner Bay of Biscay during the warm period [11]. In fact, in the Bay of Biscay the southeastern corner has shown the strongest warming trend [53], and our two neritic sites have shown progressively lower rates westward, but lower than the 0.23 °C decade^-1 SST rate of change obtained further east along the Basque coast from 1980 to 2019 [54]. However, as mentioned above, this comparison must be taken with caution due to the lack of data from the 80’s and most years of the 90`s in our series. The surface layer of neritic water masses in the estuary of Bilbao showed a thermal behaviour different from that of the subsurface layers. Similarly, in the more vertically homogeneous shelf waters, spatial and temporal differences in warming have been identified at contrasting depths [22].

4.3. Trends and Change Rates at Intra-Annual Scales

In our systems, the interannual trends and rates of change of seasonal and monthly water temperatures showed strongest warming in spring (April to June), lower in autumn (November/December), and lowest or even cooling trends in summer (August) and winter (February/March). This pattern is consistent with observations in the Bay of Biscay from 1982 to 2014 [11]. This latter study reported higher warming trends from September to November and from April to June and lower SST increases from January to March and from July to August. However, the seasonal and monthly rates of change they recorded (range: > 0.3 °C decade^-1 to < 0.2 °C decade^-1) were smaller than those in our two systems indicating that, while intra-annual variations in water temperature trends at our local scale are in accordance with the regional patterns, in the nearshore-inshore waters rates of change are higher.

Our results also confirm that warming was mainly driven by an extended warm season rather than milder winters or warmer summers [11]. This pattern appears to characterize the first quarter of the 21st century in the Bay of Biscay, but it does not match observations during the last decades of the 20th century. For example, during the 1965-2004 period, [55] reported a clear seasonal dependence of warming rates, with summer temperature trends being twice as strong as those of winter in surface waters. In the southeastern Bay of Biscay, from 1972 to 1993, both winter and summer showed warming trends, but the winter increase was slightly higher [53]. Likewise, along the Cantabrian coast, spring and summer showed similar warming trends from 1985 to 2005 [56]. Results from other marine (e.g. NW Mediterranean [57]; northeastern continental shelf of North America [13]) and estuarine areas (e.g. Yangtze River Estuary [21]; Hudson River Estuary [58]; South Florida estuaries [16]) further confirm that seasonal variations in water temperature trends and rates of change can vary greatly as a function of the time-window and sites/areas analysed. Overall, these differences highlight the need for joint analysis of water temperature series from different regions and time periods to better understand the spatial and temporal effects of climate change in coastal and transitional systems.

4.4. Interannual Patterns

The alternation of warming and cooling periods identified in our systems, with inflexion points in 2003-2006 and 2013-2015, matched the large scale subpolar North Atlantic SST trend reversal from warming (1994-2004) to cooling (2005-2015) [59]. Earlier and longer (from 1965 to 2004) water temperature data from 0 to 100 m depth in the Bay of Biscay also showed an overall, though non steady, warming during 1965-2004, when an initial cooling period occurred until the early 1970s [55]. A six-decade (1947-2007) SST analysis from a coastal location close to our study area identified cycles of approximately 8, 11, and 18 years, linked to climate cycles [52]. These alternating cooling-warming trends in the Bay of Biscay were more patent during the first decade of the 2000s (2002-2014) under atmospheric conditions favourable for plume enhancement [60]. Our results support the relationship between cooling and increased river discharge, but overall, the cooling-warming cycles also coincide with broader spatial scale cycles that affect the Bay of Biscay and the temperate Atlantic region. Regarding extreme temperatures, in the cooling-warming cycles in neritic waters we identified the warm summer values of 2020 and 2023 and the cold winter values of 2015, in addition to the warm summer values of 2003 and 2006 and the cold winter values of 2005. This timing for the maximum and minimum values was also observed for the SST series from a more eastward site on the Basque coast [52].

4.5. Seasonal Patterns

The largest difference in seasonal patterns was observed between neritic subsurface waters of the estuary of Bilbao and surface waters of the inner estuary of Bilbao, and it seemed to be mostly due to the dominant oceanic and continental influence respectively. This difference involved larger temperature variations in the shallower waters located landwards from the mouth, a pattern observed in other estuaries too e.g., [61]. Annual water temperature minima were lowest and most skewed towards January in above-halocline waters of the estuary of Bilbao and transitional shallow waters of the estuary of Urdaibai. In contrast, in neritic and below-halocline waters of the estuary of Bilbao, minima were milder and skewed towards March, as in the adjacent coastal/shelf waters of the Cantabrian Sea [56]. However, no such differences were observed for the annual maximum values and their timing. A similar pattern has been reported in other European marine areas, including the Baltic Sea, where summer SST peaks in August, but winter minima occur in February in shallow waters and in March in deeper areas [62]. Seasonal patterns also showed interannual differences. Throughout the study period, the seasonal warming phase shortened due to a delay in winter minima, an advance in summer maxima, and an earlier start of an extended warm period. This pattern agrees with observations of trends made at a broader spatial scale across the North Atlantic in recent decades. These data have shown that despite the timing of the spring warming transition has remained constant from 1998 to 2022, there are small areas that have experienced either delayed or advanced transitions, such as the earlier spring warming transition registered in the Bay of Biscay [63]. Similarly, SST trends in the North American northeastern continental shelf showed earlier summer onset and longer durations between 1982 and 2014 [13].

4.6. Relationships with Local Climate and Hydrographic Factors

In our two systems, at the interannual scale air temperature was the main driver of water temperature variability. In addition, warming rates in transitional waters (0.25-0.24 °C decade^-1) were comparable to trends in offshore air temperatures from the Copernicus database. However, overall, the waters (neritic and transitional) of both estuaries warmed more slowly than air temperature. In contrast, Chesapeake Bay showed faster warming of some estuarine waters than of air temperature, particularly in the main stem [50], while the Hudson estuary showed no significant difference between air and water temperature increases [58]. However, the seasonal patterns of air and water temperature from our two systems showed differences, which suggests that local factors modulate estuarine water warming responses. Unlike air temperature, spring (April-June) water temperature showed strong warming trends. A similar pattern was reported for the Hudson River estuary, where water warming maximum occurred in April-August but air temperature one in December and February-March [58].

Despite differences in river flow between estuaries, the overall trend of increase in freshwater discharge had a cooling effect on water temperature in both systems, which was reflected in a decline until 2013–2015, followed by a rise as river discharges first increased and then decreased. A similar cooling influence of freshwater discharge has been reported for the Bay of Biscay during 1982-2014, where river-influenced areas cooled while offshore waters warmed [60]. In the present study, the highest cooling effect was observed in above-halocline water masses of the inner estuary of Bilbao. However, responses of water temperature to river discharge and air temperature vary across estuaries. For example, in the Hudson River estuary increased freshwater input did not counteract water warming related to the increase in air temperature [58], while in the San Francisco Bay-Delta system, the combined effects of atmospheric and riverine influences differed spatially [34]. In our study, water temperature correlated strongly with river flow in spring and autumn, but no relationship was found in summer. Similarly, in the upper San Francisco Estuary, the negative relationships between water temperature and river inflow were most common, although positive relationships were found, in winter [35].

Significant positive correlations between water temperature and upwelling index were found across all water-masses in the estuaries of Bilbao and Urdaibai, but only in spring. Spatial differences in SST show that regions in the SW Bay of Biscay experiencing intense upwelling show lower SST values [53] and coastal upwelling mitigates warming in nearshore waters in comparison to adjacent ocean waters along the western Iberian Peninsula [12,64]. Upwelling and downwelling events also shape the thermohaline properties of the northwestern Iberian estuaries [32,65]. In our study, downwelling was the dominant process, and it showed a weakening trend (upwelling index increased), that was most intense in April; accordingly, as stated above, we only found a significant positive correlation between upwelling index and temperature in spring. This interannual pattern is opposite to what we would expect if upwelling/downwelling had a dominant effect on SST e.g., [66,67,68]. Coastal upwelling and downwelling processes are driven by local atmospheric forces, and on the Basque coast, the seasonal pattern shows that spring downwelling relaxation is related to reduced wind stress and changes in wind direction, mainly of the prevailing northwesterly winds [40]. This suggest that water temperature trends in our estuaries may be influenced more by changes in more local winds than by shelf water advective processes. Surface water warming can be influenced by both atmospheric surface heat flux and oceanic horizontal thermal advection, with winds playing a significant role in air-water heat flux [69] (Pei et al., 2017). However, studies report positive, negative and uncertain relationships between SST and surface wind speed in coastal waters e.g., [69,70,71,72]. Therefore, the effect of winds on water temperature in our study area remains unclear, highlighting the need to include wind data in future analyses of water temperature variability.

4.7. Relationships with Climate Teleconnections

The AMO index showed the most consistent relationship with interannual temperature variations across all water masses and seasons in both estuaries. According to [73], the AMO has been the dominant teleconnection pattern in the Bay of Biscay, i.e. it has explained the highest percentage of interannual SST variability, over a recent period of 150 years, while the influence of the NAO and EA patterns has been smaller and limited to specific years or shorter periods. In our study, the AMO was always positively correlated with water temperature in spring and autumn, the two seasons with the highest warming rates, while the NAO was correlated inversely only with autumn temperature. This contrasts with observations from coastal areas of the Bay of Biscay located northwards from our study sites and influenced by river plumes, where SST was positively correlated with the NAO index between 1982 and 2014 [60]. However, it agrees with the sequential negative correlation of the NAO index with air temperature and the positive correlation between air temperature and water temperature obtained in our study area from 1997 to 2005 [74]. Such change of sign in the relationship between NAO and SST across the Bay of Biscay might be due to the transitional position of this marine region between the northern and southern European centres, which exhibit opposite temperature responses to the NAO index [75]. On the North American northeastern continental shelf, SST change trends also show regionally variable links to the NAO index [13]. In our sites the EA and OA patterns showed weaker relationships with water temperature overall and seasonally. However, the positive relationships obtained in our study area are consistent with findings from nearby areas, where the positive phase of the EA was linked to increased incidence, duration, and intensity of marine heatwaves in a central coastal area of the Cantabrian Sea (SW Bay of Biscay) between 1998 and 2019 [76].

5. Conclusions

Our results reveal common patterns in water temperature variability across both estuaries, which agree with those observed over wider areas in the Cantabrian Sea, the Bay of Biscay, and the North Atlantic. However, differences between and within estuaries highlight the sensitivity of estuaries to climate warming in a global context, as well as, to the specific features that characterize the Basque coast. Based on temperature ranges and seasonal and interannual variability, two main thermal behaviours were identified. One is associated with water masses from the estuary mouths (neritic) and below the halocline in the stratified estuary of Bilbao, which was characterized by lower temperature variability, later annual minima and an interannual delay in the alternation of warmest and coldest events influenced by shelf waters thermal properties and dynamics. The other one is linked to water masses in the intermediate-inner reaches of the estuary of Urdaibai and surface waters of the estuary of Bilbao, which showed a higher variability and a stronger influence of local land-atmosphere conditions.

Temperature trends showed similar maximum warming rates in transitional water masses of both estuaries, but in different water masses, because these maxima were observed in below-halocline higher salinity waters (euhaline) in the estuary of Bilbao and in shallow lower salinity waters (polyhaline) in the estuary of Urdaibai. Warming rates were higher in the neritic waters of the estuary of Urdaibai than in the estuary of Bilbao, likely due to between-estuary differences in depth and geographic location. Overall, decreasing trends in above-halocline inner waters of the estuary of Bilbao were linked to river discharge and over-land climate factors.

Highest warming rates occurred in spring, particularly in May, except for the above-halocline waters of the estuary of Bilbao, likely due to the higher influence of river discharge in these waters. The strongest cooling in August sharpened summer negative trends in the neritic waters of both estuaries. A similar pattern was observed in winter in the above-halocline waters of the estuary of Bilbao. These cooling trends were attributed to external factors of marine and continental origin, respectively. A delay in the coldest annual temperatures and an advance in the warmest ones, which shortened the annual warming period, were the main water temperature phenological changes. Additionally, a sharp increase in May extended the warm period, reflecting a transition from cooler to warmer conditions over the study period.

Water temperature was primarily driven by air temperature but modulated by river flow, which had a cooling impact, particularly in spring and autumn. The relationship between upwelling index and temperature was complex, suggesting the influence of additional factors.

Teleconnection patterns showed the AMO index had the most consistent relationship with temperature, especially in waters with higher influence of marine waters (neritic waters and below-halocline waters of the estuary of Bilbao), while the NAO was mainly linked to autumn temperatures. The EA and AO had weaker associations with temperature changes.