Intensified upwelling: normalized sea surface temperature trends expose climate change in coastal areas

Eastern boundary upwelling systems (EBUSs) provide valuable natural resources due to their high primary production. However, there is significant uncertainty in how climate change may affect the mechanisms that sustain these ecosystems in the future. Therefore, assessing the effects of climate change on EBUSs under the current global warming scenario is crucial for efficient ecosystem management. In 1990, Andrew Bakun suggested an increase in the upwelling intensity due to the rise of the ocean–land pressure gradient. Since there is a significant link between thermal gradients and offshore Ekman transport, we use sea level pressure (SLP) and deseasonalized sea surface temperature (SST) data from remote sensing to elucidate this hypothesis and validate it using in situ observations. SST is an indicator of coastal upwelling, and our long-term analysis of monthly and deseasonalized SST records shows that the seasonal and synoptic processes have minimal influence on the SST–upwelling intensity relationship. Upwelling within the same EBUS is not usually evenly distributed along coastlines, leading to upwelling in specific areas, referred to as upwelling centers. We compare the SST trends in the main upwelling centers of the four EBUSs with those in open ocean waters through a new index, α_UI, designed to characterize upwelling changes in the EBUS. An a-dimensional number allows us to normalize the trends independently of the upwelling system and compare all of them. Furthermore, we have complemented the SST index with sea level pressure gradient data. This new index (supported by SLP gradient trends) indicates intensification in all the EBUSs, revealing a coherent pattern within EBUSs in the same ocean (i.e., Canarian and Benguela or Californian and Humboldt upwelling systems).

1 Introduction

The world's major coastal upwelling areas exist along the eastern margins of the Pacific and Atlantic oceans. These extended coastal upwelling systems are known as eastern boundary upwelling systems (EBUSs) and sustain the most important fisheries in the world (Pauly and Christensen, 1995). The rise of cool, nutrient-rich waters supports the high primary production needed to maintain these complex ecosystems (Ekman, 1905). Their economic and ecological relevance explains the association of the EBUS with the world's major large marine ecosystems (LMEs) (Fig. 1). LMEs are a globalized approach to a management framework that defines and ranks marine regions based on their gross primary production (Sherman and Hempel, 2008), and four of the largest LMEs are embedded in the EBUS worldwide (Kämpf and Chapman, 2016).

Figure 1 Location of the EBUSs (enclosed and labeled areas) in the world, associated with eastern boundary currents.

Bakun (1990) hypothesized that the sea–land temperature gradient will increase under climate change, and, therefore, it should increase the upwelling intensity. The hypothesized increase in the temperature gradient arises from an increased atmospheric pressure gradient between the low-pressure cell that develops over the heated land mass and the high-pressure cell existing over the colder ocean. Therefore, as land warms faster than the ocean, it enhances the low-pressure cells. Thus, the increase in the pressure gradient drives more intense upwelling-favorable winds, intensifying the cold imprint in the sea surface temperature (SST) near the shore. However, previous studies that have tested Bakun's hypothesis with in situ data have found contradictory long-term trends (Barton et al., 2013; Belkin, 2009; McGregor et al., 2007; Rykaczewski et al., 2015; Sambe et al., 2016).

Sydeman et al. (2014) performed a meta-analysis of 18 trends obtained from independent studies of wind stress, both from observational and model data. Observational data were more likely to report an increase in wind stress in the four EBUSs than model data (excluding Benguela, where no observations were considered). However, the model data were less consistent between EBUSs, showing agreement with observations only in the case of the California system (both supporting intensification) and the Iberian system (with a consistent weakening of the wind stress). García-Reyes et al. (2015) also remarked that in climate change models, the upwelling-related cooling trends were difficult to reproduce due to the small spatial scale of the coastal upwelling process. Such a controversy between the model and different observations reflects the complexity of EBUS dynamics.

Upwelling within the same EBUS is not usually evenly distributed due to irregular coastlines and seafloors, resulting in more pronounced upwelling in specific upwelling centers. In these areas, the sea surface temperature (SST) drops significantly as cold subsurface water rises, leading to a stronger relationship between SST and upwelling intensity (Kämpf and Chapman, 2016). Consequently, these upwelling centers exhibit a stronger signal-to-noise ratio between SST and Ekman transport, making them ideal for studying long-term upwelling trends.

The main motivation of this study is to assess the impacts of climate change on the four EBUSs. This is pursued by using satellite-derived SST trends as a proxy for changes and a new index that normalizes the upwelling trend of the coastal upwelling with the oceanic background trends. Specifically, we have chosen points representative of each dynamical regime at each EBUS: offshore oceanic waters (OC1), non-upwelling (DW1) areas nearshore, and upwelling centers (UP1 and UP2). These points were chosen based on the consistency of the year-round upwelling centers deduced from the mean SST field and the relevance of the area as spawning and nursery emplacements for the pelagic fisheries associated with upwelling centers. Then, we compare the SST trends in the main upwelling centers of the four EBUSs with those in open ocean waters through a new index, α_UI, designed to characterize upwelling changes in the EBUS, following Bakun's (1990) hypothesis. The paper is organized as follows: Sects. 2 and 3 describe the dataset and the analysis carried out in the study. Section 4 describes the relevant results, which are discussed and contrasted with other studies, and, finally, Sect. 5 summarizes and presents the conclusions.

2 Data

We based our study on the SST blended analyses for sea surface temperature of the National Oceanic and Atmospheric Administration (NOAA) (Reynolds et al., 2007), which combines SST satellite retrievals with in situ measures from ships and buoys. This dataset has a spatial resolution of 0.25° and covers nearly 40 years (from 1982 to 2021). Following Barton et al. (2013), the 40 years used in this study allow us to estimate significant trends. Their analysis involved a comprehensive examination of both wind stress and SST. They segmented these datasets into various subsets of different lengths. Within this analysis, the trends derived from wind datasets are not significantly different from zero for all considered subset periods.

