Introduction

The Somali coastal region is unique in the global ocean, with the world’s strongest coastal seasonal upwelling, the Somali upwelling, and a seasonal reversal in both atmospheric and oceanic circulation due to the monsoon winds1,2. The regional ocean circulation is dominated by the Somali Current (SC), the only major western boundary current in the world that reverses seasonally, which becomes one of the fastest open-ocean currents when flowing northward during the Southwest monsoon (May to September)1,3,4. The Somali upwelling is the largest upwelling system in the Indian Ocean and the fifth largest coastal upwelling in the world4,5. This exceptional upwelling occurs during the Southwest monsoon partly due to the cross-equatorial atmospheric flow between coastal East Africa and India, the south-westerly monsoon winds, which emerges in May and strengthens into the Findlater Jet (a low-level atmospheric jet), from June to September2,6,7,8. The southwest winds lead to a positive alongshore wind stress causing offshore Ekman transport9,10. Additionally, as the Findlater Jet extends across the Arabian Sea, the shear zone north of the jet exerts strong positive Wind Stress Curl (WSC) driving an Ekman pumping (open ocean upwelling) component along the Somali coast1,4,8,11,12. The wind-induced Somali upwelling component is estimated to be six times stronger than that in the most intense of the eastern coastal upwelling zones1.

The Somali upwelling brings up large amounts of nutrients ( > 15 μM) that approach regional deep-ocean values (20 μM), leading to an enhanced phytoplankton bloom (bloom hereafter), in the ocean surface layer, and hence important productivity off the Somali coast3,13,14,15, ultimately sustaining rich fisheries1. Somali fishery resources are considered among the richest in Africa16,17. The potential yield of the small pelagic fish within the Somali Exclusive Economic Zone, evaluated at 75,000 metric tons, is one of the largest in the western Indian Ocean18. The large pelagics, including tuna and tuna-like species (e.g., yellow-fin, big-eye, skipjack) have high catch rates and are more lucrative despite being highly migratory15,17,19. The traditional fishing grounds are located off the Somali coast between 5°N and 10°N due to the dominant upwelling during the Southwest monsoon15,17,19.

Adding to the strong winds and rough sea state (high waves) during the Southwest monsoon are the fast-flowing SC and a series of associated anticyclonic eddies. In typical conditions, the northward flowing SC deviates from the coast at around 4°N and ~10°N to flow eastward at the northern parts of the clockwise Southern Gyre (SG) and the Great Whirl (GW), a strong and large anticyclonic eddy3,4,14,20,21,22. The SC interacts also with other anticyclonic eddies such as the Socotra Eddy (SE), which develops to the east of the Socotra Island and is usually seen during mid-Southwest monsoon21,23. The GW is a major physical phenomenon during the Somali upwelling season. As the GW intensifies, it can reach 1000 m depth and transport an estimated 10–60 Sv of water, with speeds of 0.1 m s−1 on average21,24,25. The GW position is denoted by the position of the zero WSC contour line off the Somali coast2,25,26. The GW is generated to some degree by the remote forcing (downwelling Rossby waves) prior to the onset of the Southwest monsoon21,25. The westward propagation of downwelling Rossby and the negative WSC results also in an overall downwelling effect that includes other anticyclonic eddies as well as the GW27.

The GW has been extensively studied from a physics point of view21,22 including its impact on the atmospheric boundary layer26. The GW has been mostly studied using satellite observations and models as it is difficult to study in-situ due to the piracy off the Somali coast28. The GW seasonal behaviour is thought to affect the ecosystem of the Arabian Sea29. However, it is still unclear what role the GW plays in the Somali upwelling productivity (i.e., resulting bloom). Hitchcock et al. 30 investigated the fate of the upwelled waters along the northern edge of the GW from in-water surveys in 1995. They found that maximum Chlorophyll-a (Chl-a, a proxy for phytoplankton biomass) levels located at the northern edge of the GW are mainly the result of the rapid offshore advection of nutrients rich waters from the coastal upwelling. More recently, Dai et al. 31 used remote sensing and modelling approaches to examine the GW influence on Chl-a off the Somalia coast but only for 2017 and in terms of how phytoplankton can be imported into the interior of the GW. In this study, we propose to examine the role of the GW in the areal extent of the bloom along the Somali coast, focusing on the months of active upwelling (May–September), and using ~25 years of satellite-derived ocean colour observations (1997–2021). Understanding the variations of the GW, including sudden changes or regime shifts, and how it interacts with the Somali upwelling productivity can help in the better management of fisheries.

Our study provides observational evidence indicating that the absence or smaller coverage of a giant ocean vortex, the GW, paves the way for other less dominant anticyclonic eddies to bring a stronger downwelling effect, decreasing the areal extent of the bloom for the past two decades. Moreover, we found that an abrupt change in the Indian Ocean Dipole (IOD) in 2006 has likely contributed to the regime changes seen in the extents of the bloom and the downwelling effect. We also show a decline in the bloom intensity (Chl-a levels) over the upwelling zone for the past two decades, following decreases in upwelling strength.

Results

Variability in upwelling and productivity

In climatological Southeast monsoon conditions, both the Alongshore Wind Stress (AWS) magnitude and positive WSC (WSC + ) intensify (reaching values > 0.16 N m−2 and >2.5 10−7 N m−2 per 104 km respectively) along the Somali coast (Fig. 1a and Supplementary Fig. S1a), compared to the Northeast monsoon, inducing upwelling4,9,10. The surface signature of the wind-driven Somali upwelling is visible in the satellite derived fields as cooler Sea Surface Temperature [SST] ( < 26 °C) associated with higher phytoplankton biomass ( > 0.55 mg m−3) and lower negative Sea Level Anomaly [SLA-] (<-0.05 m) relative to the surrounding waters (Fig. 1b–d). The productive upwelling signal is also reflected in the satellite derived elevated primary production ( > 1.4 g C per m2 per d) compared to offshore areas (Fig. 1e). This upwelling signature is more marked at the SC ~ 10°N wedge than that at ~4°N which matches the SC behaviour (of stronger deviation at ~10°N than ~4°N) (Fig. 1b–e, Supplementary Fig. S1b), further indicating its advective impact on the upwelling spatial extent4,15. The overall agreement of the climatological spatial distributions of both, the wind variables (i.e., high AWS and WSC + ) and the upwelling effects (i.e., low SST, low SLA-, high Chl-a and primary production) along the coast (Fig. 1a–e, Supplementary Fig. S1a), strengthens the idea that the wind is an important forcing for the upwelling4,9,26,32.

