The impacts on marine species from secular warming and heatwaves are well demonstrated, however the impacts of extreme cold events are poorly understood. Here, we link the death of organisms from 81 species to an intense cold upwelling event in the Agulhas Current, and show trends of increasing frequency and intensification of upwelling in the Agulhas Current and East Australian Current. Using electronic tagging, we illustrate potential impacts of upwelling events on the movement behaviour of bull sharks, Carcharhinus leucas including alterations of migratory patterns and maintenance of shallower dive profiles when transiting through upwelling cells. Increasing upwelling could result in “bait-and-switch” situations, where climate change expands subtropical species’ distribution, while simultaneously exposing climate-migrants to increased risk of cold-stun/mortality events at poleward distributional limits. This shows the potential impacts of increased cold events, an understudied aspect of climate-change research, and highlights the complexities of climate change effects on marine ecosystems.

Study areas

Both study areas are located along Western Boundary Currents (WBCs): the Agulhas Current in southern Africa and the East Australian Current in eastern Australia. In southern Africa, the area spans from the warm-temperate Breede River estuary on the south coast of South Africa (34.5°S) to northern Mozambique (16°S) (Figs. 1, 5a main MS). In eastern Australia, the area extends from the warm-temperate Sydney Harbour estuary (33.5˚S) to tropical Queensland (16.5˚S) (Fig. 5b main MS). On both continents, the study areas were divided into two distinctive zones (Fig. 5a, b main MS). First, an upwelling zone, comprising multiple coastal upwelling cells inshore of each WBC, extending from the southern limit equatorward to 32.4°S along the Agulhas Current and 28.6°S along the East Australian Current. This zone includes the Breede River and Sydney Harbour estuaries in South Africa and Australia, respectively (Fig. 5a, b main MS). The second zone, in each case, was the subtropical/tropical section of coastline equatorward of the upwelling zones.

Both the Agulhas and East Australian Current closely follow the continental shelf, with their warm cores usually situated between 50-100 km offshore from the shelf edge16, 56. Note that although upwelling can occur in the subtropical/tropical areas to the north, it is most prominent and intense in both upwelling zones16, 30.

In the case of the Agulhas Current, coastal upwelling can be caused by current forcing as well as by meanders in the current^{16, 19, 31}. Such meanders can be associated with cold-core inshore eddies that travel south-ward along the shelf, often called “Natal Pulses” ¹⁷. Current-driven coastal upwelling is most common along the south-east coast of South Africa, due the proximity of the shelf edge to the coast^{16, 19, 31}. Upwelling-favourable easterly winds can also contribute to coastal upwelling here¹⁶. As the shelf widens south of Port Elizabeth, this becomes the main mechanism of coastal upwelling, with a weakening influence of the current. Similarly, along the East Australia Current, current and eddy-driven upwelling mostly occurs along the narrow shelf from 24˚S to 32.5˚S, after which the current separates and wind-driven upwelling becomes the most prominent mechanism^{31, 32, 56}.

Collection of environmental variables

On March 2nd 2021, dead marine organisms began washing up on beaches inshore of the Agulhas Current,  between Port Elizabeth to ca. 100 km north of East London (Fig. 1a main MS). Specimens were collected and cataloged by Greg Brett in association with the East London Museum (Supplemental Table 2). To determine whether rapid-temperature drops associated with upwelling could be responsible for the deaths and, if so, to qualify characteristics of lethal upwelling events, we investigated sea-surface temperature (SST) from satellite imagery, in-situ temperatures at 30 m depth and wind speed and direction.