In contrast, SST trends demonstrated statistical significance with a 40-year length dataset. In addition, we incorporated reliable in situ data in the North Pacific and the North Atlantic oceans for validation. These records were obtained from the National Buoy Data Center (NBDC) and the California Cooperative Oceanic Fisheries Investigations (CalCOFI) program for the Pacific. For the Atlantic, the data were gathered from Puertos del Estado and the Radial Profunda de Canarias (RaProCan) observational program of the Spanish Institute of Oceanography in the Canary Islands (Tel et al., 2016). In situ data, limited to the Northern Hemisphere, are used to validate the satellite observations. Given the significantly greater density existing in the Pacific Ocean compared to the Atlantic Ocean, a lower error is expected in the reanalysis for the Pacific Ocean. We have also used the ENhanced ocean data assimilation and climate Prediction (EN4) dataset from the Met Office Hadley Centre (Good et al., 2013), a collection of global observations from diverse sources interpolated into a monthly product and a spatial resolution of 1°. Due to the limited sampling in the cruise data, we cannot align them with the monthly resolution of the EN4 product. To maintain statistical rigor, we choose not to employ the cruise data for the validation of EN4.

Additionally, two sea level pressure (SLP) datasets are used: the NCEP and ERA5 reanalysis data. ERA5, from the European Centre for Medium-Range Weather Forecasts, offers a high resolution of ∼0.25° and uses advanced data assimilation methods, covering the period from 1950 to the present, while the NCEP/NCAR reanalysis, from the National Centers for Environmental Prediction and the National Center for Atmospheric Research, offers a resolution of 2.5° and covers the period from 1948 to the present.

Figure 2 Average SST (°C) maps for each of the upwelling systems: CalUS (a), CanUS (b), HuUS (c), and BeUS (d). Overlaid are the locations of the moorings (red diamond) and cruise data (black dots) used for the satellite validation with their corresponding names. The cruise stations were divided into upwelling (blue dots) and open ocean areas (red dots). The representative points of UP, DW, and OC in each basin are shown as black dots. The color scales are presented at the right margin of each graph.

3 Method

3.1 Selection of representative dynamical regimes in EBUS

We have selected areas representative of the different dynamical regimes for further analysis to avoid mixing observations in dynamically different areas. Areas UP1 and UP2 are year-round upwelling centers. DW1 areas are areas of convergence where the upwelling does not dominate on an annual average, and these areas are characterized by higher SST averages than those in the upwelling centers. OC1 is representative of open ocean (>100 km offshore) areas with trends driven by global warming. Given the distinctive and unique features that characterized each EBUS, as mentioned in Sect. 1, we based the selection of the representative locations on the literature and the SST mean field (Fig. 2).

Californian Upwelling System (CalUS)

In the Californian system, the strongest wind stress takes place in spring in the southern portion of the CalUS. In summer the strongest winds occur offshore of northern California around 38° N. This wind stress pattern diminishes as we move away from this latitude in both directions (Bakun and Nelson, 1991). Among the four regions depicted by Kämpf and Chapman (2016), the strongest upwelling occurs in two well-known upwelling centers (see Fig. 2a) that we have selected for the analysis in this study: Cape Mendocino (UP1) (Abbott and Zion, 1987) and north of Point Conception (UP2) (Dugdale and Wilkerson, 1989); south of it, the wind's seasonal variability is different.

Canarian Upwelling System (CanUS)

Following Kämpf and Chapman (2016), the CanUS is divided into two distinct upwelling areas, which experience limited continuity of flows between them. This division arises due to the coastline interruption in the Strait of Gibraltar that leads to two upwelling areas: the Iberian Upwelling System and the Canarian Upwelling System. Given the lack of permanent upwelling centers in the Iberian Upwelling System, we have focused on the Canarian Upwelling System. Following Cropper et al. (2014), we selected the two upwelling centers (UP1 and UP2) on the cold SST cores (Fig. 2b) within the permanent upwelling area.

Humboldt Upwelling System (HuUS)

The HuUS is characterized by prominently strong upwelling zones along the Peruvian–Chilean coastline due to the topographical influence of headlands, as described by Figueroa and Moffat (2000) and Mesias et al. (2003). Notable regional upwelling centers encompass the continental shelf near Pisco (13.7° S), Antofagasta (21–25° S), and the Mejillones Peninsula (23° S) and extend further south to Coquimbo Bay (30° S), Valparaíso (33° S), and the Bay of Concepción (37° S). Moreover, in the Mejillones Peninsula (23° S), both observational (e.g., Marín, 2003) and modeling studies (Escribano et al., 2004) have revealed that the dynamics of coastal ecosystems in this area rely heavily on the generation of upwelling filaments. Additionally, the continental shelf near Pisco (13.7° S) emerges as an exceptionally productive and distinctive upwelling center. Hence, the two upwelling centers (Fig. 2c) selected are Pisco (UP1) and the Mejillones Peninsula (UP2).

Benguela Upwelling System (BeUS)

As described by Kämpf and Chapman (2016), within the Benguela Upwelling System, numerous upwelling centers extend along the shelf area of the Benguela Current region. These centers include Cape Frio (18.5° S), Walvis Bay (22.95° S), Lüderitz (26.45° S), Namaqualand (28.55° S), Cape Columbine (32.85° S), and Cape Town (33.95° S). Lüderitz stands out as a particularly noteworthy upwelling center within this system (Andrews and Hutchings, 1980; Lutjeharms and Meeuwis, 1987; Peard, 2007). As defined by Hutchings et al. (2009), Lüderitz represents an intensive perennial upwelling center characterized by intense winds, high turbulence, and robust offshore transport. Another significant area of interest lies in Cape Columbine, primarily due to its biological importance (Andrews and Hutchings, 1980; Bang and Andrews, 1974; Andrews and Cram, 1969). Given these distinctive features, we have selected Lüderitz (UP1) and Cape Columbine (UP2) as the upwelling centers for our analysis (Fig. 2d).

3.2 Trend analysis