Fig. 1: Seasonal and interannual variations of upwelling parameters along the Somali coast.
Fig. 1: Seasonal and interannual variations of upwelling parameters along the Somali coast.
Full size image

Climatological spatial distribution for June (representative of the Southwest monsoon) 1998-2020 of: a wind stress, b SLA c SST, d Chl-a, e PP. The dashed black contour line denotes the upwelling site area. f–j Seasonal cycles and k–o interannual (as averages of Southwest monsoon months) timeseries of upwelling parameters (AWS, SLA, SST, Chl-a, PP) over the upwelling site for 1998 – 2020. The light shaded areas in f–j represent ± 2 standard deviations from the climatology. The vertical lines in f–j mark the time of onset (June) and end of the Southwest monsoon (September). The dashed thin lines in k–o indicate the timeseries linear trend for 1998-2020. The dotted thin line in n indicate Chl-a linear trend for 2002-2020. The trend significance was calculated using Mann-Kendall test and t-test.

The temporal variations of the upwelling parameters over the immediate “upwelling site” (Fig. 1a–e, black dashed contour, see “Methods” for details on the choice of region) reveal a strong seasonality (Fig. 1f–j, Supplementary Fig. S1c). From May to September (Southwest monsoon), increases in the ASW, WSC + , Chl-a and primary production (by up to ~0.2 N m−2, ~4*10−7 N m−2 per 104 km, 0.7 mg m−3 and 0.6 g C per m2 per d respectively) are observed associated with a decrease in SLA, becoming negative, and SST (by down to −0.2 m and 4.5 °C respectively) compared to the rest of the year (Fig. 1f–j, Supplementary Fig. S1c). The opposite effects are seen from around December to March during the Northeast monsoon (Fig. 1e-h, Supplementary Fig. S1c). Interannually, the Chl-a Southwest monsoon (May-September) averages for 1998-2020 show some variations with more pronounced troughs (such as in 2019 and 2012) than peaks, but more importantly a continuous decreasing trend becoming significant from 2002 (Fig. 1n). Satellite-derived monthly primary production shows also significant decline from 1998 to 2020 (Fig. 1o), further demonstrating the low productivity trend over the upwelling site. The negative productivity trend is associated with significant increasing trends in SLA and SST and a negative significant trend in the AWS over the same period (Fig. 1k–o). Hence, it is likely that the Chl-a negative trend depicted for 2002-2020 (Fig. 1n) is on-going from 1998, especially that the Chl-a is highly negatively (positively) correlated to SLA (AWS) (Table 1) and the number of missing data is much higher in these first years of the continuous satellite ocean colour data record (cf. Supplementary Fig. S2b). No significant trend is seen in the WSC for Southwest monsoon averages of 1998–2020 (Supplementary Fig. S1d). Note that the rest of the analysis is carried out using Chl-a as it shows similar variability and trend to the primary production (Fig. 1d, e, i–j, n, o) in addition to both variables being highly correlated (0.8 with p-value < 0.05), which suggests that Chl-a variations parallel those of primary production in the region (Supplementary Fig. S2a).

Table 1 Pearson correlation coefficient (Γ) significant at the 95% level between upwelling parameters over the upwelling site for September 1997 – June 2021.
Full size table

Chl-a shows high and significant positive (negative) correlations (of the order of 0.55-0.8 with p-values < 0.05) with the AWS and WSC (SST and SLA) in the climatology, monthly means and Mays (the onset of the Southwest monsoons) for the nearly 24 years of observations (Table 1). This is in agreement with the wind driven upwelling along the coast. Although a good agreement is found for the seasonal cycle between the Chl-a, SST, wind forcing (AWS, WSC) and SLA (Fig. 1f–j, Supplementary Fig. S1c), at the interannual scale, considering all months or the Southwest monsoon months, the SLA parameter is the most closely linked to Chl-a over the upwelling site (Table 1). Previously, lowered SLA was also associated with the coastal upwelling along the Somali coast30,33,34. The fact that SST was not found to be the best upwelling indicator here is as expected. Indeed, the overall surface cooling off the Somali coast, although uniform, during the upwelling season is not entirely due to the upwelling effect and the changes in the net negative buoyancy flux and strong wind energy input contribute to the decrease in SST4. The AWS and WSC were found significantly correlated to Chl-a during the onset of the Southwest monsoon (Mays 1998–2021) (Table 1), indicating the triggering effect of the wind forcing on the upwelling biological response. Reanalysis-based findings of a decreasing tendency in upwelling-favourable winds along the Somali coast10 and in the Northwest Indian Ocean35 agree with the significant negative trend in the AWS for 1998-2020 over the upwelling site (Fig. 1k). Moreover, a longer-term trend analysis of zonal wind stress, for 1958-2021 Southwest monsoon months, show a significant decreasing trend off the Somali coast (Supplementary Fig S3a).

