Model evaluation

We compare the CESM2(WACCM6) simulated SST for the historical period (1982–2014) with OISST and ERA5 results (Fig. 1). The simulated, observed, and reanalysis SST show similar patterns, with the highest temperatures occurring in the region between the AC and ACR. The CESM2(WACCM6) SST distribution has a high correlation (r = 0.87) with both OISST and ERA5. However, the CESM2(WACCM6) SST pattern fails to reproduce the meandering pathway of the ACR, which is bounded by the 19 °C isotherm (Fig. 1a, b and e). The model also overestimates the SST north of 37°S, especially in the western coast, while underestimating it south of 37°S, particularly in the vicinity of the ACR (Fig. 1c and f). These biases occur because the model simulates the transport of Agulhas warm water too far northwestward instead of retroreflecting it southeastward. To know how well the model captures the temporal variability of the SST, we also compare the normalized standard deviation (NSD) in CESM2(WACCM6), OISST and ERA5 (Fig. 1d and g). The simulated NSD falls within 0.8 and 1.2 over most parts of the region, except over a small area offshore of the Western Coast where CESM2(WACCM6) grossly overestimates the SST variability. While the discrepancies between CESM6(WACCM6) and the observations (OISST and ERA5) in the mean SST and its variability may be due to biases in CESM2(WACCM6), they may also result from shortcomings in the satellite data caused by cloud cover23. However, the model’s coarse resolution may result in a flattened bathymetric gradient, particularly in the AB region, which may lead to an underestimation of local upwelling and mixing processes. This limitation also explains the model’s inability to reproduce mesoscale features such as AC rings and eddies, potentially leading to a systematic underestimation of small-scale thermal extremes and, consequently, MHW metrics.

We also compare the characteristics of MHWs (i.e., number, duration, and intensity) in the CESM2(WACCM6) historical simulation with the OISST dataset (Fig. 2) and ERA5 dataset (Fig. S1). The simulation gives realistic patterns of MHWs characteristics (number, duration, and intensity) with reference to OISST observations in the period 1982–2014 (Fig. 2a, d and g) and 1982–201548. The simulation agrees with observations that the most frequent and most transient MHWs in the region occur south and southwest of the AB, although the simulated MHWs are less frequent and more persistent than observed. The peak activities in the area might be attributed to the southward moving AC, which is known to transport warm subtropical water masses from the Indo-Pacific Ocean to the Atlantic via the Agulhas Leakage4. CESM2(WACCM6) and the observations agree that the most persistent MHWs activities (> 14 days duration) reside at the north-western part of the region, although CESM2(WACCM6) generally overestimates the duration by more than 8 days over this area (Fig. 2d, e, f, and Fig. S1d, e, f). CESM2(WACCM6) also shows some persistent events in the southeastern basin, which are slightly reproduced in OISST compared to ERA5, and where the duration is also overestimated. The major discrepancy between CESM2(WACCM6) and OISST is that CESM2(WACCM6) generally underestimates the intensity of MHWs with maximum bias (> 2.0 °C) residing south of 39°S (Fig. 2i and Fig. S1i). The overestimation of MHWs duration to the northwest of the study domain and the underestimation of MHWs intensity south of 39°S may be attributed to the diversion of the warm AC from south to east in the simulation (Fig. 1c). This suggests that the model overestimates the Agulhas leakage. These biases are consistent with the model cold bias south of 39°S and warm bias north of the latitude (Fig. 1c). Nevertheless, the model bias over the AB region is relatively small. Over this region, the model biases are about 1 event/day for the MHWs frequency, 5 days for the MHWs duration and 1 °C for the intensity. To highlight the differences between CESM2(WACCM6) and OISST (Fig. 2) as well as between CESM2(WACCM6) and ERA5 (Fig. S1 in the supplementary material), we analyzed the time series of MHW event evolution at various locations over the period 2010–2014 (Fig. 3, and Figs. S2 and S3 in the supplementary material). Figure 3 focused on the two contrasting locations (indicated by the green stars in Fig. 2f and i) which exhibits opposite MHW metric patterns, one with fewer MHW events, longer durations, and lower intensity (14.6°E, 33.9°S), and another with more frequent events, shorter durations, and higher intensity (16.9°E, 41.4°S). Over the five-year period, OISST, ERA5, and CESM2(WACCM6) recorded 9, 7, and 4 MHW events, respectively, at 14.6°E, 33.9°S. The corresponding mean intensities were approximately 1 °C for both OISST and ERA5, and ~ 0.5 °C for CESM2(WACCM6), with mean durations of ~ 10 days for OISST and ERA5, and ~ 35 days for CESM2(WACCM6). At 16.9°E, 41.4°S, the same datasets showed 13, 10, and 9 MHW events, respectively, with mean intensities of ~ 3.2 °C, ~ 3.1 °C, and ~ 1 °C, and mean durations of ~ 14 days, ~ 13 days, and ~ 15 days, respectively. Figures S2 and S3 present additional examples from randomly selected locations. The results reveal consistent discrepancies between the model and the observations, confirming that MHWs are generally more frequent, intense, but shorter in the observations across the study area. These differences may be attributed to the coarse resolution of the CESM2(WACCM6) model and its smoother output compared to observational data. Such model biases could affect future projections by underestimating the frequency and intensity of MHWs, while overestimating their duration within the study area.