With a minimum data length of over 30 years, a climate series can usually be described as a combination of multiple variabilities at different timescales. Since we are interested in the trend of the record, we removed the high-frequency variability (<1 year) by averaging the daily NOAA SST analyses data (1982–2021) into monthly means and removing the seasonal cycle. The seasonal cycle is a recursive signal throughout the entire record, and, therefore, it does not influence the trend but induces noise in the target scale. The monthly climatology was subtracted from the record to remove the seasonal cycle. After this pre-analysis, we calculated the trend with the ordinary least squares method. We evaluated the strength of these correlations using the Pearson correlation coefficient (PCC). The following qualitative classification is used throughout the paper: perfect (1), very high (>0.9), high (>0.7), and moderate (>0.5). Additionally, we employed the simple Mann–Kendall (MK) test to evaluate the statistical robustness (Kendall, 1975; Mann, 1945). The MK tests verify whether an n-length series holds a monotonic increase or decrease trend. In addition, we need to consider the instrumental error since, historically, a warm coastal bias is found in satellite records compared to in situ records (Smale and Wernberg, 2009). This bias was assessed in the Northern Hemisphere by validating the data using in situ observations. Additionally, to assess the drivers of change in upwelling intensity, we calculated the sea level pressure gradients for each EBUS. The gradients were calculated between the cores of the high- and low-pressure systems (exact positions provided in the Supplement, Fig. S1). To corroborate a more recent hypothesis that suggests an alternative mechanism, a poleward shift of the oceanic high-pressure system would stimulate latitude-dependent changes in the magnitude and timing of the upwelling winds (Rykaczewski et al., 2015). A displacement of the pressure systems would increase the standard deviation of the trends around their cores.

3.3 Angular index of upwelling intensification (α_UI)

To test Bakun's hypothesis, a new index, named the angular index of upwelling intensification (α_UI), is proposed. This new index uses the angle between the trend of the most robust upwelling cell at each EBUS and the trend at the corresponding open ocean area. If the upwelling intensifies, as Bakun proposed, the trends in the open ocean and the cell are expected to differ significantly, resulting in a higher angle between the trends. To calculate α_UI, two vectors (in the time–temperature space) may be constructed from the upwelling (Up) and oceanic (Oc) trends. The rotation sense (clockwise or counterclockwise) is used to calculate the relative orientation of (Up) over (Oc), requiring the consideration of an additional unit vector (n) normal to Up and Oc.

The mathematical formulation of α_UI is

$$\begin{array}{rl}{\mathit{\alpha }}_{\mathrm{UI}}& =\arctan \left({\displaystyle \frac{\left(\left(\mathit{U}\mathit{p}\times \mathit{O}\mathit{c}\right)\cdot \stackrel{\mathrm{‾}}{n}\right)}{\mathit{U}\mathit{p}\cdot \mathit{O}\mathit{c}}}\right)=\\ & \arctan \left({\displaystyle \frac{\left|\mathit{U}\mathit{p}\right|\left|\mathit{O}\mathit{c}\right|\sin \left({\mathit{\alpha }}_{\mathrm{UI}}\right)}{\left|\mathit{U}\mathit{p}\right|\left|\mathit{O}\mathit{c}\right|\cos \left({\mathit{\alpha }}_{\mathrm{UI}}\right)}}\left|\stackrel{\mathrm{‾}}{n}\right|\cos \left(\mathit{\beta }\right)\right).\end{array}$$

We used the four-quadrant arc tangent in this analysis since it allows us to determine the sign of the angle based on the signs of the arguments.

Following the right-hand rule, if the cross product of (Up) and (Oc) is counterclockwise – that is, if the open ocean trend (Oc) is greater than the upwelling trend (Up) – the resulting vector of the cross product has the same orientation as $\stackrel{\mathrm{‾}}{n}$. This implies that the dot product ($\mathit{U}\mathit{p}\times \left(\mathit{O}\mathit{c}\right)\cdot \mathit{n}$) is positive since the angle (β) between n and (Up×Oc) is 0°. In the case where the upwelling trend is greater than the open ocean trend, Up×Oc is negative since β is now 180°. Note that this methodology is susceptible to the order of the vectors, as $\mathit{U}\mathit{p}\times \mathit{O}\mathit{c}=-\mathit{U}\mathit{p}\times \mathit{O}\mathit{c}$. Since we are interested in the relative position of the upwelling trend concerning the oceanic waters, the order used is Up×Oc. It is important to note that the angles derived from trigonometric functions are not influenced by units associated with the original vectors. Therefore, this new index is independent of temperature and time units.

We also conducted a probabilistic assessment of uncertainties for α_UI, taking into account the uncertainties associated with upwelling and open ocean SST series. We performed an error estimation using the Monte Carlo method: individual data points were separately and randomly sampled 10 000 times within their respective uncertainty ranges for (Up) and (Oc). These new sampled series were then used to calculate α_UI. The standard deviation of the 10 000 simulations represents the uncertainty of α_UI.

3.4 Satellite validation

Taking advantage of the in situ observations available in the study areas defined in Fig. 2, we performed a validation through linear regression between in situ and the deseasonalized satellite data. Records with no systematic error correspond to a linear regression slope with a value of 1 (perfect). Since the Pacific Ocean is better sampled than the Atlantic, better reanalysis performance and higher correlation are expected in the Pacific Ocean.

4 Results

4.1 Satellite validation

We first carried out a regional validation, in the Northern Hemisphere, with in situ observations (Fig. 2a and b) to test the consistency of the NOAA SST analyses between oceans. A linear fit between in situ and satellite data gives insights into the reliability of the SST data in both oceans. As described in Sect. 2, the in situ dataset was divided into two categories: mooring and cruise data. Because each record has a different spatial and temporal resolution, categories are not comparable. Nonetheless, all the linear fits are statistically significant.

In the Pacific Ocean, the moorings selected were those as close as possible to the areas UP1, UP2, DW1, and OC1 described in Sect. 3.1 and 3.4. We found the strongest correlations inside the upwelling center of this region, where the UP1 buoy and UP2 buoy presented a linear regression fit of 0.99, supported by a very high correlation strength (0.94 and 0.92, respectively; Table 1). In the area where the upwelling meets the oceanographic background, the correlation of the offshore UP1 buoy with the NOAA SST data decays slightly, presenting a linear slope of 0.84 and a correlation strength of 0.86. For the south buoy, close to DW1, the lineal fit and correlation (0.97 and 0.92) are closer to the values inside of the upwelling cell than those of the offshore UP1 buoy (Table 1). The results obtained for the offshore UP1 buoy exemplify the effect of using a large averaging area surrounding the upwelling centers instead of a point, as the average introduces noise by adding points where upwelling might not be taking place, reducing the correlation in the transitional zone.