The Southwest monsoon Chl-a maxima over the upwelling site reveals “highs” and “lows” with a similar order of magnitude for elevated values (>0.55 mg m−3), including their seasonal progression, while more pronounced changes are seen in the spatial extent of these higher Chl-a levels (Supplementary Fig. S4). The extreme spatial extents in the bloom (see “Methods” and Fig. 2 for details) are associated with the upwelling signal (as indicated by low SLA- coverage), developing at the inner edges of the SC, as it veers from and reverses back to the coast, and always bounded by the GW (Supplementary Figs. S7, S8 and Note 1). The latter can sometimes take the place of the SG and/or other anticyclonic eddies such as SE (Supplementary Fig. S7, S8 and Note 1) as previously modelled23. In agreement with our observations, the GW and SC flow along its northern arm have been shown to enhance the bloom extension off the coast30,31.

Fig. 2: Bloom surface coverage along the Somali coast during the Southwest monsoon.
Fig. 2: Bloom surface coverage along the Somali coast during the Southwest monsoon.
Full size image

a Monthly timeseries of the bloom (approximated by Chl-a) surface coverage (in km2) for May–September months of the period September 1997–June 2021 (in green). The values falling on the upper side of the 253940 km2 level, which is the mean of the surface coverage maxima timeseries (dashed horizontal line), indicate extreme large blooms and are highlighted in light red dots. The values falling on the lower side of the 53333 km2 level, which is the mean of the surface coverage minima timeseries (dotted horizontal line), indicate extreme small blooms and are highlighted in light blue dots. The most extreme large bloom (Sep-03) and smallest bloom (May-12) are further highlighted with dark red and dark blue circles, resp. Example of the b Chl-a (in mg m-3) and c SLA (in m) surface distributions for Sep 06. The black contour lines in b denotes the 0.55 mg m-3 Chl-a contour line. The dashed ellipses and black stars in c represent the edge and centre of anticyclonic eddies. The dashed black box in b, c indicate the location of the bloom box.

Abrupt change in the extents of the bloom, downwelling effect and the GW

During the Southwest monsoons of 1997–2021, a significant negative correlation (-0.6 with p-value < 0.05) was observed between the extent of the bloom and positive SLA (SLA + ) surface coverage, indicative of downwelling within the bloom region (Fig. 2). The bloom extent correlates also significantly (0.85 with p-value < 0.05) with SLA- surface distribution over the bloom region (Supplementary Fig. S9a), further confirming that the bloom extent during the Southwest monsoon is mainly associated with the upwelling signal. Additionally, a strong significant and positive correlation (0.7 with p-value < 0.05) was found between the SLA+ surface coverage inside the bloom region and the surface coverage of the non-bloom area (Supplementary Fig. S9b). These relationships suggest that downwelling factors, which include the GW and other anticyclonic eddies (e.g., Supplementary Figs. S7, S8), contrast the areal extent of the bloom and the upwelling signal. Similar bloom scenarios can correspond to different conditions in size and position (closeness to the coast, north/south migration) of the GW which vary considerably from year-to-year and month-to-month (Fig. 2a, Supplementary Fig. S9c, d and Supplementary Note 2).

The bloom extent timeseries shows an overall decreasing tendency over the past two decades (2002–2021) similarly to the phytoplankton biomass downward trend over the upwelling site (Figs. 1n and 3a). The negative tendency in the bloom surface coverage is also associated with an increasing tendency in the SLA+ surface coverage inside the bloom box (Fig. 3). This suggests that the downwelling effect has been limiting the bloom extent over the past two decades. Visual inspection of the results in Fig. 3 suggests a possible change in regime around to 2007, so we test this by applying a changepoint analysis. The changepoint analysis reveals evidence of an abrupt shift in May 2007 after which the SLA+ (bloom) surface coverage increased above (decreased below) the overall mean with a 1.23 105 km2 (-0.08 105 km2) difference between 1998–2007 and 2007–2020 (Fig. 3). After the stepwise increase (decrease) of 2007, the extent of the SLA+ (bloom) remained at a high (low) regime without returning to its initial state (Fig. 3). The GW surface coverage contribution to this downwelling effect over the bloom region varies from at least ~4% up to ~20% and is on average 9.54% (with a standard deviation of 3.2%) for the full observational period (Fig. 4a). Although these extrema (i.e., the 4% and 20%) coincided with extremely large and small bloom extents (Figs. 2a and 4a), it does not imply that the greater (lower) the percentage of the GW in the SLA+ coverage the larger (smaller) is the bloom (i.e., their relationship is not straightforward). For example, the smallest GW percentage ( ~ 5%) in the downwelling effect excluding May (where the GW do not always form and the upwelling season is at its very start), reached in July 2013, has an average bloom extent (Figs. 2a and 4a). Despite this nonlinear aspect, the GW has the largest impact in the SLA+ coverage among the anticyclonic eddies both in terms of surface coverage and number of occurrences. Indeed, the GW surface coverage is most of the time larger than that of any other anticyclonic eddy in the region (Fig. 4b). Additionally, the GW is present 88.14% during the Southwest monsoon months (104 out of 118 months) from September 1997 to June 2021, while other anticyclonic eddies take place (with or without the GW) only 43.22 % of the Southwest monsoon months (51 out of 118 months) for the same period (Fig. 4b).

Fig. 3: Regime shift of the bloom and downwelling extent for Southwest monsoons of September 1997–June 2021.
Fig. 3: Regime shift of the bloom and downwelling extent for Southwest monsoons of September 1997–June 2021.
Full size image

Monthly variations for May-Septembers 1997-2021 of the surface coverage (in km2) of a the Chl-a within the 0.55 mg m-3 contour line versus b that of the SLA+ within the bloom box. The Pearson correlation coefficient (Γ) significant at the 95% level is given in a. The solid red and dashed vertical grey lines in a) and b) represent the significant regime shifts (most significant change from the mean) and detected changepoint in the timeseries mean.

Fig. 4: GW surface coverage variability and its contribution to the total downwelling effect over the bloom region.
Fig. 4: GW surface coverage variability and its contribution to the total downwelling effect over the bloom region.
Full size image