Fig. 2

figure 2

Comparison of marine heatwave (MHWs) metrics from CESM2(WACCM6) historical simulation and OISST over the period 1982–2014: (a–c), annual mean MHWs frequency; (d–f), annual mean MHWs duration; and (g–j), annual mean MHWs intensity. The first and second columns represent the MHWs metrics from OISST and CESM2(WACCM6) historical simulation, respectively; and the third column is the difference between CESM2(WACCM6) historical simulation and OISST (CESM2(WACCM6) minus OISST). The green stars on (f) and (i), respectively, are the locations chosen to show the evolution of the MHW events in Fig. 3. Contour lines in (h) are for the maximum intensity in CESM2(WACCM6) historical simulation (1 °C and 1.05 °C). The spatial correlation between OISST and CESM2(WACCM6) historical simulation (r) on (c), (f) and (i) are for the frequency, the duration and the intensity, respectively. All the correlations are significant.

Fig. 3

figure 3

Example of marine heatwave (MHWs) evolution in OISST (a and d), ERA5 (b and e) and CESM2(WACCM6) (c and f) at 14.6°E 33.9°S (left panels) and 16.9°E 41.4°S (right panels) over the 2010–2014 period. The thick blue and green lines indicate the climatological mean and 90th percentile threshold, respectively, and the shaded red colours represent the MHW events. The green stars on Fig. 2f and i indicate the locations considered for this analysis.

Projected changes in MHWs off the South African coast

Projected future changes as shown in Fig. 4, highlight MHWs characteristics (2050–2069) with respect to the present-day climate (2015–2034). Under the SSP2-4.5 scenario without SAI (climate change scenario) an increase in MHW frequency, duration and intensity is projected over the study region, with the maximum increase (up to 150%, 200%, and 15%, respectively) occurring along the band of the AC stretching to the Atlantic Ocean and the minimum increase occurring south of the band (Fig. 4f and j). This leads to a higher increase in cumulative intensity (Fig. S4 in the supplementary material), which reflects both the duration and intensity of MHWs, and can also serve as indicator of their potential impact on marine ecosystems36. The corresponding projected changes in the SST and ocean currents (Fig. 5b) show a significantly greater increase in SST near the coast and north of 37°S. It also shows that while the AC weakens (opposite to its usual direction), the westward component of the surface currents strengthens in the region of the highest SST increase, including along the shelf-edge AC. This is consistent with Jury (2020), who showed an accelerating shelf-edge AC with a zonal trend of − 0.006 m s⁻1 yr⁻1 over the 1980–2015 period at longitudes 21°E–28°E. It suggests that the band of maximum increase in MHW activities may be driven by the intrusion of warmer water from the Indian Ocean into the South Africa coasts and the Atlantic Ocean north of 37°S via the AC (i.e. enhanced Agulhas leakage) while the dynamic subtropical front (DSTF) exhibits minimal MHW activities. The DSTF region is characterized by a strong SST and sea surface height (SSH) gradient between 37°S-40°S17, along with a strong current gradient, deepest isotherms and strong mixing processes, which may help regulate extreme SSTs, despite the weakening wind strength (Figs. 5f and 6a, e, i). The projected changes in the vertical sections of the zonal and meridional current components across the DSTF (within 20°E-24°E, 35°S-42°S box) indicate a weakening of the AC (whose core is located around 37°S) and the ACR, further confirming the weakening of the surface currents shown in Fig. 5b (Fig. 6a, b, e and f). This is accompanied by a weakening of the net downward current below 150 m depth, likely combined with strong mixing within the water column, resulting in a minimum temperature increase between 37°S and 39°S, with the lowest variations (about 0.1 °C) occurring at 150 m depth around 38°S (Fig. 6i, j, k, l). This could explain the minimum increase in MHW activities in the DSTF region, despite the increase in net heat flux of the ocean and the weakening of westerly winds, consistent with Jury (2020), who found a trend toward easterly winds across our entire study area (Fig. 5f, i, j).