The cruise observations in the Pacific Ocean (CalCOFI) were divided into two different areas, open ocean (CalCOFI OC) and upwelling (CalCOFI UP), and we avoid data at transitional areas, as seen in Fig. 2a. The results of the linear fit between the cruise and the reanalysis data are similar to those obtained with the mooring data (Table 1 and Fig. S2). The average within the upwelling cell, CalCOFI UP, has a very high linear regression value (0.91), while for the offshore stations, CalCOFI OC, the regression is 0.81. Nevertheless, the strength of the correlation is sensitive to the amount of data available from the cruise data in the ocean versus the upwelling areas. Additionally, due to the time resolution of the cruise, the correlation strength is moderate (0.68) for CalCOFI OC and high (0.71) for CalCOFI UP.

We also compared the EN4 dataset (Table 1, last two columns) against the NOAA SST analysis. As described in Sect. 2, due to the coarse temporal resolution of the cruise data, it was not possible to compare the NOAA SST analysis with the EN4 monthly. Overall, the EN4 has lower correlation and regression values than the NOAA SST analyses. The correlation has the same pattern in both datasets (NOAA and EN4): moderate in the northern locations (0.60 and 0.63) and high in the southern locations (both with 0.76), although the data length of the EN4, unlike the cruise records, is the same as the mooring.

Table 1 Values of the linear regression and Pearson's correlation strength between the in situ data (listed in the first column, with the designation within the paper and, in brackets, the official NBDC station ID). The in situ data are compared with the NOAA SST reanalysis (second and third columns) and the EN4 product (last two columns) validation for the Pacific Ocean. The last two rows are for the CalCOFI cruise data of the Pacific Ocean.

In the Atlantic Ocean, the number of in situ measurements available is more limited than in the Pacific Ocean. Therefore, the long-term records are shorter than in the Pacific Ocean and are only available in the surroundings of the Canary Islands. The results reflect these limitations since none of the linear regressions or correlation coefficients ever exceeds 0.9, unlike the results for the Pacific Ocean (Table 1). As in the Pacific Ocean, the mooring data show a better linear regression and higher correlation coefficient than the cruise data (Table 2 and Fig. S3). The correlation coefficients for the buoys of Santa Cruz, Gran Canaria, and Tenerife Sur are high, with values of 0.89 (0.84), 0.84 (0.87), and 0.83 (0.88), respectively. The result of the Las Palmas east buoy is comparable to the linear slope for the cruise data, 0.71. However, the Las Palmas east buoy has a higher correlation coefficient (0.89), which agrees with the other moorings.

For the cruise data, we follow the same approach as in the Pacific Ocean, dividing the data into open ocean areas (RaProCan OC) and upwelling centers (RaProCan UP). As observed for the results in the Pacific, the regression is better for the upwelling cell with values of 0.73 and a correlation strength of 0.71. For the open ocean areas (RaProCan OC), the linear slope is slightly lower (0.63) but has a higher correlation strength coefficient (0.77).

The EN4 dataset was also used in these areas and showed linear regressions and correlation coefficients that are similar to the cruise data, as shown in Table 1. In Table 2, the linear regressions of EN4 are close to the results of the moorings, even though EN4 performs similarly compared to Table 1. Thus, EN4 seems to be a better alternative in the Atlantic Ocean compared to the Pacific Ocean but only due to the lower performance of the NOAA SST analyses in the Atlantic Ocean. However, the correlation strength of EN4 showed lower values than the NOAA SST analyses, with three out of five locations under 0.70. Hence, using the NOAA SST analyses for long-term analysis in both oceans is a better approach.

Table 2 Values of the linear fit and correlation strength between in situ data (listed in the first column) and both the NOAA SST analyses (second and third columns) and the EN4 product (last two columns) for the Atlantic Ocean. The last two rows are for the RaProCan cruise data of the Atlantic Ocean.

Figure 3 Mapped SST trends (°C per decade) for the major EBUSs: (a) California, CalUS; (b) Canary Islands, CanUS; (c) Humboldt, HuUS; and (d) Benguela, BeUS. The color scale indicates the trend values at the bottom of the figure. The black dots in the shaded areas indicate non-significant trends. The locations of the areas selected (UP1, UP2, OC1, and DW1) are marked with solid black circles and labeled with their names. Black contours enclose the isotrends −0.20, 0, and 0.20 °C per decade.

4.2 Trend patterns in the EBUS

We begin by identifying the overall pattern of the long-term trends in each EBUS, which we show in Fig. 3. A general pattern of cooling (negative) trends within the upwelling centers and warming (positive) trends offshore is found in most of the regions. An exception is found in the HuUS, where there is also a cooling trend offshore. This cooling mode, however, is not as pronounced as the trends onshore driven by the upwelling process. On top of these general patterns, each region presents unique features that are described in the following paragraphs.

The mean trend for the CalUS (Fig. 3a) is 0.10 °C per decade (SD =0.06). The minimum values, −0.17 °C per decade, are located near Cape Mendocino around 43–39° N and 32° N, corresponding to the permanent upwelling centers. On the other hand, maximum values (excluding the shallower Gulf of California), reaching up to 0.20 °C per decade, are located south, at 30° N, where the winds are non-favorable for year-round upwelling, and there is convergence (Kämpf and Chapman, 2016). At around 22° N, there is an offshore negative trend area with mean values of −0.02 °C per decade. However, the MK test revealed that these trends in the offshore area are not significant. The extensive non-significant regions in the MK test support the idea that the trends in the CalUS coastal upwelling are only statistically distinguishable from zero within the upwelling centers.