Monthly variations of the a the percentage of the GW in the SLA+ over the bloom region and b GW surface coverage (in km2) for May-Septembers 1998-2021. The blue dots in b represent the surface coverage of other anticyclonic (A) eddies occurring with or without the GW over the bloom region. The solid red and dashed vertical grey lines in a represent the significant regime shifts (most significant change from the mean) and detected changepoint in the timeseries mean.

The GW percentage in the SLA+ coverage over the bloom box exhibits not only an overall decrease for 1998–2020 but also a significant changepoint (most significant change in the mean) from August 2006, after which a stepwise decrease by -2.51 % is evident without a return to the initial state (Fig. 4a). This might suggest that the GW is becoming less important, however, the GW surface coverage itself does not show a significant changepoint nor a trend for 1998–2020 indicating an overall stable regime (Fig. 4b). Although the GW surface coverage regime itself is stable, extremes of large or small GW surface coverage are present (Fig. 4b). The extremes and overall regime of the GW surface coverage could be related to changes in the remote forcing (Rossby waves) prior to the start of the Southwest monsoon as this mechanism initiates the GW21,25,26. The GW changes could also be due to the alterations in the WSC and its zero line, which denotes the positions of the Findlater Jet and the GW, ultimately controlling the negative versus positive WSC distribution in the Arabian Sea2,25.

Controls on the bloom extent regime

By becoming increasingly important, the downwelling effect that is not due to the GW (e.g., other anticyclonic eddies) has been driving the increasing tendency in the downwelling effect (SLA + ) over the bloom region for the past two decades (Fig. 5a). The non-GW downwelling effect surface coverage abruptly increased from its mean by 1.34 105 km2 (significant changepoint detected in June 2007) between 1998–2007 and 2007–2020 (Fig. 5a). This suggests that the less dominant eddies (in terms of surface coverage and occurrence) are getting a more important effect and leading the SLA+ extent regime over the bloom region. Indeed, when the GW is absent or has a smaller coverage than other anticyclonic eddies, the generated downwelling effect is overall stronger (i.e., of a wider coverage) than when the GW is present with a high coverage (Fig. 5b). Thus, it is the absence of the GW (or its smaller coverage) that is more important for the bloom extent regime than its presence (which may enhance or bound the bloom extent (e.g., Supplementary Figs. S7, S8 and Note 1)). Both the SLA+ with absent/smaller GW coverage and that with a GW of higher coverage than other anticyclonic eddies have increased abruptly above their mean after 2007 by 1.47 105 km2 and 1.33 105 km2 respectively (significant changepoint detected in May 2007) without moving back to their initial state (Fig. 5b).

Fig. 5: Changes in the non-GW versus the GW downwelling parts and Link to the IOD regime change.
Fig. 5: Changes in the non-GW versus the GW downwelling parts and Link to the IOD regime change.
Full size image

Monthly variations for May-Septembers 1997-2021 surface coverage (in km2) of a the SLA+ of the non-GW downwelling effect over the bloom region; and of b the SLA+ with absent or lesser GW coverage than other anticyclonic eddies (A) (in black) versus that of the SLA+ with present GW of higher coverage than other A (in magenta). The part of the timeseries that are dashed in b) indicate non available data. c Monthly IOD Dipole Mode Index (DMI) timeseries for 1979–2020. The solid red and dashed vertical grey lines in a–c represent the significant regime shifts (most significant change from the mean) and detected changepoint in the timeseries mean. d Cumulative sums of monthly IOD anomalies for 1979–2020. The Cumulative sums summarizes major changes by smoothing high frequency noise along the time-series (cf. Methods for more details). Increasing slopes represent positive temporal trends (periods with above average values) while decreasing slopes highlight negative temporal trends (periods with below average values), and horizontal segments indicate average conditions.

The abrupt shifts observed in the extents of the bloom and downwelling effect are all occurring around late Southwest monsoon 2006—early Southwest monsoon of 2007 (Figs. 3, 4a, 5a, b). It might be argued that this shift can be attributed to decadal variability since it is important in the tropical Indian Ocean36. Although the Chl-a record is not long enough, extending the SLA-based timeseries in Figs. 2b, 3a, 4a, b back to 1993 demonstrates further that the abrupt shift is indeed after 2006 and not just a phase of decadal variability (Supplementary Figs. S5 and S6). Given 2006–2007 is a strong positive Indian Ocean Dipole (IOD) year (e.g., Fig. 5c) and that the GW start time and lifetime were recently found correlated to the IOD, generating later, and lasting less than normal during strong positive IOD years37, we hypothesised the potential occurrence of an abrupt shift in the IOD possibly leading to the changes in the SLA+ and bloom coverages. The changepoint detection algorithm reveals a significant changepoint in the IOD 1979–2020 timeseries, indicative of a regime shift (most significant change in the mean), by 0.236 difference pre and post shift, from July 2006 (Fig. 5c). This is 1–11 months prior to the abrupt shifts observed in the bloom and downwelling effect extents displayed in Figs. 2a, 3, 5a, b. The application of normalized cumulative sums of anomalies (see “Method” for details) to the IOD timeseries confirms the timing of this abrupt shift (Fig. 5d). The cumulative sums of IOD anomalies show values consistently below the mean, producing an overall decreasing trend, prior to July 2006 and an opposite pattern after July 2006, mirroring the abrupt change detected from the changepoint analysis (Fig. 5c, d). Additionally, significant positive correlations (of the order of ~0.4 with p-value < 0.05) are found between the IOD and the total downwelling effect and non-GW downwelling effect over the bloom region (Figs. 3b, 5a) during the Southwest monsoons (May-Septembers) of 1998–2020. These results suggest that the increase in the downwelling effect (leading to a decrease in the bloom extent regime) over the bloom region around 2006–2007 is likely due to the abrupt shift in the IOD after July 2006. The 2006 shift observed here in the 1979–2020 IOD timeseries is about three decades long. This is not totally surprising as, although the IOD is an interannual variability mode38, it has an important decadal to multi-decadal variability39,40. The changes in anomaly are also unlikely to be due to the warming of the Arabian Sea as it has been continuous since the early 1990s41.