We used the Aqua Modis sea-surface temperature data set from the Moderate Resolution Spectroradiometer on the AQUA NASA satellite, a daily, gridded level-3 data set at 0.04166° resolution, with data availability from 2003-2022. To analyse the duration as well as maximum spatial extent and rate of temperature declines for potential upwelling before, during and after the dead specimens washed up on local beaches, we generated mapped images in ArcGIS based on NetCDF-files downloaded from the OceanColor webserver (https://oceancolor.gsfc.nasa.gov/l3/). Satellite-derived temperature data was complemented with in-situ temperatures collected by a temperature logger deployed 3 km offshore from Port Alfred (Fig. 1a main MS) at 30 m,, which recorded temperature at hourly intervals.

Finally, wind data collected by the South African Weather Service at Port Alfred was used to determine whether upwelling-favourable easterly winds (45°-135°) occurred before and during the period dead marine organisms washed up on local beaches. Wind speed in m/s and wind direction in degrees was collected at 6-hour intervals daily between 8 am and 8 pm local time.

Modelling trends in upwelling characteristics over time

To investigate whether there has been an increasing trend in the proportion of upwelling-favourable winds on the south and south-east coast of South Africa, we used wind data collected by the South African Weather Service from three stations within the upwelling zone. Note that a similar analysis has previously been done for eastern Australia³⁵. Daily wind speed in m/s and wind direction in degrees was collected at 6-hour intervals between 8 am and 8 pm local time. Data was available from November 2001 until March 2022 at Port Alfred, from October 1992 until March 2022 at Port Elizabeth, and from October 1988 until March 2022 at Plettenberg Bay.

Based on the orientation of the coastline in the upwelling zone and wind direction in degrees, the data set was divided into upwelling-favourable winds (north-east to south-east (45°-135°) and non-upwelling favourable winds (all other directions)⁵⁷. To estimate the daily prevailing wind speed and direction, data was averaged per day.

Since we aimed to investigate potential impacts of increasing trends in upwelling on mobile marine species at their distributional limit, only wind data between October and April was included in the models. Based on bull shark tracking data from this study, this was the period during which animals tagged in the Breede River and Sydney Harbour estuaries were mainly present in the upwelling zone. Furthermore, this period also encompasses the austral summer months which are the main period of upwelling activity inshore of the Agulhas and East Australian Currents^{16, 19, 30, 34}. Thus, annual data from October to April were combined to represent corresponding upwelling seasons. As a result, we had data for 21 upwelling seasons at Port Alfred, 30 at Port Elizabeth and 34 at Plettenberg Bay. Although timeseries at Port Alfred were shorter than 30 years, which is the recommended length for trend analysis²¹, we included this site to inspect a potential trend and compare it with Port Elizabeth and Plettenberg Bay sites.

We then calculated the proportion of days during which upwelling-favourable winds occurred in each upwelling season. The resulting values were modelled against time (upwelling season) using generalised linear models (GLM) with binomial error distribution. Model fit and diagnostics, including tests for temporal autocorrelation were evaluated using the Dharma package in R. Only trends in the South African study areas were modelled as those in the east Australian coastline has previously been examined³⁵.

In order to investigate whether the number, mean intensity (in °C) and the mean rate of onset (in °C/day) of potential upwelling events has increased in upwelling cells inshore of both WBCs, we used the long-term, interpolated, reanalysed, daily Global Ocean OSTIA Sea Surface Temperature data set⁵⁸. This gridded, level-4 data set covers the period from 10/01/1981 until present at a resolution of 0.05°, combining satellite data as well as in-situ data. We chose this data set because of long-term data availability of 41 years, the high standard of quality control and overall accuracy and applicability⁵⁹. Although level-4 data sets can underestimate thermal signatures of upwelling zones due to warm biases, the OSTIA data set has been shown to perform better than other long-term data sets³⁶. Furthermore, this bias is more pronounced in Eastern Boundary Current compared to Western Boundary Current upwelling systems^{33, 59}.

Next, we extracted sea-surface temperature data for prominent upwelling cells inshore of each WBC. This was off Port Alfred, Port Elizabeth and Plettenberg Bay for the Agulhas Current and Port Stephens, Diamond Head and Ballina for the East Australian Current^{16, 19, 30, 34, 37, 47}. We then divided the shelf into inshore (0-15 km from shore) and midshelf (15-30 km from shore) sections. For each station, a pixel was randomly selected from inshore and midshelf sections from the gridded OSTIA SST data set along a transect (Fig. 7 main MS). Each transect extended from the coastline, across the shelf in a direction perpendicular to the shelf edge.

Choosing inshore and midshelf pixel allowed us to compare trends across areas of varying shelf width where each WBC would influence upwelling differentially. Additionally, as a control, we selected an offshore pixel within the warm core of each current along the same transects for each station. No upwelling occurs within the offshore core of the currents and possible cold events would be due to temperature anomalies other than upwelling, hence no trends should be evident.

We used the R package heatwaveR to detect cold events in each randomly selected pixel over time⁶⁰. The algorithm implemented in this package is based on the hierarchical definition of extreme temperature events in21 and follows a similar methodology to characterising such events. In short, marine extreme temperature events are defined as discrete, anomalously warm or cold events with characteristics such as intensity, rate of onset, and spatial extent²¹. For extreme cold events, these characteristics are calculated in cases when temperatures fall below the 10th percentile threshold, which is based on a baseline climatology at least 30 years in duration to avoid impacts of decadal variability. Thus, data sets with long-term data availability, such as the OSTIA SST data set, available from 01/01/1981-present, are preferred²¹. Before calculating the baseline climatology at each pixel, we detrended the daily OSTIA sea surface temperature data as long-term trends in the climatology can affect detection of extreme temperature events using the detrend() function61.

Within the detect_event() function of the heatwaveR package, we set the minimum duration for cold events to one day, in order to include short-term upwelling events which are common along both WBCs57. From the output of this function, we summed all discrete cold events each year from October to April, the time period tagged sharks were mainly present in both upwelling zones. We also calculated the annual mean intensity (in °C) and mean rate of onset (°C/day) for events from October to April. Visual inspection of the output suggested that trends in either variable remained similar using detrended and non-detrended data.

We used generalised linear models to investigate trends in the number, mean intensity and mean rate of onset of upwelling events at each pixel between 1981 and 2022, with upwelling season (aggregated annually between October and April) as the independent variable. The number of upwelling events was modelled assuming a Poisson error distribution with a log-link function, while mean intensity and mean rate of onset were modelled assuming a Gaussian error distribution. Model fit and diagnostics were inspected using the Dharma package. When significant non-linearity and/or quantile deviations were detected, we fitted generalised additive models (GAMS) (Supplemental Tables 6 and 7), which were also inspected using Dharma and the gam.check() function in the mgcv package.

Finally, as decadal variability and biases associated with L4-satellite SST products can confound trends in upwelling signals, we utilised temperature logger data from four locations within the upwelling cells of the Agulhas current (deployments from 1991 to 2022) and three different satellite SST products to validate our analyses above. Loggers were Star OddiTM temperature loggers with accuracy of ± 0.05°C located at Ystervarkpunt (78 km east of the Breede River), Mossel Bay (127 km east of the Breede River), Tsitsikamma (37 km east of Plettenberg Bay) and Mosterts Hoek (85 km west of Port Elizabeth) and serviced every 1-2 years (for deployment details see Supplemental Table 5). Satellite SST products used were the L4 OSTIA SST dataset from our original analyses (5 km resolution), the L4 CCI SST product (5 km resolution) and the L3 Pathfinder AVHRR night and day products combined into a single product (5 km resolution). Warming/cooling trends in SST products and logger data were estimated using the Theil-Sen estimator. The estimator is nonparametric and hence more robust than the parametric least-squares regression as it is resistant to outliers and does not draw from any probability distribution. Theil-Sen uses the median to estimate the slope as opposed to least-squares which makes use of a weighted mean for slope estimation. Spatial SST trends were calculated using each grid point within the designated area of the gridded SST product. Additionally, the SST and logger data trends were calculated using the anomalies to avoid the complicating effects of seasonality on trend values. Trend analyses for all logger locations and SST products ranged from June 1991 to August 2022.

Animal tagging and tracking

Tagged bull sharks were part of previously published data sets from southern Africa and eastern Australia and were acoustically tracked in coastal receiver arrays. For details about tagging methods and receiver arrays see27, 39, 62 , Supplemental Table 8). In short, in Australia, 25 bull sharks were internally tagged in Sydney Harbour with InnovaSeaTM V16TP tags that also recorded depth and temperature at the time of detection (tag nominal delay 30-90 s). Out of 41 bull sharks acoustically tagged in southern Africa as part of a separate study 27, 62, 15 individuals were also equipped with temperature and depth recording pop-up satellite archival transmitters (Wildlife ComputersTM Mk10-PSAT and miniPAT tags) and included again in this study (Supplemental Table 8). Five of these sharks were tagged inside the warm-temperate Breede River Estuary, and 10 at the subtropical Pinnacle Reef. PSAT-tags archived depth, temperature, and light level data at pre-programmed intervals of three seconds (n=1), five seconds (n=11) and 25 seconds (n=3). Additionally, binned data at 24-hour intervals was collected for temperature (9-11°C, 11-13°C, 13-15°C, 15-17°C, 17-19°C, 19-21°C, 21-23°C, 23-25°C, 25-27°C, 27-29°C, 29-31°C, 31-33°C) and depth (0-5m, 5-10m, 10-20m, 20-30m, 30-40m, 40-50m, 50-100m, 100-200m, 200-300m, 300-400m) for all PSAT-tags. Tags were released from animals after pre-programmed deployment durations (90-360 days). Tags then floated to the surface and transmitted data to the Argos satellites. Out of the 15 PSAT-tags deployed, 10 were physically recovered by the authors via boat, providing full archived temperature and depth data (one from a bull shark tagged at the Breede River and 9 from Pinnacle Reef). Transmitted data was initially decoded using the manufacturer’s software (WC-DAP 3.0, Wildlife Computers, Inc., Redmond, WA). The Wildlife ComputersTM GPE3 State-Space model was used to generate estimates of most likely locations from acoustic detections and PSAT-tag light level data using a swimming speed of 0.7 m/s64 and to create movement tracks for each shark (Supplemental Figs. 19-23).

Then, in the statistical software R63 each shark location was matched with depth and temperature data collected by the PSAT-tags to presence in the two regions of interest: The upwelling zone (distinguishing locations inside the Breede River and outside), and the subtropical/tropical areas past the upwelling zone (Supplemental Figs. 19-23). For example, if an acoustic detection or geolocation estimate occurred in the upwelling zone (Fig. 5 main MS), the corresponding depth and temperature value at this timestamp was assigned to the upwelling zone. Note, this was only done for sharks tagged in the Breede River, as animals tagged at Pinnacle Reef remained  in subtropical/tropical areas with limited upwelling. Similarly, for sharks acoustically tagged in Sydney Harbour, detections at receivers were assigned to either being in the upwelling zone, (distinguishing detections inside Sydney Harbour and outside), or past the upwelling zone (subtropical-tropical Queensland) (Fig. 5b main MS). 

Studies from multiple geographical locations have shown that bull shark abundance declines significantly when seasonal temperatures fall below a long-term thermal affinity threshold of 19°C38, 39, 40, 41. To further test this threshold and investigate thermal preferences for bull sharks, we established a tentative thermal range, using bull shark occurrence records in southern Africa and eastern Australia from the Global Biodiversty Information Facility (GBIF)65in combination with our acoustic detections. For this, we filtered combined bull shark occurrences to daily, verifed human and machine observations and spatio-temporally matched each occurrence with a pixel from the AQUA Modis L3-satellite temperature data set. Similar analysis using occurrence data was also performed for three megafauna species that have had fatalities in severe upwelling events, the whale shark Rhincodon typus manta rays, Mobula birostris and Mobula alfredi.

For South African sharks for which PSATs were retrieved and full archival data was available (Supplemental Table 15), we plotted continuous temperature and depth over time, while indicating the occupied zone. Binned data from all Breede River bull sharks (full archival and partial data) was used to calculate percent of time at depth and temperatures to compare this in both upwelling and non-upwelling zones Pinnacle Reef sharks were detected only outside of upwelling cells). The maximum duration in consecutive hours, for which the Breede sharks were experiencing temperatures below 19°C for each zone was also calculated. Similarly, we binned temperature and depth data recorded from acoustic detections of bull sharks tagged in Sydney Harbour to determine percent of time at depth/temperature. Acoustic detections at depth/temperature were also plotted over time for sharks tagged in Sydney Harbour and indicated in which zone they occurred. Lastly, data from temperature loggers deployed outside the Breede estuary at 20 m and off Port Alfred at 30 m, as well as the loggers used for Theil-Sen estimation (see previous section)  were used to spatio-temporally match temperature and depth recorded by PSAT-tags of the Breede River sharks with upwelling events at nearby loggers based on shark acoustic detections/PSAT-geolocations. This was plotted in combination with a detection timeline from combined acoustic and PSAT-tag locations to indicate when and where PSAT-tagged sharks encountered active upwelling cells. Temperature data for outside Sydney Harbour was recorded by a mooring located 14 km offshore, between 12-20 m depth, which was used to plot upwelling events with acoustic shark detections.