In the CanUS (Fig. 3b), the mean trend in the region is warmer than in the CalUS, with a value of 0.20 °C per decade (SD =0.04). The minimum values (−0.20 °C per decade) are confined to two upwelling centers in the permanent annual upwelling zones, located north of Cape Ghir and south of Cape Bojador. The maximum warming trend is 0.60 °C per decade, located west of Cape Timiris in the Mauritania–Senegalese convergence zone. Because of the low spatial variability of the trend in the CanUS, the mean offshore value is the same as the average value of the entire CanUS.

For the HuUS (Fig. 3c), the mean value of the region is 0.007 ° per decade (SD =0.09). A cooling signal stands out as the general pattern in the tropic, with a mean value of −0.15 °C per decade. Despite this cooling of the overall trend, there are two clear upwelling centers, at 13 and 24° N, with minimum values of −0.36 °C per decade. In these upwelling centers, the negative trends are stronger than in the open ocean at the same latitudes.

Finally, in the BeUS (Fig. 3d), the mean trend is 0.19 ° per decade (SD =0.09), closer to the average of the CanUS. Two warm fronts with trends over 0.40 °C per decade can be found in the north and south ends. In the year-round upwelling area, the Lüderitz cell (Andrews and Hutchings, 1980; Lutjeharms and Meeuwis, 1987; Peard, 2007), the minimum value is −0.25 °C per decade. In contrast, for the open ocean area between the two warm fronts, the values are similar to those found in the CanUS (Fig. 3b), around 0.20 °C per decade.

The average trend for each region reveals stronger open ocean warmings for the Atlantic EBUS (0.20 and 0.19 °C per decade for the CanUS and BeUS, respectively) than for the Pacific Ocean (0.10 and 0.007 °C per decade for the CalUS and HuUS, respectively), with the weaker trends observed in the HuUS. However, the trends along the coast are rather heterogeneous, responding to the variability of the local upwelling dynamics, as seen in Fig. 3, where several permanent upwelling centers (negative trends) exist along the coast on both continents.

4.3 Trends in the upwelling cells and open ocean areas

Although we gave an overview of the long-term trends for each EBUS in the previous section, the validation of Bakun's hypothesis would require a finer resolution to describe the local dynamics of the upwelling and identify areas where the coastal upwelling is the main forcing. In this sense, we have chosen representative points (see criteria in Sect. 3.4 and Figs. 2 and 3) of both EBUS and non-EBUS areas instead of using the large-averaged regions shown in Fig. 3. These points include the upwelling centers (UP1, UP2), the nearshore areas where the upwelling is not the primary process (DW1), and the open ocean areas (OC1).

Figure 4 SST time series and trends (°C per decade) for the selected areas for each of the EBUSs (DW1, UP1, and UP2) and the areas representative of the open ocean (OC1) in each EBUS: (a) California, (b) Canary Islands, (c) Humboldt, and (d) Benguela. The blue lines are for negative (cooling) trends, and the red ones are for positive (warming) trends. On the right side of the y axis, the trends are shown in degrees Celsius per decade, along with the area's label and confidence interval. Mann–Kendall significant trends (p value <0.05) are marked with an asterisk next to their value.

The CalUS presents the weakest trends of all the EBUSs in both year-round upwelling centers, being −0.06 °C per decade for UP1 and non-significant in UP2 (−0.03 °C per decade, Fig. 4a). For the open ocean area, the OC1 trend is positive with a value of 0.14 °C per decade, a lower value than the trend of 0.22 °C per decade of the DW1. In the Atlantic, the CanUS (Fig. 4b) possesses the only positive trend (0.07 °C per decade) of all the EBUSs in an upwelling area (UP1). Despite UP1 being positive, the trend is closer to zero than in OC1 (0.20 °C per decade). Furthermore, the UP2 cell in the CanUS shows a trend twice as negative (−0.15 °C per decade) as the one in the CalUS upwelling UP1. The OC1 and DW1 areas show a warmer trend in the CanUS than in the North Pacific region (Fig. 4). All of the above suggests that in the Northern Hemisphere EBUS, Bakun's hypothesis is fulfilled, and this is even more significant in the CanUS despite having been dismissed in previous studies (Sydeman et al., 2014).

In the HuUS (Fig. 4c), we find a different behavior than the one seen in the other EBUSs, showing negative trends in all the representative locations. The upwelling centers of the HuUS present the greatest cooling trend of all the EBUSs, −0.30 °C per decade at UP1 and −0.26 °C per decade at UP2. As observed in the CanUS, the values for the HuUS in OC1 (0.06 °C per decade) are similar to the ones of DW1, although DW1 is non-significant in the HuUS. Its counterpart in the Atlantic Ocean, the BeUS (Fig. 3d), presents significant large negative trends in the year-round upwelling areas (−0.23 °C per decade for UP1 and −0.10 °C per decade for the UP2) and positive trends of 0.14 °C per decade in OC1. The trend of DW1 is the warmest (0.57 °C per decade) found in all the EBUSs, and it is related to the warm inflow from the Indian Ocean.

Overall, the trends show warming in the OC1 areas and cooling in the upwelling areas, except for the HuUS, where the trend at OC1 is also slightly negative. The contrast between the trends of upwelling and open ocean areas found throughout the EBUS indicates upwelling intensification. To quantify this intensification, it is necessary to have an index that compares the intensification of the upwelling with the global warming trend in the open ocean area, which can be compared for all the EBUSs. To compare the intensification between upwellings and to further understand the impact of the oceanic background on these trends, we use the index α_UI described before to normalize the trends of each EBUS.

Figure 5 EBUS SST trends (°C per decade) against non-EBUS area trends for each region. The EBUS annual series (continuous line) and trend (dotted line) are shown in blue; the same is true for the non-EBUS annual series but in red. On the opposite side of the y axis is α_UI with its corresponding EBUS.

4.4 The relation between oceanic and EBUS trends

Global warming induces the increase in oceanic SST, and under Bakun's hypothesis, it also enhances the upwelling-favorable winds responsible for intensifying the upwelling areas. We define an angle between the upwelling and oceanic trends, as described in the Methodology section, to discern Bakun's hypothesis from the global increase in SST. As described previously (Fig. 4), the upwelling centers UP1 (except for the CanUS, where UP2 is used due to the warming detected in UP1 in Sect. 4.3) have the strongest cooling, and we use these trends (hereafter in this section, UP) to create the angle contrasts with the positive warming trend of the open ocean (OC1) area for each EBUS (Fig. 4).

