Most eastern boundary upwelling regions represent thermal refugia in the age of climate change

Introduction

The four major Eastern Boundary Upwelling Systems (EBUS) - Humboldt, Benguela, California, and Iberian/Canary - are characterized by a narrow coastal band (~50–150 km) of surface temperatures that are cooler than temperatures offshore (Bakun and Parrish, 1982; Chavez and Messié, 2009). This is caused by wind-driven coastal upwelling that brings deep, cold water to the surface at the coast. Cool, nutrient-rich upwelled waters fuel large phytoplankton blooms that sustain the rich marine ecosystems of EBUS (Kämpf and Chapman, 2016), including 20% of world fisheries catch (Pauly and Christensen, 1995; Bindoff et al., 2019) and high biodiversity of upper trophic level species (Chavez and Messié, 2009; Watermeyer et al., 2020]. Furthermore, the cool coastal sea surface temperature (SST) footprint is important to adjacent terrestrial ecosystems in driving low clouds (fog), low temperature, and high humidity (Johnstone and Dawson, 2010) that are favorable to forested and agricultural areas (Black et al., 2014). Overall, EBUS provide services to ~80 million people (García-Reyes et al., 2015). On a global scale, the cool EBUS SST footprint affects climate through air-sea interactions (Large and Danabasoglu, 2006; Curchitser et al., 2011; Bakun et al., 2015), as well as CO₂ outgassing (as upwelled water is high in CO₂ content) and capture (through enhanced primary productivity) (Mackas et al., 2006; Ikawa et al., 2013; Brady et al., 2019).

Changes in the cold footprint of upwelling have not been studied extensively as have changes in temperatures along EBUS, although they are clearly related. For example, Santora et al. (2020) quantified the upwelling footprint as the coastal area below a temperature threshold associated with upwelled water, calling it the habitat compression index (HCI). Impacts of warm conditions along EBUS are well documented, for example, during El Niño events along the Pacific EBUS (Chavez et al., 2002; Escribano et al., 2004; Jacox et al., 2016; Diaz-Astudillo et al., 2022), or a marine heatwave (MHW) in the Northeast Pacific in 2014–2016 (Gentemann et al., 2017; Santora et al., 2020); both types of events caused large negative consequences for coastal socio-ecological systems. In the California, variability in the HCI has been associated with changes in zooplankton and forage indices (Schroeder et al., 2022), as well as the risk of humpback whale entanglement in fishing gear due to reduction of cold favorable habitat, which leads to redistribution and overlap of whale prey into crab fishing grounds (Santora et al., 2020). The HCI has demonstrated that changes in upwelling footprint is a good indicator of suitable habitat.

Anthropogenic global climate change is impacting EBUS similarly to other oceanic regions by increasing global ocean temperatures (Cooley et al., 2022; Bograd et al., 2023), stratification (Bakun et al., 2015; Li et al., 2020), and frequency of extreme events such as MHWs (Gentemann et al., 2017; Varela et al., 2021; Izquierdo et al., 2022; Thoral et al., 2022). On the other hand, upwelling-favorable winds are predicted to strengthen in poleward regions of EBUS (Sydeman et al., 2014; Bakun et al., 2015; Rykaczewski et al., 2015; Wang et al., 2015; Abrahams et al., 2021; Bograd et al., 2023), potentially counteracting increasing temperature trends. However, how temperatures in coastal upwelling regions are and will respond to these counteracting drivers has been a subject of debate, with unclear or contradicting trends (Bograd et al., 2023) depending on dataset, period of study, and spatial scale (Sydeman et al., 2014; Bakun et al., 2015; García-Reyes et al., 2015). Few studies have examined SST trends at the spatial scales at which coastal upwelling occurs. Varela et al. (2018) and Seabra et al. (2019) showed that decadal trends in SST along EBUS have changed little, or even cooled, at nearshore locations (within ~50 km) due to the cooling effects of upwelling. For these reasons, EBUS have been hypothesized to be “thermal refugia” from global warming (Bakun et al., 2015; Hu and Guillemin, 2016; Lourenço et al., 2016; Barceló et al., 2018; Varela et al., 2018; Lourenço et al., 2020), but observed regional and seasonal changes in the spatial extent of the surface expression (i.e., “footprint”) of coastal upwelling has yet to be investigated.

A general question of interest is, therefore: as global ocean temperatures rise and marine heatwaves become more frequent (Bindoff et al., 2019), are the upwelling-derived cool-water habitats in EBUS coastal areas changing in a manner that may be detrimental to marine life and coastal societies? This question has yet to be asked across all EBUS. In particular, trends in the upwelling footprint would depend on the relative importance and rate of change of different mechanisms driving temperature change in different regions (see review in Bograd et al., 2023 and in the IPCC’s Cooley et al., 2022). In this study, we address this question by investigating changes in spatial extent of cool SST footprints in EBUS using high-resolution observations collected over the last two decades. We test the null hypothesis of no significant change (expansion or contraction) in the cool-water SST footprint over time, consistent with the concept that EBUS coastal habitat represents thermal refugia to global warming due to the effects of upwelling.