Method References

56.          Archer MR, Roughan M, Keating SR, Schaeffer A. On the variability of the East Australian Current: Jet structure, meandering, and influence on shelf circulation. Journal of Geophysical Research: Oceans 2017, 122(11): 8464-8481.

57.         Abrahams A, Schlegel RW, Smit AJ. A novel approach to quantify metrics of upwelling intensity, frequency, and duration. PLoS One 2021, 16(7): e0254026.

58.          Good S, Fiedler E, Mao C, Martin MJ, Maycock A, Reid R, et al. The current configuration of the OSTIA system for operational production of foundation sea surface temperature and ice concentration analyses. Remote Sensing 2020, 12(4): 720.

59.          Yang C, Leonelli FE, Marullo S, Artale V, Beggs H, Nardelli BB, et al. Sea surface temperature intercomparison in the framework of the Copernicus Climate Change Service (C3S). Journal of Climate 2021, 34(13): 5257-5283.

60.          Schlegel RW, Smit AJ. heatwaveR: A central algorithm for the detection of heatwaves and cold-spells. Journal of Open Source Software 2018, 3(27): 821.

61.          Schlegel RW, Oliver EC, Hobday AJ, Smit AJ. Detecting marine heatwaves with sub-optimal data. Frontiers in Marine Science 2019, 6: 737.

62.          Daly R, Smale MJ, Cowley PD, Froneman PW. Residency patterns and migration dynamics of adult bull sharks (Carcharhinus leucas) on the east coast of southern Africa. PloS one 2014, 9(10): e109357.

63.        R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/

64.          Brunnschweiler JM. Sharksucker–shark interaction in two carcharhinid species. Marine Ecology 2006, 27(1): 89-94.

65.         Stuart-Smith, R. D., Edgar, G. J., Barrett, N. S., Kininmonth, S. J., & Bates, A. E. (2015). Thermal biases and vulnerability to warming in the world’s marine fauna. Nature, 528(7580), 88-92.