For CanUS, the α_UI obtained between the UP and OC1 trends is 20°±2°, whereas it is 11°±3° in the CalUS, 21°±2° for the BeUS, and 14°±3° for the HuUS. The smallest angle is found in the CalUS because of the low cooling and warming trends described in Sect. 4.3. This result for the CalUS is in agreement with the overall non-significant trends in the Mann–Kendall test (Fig. 3a). The BeUS and CanUS have the highest contrast between the two regression lines (recall that we observed similar UP and OC1 trends in the Atlantic Ocean in Sect. 4.3), presenting similar α_UI.

The BeUS and CanUS show a weaker negative trend than the HuUS, but the oceanic background at the HuUS leads to a smaller angle. At the HuUS, a negative SST trend is observed for the whole study area. The existing hypotheses suggest that this trend is either led by a stronger cool phase of the El Niño–Southern Oscillation (ENSO) or related to the Southern Ocean SST changes (Meehl et al., 2016; Kang et al., 2023b, a). Nevertheless, at the HuUS, we have the most prominent upwelling negative trend. However, when normalized with the open ocean trend, Bakun's effect is reduced. In the annual upwelling series of the CalUS and HuUS, two prominent peaks associated with the warm phase of ENSO around 1983 and 1997 are observed.

In general, we found positive α_UI for all the EBUSs, supporting the intensification of the upwelling–oceanic gradient, as expected from Bakun's hypothesis.

Figure 6 Spatial distribution of α_UI (over the period of 1982–2021) along the coast for the CalUS (a), CanUS (b), HuUS (c), and BeUS (d). α_UI is calculated between grid points along the coast and OC1.

4.5 Latitudinal distribution of α_UI

Many authors have previously tested the Bakun hypothesis, providing little consensus on both historical and projected records (Barton et al., 2013; Belkin, 2009; McGregor et al., 2007; Sambe et al., 2016; and Sydeman et al, 2014). Such controversy has yielded alternative hypotheses to explain changes observed in the magnitude and timing of upwelling processes. Rykaczewski et al. (2015) suggest an alternative mechanism to the intensification of the upwelling process. They suggest a poleward shift of the oceanic high-pressure system which would stimulate latitude-dependent changes in the upwelling winds. To address this, we have calculated the latitudinal distribution of α_UI (see Fig. 6) in each EBUS. The spatial variability of the upwelling intensity proxy, α_UI, reveals distinct patterns and regional differences. In the CalUS, upwelling intensification demonstrates consistent upwelling activity between 35 and 45° N, with α_UI values reaching up to approximately 10° (Fig. 6a). Conversely, in the CanUS, significant upwelling intensification is observed between 20 and 30° N, with α_UI values peaking at 20° and in locations consistent with our dynamical analysis based on the literature review (Fig. 2b, UP1 and UP2). Similarly, in the HuUS upwelling intensification is confined to low latitudes (10–20° S, Fig. 6c), and the values are close to those of the CalUS (index values around 10°), as seen in the previous section. In contrast with the other regions, the BeUS shows intensification at high latitudes, with maximum values of α_UI (20°) in the upwelling center of this region – the Lüderitz upwelling center at 25° S and Cape Columbine (around 32° S). While results of the BeUS and CalUS appear to be consistent with the findings by Rykaczewski et al. (2015), there is no supporting evidence in the other regions. To elucidate the possible mechanism responsible for such differences, we attend to the driver of the upwelling-favorable wind, the sea level pressure gradient.

Figure 7 EBUS SLP gradient trends and temporal series for NCEP (blue lines) and ERA5 (red lines) datasets over the period of 1982–2023.

4.6 SLP gradients

The coastal upwelling intensification postulated by the Bakun mechanism in 1990 would involve a stronger increase in near-surface temperature over land than over the ocean, which would lead to an intensification of the continental thermal low-pressure system relative to the ocean. To test this driver mechanism, we calculated the trends (Fig. 7) in the gradient between the continental thermal low and the oceanic high pressure.

ERA5 data show positive and significant trends across all EBUSs (see Table 3), while NCEP data indicate negative trends in the BeUS. Despite these differences, both datasets show good overall agreement. The strongest SLP gradient trends are found in the HuUS region, whereas the weakest trends occur in the BeUS. Given its coarser resolution (2°) compared to ERA5 (0.25°), NCEP data are considered less reliable. Despite these findings, both datasets support an intensification of the pressure gradient.

Table 3 Values of the trend, over the period of 1982–2023, for ERA5 (first row) and NCEP (second row) for all the EBUSs. Parentheses enclose the spatial standard deviation.

5 Discussion and conclusions

Bakun proposes an intensification of the upwelling due to the increase in the continental low-pressure system driven by global warming. However, controversies arise from discrepancies between wind stress datasets and differences in the methodologies used.

On the one hand, Barton et al. (2013) highlighted a lack of consensus among various wind datasets, since they did not observe statistically significant changes in the meridional (upwelling-favorable) wind. These discrepancies in wind data are consistent with those noted by Narayan et al. (2010), who, despite finding significant increases in coastal upwelling areas when using the Comprehensive Ocean-Atmosphere Data Set (COADS) wind stress, also found that the NCEP/NCAR wind stress indicated a significant decrease in the upwelling off NW Africa and a non-statistically significant trend for Lüderitz, California, and the Peruvian upwelling areas. Furthermore, they also observed that the ERA-40 dataset showed an increasing coastal upwelling in the NW African and Peruvian upwelling areas but a decrease in the California upwelling areas, with a non-statistically significant trend in the Lüderitz upwelling areas. Therefore, using wind data as a proxy for upwelling leads to a wide spread of results as it strongly depends on the data product used.