Data & methods

In this study we used an area of 300 km from shore (approximately the extent of national exclusive economic zones (EEZ); see Figures 1, 2) along the EBUS to fully capture variability in the SST associated with upwelling coastal process. This occurs mostly in a narrow band of ~30 km from shore, corresponding to the Rossby radius of deformation (Allen, 1973; Carr and Kearns, 2003), but the cooler SST upwelled can extend to 150 km or more due to offshore transport, especially downstream of upwelling centers like promontories or capes (Chavez and Messié, 2009; Hutchings et al., 2009). These 300-km subregions are embedded in the equatorward eastern boundary currents that occur year-round and can extend hundreds of kilometers offshore (Thiel et al., 2007; Checkley and Barth, 2009; Mason et al., 2011). We selected three subregions within each EBUS (Figure 1, Table 1) based on subregions identified in the literature as dynamically coherent (Tarazona and Arntz, 2001; Thiel et al., 2007; Arístegui et al., 2009; Checkley and Barth, 2009; Hutchings et al., 2009). The subregions were identified by their distinctive seasonality and magnitude in upwelling-favorable winds and sea temperature, with equatorial subregions exhibiting larger temperatures and year-round winds, and poleward subregions with lower temperatures and stronger seasonality in winds. Central subregions in most EBUS exhibit a strong center of seasonal upwelling.

Figure 1 Mean SST (2002-2022, MURSST) in the four major EBUS ((A) California, (B) Iberian/Canary, (C) Humboldt, (D) Benguela), within an area extending 300 km from the shore. Selected subregions within each EBUS are delineated. The black line delimits the threshold SST_T for each subregion that identifies the cool SST footprint. To the right of each map are the trends in annual SST footprint from the GLS model (diamond) and the GLM (circle) for each subregion represented as the percent of the total area per year. Statistically significant trends, based on the two-tail Monte Carlo test, are colored in red and denoted with an asterisk (See Table 1). 95% confidence intervals are shown from the GLS model and indicated by black lines.

Figure 2 (A) Map of the poleward Benguela subregion illustrating the method of study area delineation. Subregion broad boundaries are shown as a black box, the total area of study is delimited by the ‘blue’ area from the coast to 300 km offshore (this area represents 100% of the possible area of the SST footprint), and the 30-km ‘red’ coastal band is where upwelling occurs and where the SST threshold (SST_T) was calculated as the mean SST within it. (B) Cool SST footprint for September 2016, with color representing monthly SST and the green contour line indicating the subregion SST threshold (SST_T = 15.6°C); the portion of the total area with temperature below the SST_T (green line) is defined as the SST footprint for this month.

Table 1 Latitudinal range, SST threshold (SST_T) for each subregion (defined mean SST of the coastal band within 30 km from shore for each EBUS), and for comparison, the entire subregion mean SST (from shore to 300 km away from the coast).

We used Multi-scale Ultra-high Resolution (MUR) SST V4.1 data (JPL Mur MEaSUREs Project, 2015; Chin et al., 2017) for a 20-year period (June 1, 2002 to May 31, 2022). MUR SST data are gap-free at daily temporal scales and have a spatial resolution of 1x1 km, which is ideal for this analysis; the fine spatial resolution can most accurately capture the SST footprint extension and its changes. We aggregated the SST data monthly to increase the signal to noise ratio, as we were interested in interannual and longer time scale changes. Another advantage of using MUR SST, despite its relatively short temporal reach, was that it is available worldwide, allowing a comparison across EBUS in recent times. Subregions of SST were then delimited (masked) by the latitudinal range shown in Table S1, the coast, and a longitudinal line 300 km away from the shore. This masked area (colored in Figures 1, 2) constituted the 100% of the area for each subregion.

Inspired by the HCI Index definition from Santora et al. (2020), we defined an indicator of the EBUS cool SST footprint: an area with SST below a given threshold proxying a coastal cool habitat associated with upwelling and mixing of cool and nutrient-rich waters, and therefore biological productivity (Schroeder et al., 2022). This threshold (SST_T) was calculated as the long-term mean SST value of the 30-km coastal band within each subregion (Table 1, red band in Figure 2A) where coastal upwelling occurs. The SST_T used in this research might differ from one more specifically derived based on local upwelling measurements that factor in wind stress or nutrient content, but this one-size-fits-all definition is an objective way to compare SST footprint variability for all subregions, and it was necessary given the heterogeneity in coastal upwelling among and within EBUS. Coastal upwelling shows heterogeneity in phenology, magnitude, and sign of winds and SST through the year (Figure 3) (Bakun and Nelson, 1991; Carr and Kearns, 2003; Escribano and Schneider, 2007; Bograd et al., 2009). However, we used the same subregional threshold SST_T value for the entire subregion, as it is representative of the local upwelling along the subregion coast.