Although most past studies failed to demonstrate strong links between the IOD, as a climate variability mode intrinsic to the Indian Ocean, and the GW development or intensity21,24, the recent work by Dai et al. 37 showed that the GW start time (lifetime) is significantly positively (negatively) correlated to the IOD during September-November months. The authors demonstrated also that anomalous northeasterly wind and Rossby waves during strong positive IOD years result in a late generation time and short lifetime of the GW37. Here the timing of the abrupt shift found in the IOD is few months prior to the shifts detected in the GW percentage in the SLA+ over the bloom region and the other downwelling effect metrics (cf. Figs. 3, 4a, 5). This delayed effect is coherent with Beal and Donhue21 suggestion of the formation of the GW being the result of remote forcing transmitted by Rossby waves before the onset of the Southwest monsoon. Our results support further Dai et al. 37 proposal that the IOD might influence the GW variability.

It is also acknowledged that the Rossby waves propagation emanating from both the Northwest Tropical Indian Ocean (NTIO) (as reflected coastal trapped Kelvin waves) and from the South Tropical Indian Ocean (STIO) (as reflected equatorial Kelvin waves) play a key role in the formations of strong positive IODs and the GW37,42. Positive IODs were also associated with combined downwelling Rossby waves (as reflected equatorial Kelvin wave and direct wind induced), which lead to thermocline deepening and sustained anomalously warm SST in the Western equatorial Indian Ocean43. These Rossby waves are excited by changes in the anticyclonic Wind Stress Curl Anomalies (WSCA)44,45,46. A time-space diagram along the latitude line 9°N of SLA shows a westward movement, indicative of Rossby wave propagation towards the Somali region44, with more positive SLA starting to exceed 0.2 m after 2006 (Fig. 6a, b, green contours). This apparent shift in behaviour in the downwelling Rossby waves propagation in the NTIO is associated with change in the anticyclonic WSCA detected along the 7.5°N latitude line off India west coast with higher values exceeding 0.5 N m−2 per 104 km after 2006 (Fig. 7a, b). We also detect a change in the Rossby waves propagation in the STIO along 7.5°S with more frequent positive SLA higher than 0.2 m after 2006 (Fig. 6c, green contours), except for 1997–1998 which is likely due to the exceptional conditions of 1997–1998 super El Nino47,48,49. Indications of change in the anticyclonic WSCA in the STIO can be observed along 5°S off Sumatra from 2006 with more positive WSCA higher than 1.12 N m−2 per 104 km (Fig. 7c, green contours). The changes in the WSCA in the NTIO and STIO can be attributed to alterations in the Indian Ocean Walker circulation via anomalous easterlies over the equator and southeasterlies off Sumatra which represent the first atmospheric triggers of the IOD45,46,50,51. Other studies have linked the Rossby waves prorogation in both the NTIO and STIO to changes in the thermocline in the tropical Indian Ocean which in turn alter the SST leading to the IOD temporal signature43,46,51,52. Evidently, the SST changes would in turn affect the easterlies and southwesterlies, but investigating this internal feedback loop requires a sperate modelling study in future research. Our results of changes in the Rossby waves propagation and the WSCA around 2006 in the NTIO and STIO, supports the mechanistic link between the IOD, the GW and the upwelling productivity off the Somali coast. Although the 2006–2007 positive IOD is not the strongest positive IOD event on record, it was unique because it marked the start of a rare manifestation of three consecutive positive IODs53,54, followed by an escalating occurrence of positive IODs in recent years as shown by Albert et al. 55 (their figure 13). In terms of upwelling intensity, the decrease in the zonal wind stress during the Southwest monsoon (1958–2021) off the Somali coast post 2006 (Supplementary Fig. S3b) suggests that the local forcing has also changed after 2006.

Fig. 6: Westward propagation of Rossby waves as depicted from SLA signal in the NTIO and STIO over the period 1993-2020.
Fig. 6: Westward propagation of Rossby waves as depicted from SLA signal in the NTIO and STIO over the period 1993-2020.
Full size image

a SLA (in m) for January-April (JFMA) 2006 over the Indian Ocean and 9°N and 7.5 °S latitude lines along which b, c time-space diagram of SLA are derived between 67-76.5°E and 38-81°E respectively. The green lines in b, c highlight the SLA 0.2 m contours and the dashed black lines indicate the change in SLA behaviour from 2006.

Fig. 7: Changes in WSCA in the NTIO and STIO over the period 1993-2020.
Fig. 7: Changes in WSCA in the NTIO and STIO over the period 1993-2020.
Full size image

a WSCA (N m−2 per 104 km) for June 2006 over the Indian Ocean and 7.5°N and 5 °S latitude lines along which (b-c) time-space diagram of WSCA are derived between 70−73°E and 90.5−100.5°E respectively. The green lines in b highlight the WSCA 0.5 N m−2 per 104 km contours and those in c the WSCA 1.12 N m-2 per 104 km contours. The dashed black lines indicate the change in WSCA behaviour from 2006.

Discussion

