Variation and Change of Upwelling Dynamics Detected in the World’s Eastern Boundary Upwelling Systems

Global increases in temperature are altering land-sea temperature gradients. Bakun (1990) hypothesized that changes within these gradients will directly affect atmospheric pressure cells associated with the development of winds and will consequently impact upwelling patterns within ecologically important Eastern Boundary Upwelling Systems (EBUS). In this study we used daily time series of NOAA Optimally Interpolated sea surface temperature (SST) and ERA 5 reanalysis wind products to calculate a series novel of metrics related to upwelling dynamics. We then use these to objectively describe upwelling signals in terms of their frequency, intensity and duration throughout the four EBUS during summer months over the last 37 years (1982–2019). We found that a decrease (increase) in SST is associated with an increase (decrease) in the number of upwelling “events,” a decrease (increase) in the intensity of upwelling, and an increase (decrease) in the cumulative intensity of upwelling, with differences between EBUS and regions within EBUS. The Humboldt Current is the only EBUS that shows a consistent response from north to south with a general intensification of upwelling. However, we could not provide clear evidence for associated changes in the wind dynamics hypothesized to drive the upwelling dynamics.


INTRODUCTION
Coastal upwelling is a major oceanic process driven by prominent currents, of which those within Eastern Boundary Upwelling Systems (EBUS) are most important globally (Bakun and Nelson, 1991;Messié et al., 2009;Gruber et al., 2011;Pegliasco et al., 2015;Varela et al., 2015Varela et al., , 2016Varela et al., , 2018Bonino et al., 2019;Brady et al., 2019). EBUS include the California (CCS), Humboldt (HCS), Canary (CnCS), and Benguela (BCS) current systems (Figure 1), with each of these significantly impacting their associated coastal ecosystems. These systems are present along the western shores of landmasses in the Pacific and Atlantic Oceans where they comprise vast regions of coastal ocean water (Bakun, 1990;Pauly and Christensen, 1995;Bakun et al., 2010Bakun et al., , 2015Santos et al., 2012b,a;Seabra et al., 2019). EBUS encompass a multitude of coastal regions stratified across several latitudes. Accordingly, these systems vary heterogeneously in terms of their environments and atmospheric conditions (Wang et al., 2015). As a result, variability within each EBUS naturally differs as the oceanic variables are driven by and respond to the combined effects of an assortment of differing oceanic and atmospheric processes and phenomena. For example, ocean temperatures within the CCS vary due to the combined effects of both El Niño-Southern Oscillation (ENSO) and Pacific Decadal Oscillation (PDO) (Jacox et al., 2015). Similarly, North Atlantic Oscillation (NAO), Benguela Niños, and Pacific ENSO are primarily responsible for driving variability in environmental conditions in the CnCs, BCS, and HCS respectively (Minobe, 1999;Chhak and Di Lorenzo, 2007;Di Lorenzo et al., 2008;García-Reyes et al., 2015;Gómez-Letona et al., 2017). One potential source of shared variability is suggested to stem from anthropogenically-mediated climate change (Bakun, 1990); although the anticipated direction of the response remains uncertain Varela et al., 2018;Bonino et al., 2019), no study, yet report shared variability and a consistent response across all EBUS (McGregor et al., 2007;Narayan et al., 2010;Patti et al., 2010;Pardo et al., 2011;Varela et al., 2015Varela et al., , 2018Bonino et al., 2019).
Eastern Boundary Upwelling Systems are a common focus of oceanographic and biological research because the complex interplay of biotic and abiotic processes within EBUS results in them being highly productive (Bakun, 1990;Pauly and Christensen, 1995;Tretkoff, 2011;Varela et al., 2015), diverse, and abundant in marine life. EBUS provide up to 20% of the world's fishery output despite only covering <1% of global ocean area (Bakun et al., 2010;Bakun et al., 2015). Warm-season upwelling, driven by equatorial winds advected offshore by the Coriolis effect, is primarily responsible for the high levels of biological productivity present within EBUS (Huyer, 1983;Borges et al., 2003;Chavez and Messié, 2009;García-Reyes and Largier, 2010;Varela et al., 2018). These regions provide both lucrative economic (Costanza et al., 1997) as well as significant recreational services to people living along these coastlines, linking indirectly to the rest of the world. Recent studies have shown ecological changes in EBUS ecosystem structure (Fréon et al., 2009;Wang et al., 2015), hence monitoring these systems is becoming increasingly important. EBUS are formed in part by wind-driven ocean circulation, and their upwelling is dependent on wind direction and strength (Capet et al., 2004;Messié and Chavez, 2015;Steinfeldt et al., 2015). Recognizing that alongshore winds that drive upwelling are initiated by changes in atmospheric pressure gradients at the cross-shore, Bakun (1990) hypothesized that long-term changes in climatic conditions would likely intensify continental oceanic pressure gradients (Bakun, 1990;García-Reyes et al., 2015) and subsequently result in an increase in the frequencies and intensities of upwelling-favorable winds. Understanding, therefore, how these winds will change is of high importance for anticipating how upwelling might respond (Varela et al., 2015). For example, weaker upwelling may limit nutrient enrichment and potentially impact primary production (Figueiras et al., 2002;Chhak and Di Lorenzo, 2007;García-Reyes et al., 2015). In contrast, stronger upwelling may increase nutrient input and therefore offshore transport (Bakun et al., 2010(Bakun et al., , 2015. Increased wind intensity could also induce changes in water turbulence (Cury and Roy, 1989) which could affect chemical mechanisms and processes like ocean acidification and deoxygenation (Gruber et al., 2011) that may ultimately impact on productivity (García-Reyes and Largier, 2010). Elucidating the drivers of upwelling and their relationship to nutrient enrichment and productivity could therefore inform conservation and facilitate management of fishers and other EBUS dependent resources. Such understanding is particularly important toward predicting the impacts of future climate change scenarios.
Recent reports from studies investigating patterns related to changes in upwelling winds have produced conflicting results, and a consensus on current trends has not been reached (Phillips, 2005;Bakun et al., 2010;Varela et al., 2015). This is largely due to contrasting findings regarding seasonal, interannual, and decadal fluctuations of unidirectional winds because of limited time series (Bonino et al., 2019;Seabra et al., 2019). Some researchers question whether the impacts of differential heating on the pressure gradient force drives 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 Brady et al., 2017;Grise et al., 2019;Grise and Davis, 2020). Further complications preventing the agreement on findings in EBUS largely include (a) researchers deriving conclusions from non-comparative datasets, (b) data and analyses being treated in an inconsistent manner, (c) variable quality between datasets, and (d) inconsistencies in measurement techniques (Jacox et al., 2015;Bonino et al., 2019).
Upwelling has been investigated for several decades using a variety of statistical models and simulations (e.g., Bakun, 1990, Shannon et al., 1992Bakun et al., 2010;Varela et al., 2015;Wang et al., 2015;Bonino et al., 2019). However, a common theme among these attempts has been the broad temporal scale at which estimates were made. Most of studies on upwelling trends have used wind and temperature variables captured and averaged at monthly intervals. Here, we aimed to test the efficacy of a novel method for detecting upwelling signals and characterizing them in terms of intensity, frequency, and duration of upwelling "events" in a more reliable and objective manner. This is accomplished by using data at a fine temporal scale of daily intervals. In this context, the atmospheric and oceanic mechanics responsible for coastal upwelling are interdependent, and changes in one variable, such as wind speed, should directly affect variables such as upwelling intensity (Varela et al., 2015Bonino et al., 2019). The objective of this study was to quantify the variation and changes in upwelling signals over time. The phenology of many marine ecosystem processes is highly affected and dependent on upwelling events such as its duration, frequency, and intensity (Barth et al., 2007). For example, intertidal communities require intermediate disturbance on rocky shores to maintain diversity (McGregor et al., 2007;Landry et al., 2009). Upwelling events may cause substantial or minimal disturbance depending on the frequency, duration, and intensity of the events (Benoit-Bird et al., 2019). However, upwelling also affects other marine communities and changes in upwelling patterns are likely to dramatically impact them (Wang et al., 2015). We therefore set out to detect if changes occur in the (a) SST patterns, (b) intensity, duration, and frequency of upwelling-favorable wind, and (c) the frequency, mean, and cumulative intensity of upwelling signals during summer months over a period of 37 years. Understanding changes within these systems could allow for predictive management and conservation of fisheries and marine resources.