On the other hand, SST-based indexes are usually constructed from thermal differences between coastal and offshore SST areas taken at the same latitude and following the coastline (Benazzouz et al., 2014; Gómez-Gesteira et al., 2008; Santos et al., 2012). This methodology does not consider the regional upwelling dynamics and averages upwelling centers with areas without upwelling. Abrahams et al. (2021) introduced an upwelling metric based on marine heatwave detection techniques, examining upwelling-favorable winds and SST data together. Their findings revealed a strong association between a decrease in SST and an increase in upwelling intensity. Their novel methodology holds significant importance for unraveling the connection between the physical upwelling phenomenon and its ecological impacts. However, predicting ecological impacts remains challenging. While intensified upwelling could mitigate habitat warming, it may also increase ocean acidification and hypoxic events and reduce suitable food for fish larvae (Abrahams et al., 2021; Bakun et al., 2015). Nonetheless, they successfully establish a link between decreased SST and changes in upwelling intensity, even when trends in wind dynamics do not fully account for the upwelling response, reinforcing the notion that SST is a suitable proxy for upwelling intensity. However, their SST metrics exhibited inconsistencies across upwelling areas, except for the Humboldt system. These inconsistencies may be due to the averaging of data across extensive areas, mixing upwelling areas with areas without the associated cold water of upwelling. Abrahams et al. (2021) also explored these metrics in upwelling-favorable wind data, and their results indicated that decadal trends were generally not significant. As previously discussed, wind products often yield contradictory results despite their direct relevance to upwelling. Hence, our study complements Abrahams et al. (2021) since we have focused on SST to understand longer-term changes in the upwelling intensity, using areas with an optimal signal-to-noise ratio, namely the upwelling centers, revealing upwelling-related cool water in all eastern boundary upwelling systems (EBUSs).

Therefore, in this study, we assess the intensification of the upwelling from a regional perspective by using SST trends at locations representative of upwelling and of an open oceanic reference location for each EBUS.

Additionally, we tested the effects of averaging areas on the index (see Fig. S4 and Table S1). Our findings indicate that the averaged response is influenced by the dynamical regions involved rather than by the size of the region averaged. This is evidenced by the invariant results when including the three coastal areas (DW1, UP1, and UP2). In contrast, focusing on specific upwelling zones, particularly around upwelling centers, made the intensification more evident. Moreover, we verified the stability of the trend both spatially and temporally by performing the analysis of Barton et al. (2013) across all the EBUSs (Fig. S5).

Furthermore, to assess the strength of the net upwelling intensification, we proposed an index that allows for inter-basin comparisons while attending to their regional background. SST, often used as an indicator of coastal upwelling, can be influenced by various factors, such as changes in surface mixing and offshore storm activity. However, in our long-term analysis of monthly and deseasonalized SST records, the seasonal and synoptic processes have minimal influence on the SST–upwelling intensity relationship. Moreover, Wang et al. (2015) explored the connection between sea–land thermal gradients and offshore Ekman transport using the CMIP5 models. Their findings underscore the significant link between thermal gradients and offshore Ekman transport, even under greenhouse gas emission scenarios. McGregor et al. (2007) and Santos et al. (2012) also support this relationship, emphasizing significant correlations between coastal SST and offshore Ekman transport, reinforcing the utility of coastal SST as a proxy for assessing upwelling intensity.

To assess the quality of our results, we validated the NOAA SST reanalysis with in situ data from both the Atlantic Ocean and the Pacific Ocean before estimating the trends in all EBUSs. Overall, the Atlantic Ocean had lower correlations with the satellite data than the Pacific Ocean, likely due to shorter in situ records. Nevertheless, we found high and robust correlation coefficients (>0.7) that sustain the satellite SST trends in oceanic and upwelling areas. We observed negative SST trends in all the EBUSs, being stronger in the Southern Hemisphere (with the strongest located in the HuUS UP1, showing a trend of −0.30 °C per decade) than in the Northern Hemisphere (with the weakest in the CalUS UP2, with a trend of −0.06 °C per decade). Our results are consistent with the meta-analysis by Sydeman et al. (2014), who concluded, from observational and model data, that a significant intensification of upwelling exists, except for the case of the CanUS.

Other studies have investigated the SST trends in the EBUS but with an approach that did not consider the heterogeneity of the upwelling areas. For instance, in the CalUS, Seabra et al. (2019) reported a 0.06 °C per decade warming rate over the period of 1982–2018. However, their approach involved averaging a 500 km nearshore area, excluding non-significant regions. Thus, almost half of the extension was not considered, resulting in the average of different dynamical areas and the exclusion of upwelling centers. Belkin (2009) performed a similar analysis but included the entire CalUS nearshore area. They found a net change of −0.035 °C per decade over the 1982–2007 period, agreeing in sign with our study but showing a weaker trend due to the use of a large average. In contrast, Siemer et al. (2021) found negative trends of −0.14 °C per decade over the 1982–2019 period for the CanUS permanent upwelling area, like our results (−0.15 °C per decade). However, this trend fades away and becomes positive when they average the whole coastal upwelling area, highlighting the relevance of the methodology used in this study. Likewise, many studies carried out in this area present positive trends for the upwelling due to the method used (Belkin, 2009; Demarcq, 2009; Seabra et al., 2019). In line with our study, Seabra et al. (2019) found the largest cooling trends ($-\mathrm{0.07}\pm \mathrm{0.08}$ °C per decade) in the HuUS. However, like in other EBUSs, using averaged areas increased the trend values. This pattern is also observed by Belkin (2009), where the net change in the averaged nearshore area results in −0.05 °C per decade. For the BeUS, similar to the CanUS, averaging the entire coastal upwelling area results in the fading of the observed upwelling trend. Hence, a warming rate of 0.17 °C per decade is found in Seabra et al. (2019). In contrast, Santos et al. (2012) investigated trends over the period of 1982–2010 close to the shore without averaging areas and found a negative trend in the BeUS, strongly agreeing with our results (−0.13 °C per decade). Hence, a warming rate of 0.17 °C per decade is found in Seabra et al. (2019). In contrast, Santos et al. (2012) investigated trends close to the shore without averaging areas and found a negative trend in the BeUS, strongly agreeing with our results (−0.13 °C per decade).