The results presented in this study highlight the key role ocean vortices can play in regulating the productivity of one of the world’s largest upwelling systems and demonstrate how atmospheric (e.g., wind) or climate mode (e.g., IOD) mechanisms could influence it. With the decline in the phytoplankton biomass over the upwelling zone for the past two decades, following decreases in upwelling strength and wind forcing, an investigation into the wind-biological interactions, using ocean–atmosphere-biogeochemistry coupled models could provide further insights. How the upwelling strength might affect deoxygenation occurring in nearby regions of the Arabian Sea8 is also yet to be examined in future research. Our results demonstrate also that the bloom extent experienced an abrupt decrease associated with an abrupt increase in the downwelling coverage after 2007, linked to an abrupt increase in the non-GW downwelling effect in 2006-07. More specifically, when the GW is absent or has a smaller areal coverage than other less dominant anticyclonic eddies, the generated downwelling effect is overall stronger than when the GW is present with a high coverage, limiting further the bloom extent. Therefore, we conclude that the absence of the GW is more important for the bloom extent regime than its presence. Simply put, the GW absence leads to a stronger downwelling impact on the bloom extent than when it is present; and on the long term this behaviour affected the bloom extent regime. The abrupt changes in the extents of the bloom and downwelling regimes observed around 2006-07, are likely due to the abrupt shift in the IOD detected in July 2006. Changes in other remote large-scale influences such El Niño Southern Oscillation, may produce changes that affect nutrient supply to the upwelling waters56 and remain to be investigated.

The declining productivity, seen here in both the bloom intensity and extent and their controls, provides a foundation for future research on the productivity regime changes using climate model projections. This is particularly important because the Southwest (Indian summer) monsoon, which is when the Somali upwelling occurs, has been identified as one of the climate system’s tipping points57. The fate of both the GW and anticyclonic eddies interacting with the upwelling productivity need to be examined in such scenarios with a disruption of the monsoon systems in the future climate projections. This, in addition to continuous monitoring of environmental conditions using satellite observations, could provide useful information on potential low productivity scenarios, and in turn, fisheries, which are of socio-economic importance to the local communities and stakeholders in Somalia and the northwestern Indian Ocean.

Methods

Satellite derived data