Fig. 4

figure 4

Spatial distribution of the marine heatwaves characteristics (annual mean frequency, duration and intensity) off the coast of South Africa in the present-day climate (PRS: 2015–2034; panels (a), (e) and (i)), along with projected future changes for the period 2050–2069 under the SSP2-4.5 scenario without and with stratospheric aerosol injection (SAI) (CC—PRS: SSP2-4.5 without SAI minus present-day climate, and SAI—PRS: SSP2-4.5 with SAI minus present-day climate, respectively). The impact of the SAI on global warming is also shown (SAI—CC: SSP2-4.5 with SAI minus SSP2-4.5 without SAI). Projected changes are expressed in terms of the rate of change. The gray lines on the panels represent the 200-m isobath. In panels (b), (c), (f), (g), (j), and (k), the dashed contours represent negative percentages, while the solid contours indicate positive percentages. The white boxes in panels (d), (h) and (l) represent the Agulhas bank region considered in this study.

Fig. 5

figure 5

Same as Fig. 4, but for SST (a–d), wind speed at 992 hPa (e–h), and net ocean surface heat flux (i–j). The arrows in (a), (b), (c) and (d) represent the surface current and its differences according to the corresponding scenarios (CC-PRS, SAI-PRS and SAI-CC), while the arrows in (e), (f), (g) and (h) represent the wind direction and its differences according to the same scenarios. The black box on panel (k) is the box used for Fig. 6.

Fig. 6

figure 6

Same as Fig. 4, but for the zonal component of the current (a–d), the meridional component of the current (e–h), and the temperature (i–l) of the zonal mean of the box 20°E-24°E, 35°S-42°S (black box on panel (k) of Fig. 5). The green contours on (e), (f), (g) and (h) are the values of the vertical component of the currents: full and dashed contours are in the upward and the downward direction, respectively, in centimeters per day (10−2 m d−1). The arrows on (i), (j), (k), (l) show direction of the meridional and vertical current (PRS, CC—PRS, SAI—PRS, SAI -CC), amplified here for reading. The white contours on (i) represent the 17 °C and 21 °C isotherms. AC and ACR are the Agulhas Current and retroflected branch, respectively.

The implementation of SAI mitigates the impacts of climate change on MHWs (Fig. 4). With respect to PRS, SAI reduces the increases in MHWs frequency to ± 20% over the entire study domain, the duration to 4h and l). The analysis of cumulative intensity reveals similar patterns to those of duration and intensity, as it represents the combined effect of both (Figure S4 in the supplementary material). The reduction occurs because the SAI apparently offsets the impacts of climate change on SST over the study domain. However, it overcompensates the impacts (by about 0.4 °C) in the DSTF region (Fig. 5c). The overcompensation arises because SAI partially reverses the impacts of climate change on westerly winds and current south of 39°S (i.e. winds become stronger and current become weaker with respect to CC: Fig. 5d and h). The cross-section of temperature between 20°E and 24°S also shows the overcompensation in depth between 37°S and 39°S (Fig. 6k). This is explained by the strong mixing taking place in the water column due to both the weakening of the zonal and meridional current, and the downward current below 100 m depth with respect to CC (Fig. 6c, d, g and h).

Impacts of climate change and SAI on the dominant patterns of persistent sea surface temperature anomalies associated marine heatwaves events

Figure 7 shows the SOM classification of SST anomaly patterns (for days with extreme events occurring in the Agulhas Bank area) in the three datasets (PRS, CC, and SAI) into nine major patterns (3 × 3 nodes). The nodes at edges feature the most extreme patterns while the inner nodes provide smooth transitions between the extreme patterns. A closer look at the SOM classification can be broadly divided into two groups (Figs. 7 and 8).

Fig. 7

figure 7

Dominant patterns of sea surface temperature anomalies associated with marine heatwaves (MHWs) events off the South African coast in the present-day climate (PRS) and future climate under the SSP2-4.5 scenario without and with SAI (CC and SAI, respectively). The relative frequency of each pattern (node) is shown in the top center of each pattern. The percentage contribution of each dataset is indicated. The nodes are classified in two groups: SWC for South and West Cooling and ACC for Agulhas Current Cooling.

Fig. 8

figure 8

Seasonal variability associated with each node of the 3 × 3 SOM in percentage.