While all the upwelling trends are negative and support Bakun's hypothesis, the oceanic trends behave differently across basins. We observed warming in all the open ocean areas except in the HuUS, where a cooling of −0.06 °C per decade is observed. Dong and Zhou (2014) studied the influence of the Interdecadal Pacific Oscillation (IPO) on global warming trends. Their empirical orthogonal function (EOF) analysis results indicate that the transition to the negative phase of the IPO modes is responsible for the cooling trends observed in the Pacific.

The warming in the CalUS and BeUS is 0.14 °C per decade (over the period of 1982–2023), while this trend is slightly more prominent in the CanUS. Seabra et al. (2019) revealed oceanic warming rates (0.06°C per decade) over the period of 1982–2018 on the averaged upwelling in the CalUS, which is lower than the OC1 trend (0.14 °C per decade). The open ocean positive trend of the CanUS is identical to the one in Siemer et al. (2021) and further agrees with other studies (Belkin, 2009; Good et al., 2007; Signorini et al., 2015). The result of Seabra et al. (2019) in the HuUS also showed a very similar trend (−0.07 °C) compared with our OC1 trend. Finally, in the BeUS, good agreement is found with the average warming rate of Seabra et al. (2019). Our study demonstrates good agreement with existing literature on oceanic trends despite the differences in methodologies employed.

Although long-term changes, such as the North Atlantic Oscillation or the Pacific Oscillation, can impact the SST gradient, their effect would not surpass the ability of our analysis to support Bakun's hypothesis. In that sense, Nayaran et al. (2010) found that correlations between upwelling indices and climate indices like the Atlantic Multidecadal Oscillation Index (AMOI) lack significance. Similarly, the North Atlantic Oscillation Index (NAOI) exhibits a notable negative correlation with meridional wind stress off NW Africa, yet its correlation with the SST index remains insignificant. In the case of the CalUS, the Pacific Decadal Oscillation Index (PDOI) shows a weak but statistically significant correlation with the coastal upwelling SST index off California. However, no substantial correlation is found with alongshore wind stress. Cross-correlation analyses also reveal a lack of significant correlations across various time lags. On the other hand, Bonino et al. (2019) found that local drivers and trends favoring upwelling (e.g., equatorward wind stress, cyclonic wind stress curl, and thermocline depth variation) explain the low-frequency modulation of upwelling. Bonino et al. (2019) also explored the link between wind-based upwelling indices and climate modes. They found that Atlantic and Pacific upwelling variabilities are mainly independent, while intra-basin domain variabilities present some coherency, which is consistent with our results. This intra-basin covariability is especially marked in the Pacific Ocean, where the shared variability is majorly due to the ENSO mode. In contrast, in the Atlantic Ocean, coherent variability is associated with upwelling trends, whereas only in the CanUS is it linked to the AMO. These results suggest that long-term climate indices may influence coastal upwelling dynamics, which is especially important in the Pacific. However, by normalizing the trend for the oceanic background of our index, α_UI, our results should account for the effects of local climate indices.

To assess Bakun's hypothesis and, thus, the upwelling capacity to overcome the oceanic warming effect, we define the angle (α_UI, readers are referred to Sect. 3.3) between oceanic water and upwelling trends. Because this new index is directly based on trends, it captures only the low-frequency variability. Additionally, we verified the method's robustness using a probabilistic assessment of the uncertainties that showed consistent intensifications for all EBUSs (Fig. 5). This new approach differs from the traditional trend analysis since it normalizes the upwelling trends by comparing them with open ocean changes.

The EBUS in the Pacific Ocean yields minimum α_UI (10°±3° and 14°±3° for the CalUS and HuUS, respectively), which is consistent with the low signal-to-noise ratio of global warming in this ocean, given its natural variability. The overall cooling signal caused by the IPO enhances the HuUS open ocean negative trends. Still, our index normalizes the upwelling trend to the full basin variability, suggesting the possibility of a mild Bakun effect even at the HuUS. In the Atlantic Ocean, α_UI for the CanUS and BeUS is 20°±2° and 21°±2°, respectively, twice as large as in the Pacific Ocean. α_UI presents wider angles at the Southern Hemisphere EBUS than in the Northern Hemisphere EBUS. Nevertheless, our results show a significant difference between oceanic and coastal trends reflected in positive α_UI in all EBUSs (Fig. 5).

The SST changes in the EBUSs respond mainly to changes in the upwelling processes, which are ultimately driven by the pressure gradients. We analyzed the pressure gradient trends in all four EBUSs. Our findings further support the intensification of the pressure gradients driven by climate change, as stated by Bakun (1990). However, there are probably other contributors to the intensification of the upwellings. Some researchers question whether the impacts of differential heating on the pressure gradient force drive intensification of coastal upwelling. Rather, a complementary hypothesis proposes that evidence of an intensifying pressure gradient force is limited to poleward migration of the Hadley cell (Arellano and Rivas, 2019; Rykaczewski et al., 2015; Wang et al., 2015). Nevertheless, these projections are only supported by observational records in the Humboldt and Benguela systems (Sydeman et al., 2014). In contrast, we have tested this hypothesis on the historical record by computing the latitudinal distribution of α_UI. The results shown in Fig. 6 partially agree with Rykaczewski et al. (2015), as only the CalUS and BeUS presented a poleward intensification of α_UI. To further understand the drivers of these changes, we examined the spatial stability of the trends in the SLP continental–oceanic gradient through Monte Carlo simulation. The discrepancy between the latitudinal distribution of α_UI and the small standard deviation of trends around the cores of the pressure systems suggests that the hypothesis of poleward displacement of the high-pressure systems remains inconclusive.

In summary, in this study, we use SST at discrete locations and the pressure gradient to explore Bakun's hypothesis in the four major EBUSs. Cooling trends are observed for all upwelling areas (with the strongest in the HuUS and the weakest in the CalUS) and mainly warming trends offshore except for the HuUS. In addition, a novel index α_UI that normalizes the upwelling trends to their background open ocean trend is proposed. This index is easy to estimate, allows for inter-basin trend comparisons, and helps us to understand the role of changing upwellings in a changing climate. The index reveals that the Bakun hypothesis remains a possible mechanism for upwelling intensification in all four EBUSs, although the Atlantic Basin shows a stronger intensification effect than the Pacific Ocean.