We use monthly reprocessed satellite-derived Chl-a data available from the Ocean-Colour Climate-Change Initiative project (OC-CCI version 558, http://www.esa-oceancolour-cci.org/), at a spatial resolution of 4 km. This dataset is the most consistent timeseries of multi-satellite (MODIS-Aqua, SeaWiFS, MERIS and VIIRS) global ocean colour data56. The Chl-a monthly composites span the period September 1997 to June 2021. They are used to examine the Chl-a variability seasonally and interannually off the Somali coast. Although the OC-CCI processing applies a multi-chlorophyll algorithm approach59, satellite-derived Chl-a could still be overestimated in coastal and/or shallow waters (generally shallower than 30 m60). The bottom reflectance issues, and suspended material (such as sediments and dissolved organic matter) present in shallow optically complex waters which may lead to high water leaving radiance, and in turn, an overestimation of the correction term61. However, these high Chl-a values in coastal/shallow areas are not necessarily erroneous. They could reflect the highly productive coral reef zones, or chlorophyll-rich detritus regions enhancing phytoplankton growth near the coast62,63. Regardless, the Somali region consists mainly of deep oceanic waters including the “upwelling site” region (see Section 4.4 for details and Fig. 1a–e, black dashed contour) where 98.15% (87%) of the data represent waters deeper than 30m (100 m).

The satellite-derived primary production data used here are processed by ACRI-ST based on the Copernicus-GlobColour processor and retrieved from the Copernicus Marine Environment Monitoring Service (CMEMS) (http://marine.copernicus.eu/services-portfolio/access-to-products/). This is a global multi-satellite product available monthly with a horizontal resolution of 4 km and spanning the period September 1997 to January 2023. Primary production was calculated using Antoine and Morel Algorithm64, and uses satellite-derived ocean colour products (merged Chl-a, Photosynthetically Active Radiation [PAR]65 and Kd(PAR)), SST and a mixed layer depth climatology, estimated based on de Boyer Montegut66 definition. The primary production algorithm is validated globally against in situ data67.

We utilize altimetry-derived SLA and absolute geostrophic currents processed by CLS (previously by AVISO (Archiving, Validation and Interpretation of Satellite Oceanographic Data)) and distributed by CMEMS (http://marine.copernicus.eu/services-portfolio/access-to-products/). These multi-satellite (Jason-3, Sentinel- 3 A & 3B, HY-2A, Saral/AltiKa, Cryosat-2, Jason-2, Jason-1, T/P, ENVISAT, GFO, ERS1/2) dataset are daily composites gridded at 25 km spatial resolution from the delayed time DUACS_DT2018 version and spanning the period 1993-2021. The SLA product is referenced to a twenty-year (1993-2012) mean68. Only the period September 1997- June 2021 is considered, to match the satellite Chl-a data. The SLA data is used to characterise the upwelling surface signature and eddies. The geostrophic currents are used to examine the surface circulation. To overcome the geostrophic limitation near the equator, the geostrophic current field product applies the β plane approximation69. Satellite altimetry data can be impacted by the sensor contamination near the coast70. However, our region (Fig. 1a-e) comprises mainly open waters. The product version has also a reduction in error on geostrophic currents in coastal zone by >15%71.

The reprocessed SST Met Office Operational-Sea-Surface-Temperature-and-Sea-Ice-Analysis (OSTIA) product, acquired from CMEMS (http://marine.copernicus.eu/services-portfolio/access-to-products/) is used. This multi-satellite global SST dataset is provided daily at 5 km spatial resolution from 1981 to 2021. Monthly SSTs are computed for the same period as Chl-a (September 1997–June 2021).

Wind parameters

We use the ERA-5 reanalysis 10 m zonal and meridional winds produced by the European Centre for Medium-Range Weather Forecasts (ECMWF). This dataset is made available by the Climate Data Store (CDS) (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels-monthly-means?tab=overview; Hersbach et al. 72) at a 25 km spatial resolution and as monthly means from 1979 to present. However, only the period September 1997 to June 2021 is used to match the Chl-a time coverage. We calculate the AWS along the Somali coast using the zonal and meridional wind components. The wind stress fields are computed using the Gill73 equation:

$${{{{{\rm{\tau }}}}}}={{\rho }}_{a}\,{C}_{d}{{{{{\rm{|}}}}}}v{{{{{\rm{|}}}}}}v$$
(1)

where τ and v are the wind stress and wind vectors, respectively. ρa is the air density (1.2 kg m−3) and Cd is the drag coefficient determined following Large and Pond (1981) equation:

$$\begin{array}{cc}{C}_{d}=(0.49+(0.065\,{{{{{{\rm{U}}}}}}}_{{{{{{\rm{wind}}}}}}})){10}^{-3} & {{{{{\rm{if}}}}}}\,{{{{{{\rm{U}}}}}}}_{{{{{{\rm{wind}}}}}}} > ={11{{{{{\rm{ms}}}}}}}^{-1}\\ {C}_{d}=0.0013 \hfill & {{{{{\rm{if}}}}}}\,{{{{{{\rm{U}}}}}}}_{{{{{{\rm{wind}}}}}}} < {11{{{{{\rm{ms}}}}}}}^{-1}\hfill\end{array}$$
(2)

where Uwind is the winds at 10 m. The coastal AWS provides information on the Ekman transport wind-driven upwelling, which bring nutrient-rich waters towards the surface. This variable has been previously used to provide estimates on upwelling intensity, especially in regions around the Equator such as the Somali region where Ekman transports diverge (e.g.,9). We also examine the WSC which is a sign of vertical pumping and proportional to the vertical velocity. WSC induced upwelling along the Somali coast is known to be important during the start of the Southwest monsoon4,26. WSC+ indicates ocean surface divergence forcing waters upward (upwelling - Ekman pumping) in the northern hemisphere10,32 and vice versa for negative WSC. The WSC is computed using the wind stress fields in the following equation:

$${{{{{\rm{WSC}}}}}}=\frac{\partial \tau {{\mbox{y}}}}{\partial x}-\frac{\partial \tau {{\mbox{x}}}}{\partial y}$$
(3)

where τx and τy are the zonal and meridional components of the wind stress (cf. Equation (1)) and x and y are the zonal and meridional dimensions.

We also use monthly zonal wind stress for 1958–2021 calculated from the Japanese 55-year atmospheric reanalysis74 based surface dataset (JRA55-do)75 to investigate the local wind forcing affecting the upwelling intensity at a longer term. This dataset is available at a spatial resolution of 55 km and on a 3 h basis.

IOD data

Monthly IOD Dipole Mode Index (DMI) timeseries for 1979−2020 is used here. This IOD timeseries is the HadISST v1.1 version downloaded from https://psl.noaa/gov/data/timeseries/DMI.

Eddy detection algorithm and GW identification

We use an automated eddy detection algorithm to determine the position and surface coverage of the GW off the Somali coast. Automated eddy detection algorithms can be based on physical quantity, flow geometry characteristics or they can be hybrid, which is a combination of the two former categories76. Here, we apply a hybrid method where SLA (as a physical parameter) is used to spot eddies and streamlines (as flow geometry characteristics) are utilized to determine eddy edges. SLA is conventionally used in identifying and tracking mesoscale eddies in ocean regions deeper than 200m77,78. Since streamlines associated with eddies are parallel to SLA contours under the geostrophic assumption78,79, closed SLA contours are considered here as a proxy for closed streamlines.

The hybrid eddy-identifying method used in this study was employed by Chaigneau et al. 80,81, Zhang et al. 78, Chelton et al. 81 and Xu et al. 82,83. The algorithm identifies eddies by finding their centre and edges (i.e., outermost closed contour of SLA which defines a compact structure)78,83. The centre of a possible eddy is determined by looking for the mass centre of an innermost closed SLA contour. All closed contours, surrounding the centre, with their SLA values changing monotonously outward from the centre are identified. The outermost closed SLA contour defines the eddy edge78. The following eddies-identifying criteria, adapted from Xu et al. 82,83 and Zhang et al. 78 are considered: (1) Eddies amplitude, defined as the SLA difference between the centre (i.e., SLA extremum in the eddy interior) and the edge (i.e., value of the outermost closed SLA contour), is at least 3 cm. The choice of this ±3 cm threshold is due to the altimeter observation error (2-3 cm), so that eddies with amplitudes less than 3 cm are excluded. (2) The zonal and longitudinal extent of the eddy area (i.e., area enclosed by the eddy edge) are at least 0.5° long. (3) The distance between two arbitrary points on the SLA contour should not exceed 400 km. (4) The eddy centre is taken as the midpoint between the centroid of the area within the eddy edge and the location of the SLA extremum, in order to minimize the error. (5) Eddies radius, defined as the radius of a circle with the same area as that of the eddy, should be larger than 45 km. This threshold is established due to the spatial resolving capability of the altimetry data. (6) Data in the equatorial band (i.e., with latitudes in ±3.5 °N range) are disregarded as the Coriolis parameter approaches zero there, lowering the performance of the eddy-identifying method.

The eddy detection algorithm is an adapted version of the MATLAB subroutine SimpleEddyDetection.m that is available freely online84. The algorithm takes as inputs the longitude, latitude, SLA (in m) and date (year, month, day) of the SLA measurements. The output consists of a structure array with the properties of detected eddies: polarity (i.e., cyclonic or anticyclonic), centre, amplitude, radius and edge (defined by the largest closed sla contour). To distinguish the GW from the rest of eddies returned by the algorithm, only the anticyclonic eddies that fall within the GW centre known latitudinal and longitudinal extent (i.e., from 4 to 12°N and from 45 to 56°E)22 are kept. Then the GW locations (i.e., centre and edge) found by the algorithm are plotted superimposed to the SLA monthly maps and to the WSC zero contour line over the full observation record (e.g., Supplementary Figs. S7, S8). In fact, the WSC zero contour line is always tangent to the GW edge21,26. This condition is used to validate the GW location. It is also visually tracked from the start to the end of the Southwest monsoons of the full observational record to further check the identified location.

There have been different approaches and algorithms to detect and track the GW based on streamlines (approximated by SLA), surface currents or Absolute Dynamic Topography21,22,26,81,85,86,87. Although the choice of the method impacts the GW metrics, including its surface coverage (e.g., Melzer et al. 22 method led to the largest estimation of the GW size), the GW characteristics obtained here agree with a number of these previous investigations21,22,87. Our approach was also verified using the WSC zero contour line and a visual comparison of the detected GW and its signature in other ocean fields (e.g., surface currents and Chl-a).

Choice of region for analysis of temporal upwelling parameters

We select an area representative of the immediate “upwelling site” (Fig. 1a–e, black dashed contour), based on the climatological spatial distribution of the upwelling triggers (i.e., high WSC+ ( > 2.5 N m-2 per 104 km) and high wind speed, Fig. 1a, Supplementary Fig. S1a) and effects (i.e., cool SST ( < 26 °C), low SLA- (<−0.05), high Chl-a (>0.55 mg m−3) and primary production ( > 1.4 g C per m2 per d), Fig. 1b–e) along the Somali coastline.

Bloom surface coverage metric

To characterize the surface coverage of the Somali bloom (Somali upwelling productive signal), we used a box, called hereafter “the bloom box” (Fig. 2), in a similar fashion to how other large blooms in the Indian Ocean were described such as the Southeast Madagascar bloom88,89. Although the Somali upwelling productive signal has a well-defined shape in the climatological mean with two wedges at ~4°N and 10°N (Fig. 1b–e, Supplementary Fig. 1b), the interannual variations of its geographical extent can differ considerably from the climatology (e.g., Supplementary Figs. S7, S8). Due to this substantial spatial variability, we tested different Chl-a contour lines (0.45, 0.55, 0.6 and 0.65 mg m−3) comparing them subjectively by visual analysis to the spatial distribution SLA-, which is indicative of upwelling26,30. The pixels with higher Chl-a than a threshold value should match overall the spatial extent of the SLA- (e.g., Fig. 2b, c). The 0.55 mg m−3 Chl-a contour is selected as it was found to be the most representative of the spatial extent of the Somali phytoplankton bloom for each month of the Southwest monsoon over the ~24 years of observations. The bloom surface coverage is then computed by summing the surface coverage (in km2) of all the pixels (each representing 4*4 km2) within this threshold value (0.55 mg m−3 Chl-a contour) and inside the bloom box.

This bloom spatial distribution metric is based on a statistical method with empirical thresholds, similar to upwelling productivity delimitation approaches used in the WIO90. Less empirical methods such as the recent automatic detection of upwelling in the WIO developed by Hammond et al. 91, so called “upwelling watch” is aimed for operational purposes and upwelling detection only. Here, our goal was to delineate the spatial extent of both upwelling productivity and the downwelling effect.

Using the bloom surface coverage metric, extremely “Larger” (“smaller”) extension events can be identified as Chl-a surface coverage values falling on the upper (lower) side of a threshold defined by the mean level of their maxima (minima) timeseries, here corresponding to 2.54 105 km2 (0.53 105 km2) (Fig. 2a, dashed and dotted horizontal lines).

Changepoint detection analysis

To confirm the existence of statistically significant abrupt changes in the SLA+ and Chl-a surface coverage timeseries, we apply a changepoint detection algorithm, the findchangepts MATLAB sub-routine. This algorithm is based on the Pruned Exact Linear Time (PELT) method, an optimal segmentation algorithm, which finds the optimal number of changepoints and has a linear computational cost, leading to accurate segmentation of the data92,93. The PELT method is combined with the modified Bayesian information criterion (MBIC) as the penalty function to select an optimal number of changepoints93. The approach is effective in differentiating abrupt changes from trends and memory in timeseries93. Full details of the algorithm can be found in Beaulieu & Killick93 and references therein. This method was proven to perform better than some commonly used approaches for the detection of abrupt changes in climate timeseries93, such as, the sequential t-test analysis of regime shifts approach by Rodionov94 and the Bayesian linear regression multiple changepoint-detection method by Ruggieri95. Additionally, Rodionov94 t-test method works on the assumption that the data are independent and normally distributed96, which is only valid for some variables.

In findchangepts, a changepoint is a time instant at which a statistical property (e.g., the mean, variance, spectral characteristic, etc.) of a given timeseries changes abruptly. The function partitions a timeseries into parts that minimizes the sum of residual squared error of each part from its local statistical property (e.g., mean). The function allows to set a proportionality constant, MinThreshold, corresponding to a fixed penalty added for each changepoint. MinThreshold specifies the minimum improvement in total residual error for each changepoint. The penalty grows linearly with the number of changepoints. Hence, setting a large MinThreshold reduces the number of change points detected due to noise. Not specifying a MinThreshold or a maximum number of changes, can lead in an extreme case to identify a changepoint between for every point and vanishing of the residual error. The findchangepts algorithm was applied in several oceanography and geophysics studies97,98,99,100.

Here, the changepoint analysis is applied to the gap-filled SLA+ and Chl-a surface coverage timeseries (Figs. 3, 4a, b, 5a–c) using the mean as the chosen statistical property (i.e., detecting most significant abrupt changes in the mean). Changepoints detected in the first or last two years of the timeseries are excluded, as there is not enough information before or after the candidate changepoint to confirm that it is a real abrupt change.

Cumulative sums of anomalies

Cumulative sums are an efficient statistical approach that reduce the high frequency noise in the timeseries while objectively detecting changes in slope/trends in datasets highlighting timing and duration of major (e.g., decadal) changes101,102. Here, we use a variant of cumulative sums where original values are normalized to a mean of zero and unit standard deviation (z-scores) before computing the running sum103. Normalized cumulative sums segments with increasing slopes represent positive temporal trends (positive anomalies) while segments with decreasing slopes highlight negative temporal trends (negative anomalies), and horizontal segments indicate average conditions101. The cumulative sums technique has increasingly been used to determine underlying features of climate variability modes104 and environmental timeseries105,106.