The first group, the south and west cooling (SWC; Nodes 1, 4, 5, 7, and 8) associates MHW activities in the AB region with weak SST anomalies south of 39°S. This suggests that the driver of the MHW activities may be confined to the local dynamics of the AC (Fig. 7). These include the meandering of the AC entering the AB region and the leakage of its warm water into the Atlantic Ocean5,12,26. In this group, the SST anomalies associated with MHW activities are weaker in the AB region and might be due to the northward displacement of cold water from the DSTF region towards the AB which weakens the AC (Fig. 6b, c, f and g). Because of the direction of winds, especially in the SAI scenario, the cooling can be also due to the Ekman drift of cold water from south and west towards the AB region (Fig. 5g and h). The group accommodates the majority of SST anomaly patterns in the three climate datasets (PRS: 81%; CC: 54%; SAI: 78%) but with the lowest percentage in the CC dataset. Global warming does not eliminate the occurrence of the patterns of SST anomalies associated with MHW activities in this group but favours the occurrence of some nodes within the group. For instance, while the PRS and SAI climates feature all nodes in the group, the CC climate rarely features nodes with the weakest SST anomalies (i.e., Nodes 1, 4 and 7) and show nodes with the strongest SST anomalies (i.e., Node 5). In addition, the PRS and SAI climates exhibit the patterns of this group across all seasons, whereas the CC climate rarely features them during summer (December to February; Fig. 8). Nevertheless, the results suggested that SAI intervention is able to offset the impacts of global warming on the occurrence of the SWC patterns.

The second group, the AC cooling (ACC; Nodes 2, 3, 6 and 9) links MHW activities in the AB region with stronger SST anomalies south of 39°S. Compared to the first group, this group shows stronger SST anomalies associated with MHW activities in the AB region, and the weaker SST anomalies associated with MHWs are shown along the AC pathway and its leakage. This suggests that MHW activities in this case are rooted in the interaction between the large-scale circulations, the AC, and/or to ocean atmosphere interaction (Fig. 5). The cooling along the AC pathway and its leakage suggest that the AC would be transporting relatively cold water into the Atlantic Ocean and in the AB region26. The group accommodates the minority of SST anomaly patterns in the three climate datasets (PRS: 19%; CC: 46%; SAI: 22%) but with the highest percentage in the CC dataset. In the CC scenario, the weakening of winds south of 36°S and the heat gain from ocean in the entire study area would also explain the higher SST anomalies in this group (Fig. 5f and j). The three climate datasets agree that these patterns predominantly occur from October to May, except for Node 2, which is found in all the season but rarely occurs in the CC climate.

Hence, apart from increasing the frequency, intensity and duration of MHW activities in the AB region (Fig. 4), the climate change alters the synoptic patterns of SST anomalies during the MHW events in AB region (Table S1 in the supplemetntary material). Although, with reference to PRS, the CC climate increases the number of occurrences of both SWC and ACC patterns, the increases favours ACC more than SWC, because the percentage contribution of SWC decreases by 25.5% while that of ACC increases by 26.5%. The SAI restores the percentage contribution of SWC and ACC patterns to the PRS values. Climate change eliminates the Node 4 pattern, in which SST anomalies in the AB region are connected to the weakest anomalies in the northwestern Atlantic portion of the study area. This node accounts for 20% of SST anomalies in the PRS and the SAI restores its occurrence to 16%.

Potential implications for the Agulhas Bank marine ecosystem and surrounding populations

The present results indicate that climate change would increase the intensity, frequency and duration of MHWs off the coast of South Africa. The increase may have huge societal impacts in the region. Although the intensity of MHWs off the South African coast may be weaker compared to areas such as the eastern tropical Pacific or the Australian coastline, Roy, et al.47, shows that an increase of even 0.5 °C can disrupt fish habitats. Therefore, the projected changes in MHWs under climate change and the persistent SST anomalies higher than 0.5 °C associated with the MHW events in the AB region are significant enough to potentially harm marine ecosystems and fisheries off the South African coast. While some studies suggest that certain MHW events can benefit specific marine species in particular regions, the majority agrees that the increasing frequency of these events can lead to habitat loss due to extreme temperatures, reproductive failures, range shifts, and even species mortality15,32,52,58. Additionally, MHWs can trigger harmful algal blooms, deoxygenate ocean waters, and contribute to disease outbreaks15,38,55. For instance, the 2011 MHW event in western Australia caused catastrophic fish mortality due to elevated temperatures that led to water deoxygenation11,60. SAI may potentially benefit the AB shelf sea by helping stabilize this region, which is known as a nursery ground and spawning site for commercially important species such as Agulhas sole, Cape anchovy, chokka squid, hake, and kingklip44. However, given that SAI could also alter rainfall patterns, potentially leading to detrimental consequences in sectors such as agriculture in certain regions1, this study underscores the need for more regionally focused research to comprehensively evaluate both the potential benefits and risks of SAI, in order to support informed policy decisions.