Data
To evaluate if there are changes in upwelling dynamics, this study used the gridded data of the global 0.25 • National Oceanic and Atmospheric Administration (NOAA) daily Optimally-Interpolated Sea Surface Temperature (dOISST, v.2.1) (Reynolds et al., 2007;Banzon et al., 2016). This remotely-sensed satellite product has any gaps in spatial coverage filled through the interpolation of data collected from ships, and buoys (Reynolds et al., 2007;Banzon et al., 2016). The OISST data product has been collected for nearly four decades, hence providing us with a long (>30 years) time series from which upwelling trends and their rate of change can be calculated. High resolution SST data sets are useful but often have limited time series. Most highresolution data products (i.e., those with a resolution of ∼0.01 • ) do not exceed 30 years in duration; for example, MUR (2002MUR ( -2020 and G1SST (2010G1SST ( -2019. To detect variations in upwelling signals among the four EBUS (Figure 1; Bakun, 1990;Bakun et al., 2015;García-Reyes et al., 2015;Varela et al., 2015Varela et al., , 2016Sousa et al., 2017), wind speed and direction were also necessary. Wind speed and direction variables were downloaded from the ERA5 climate reanalysis produced by ECMWF, providing daily data on regular latitude-longitude grid at 0.25 • × 0.25 • resolution (C3S, 2017). Upwelling zones representative of each EBUS were selected according to previous studies Varela et al., 2015Varela et al., , 2018Varela et al., , 2020Seabra et al., 2019). Two upwelling systems were analyzed along the Eastern Atlantic Ocean, these

Upwelling Identification
To identify whether changes exist in the frequency and intensity of upwelling signals across EBUS it must first be determined when upwelling occurs. To do this, a set of upwelling threshold values needed to be established. Given that upwelling is primarily caused by alongshore, equatorward winds, both SST and wind data were used. Wind data were used to calculate the upwelling index using the formula presented in Fielding and Davis (1989): Where µ represents the wind speed (ms −1 ), θ represents the wind direction in degrees, and τ is the orientation of the coastline (Jury 1980). This index relies on wind speed and direction to identify upwelling-favorable conditions. When the upwelling index is >0, SST usually drop, as expected, suggesting that upwelling is occurring. Here we established that the drop in SST that coincided with a positive upwelling index was close to the seasonally varying 25th percentile threshold for SST. This threshold temperature was subsequently used in combination with the upwelling index to identify upwelling signals. With the upwelling index and SST threshold set, the detect_event() function from the heatwaveR package (Schlegel and Smit, 2018) was used to calculate metrics for the upwelling signals. By using this method of detecting signals, we were able to obtain upwelling metrics (Table 1) such as the frequency of occurrence, mean intensity, and cumulative intensity. Here, we defined the mean intensity ( • C) as the mean temperature exceedance (i.e., duration below the 25th percentile) during the upwelling signal. Cumulative intensity ( • C days) is defined as the sum of daily intensity exceedances across the duration of an event. Because upwelling signals were calculated relative to seasonal climatologies of percentile exceedances, rather than an absolute definition such as temperatures below a fixed temperature threshold, these signals could occur any time of the year; however, upwelling was proven to be more dominant during summer months, as expected. To quantify metrics of upwelling-favorable winds [i.e., southeasterly (SE) winds in the Southern Hemisphere and northeasterly (NE) winds in the Northern Hemisphere], we again made use of the detect_event() function. For these analyses, by setting the seasonality and threshold to a value of 0 allowed the statistical functions of the package to compensate for using wind data rather than SST data. This allowed us to estimate wind metrics such as the duration, intensity, and count of upwelling-favorable wind events within each of the regions.
With upwelling being influenced by wind patterns, Bakun (1990) predicted that changes in wind patterns would influence the intensity of upwelling signals. As a first step, we calculated trends in SST to understand how temperatures within these cold-water regions are changing over time. We selected all coastal SST pixels along the latitudinal bands shown for each EBUS (Figure 1), as the coastal region is the most utilized for human habitation, leisure, recreational activities, or tourism. Given that the OISST time series length is greater than 30 years, it should be possible to discern long-term trends within the data from interannual noise (Hobday et al., 2018). Using linear regression analyses, we observed if there were changes in the number of times wind blew in an upwelling favorable direction, the duration of those winds, and their intensity, with an emphasis on austral (DJF) and boreal (JJA) summer months. The month variable was selected to observe which month expressed the most intense signal. To establish whether differences existed between currents and upwelling metrics over time, we assessed the upwelling metrics as a function year or month. Thereafter, using regression analysis, we compared the mean numbers of upwelling signals detected in each of the currents.

RESULTS
Our investigations of trends in upwelling metrics across individual time series for each of the EBUS currents revealed the presence of noticeable shifts in SST patterns at each EBUS over the past 37 years (Figures 2A,B). We found that SST dramatically Sum of the daily intensity anomalies over the duration of the signal differed across each EBUS (two-way ANOVA: F = 1878.76, SS = 6966, p < 0.05). This was in part driven by variable increases in SST in the BCS North (Table 2). Conversely, the BCS South and HCS showed a significant decrease in SST ( Table 2). We further found a significant negative trend in monthly SST within the HCS, specifically during the months of January near Chile (R 2 = 0.01, slope = −0.2 • C.dec −1 , and p < 0.05) and during December near Peru (R 2 = 0.03, slope = −0.4 • C.dec −1 , and p < 0.05). The number of upwelling signals decreased over time in the BCS North and CnCS, but slightly increased in HCS (Figure 2B). Increases in detections were particularly prominent when SST decreased (Figure 2A). A regression analysis showed a significant negative trend in the number of upwelling signals detected in the BCS North and CnCS ( Figure 2C). However, a positive trend is present in HCS Chile, HCS Peru. Conversely, a significant positive trend was detected in the mean intensity of upwelling in the BCS North and CnCS ( Table 2). Results of a linear regression showed significant negative trends in the mean intensity of upwelling signals in the BCS South, HCS Chile and Peru and CCS North. We also found that the cumulative intensity of upwelling signals significantly differed between the EBUS (ANOVA: F = 371.15, SS = 178853, and p < 0.05) and between month (ANOVA: F = 139.78, SS = 123484, and p < 0.05; Figure 2D). Results of a linear regression showed significant negative trends in cumulative intensity within the BCS South, HCS Chile and Peru and CCS North.
Counts of discrete upwelling-favorable wind events were not uniform among the four EBUS ( Figure 3A). A two-way ANOVA analysis showed a significant difference in upwelling-favorable wind events between currents (F = 1878.76, SS = 6966, and p < 0.05). However, a regression analysis showed no significant changes in wind events detected over the past 37 years; the CCS North was the only exception showing a significant negative trend in upwelling-favorable winds ( Table 2). A significant positive trend was detected in the number of monthly upwelling-favorable winds during July in the CCS South (R 2 = 0.25, slope = 0.05 count dec −1 , and p < 0.05) and during February in the HCS Chile (R 2 = 0.14, slope = −1.1 count dec −1 , p < 0.05) within the HCS. A significant negative trend was present during December within the HCS (R 2 = 0.10, slope = −1.1 count dec −1 , p < 0.05).
The typical duration of south-easterly winds in the BCS and HCS Chile ranged between 3-6 days with the HCS Peru being the only exception with a duration >10 days ( Figure 3B). Upwelling-favorable winds in the CCS often did not exceed 2 days, and in the CnCS the average duration was 14 days (Figure 3C). Results of an ANOVA analysis showed a significant difference in the duration of upwellingfavorable winds between the EBUS regions (F = 431.29, SS = 16,113, and p < 0.005) ( Table 2). A significantly negative monthly trends in wind duration during July months in the CCS North (R 2 = 0.004, slope = −0.2 days dec −1 , and p < 0.05) and during January in the HCS Peru (R 2 = 0.008, slope = −1.1 days dec −1 , and p < 0.05) was detected. Upwellingfavorable winds appeared to be most intense in the CCS North, often exceeding 7 ms −1 (Figure 3C). Winds were least intense in the CCS South with a speed of approximately 4 ms −1 being commonplace ( Figure 3C). There were no significant changes in the overall wind intensity over time within the EBUS regions ( Table 2). A regression analysis comparing wind intensity per month, however, showed a significant change during July in the CnCS (R 2 = 0.01, slope = 0.10 ms −1 dec −1 , and p < 0.05). We also found no significant change in the wind intensity of EBUS over time ( Table 2). Bakun (1990) hypothesized that an increase in greenhouse gases will result in considerable changes in land-sea pressure gradients that will affect global wind patterns and ultimately result in an increase in the intensity of upwelling across the world's oceans. We tested this hypothesis by analyzing upwelling trends at four prominent coastal upwelling regions, with an emphasis on austral (DJF) and boreal (JJA) summer months as this is   when upwelling is most prevalent. By combining wind and SST data we were able to observe trends in upwelling responses. Specifically, change in upwelling metrics-as established by our novel methodology-was ubiquitous across all four EBUS. Trends in the metrics are coupled such that we see a decrease (increase) in SST combined with a concurrent increase (decrease) in the number of upwelling "events" and an increase (decrease) in the mean and cumulative intensity of upwelling. However, the change was inconsistent across the upwelling regions. In terms of the SST response, the Humboldt Current off Peru and Chile and the southern Benguela Current displayed decreases of ∼0.37 • C and 0.44 • C during past 37 years, while in the northern Benguela Current region it increased by ∼0.52 • C during this time. No significant changes were observed for the Canary and California Currents in the northern hemisphere. Surprisingly, it was interesting that anticipated changes in the decadal trend in mean intensity, duration, and count of upwelling-favorable winds were generally not detected. The significant cooling trend of the southern region of the Benguela Current is in agreement with previous research (Lima and Wethey, 2012;Santos et al., 2012b) that suggested that these decreases were related to increased upwelling. Further, at the southern end of the Benguela Current the highest upwelling intensities were observed during summer seasons (Narayan et al., 2010;Patti et al., 2010). Other studies (for example Cropper et al., 2014;Benazzouz et al., 2015;Varela et al., 2015Varela et al., , 2020Santos et al., 2016) demonstrated that a cooler SST associated with an increase in upwelling in coastal areas, as was also observed in Chile and Peru in this study. Previous studies demonstrated that trends in upwelling along the Peru coastline were variable. Here we report no significant change in upwelling-favorable winds along the Peru coastline, which contrasts with the studies reported by Gutiérrez et al. (2011) and Varela et al. (2015). The California Current displayed a positive trend in the mean intensity of upwelling signals over time, which is in disagreement with the study conducted by Varela et al. (2015). Here we also report no significant trend in wind speed, which contrasts with the studies reported by Mendelssohn and Schwing (2002), Narayan et al. (2010) and García-Reyes et al. (2014). The Canary Current displayed increases over time in metrics indicative of upwelling-favorable wind (Narayan et al., 2010;Patti et al., 2010), but apparently without an associated change in upwelling. Our findings here, and the findings of the aforementioned studies, suggest that warming rates are generally depressed along the coast within upwelling regions (Santos et al., 2012b,a), but we show that this seems to be the case for the southern hemisphere upwelling systems only (excluding the northern Benguela Current where it is increasing).

DISCUSSION
Using our novel method of determining upwelling signals, our results do not support the hypothesis that intensified upwelling will result from the increased land-sea temperature difference associated with climate change. However, it is possible that natural multi-decadal scale climate variability impacted the trends discussed. For example, the occurrence of ENSO at the end of the time series could also have initiated an anomalously warm SST in EBUS and potentially affected the trends observed. These results also lead to discussion about the potential connection between ENSO within EBUS and the Benguela Niño in the Benguela Current (Peterson and Schwing, 2003;Blamey et al., 2012). ENSO represents a weakening of the Walker cell circulation (Wang, 2004). During normal Walker cell conditions there is consistent upwelling, and this upwelling contributes to cool SST (Bakun et al., 2010). Since the Canary Current shows an increase in SST, it may suggest that there has been a weakening or reversal of the trade winds as has been indicated by earlier studies (Vallis, 1986;Saji et al., 1999). ENSO drives atmospheric and oceanic Rossby waves that can influence climate processes (Battisti, 1989;Holbrook et al., 2011). The ENSO's counterpart, La Niña, produces a teleconnection that alters the mean climatic states (Diaz et al., 2001;Fogt and Bromwich, 2006;Yeh et al., 2018). As a result, it can also be reasoned that the fewer upwelling signals are detected in Chile and Peru is potentially a consequence of the effects of La Niña. During La Niña, the Walker cell circulation is enhanced (Sohn et al., 2013), meaning that pre-existing trade winds and SST strengthen. The low SSTs created by the presence of La Niña are accompanied by anticyclones. The anticyclones rotate counterclockwise in the southern hemisphere, and therefore its winds travel equatorward where they move surface water away from the shore, resulting in a net Ekman veering that encourages coastal upwelling and lowers SST (Alford, 2001). Given these results, the need to consider changes in both thermal and in-shore hydrodynamic processes such as wind driven currents when interpreting the full dynamical response such as in conditions of changing wind magnitude and direction of the coupled atmospheric-ocean system to climate warming is important. In this regard, our conclusion differs in respect to recent studies (e.g., Wang et al., 2015) that have examined the relationship between increasing summertime landsea temperature differences and intensified upwelling in response to climate change. Rykaczewski et al. (2015) recommended that expanded land-sea temperature contrasts are omnipresent in projections of future conditions, but that summertime upwelling is limited to the polar extremes of upwelling zones. If true, this would indicate that increased land-sea temperature differences do not have a dominant influence on upwelling intensity (Tim et al., 2015(Tim et al., , 2016. Additionally, changes in upwelling-favorable winds are not always directly related to broad increases in land-sea temperature differences associated with climate change . Research suggests that changes in sea level pressure (SLP) fields are expected in response to increased greenhouse gas concentration (Gillett et al., 2013). These shifts will initiate changes in the magnitude, location and timing of upwelling-favorable winds that are more consequential than increased land-sea temperature differences.
Simulations of ocean patterns in the Benguela Current using high resolution data over the past few decades have thus far failed to detect predicted intensification of upwelling. For instance, Tim et al. (2015) demonstrate that upwelling has not dramatically intensified over at least the past 50 years. Moreover, broader investigations of global climate changes have also failed to detect meaningful intensification patterns of upwelling and in addition have not forecasted that such intensification will occur in the future (Wang et al., 2015;Tim et al., 2016). One potential reason for our results contrasting with the predictions by Bakun (1990) could be that the differences in quality of wind data employed were too severe. Bakun (1990) used Wave and Anemometerbased Sea surface Wind (WASWind) data whereas we made use of ERA5 climate reanalysis. Both datasets are biased to some degree. WASWind data are based on ship-based measurements and were fund to produce an artificial upward trend in sea surface wind speeds due to variable increases in the heights of anemometers (Tokinaga and Xie, 2011). These data have since been corrected, and the trends expressed by Bakun (1990) are no longer applicable with the corrected wind data, at least in some regions such as Peru (Belmadani et al., 2014). ERA5 climate reanalysis data are also susceptible to several biases (see Astudillo et al., 2017;Taboada et al., 2019) but are considerably more accurate and reliable given the spatial and temporal resolution at which winds are measured (Graham et al., 2019;Mayer et al., 2019;Tetzner et al., 2019). For comparison, ERA5 products are provided daily at a 0.25 • × 0.25 • latitude-longitude and are available from 1979 to present, whereas WASWind data are available at monthly intervals at 4 • × 4 • latitude-longitude resolution and cover the years from 1950-2008. One potential caveat of our analysis and the interpretations thereof lays in the limited nature of the time series we assessed. Because we sought to detect and quantify patterns of upwelling at a fine spatio-temporal scale, our choice in datasets was limited to highresolution data that only cover the past 37 years. In the context of broader climate change, 37 years is likely too short a duration of time from which to inform long-term shifts in patterns of global phenomena like upwelling and so our findings here should be regarded with some caution. While we cannot provide evidence in support of Bakun's hypothesis regarding intensification of upwelling (Bakun, 1990), we cannot rule out the possibility that an analysis with longer time series could produce a different outcome. Unfortunately, high-resolution SST datasets that cover longer time periods do not exist. Some longer reanalysis products like ERA20C or 20CR provide data over a time series of more than 100 years but these data products are only available at a considerably coarser resolution and may not be suitably sufficient to reliably estimate upwelling.
There is not yet an adequate amount of data to effectively investigate the effects of changes in climatic conditions on coastal SSTs and biogeochemistry in the EBUS . Despite this, consideration is now given to changes in terms of the location and intensity of upwelling due to its necessity for the identification of high-risk regions of the coast that are prone to processes such as ocean acidification, increased hypoxia, and eutrophication which have been projected in accordance with future warming . These patterns are affected by both complex global and local processes; however, the availability of the data and their corresponding resolution in conjunction with a wide range of decadal temperature variability and biogeochemical properties also play a significant role. Most of the remotely sensed SST and wind products do not have sufficient resolution to accurately detect nuances within the upwelling process. Despite this, data providers have been yielding promising advancements in higher resolution products by using land and air-sea interactions, cloud formation, and oceanic mesoscale processes to correct the biases presented in older sources. It has since become necessary to take into consideration changes in local, alongshore winds when investigating the relationship between sensitivity of upwelling and climate change (Bakun et al., 2015). Environmental variables that drive upwelling like temperature and wind measures are highly variable in EBUS regions, both spatially and temporally. Changes in these variables are shown here to directly affect upwelling patterns. As a result, major shifts in these drivers can greatly influence the frequencies, intensities, and the durations of upwelling events. Increases and decreases in SST are negatively associated with corresponding shifts in detected upwelling. Similarly, heterogeneous winds can dramatically influence upwelling responses. Accordingly, it is vital to take into consideration both the geostrophic and ageostrophic processes that affect the boundary layer of the coastal ocean when investigating upwelling. By analyzing these coupled processes, researchers can further elucidate upon the drivers that regulate changes in the characteristics of upwelling events. This is especially true since changes in water-column stratification can distort upwelling signals (Roemmich and McGowan, 1995;Chhak and Di Lorenzo, 2007) and so external processes like wind direction are important to consider. However, we found that the dynamics of upwelling-favorable winds have not changed significantly over the past 37 years. Further analyses may be required to fully interrogate the effects of shifts in climate on upwelling dynamics in EBUS.
In this study we produce metrics that track the duration and intensity of upwelling signals over the past 37 years. These metrics are based in part on those used to quantify marine heatwaves (Hobday et al., 2018;Schlegel and Smit, 2018) and they offer a valuable new approach that can be used to determine the changing dynamics of upwelling consistently and objectively in daily time series of SST and winds. While all four EBUS experienced some change in upwelling metrics, only those in the southern hemisphere generally responded with a decrease in SST (the northern Benguela Current region trended upward). Although decreases in SST were associated with corresponding changes in upwelling intensity, a trend in wind dynamics apparently did not explain the upwelling dynamical response. We thus follow in the long tradition of not being able to conclusively support Bakun (1990) hypothesis. Continuing shifts in climate may indeed facilitate increased upwelling in certain coastal areas and thus buffer ecosystems from climate change (Bakun, 1990, Bakun et al., 2010García-Reyes et al., 2015), but it seems that this might not be seen consistently across all EBUS. The degree to which all EBUS can be characterized as resilient and robust to natural climate variability remains uncertain, as is the universal benefit of upwelling intensification to maintaining an abundance of species of high commercial and conservation value (Bell et al., 2015;Wilson and Forsyth, 2018). However, given the coarse resolution of current climate models for ocean variables, it is difficult to reproduce the relatively fine-scale upwelling features necessary for such research at present . Continuing development of geophysical time series in the future will likely yield data series that are better able to yield trends indicative of climatic change. The knowledge of future changes in upwelling systems is a key factor for estimating changes in economic and biological trends. The metrics of intensity and duration of upwelling (and of frequency of upwelling "events" that can be also established using our approach) will more than likely also relate to a diversity of measures of ecological functioning and well-being and provide researchers with another tool to track environmental drivers and responses.

DATA AVAILABILITY STATEMENT
The data and analyses used in this paper may be downloaded at: https://github.com/AmierohAbrahams/EBUS.

AUTHOR CONTRIBUTIONS
AS conceptualized the scope of the research reported here. AA undertook the statistical analyses, made the first round of interpretation, and did the majority of the writing. AS and RS provided guidance in terms of writing, statistical analyses, and the final draft of the manuscript. All authors contributed to the article and approved the submitted version.