Figure 3 Seasonal cycle of upwelling-favorable meridional wind speed (red diamonds: ERA5 data), SST from the coastal 30-km band (blue dots: MURSST), and cool SST footprint percentage (black triangles; MURSST). Values were based on monthly means for the period 2002-2022, for (A) California, (B) Iberian/Canary, (C) Humboldt, and (D) Benguela EBUS.

The monthly cool SST footprint extension was then calculated as the number of 1x1-km pixels in a given subregion that were below its threshold SST_T (see Figure 2B), and was expressed as a percent of the subregion area (masked) to standardize it. Because the threshold SST_T was based on the SST of the cool 30-km coastal band, the total subregion area was seldom below that value; however, winter temperatures, not necessarily associated with upwelling but not distinguishable from them, could be below the threshold SST_T. Finally, we calculated seasonal and annual averages of the SST footprint. The annual average spanned winter to fall for both hemispheres, covering an entire upwelling season. Winter was defined as December to February for the northern hemisphere and June to August for the southern hemisphere, and so on.

We tested our hypothesis through two methods. First, we used a generalized least square model (GLS) to calculate linear trends for the subregions within each EBUS. To build the GLS, the subregion’s annual SST footprint time series acts as the response variable, year as the time covariate, and EBUS and subregions as categorical values (See supplementary material). Following the methodology laid out by Zuur et al. (2009), the GLS included a weight parameter to help control for heteroskedasticity in the residuals. For this model, a constant variance function was used in which each EBUS was allowed to vary independently. The GLS also included a correlation term that corrected for the non-independence of the residuals, necessary given the relatively short span of the time series. A first order autoregressive structure was fit to the residuals of the model with year as the continuous covariate, and the significance is represented in the results by the confidence intervals. This method allowed us to estimate the linear trends with larger statistical power by increasing the degrees of freedom when combining data together from all subregions. The second approach used generalized linear models (GLM) to individually calculate the trend for each EBUS/subregion time series, and therefore had fewer degrees of freedom given by the length of each time series. These GLMs had a similar formula as the GLS, with percent of the area below SST_T as the response variable and year as a continuous covariate. For the GLM trends a two-tail Monte-Carlo test was used to calculate p-values from a randomized set of similarly distributed data to avoid overestimating statistical significance (p < 0.025 for the two-tail test) due to autocorrelation and length of time series. Both model methods provide similar information, but have different statistical capabilities and we consider them complementary to address the length of the time series and the large interannual variability in some regions. Furthermore, each subregion is large enough to be considered independent from each other (given the smaller Rossby radius of deformation (Huyer, 1983) and the offshore extent of upwelling compared to the distance to the next subregion), therefore, we are not concerned with inflation of p-values by the number of linear regressions calculated.

To investigate changes in the timing of the SST footprint, we calculated linear trends in seasonal means in each subregion’s footprint. Seasonal means were chosen over monthly to avoid over-calculating trends on data that did not necessarily show large changes month to month. In addition, we calculated linear trends of seasonal SST within each EBUS to further describe the spatial distribution of cool SST footprint changes. We used the GLM and Monte-Carlo tests to calculate trends and p-values for these time series. The Python and R code used to generate the indices and analyze the data can be found in a Github repository [https://github.com/farallon-institute/Garcia-Reyes_etal_2023_EBUSFootprint].

Additionally, to illustrate the seasonality of upwelling-favorable wind and SST, we used the monthly seasonality of SST in the 30-km coastal bands and the meridional wind (V) over each EBUS subregion (Figure 3). Meridional wind speed (m/s) was extracted from the ERA5 dataset (Hersbach et al., 2020) for the period 2002–2022. Wind data were provided at 10 m over the sea surface with an hourly temporal resolution and a 31-km spatial resolution. SST was extracted from MUR SST for the same period.

Results

The extent of cool SST footprints in all EBUS subregions show differences in timing and magnitude through their seasonal cycles (Figure 3), with magnitudes varying from 0% to 40–100% (area below the threshold SST_T). The largest cool footprints occur between February and April in the northern hemisphere, and between July and September in the southern hemisphere. The meridional wind and SST seasonalities shown in Figure 3 are consistent with those reported in the literature (Shannon, 2001; Carr and Kearns, 2003; Chavez and Messié, 2009; Montecino and Lange, 2009), but note that the peaks of the SST footprint do not occur simultaneously with the peaks in meridional wind speed (upwelling-favorable) in most subregions (Figure 3). The difference in peak timing indicates that the strong seasonality in temperature due to the solar cycle dominates the timing of the coastal SST annual cycle, rather than upwelling-induced cooling, independently of the cooler SST magnitudes at the coast due to upwelling. In winter, temperatures are below the mean SST threshold, and the SST footprint is large. This is not necessarily due to upwelling (although in equatorial subregions upwelling occurs year-round), but SST values are below those typical of the upwelling season anyway, so they still represent a cool habitat.

The seasonality of the SST footprint unfolds, in general, as follows. During the peak of upwelling-favorable winds, nearshore waters are cooler than offshore. As summer unfolds, the SST footprint contracts, even if upwelling-favorable winds are still present, as deep cool water no longer reaches the surface due to the increasing stratification of the water column caused by the seasonal warming of upper layers of the ocean (Jacox and Edwards, 2011). The poleward Benguela subregion shows the narrowest SST footprint due to the warm influence of the portion of the Agulhas Current that enters the South Atlantic as mesoscale eddies (Hutchings et al., 2009), which leads to a strong cross-shore SST gradient year-round (Figure 1D). In any season, SST footprints can contract for other reasons, such as the occurrence of a marine heatwave, as occurred in 2014–2016 in the California EBUS (Figure S1) (Gentemann et al., 2017).

Linear trends for the annual SST footprint show similar results for the GLS and GLM (Table 2, Figure 1), with no significant trends for 10 of 12 EBUS subregions (p < 0.025 in the two-tail Monte Carlo test in the GLM or p < 0.05 in the GLS); the exceptions are the equatorward Benguela and California subregions, which show a contraction of their annual SST footprints. The equatorward California subregion shows a significant negative large trend in the GLM, but not in the GLS due to the autocorrelation and large interannual variability of the time series (Figure S1). The central California subregion shows a significant negative trend in the GLM, but the Monte Carlo test removes the significance and the trend is small and largely not significant in the GLS (Figure 1) due to autocorrelation and large interannual variability. It is worth noting that despite the lack of statistical significance (p > 0.025 on the Monte Carlo test), trends in other subregions show negative trends of significant magnitude, like central and poleward California in the GLM, and a positive trend in central Iberian/Canary in both models.

Table 2 Linear trends of annual SST footprint area in %/year for generalized least squared (GLS) model and generalized linear model (GLM).

Seasonal linear trends in SST footprints for each EBUS subregion from the GLM are shown in Figure 4 and Table S1, while the time series are plotted in Figure S2. Trends vary in magnitude, and in some subregions in sign, across seasons. Humboldt and Iberian/Canary show no significant trends in any season (based on the Monte Carlo test), although they show not significant trend in all seasons, especially a positive trend in the fall in central Iberian/Canary and negatives ones in winter and summer in central Humboldt. The Benguela EBUS shows negative trends in winter (significant only in the equatorward subregion), and is more neutral in other seasons for the poleward and central subregions, although the equatorward subregion also shows a negative trend in spring. The California EBUS exhibits the largest negative trends among EBUS, in both its central and equatorward subregions. In the fall, all subregions show contracting (negative) trends in the Monte Carlo test, with the equatorward subregion having the largest trend, and the poleward the smallest. Central and equatorward California regions also show significant negative (contracting) trends in the summer, again, with the equatorward subregion exhibiting the largest negative trend, and the largest interannual variability in all seasons. It is worth noting that the central California spring season has a very small positive trend, the only one in this EBUS for any season.

Figure 4 Linear trends of seasonal means of SST footprint (%/year) for each EBUS subregion with 95% confidence intervals from the GLM. Colors indicate seasons; winter is defined as December to February for the northern hemisphere and June to August for the southern hemisphere. Significant trends (p < 0.025 based on two-tail Monte Carlo test) are represented with an asterisk, and confidence intervals are based on the linear regression model. (A) California, (B) Iberian/Canary, (C) Humboldt, and (D) Benguela EBUS.

Seasonal SST trends show narrow coastal regions of decreased temperatures in the central and poleward Benguela, Humboldt, and central Canary subregions, and warming trends in equatorward Benguela and California EBUS (Figures S3–S6). While the trends in most areas were not statistically significant (p>0.05, significance directly from the linear trend calculation at each grid point), the Iberian-Canary and Humboldt EBUS had a mix of small negative and positive trends offshore of the coastal band, while the Benguela and California EBUS had more ubiquitous warming trends, except during the spring season, which were negative. Any negative trend in the EBUS areas is opposed to the warming global SST trends (see Table S2 for global and hemispheric SST trend magnitudes), however, EBUS warming trends, when significant, can be larger than global or hemispherical trends by up to an order of magnitude. For example, the central and equatorward California subregions in the fall and the equatorial Benguela subregions in winter show large patches of significant warming up to 1°C/decade (from 2002 to 2022), while the global SST trend is 0.14°C/decade and the trend for the northern hemisphere is 0.21°C/decade (See Table S2) change into (NOAA Climate at a Glance). Note, however, that SST trends in Figures S3–S6 are representative of the change in temperatures in the region, but their magnitudes are based on pixels of 1 km x 1 km, as per the resolution of the data; smaller trends (and with different statistical significance) would be obtained if spatial averages were considered (as done for global or California EBUS scale analysis).

Discussion

In this paper, we quantify the changes in the coastal cool SST footprint that characterize EBUS. We test if 20-year linear trends in the extension of these footprints are contracting at annual and seasonal scales across EBUS (Table 3). We found that, taken into account the autocorrelation and variability of the time series, 10 out of 12 EBUS subregions support the hypothesis that the extension of the footprint, on annual scales, has not contracted in the last 20 years, despite the continued warming of the global oceans (Fox-Kemper, 2021). This result is in agreement with the idea that EBUS could be considered thermal refugia, and consistent with results from Seabra et al. (2019) showing that coastal areas in EBUS are responding differently to global warming, mostly warming at lower rates than offshore areas. Two subregions, equatorward California and Benguela, however, exhibit a contraction of their annual SST footprint in the period of study, with the former being the largest trend across all subregions. The large interannual variability in the equatorward California, plus the autocorrelation in the time series, removes the statistical significance in the trend in the GLS model, however, as we discuss later, diverging seasonal trends play an important role. Although not significant, the equatorward Humboldt subregion exhibits a contracting trend (negative), but on the contrary, equatorward Iberian/Canary shows an expanding trend. Interestingly, despite the lack of statistical significance, there is a pattern in the trends in which the equatorward regions exhibit the largest contraction trends (negative) and the poleward the smallest or zero (Figure 1), except in the Iberian/Canary EBUS.

Table 3 Were used in. Summary of locations and timing of statistically significant trends in the annual and seasonal cool SST footprint (downward arrow: contraction; upward arrow: expansion).

This latitudinal pattern appears to be consistent with the weakening of upwelling, stronger water column stratification, and warmer temperatures in equatorward regions of EBUS predicted under climate change (Bakun et al., 2015; García-Reyes et al., 2015; Rykaczewski et al., 2015; Wang et al., 2015). The reported strengthening of upwelling-favorable winds in recent decades (Sydeman et al., 2014; Varela et al., 2015; Bindoff et al., 2019; Bograd et al., 2023) and predicted under climate change (Bakun, 1990; Rykaczewski et al., 2015; Wang et al., 2015; Bindoff et al., 2019) is not reflected in an expansion of cool SST footprints, although we observe cooling of coastal waters during the study period across EBUS, except California. It is possible, however, that the lack of expansion is caused by the increasing offshore temperatures (Seabra et al., 2019; Fox-Kemper, 2021), pressing on EBUS coastal footprints.

The lack of contraction of annual SST footprints in most EBUS, however, is encouraging, and it is supported from previous studies showing that in coastal areas of the EBUS, upwelling appears to mitigate or counteract the global ocean warming trend observed elsewhere (Santos et al., 2012; Varela et al., 2018; Seabra et al., 2019; Abrahams et al., 2021; Izquierdo et al., 2022). This suggests that, for the period 2002–2022, local processes like coastal upwelling, winter mixing, and longwave cooling have been strong or stable enough to counteract global warming trends and/or regional climate processes that drive low-frequency variability in EBUS conditions (Bonino et al., 2019; Marin et al., 2021). Notably, for the period of this study, the Humboldt EBUS does not show the reported cooling (corresponding to contraction of the SST footprint) trends observed in other datasets with longer time series (Varela et al., 2018; Seabra et al., 2019; Abrahams et al., 2021), while California shows significant contraction and warming at seasonal scales (see below).

Given the short length of the time series in climate time-scales, the observed trends might not be a reflection of climate change only, but also of basin-scale decadal climate oscillations. This is particularly the case of the Iberian/Canary poleward subregion, which exhibits opposite trends (positive, expansion) in annual SST footprint to other EBUS, although not statistically significant, and the California EBUS where the largest variability is observed (Figures S1, S2). However, the Iberian/Canary trends are consistent with previous results showing trends in upwelling-associated indicators divergent from other EBUS (Sydeman et al., 2014; Varela et al., 2015; Abrahams et al., 2021; Bograd et al., 2023). These differences among EBUS trends (in sign and magnitude) are addressed in a modeling study by Bonino et al. (2019) which shows that different regional mechanisms drive SST variability in each basin: California and Humboldt systems respond largely and synchronously to El Niño variability, Benguela responds to long-term anthropogenic climate change, and Iberian/Canary responds to North Atlantic multi-decadal variability. While results presented here are mostly consistent with Bonino et al. (2019), (seen for example in the large interannual variability in the California EBUS, associated with El Niño in 2016 and the preceding large marine heatwave in 2014–2015 that impacted the trends and their significance), in the poleward Humboldt cooling has been associated with strengthening upwelling-favorable winds and regional cooling temperatures (Varela et al., 2015; Abrahams et al., 2021), rather than the large interannual variability associated with El Niño observed in the equatorial subregion and in the California EBUS. As global ocean warming is expected to continue (Fox-Kemper, 2021), it is likely, however, that warming trends might become significant in all subregions, despite upwelling-favorable wind trends or regional multi-decadal oscillations, due to increased offshore SST, increased ocean temperature at depth, and increased stratification that limits the efficacy of upwelling to bring deeper, cooler waters to the surface (Bakun et al., 2015; Li et al., 2020).

Since coastal upwelling is mostly a seasonal process in EBUS, we further investigate changes in SST footprints by calculating linear trends of seasonal means of SST footprint for each EBUS subregion (Figure 4, Table S1) and of spatially explicit seasonal SST along EBUS (Figures S3–S6). Corroborating the results for the annual SST footprint analyses, no statistically significant trends in seasonal SST footprint were found in the Humboldt and Iberian/Canary subregions in any season. However, in the Humboldt EBUS winter footprint trends are negative and larger in magnitude than in the other seasons, particularly in the equatorward subregion, which correspond to larger warming in this region in Figure S5. The lack of significance in the Monte Carlo test, despite the large magnitude of the winter footprint trend, is likely caused by the large interannual variability and magnitude of the SST footprint in this season. Winter also exhibits patches of significant warming in all subregions. In the Iberian/Canary EBUS, trends in all seasons in the poleward and equatorward subregions are close to zero, while in the central subregion, the fall shows a small expanding trend, also represented in a coastal band of significant cooling SST (Figure S4). Winter also shows an expanding trend that is larger but not significant, which corresponds to more widespread but not significant cooling. In the Benguela, the central and poleward subregions show no significant trends in all seasons, although winter trends are negative. Figure S6 also shows how for these subregions, a coastal cooling occurs in spring in both regions, and in the poleward regions also in summer and fall, while winter show more warming patches offshore. Therefore, the SST footprints in the Humboldt and Iberian/Canary EBUS and the central and poleward Benguela subregions are not significantly contracting at seasonal scales during the period of study, supporting the hypothesis, and consistent with the idea that these EBUS coastal areas can act as thermal refugia, as offshore areas exhibit warming.

On the other hand, the equatorial Benguela shows increasing trends in SST in winter and spring (June–November, Figure 4, Table S1), consistent with the contracting trends in the annual footprint (Figure S6). This warming trend is in agreement with the large warming (largest among EBUS) reported by Seabra et al. (2019) on a longer time scale in this region, and with a study (Bonino et al., 2019) showing that the driver of low-frequency variability in the Benguela EBUS is related to anthropogenic global change. Moreover, in the other Benguela subregions, non-significant negative footprint and positive SST trends occur in the same season, which are consistent with the findings of Bonino et al. (2019) as well. It is worth noting that the warming SST trends occur during seasons when the cool SST footprint is largest (austral winter, Figure 3), thus, the contracting trend in footprint represents a significant loss of cool habitat area in this EBUS (Table 2).

Similarly, the California EBUS exhibits statistically significant decreasing SST footprint trends and SST warming trends in all its subregions during the fall, as well as in the summer in the equatorward and central subregions (Figures 4, S3, Table S1). This appears to contradict the lack of change in the annual SST footprint in the central and poleward subregions in the GLM and in all regions in the GLS, but this is most likely due to the fact that the SST footprint is usually fully contracted during the warm seasons (Figures 3, S2), and therefore does not have a large influence on the expanded footprint during winter and spring. In the equatorward subregion, extent of the SST footprint in summer and fall is larger than in the central subregion, likely the reason for having an impact on the annual SST footprint. Therefore, the California EBUS, despite showing no significant trend in annual SST footprint, does not qualify as a thermal refugia year-round, as it shows strong and widespread coastal warming in the fall, and to a lesser extent in the summer. The significant contracting footprint and warming SST trends in the summer/fall in the California EBUS are likely due to the large interannual variability and autocorrelation in SST associated with the strong MHW and El Niño events in this EBUS during this period (Gentemann et al., 2017) (Figures S1, S2). The large trend in the equatorward California subregion can also be due to the dominance of temperature processes and advection from the south that can overtake the weak coastal upwelling (Checkley and Barth, 2009; Di Lorenzo et al., 2005).

Finally, it is interesting to note that the contracting trends in seasonal SST footprints that are larger toward the equatorial subregions, occur in different seasons in each hemisphere—summer and fall in the northern hemisphere (California EBUS), and winter and spring in the southern hemisphere EBUS—but at the same time of year (June–November), suggesting that these trends are related to a more global response than regional/seasonal drivers.

While there has been a significant global warming trend during our study period (Table S2, Fox-Kemper, 2021), the observed trends (or lack thereof) in SST footprint and in SST could reflect decadal trends due to basin-scale oscillations (Bonino et al., 2019) such as the Pacific Decadal Oscillation (Mantua and Hare, 2002) in the California EBUS or the Atlantic Multidecadal Oscillation and the North Atlantic Oscillation in the Iberian/Canary EBUS (Pardo et al., 2011). Furthermore, large interannual variability and autocorrelation, as seen in the California EBUS due to an extreme MHW (Gentemann et al., 2017), could impact linear trend calculations on short time series. Aware of the potential issues with the limited time series, large interannual variability, and autocorrelation mentioned previously, we trade longer SST data products in this study for a shorter time period of MUR SST due to its high resolution. At this scale, we can observe change in upwelling-related SST footprints. While the use of the GLS model and Monte Carlo test aims to address these caveats, results should be carefully interpreted in the timeframe they represent. Another consideration for the interpretation of our results is the definition of subregions and threshold SST_T for the cool footprint. A latitudinal redefinition of the boundaries by a few degrees does not change the sign or statistical significance of SST footprint trends (not shown), therefore, the results are not sensitive to the chosen subregion margins. The threshold SST_T was designed to represent temperatures in coastal upwelling regions and facilitate its comparison across EBUS, but does not take into account the influence of local processes and characteristics of individual subregions. For example, the SST footprint of the equatorward California and Iberian/Canary and poleward Humboldt and Iberian/Canary subregions show the influence of upstream cold currents and latitudinal gradients in SST (Figure 1), not only the coastal signature typical of upwelling regions as in poleward California and Benguela. This study aims to evaluate the cool SST footprint as objectively as possible across EBUS and subregions, as well as across their different upwelling seasonalities, and thus it is recommended to consider local processes, and in situ data when available, in the interpretation of the results for particular locations along any EBUS.

Implications for ecosystems

Marine species that inhabit or migrate through EBUS are adapted to their variable, cold, and nutritious conditions (Kämpf and Chapman, 2016). We have shown that the Humboldt and Iberian/Canary EBUS, as well as the poleward and central Benguela and California subregions, are consistent with the idea of thermal refugia habitats for their ecosystems as their cool annual SST footprints have not contracted over the past two decades in a rapidly warming ocean [see for example Medellín-Mora et al., 2016). Species can still be impacted by changes in extent or phenology of the cool habitat, however, as they are during MHWs when warm conditions can have significant negative impacts (Sanford et al., 2019; Smale et al., 2019; Santora et al., 2020; Varela et al., 2021).

EBUS, as other oceanic systems, are subject to a wide array of anthropogenic pressures that might elicit major ecosystem modifications (Watermeyer et al., 2008; Lima et al., 2020; Rivadeneira and Nielsen, 2022). In addition, EBUS are also experiencing increased ocean acidification and hypoxic and marine heatwave events (Gruber et al., 2021) that add to the pressure of long-term overfishing and marine pollution, impacts that might be challenging to disentangle from favorable temperature conditions (see for example Yemane et al., 2014). However, EBUS are potentially significant marine thermal refugia due to the richness of the ecosystems they sustain, which is disproportionate to the world ocean’s area that they cover (<1%) (Pauly and Christensen, 1995). This remains true even in those areas showing contractions of SST footprint, as they warm at a slower rate than other coastal and offshore areas (Bakun et al., 2015; Hu and Guillemin, 2016; Varela et al., 2018; Seabra et al., 2019; Lourenço et al., 2016; Lourenço et al., 2020).

Regional differences in SST footprint trends could represent habitat changes for the species they support. Contractions of the cool SST footprint can negatively impact the abundance, health, and distribution of cold-water affinity foundational species (such as macroalgae, euphausiid crustaceans, or forage fish) as well as the higher trophic level organisms they sustain (Soto et al., 2004; Ralston et al., 2015; Pérez-Matus et al., 2017; Becker et al., 2019; Santora et al., 2020; Cárdenas-Alayza et al., 2022; Schroeder et al., 2022). SST trends in EBUS may also have direct physiological effects on cold-water species as contractions result from warm conditions (Smale et al., 2019), as well as synergistic impacts due to combined stressors (Lima et al., 2020). These conditions can ‘push’ species out of the contracting and warming regions into areas that remain cool (Poloczanska et al., 2013; Poloczanska et al., 2016). For example, warming might disrupt the coping mechanisms of forage fish beyond the stability that ecological trade-off grants in the long-term, yielding habitat loss for some species (Deutsch et al., 2015; Howard et al., 2020). Contractions of the cool SST footprint can also lead to changes in habitat use by top predators which alters trophic interactions through spatial overlap among top species and also fisheries (Santora et al., 2020; Cimino et al., 2022; Schroeder et al., 2022).

These impacts have been observed particularly during warming episodes, like the 2014–2016 MHW that drastically contracted the cool SST footprint in the California EBUS (Gentemann et al., 2017). It is worth noting that these warming events can have a large impact on the ecosystem, even when thermal refugia protect these ecosystems from larger negative impacts to some degree (Varela et al., 2021; Izquierdo et al., 2022). This is of particular importance in the California EBUS as this system exhibits large variability in SST footprint during all seasons.

As the equatorward Benguela EBUS has little interannual variability, we would expect to observe evidence of ecosystem impacts from changes in the SST footprint, however, fishing pressure in that area has largely impacted stocks and it is difficult to disentangle the role of the environmental change (Blamey et al., 2015). Furthermore, expansion in the poleward direction towards refugia areas is challenged by the Luderitz upwelling cell off the coast of Namibia that acts as an ecological barrier for fish due to its strong offshore transport (Hutchings et al., 2009). On the other hand, warming areas can attract and better support species with warm-water affinities, as seen in marine species shifts with ocean warming trends (Perry et al., 2005; Poloczanska et al., 2016; Mason et al., 2019).

An interesting result is the timing at which contraction of SST footprints happens in the equatorward California and Benguela EBUS, as it has potentially different consequences for each ecosystem. In California, little impact is seen from the annual cool SST footprint, as seasonal changes occur after the peak of upwelling and during the warm seasons (spring and summer). This might result in disparate impacts to organisms, being a thermal refugia for species that migrate into the region during the cool seasons or have life stages that are more sensitive to temperature and nutrients during the cool season (Black et al., 2011; García-Reyes et al., 2013; Black et al., 2014; Ralston et al., 2015), but a negative impact for those species that are sensitive to temperatures during the summer and fall (Doney et al., 2012; Harley et al., 2012), particularly in the equatorward subregion. On the other hand, in the equatorward Benguela subregion, contraction of the cool SST footprint occurs when the footprint is largest (winter and spring), therefore, any contraction represents a significant habitat loss for species that rely on the cool temperatures and entrainment of nutrients due to upwelling and the productivity that follows. In this subregion, species that are more resilient to summer warm conditions are bound to be less impacted than those responding to winter/spring cool conditions.

Conclusion

We used high-resolution SST data and found that the coastal cool SST footprints characterizing EBUS did not exhibit contracting annual trends in 10 of the 12 subregions analyzed, in agreement with the thermal refugia concept. However, in the Benguela and California EBUS, the equatorward subregions show contraction of the cool SST footprint (by 0.6% and 1.0% of their areas respectively). These contractions occur in the austral winter/spring for Benguela and boreal summer/fall for California, corresponding to regions with significant warming during these seasons. EBUS could be important areas of conservation for marine ecosystems, as they are areas of disproportionally high productivity for local and migrating species, especially as surrounding areas warm faster.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Author contributions

MG-R, WS, and EH conceived and designed the study. GK, WS, and MG-R conducted the statistical analyses and prepared the manuscript figures. All authors wrote and reviewed the manuscript. All authors contributed to the article and approved the submitted version.

Funding

MG-R, GK, KD, and EH were supported by The Pew Charitable Trusts Contract ID 0034562. MG-R and GK were also supported by NASA award 80NSSC20K0768. LB-R was supported by COPAS Coastal ANID FB210021.

Acknowledgments

Multi-Scale Ultra High Resolution (MUR) Sea Surface Temperature (SST) was downloaded from AWS https://registry.opendata.aws/mur from 2002–2020 [dataset accessed 08 July 2022], and from the ECMWF Reanalysis v5 (ERA5) wind data was downloaded from PO.DAAC https://doi.org/10.5067/GHGMR-4FJ04 from 2020–2022 [dataset accessed 10 August 2022]. The authors thank Chelle Gentemann for making the MUR SST dataset available in Zarr format on Amazon Web Services (AWS), Bian Hoover and Helen Killeen for the helpful discussions on statistics, and Sarah Ann Thompson for editing assistance. We thank the reviewers for their comments that helped us improve the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2023.1158472/full#supplementary-material

References