The study area consists of six coastal regions located in the tropical and south Atlantic, three along the South American eastern coast and three along the African western coast (Figure 1). Regions 1 to 4 are located within the tropical band, while Regions 5 and 6 are located roughly south of the Tropic of Capricorn. These regions are mostly characterized by important levels of CHL (Figure 1) and they have a high ecological and socioeconomical (fisheries, tourism, etc.) importance. As will be shown below, these regions present a vigorous interannual variability of CHL affected by diverse oceanographic processes, either at local- or basin-scale, which are different between each other.

In each region, monthly mean time series of spatially-averaged satellite-derived CHL is correlated with 1) 23 spatially-averaged physical variables obtained from oceanic and atmospheric gridded products (Table 1), 2) 14 teleconnection climatic indices (Table 2), to elucidate the physical drivers of the regional biogeochemical characteristics. The analysis period is from January 1998 to December 2021 (24 years). A monthly climatology (1998–2021 period) was subtracted from each data series to produce monthly anomalies and hence focus on interannual variations.

The selected CHL dataset is a Level-4 gridded product by Copernicus Marine Environmental Monitoring Service (CMEMS, https://marine.copernicus.eu/). The other oceanographic variables were taken from different data sources, some of them taken directly from the source while the other resulted from a post-processing carried out for this analysis (Table 1). The large-scale teleconnection indices were mostly provided by the National Oceanic and Atmospheric Administration (NOAA) - Physical Science Laboratory (PSL; https://psl.noaa.gov/data/climateindices/list/); only those indices directly associated with, or with known effects on, the Atlantic were selected. For details about all these variables and indices, see Supplementary Section S1.

Then, in order to identify the dynamical factors that affect the CHL variability, a mechanistic-oriented correlation analysis was carried out for each region (for details, see Supplementary Section S2). This analysis consisted of 3 steps (Figure 2A): 1) correlating the CHL series whith those of the potential drivers, either the oceanographic variables (Table 1) with lag = 0 or the teleconnection indices (Table 2) with 0≤lag≤9 months; 2) ranking the significantly-correlated (p<0.05) drivers’ series according to the magnitude (absolute value) of their correlations; 3) selecting the ranked drivers’ series with no covariability between each other. This selection (Step 3) is carried out by a cross-correlation of the ranked drivers’ series (from Step 2), reducing the number of independent drivers’ series (Figure 2B). The selected drivers’ series are the focus of our mechanistic analysis and are used as predictors to reconstruct the CHL (Section 2.3).

To elucidate the link between the oceanographic variables and the teleconnection indices, i.e. the large-scale dynamical patterns associated with each index, regression maps of the oceanographic variables were calculated as functions of individual climate indices. In preparation for these regression maps, each index was standardized (1998–2021 period), i.e. its mean was subtracted and the result was divided by its standard deviation, to ensure a mean of zero and for not dealing with the indices’ units. Then, a no-intercept linear regression-model was fitted (by a least-square technique) to every point of each variable’s gridded anomaly (1998–2021 period) as function of each standardized climate index; the regression map consists of the fitted coefficients.

Given the relatively short temporal coverage of the CHL datasets, barely long enough for interannual- but insufficient for decadal-variability analyses, it is desirable to extend its temporal extension. Then, herein we implement empirical models to estimate the CHL in each region (Figure 1) using the selected drivers’ series (Section 2.2) as predictors. These models consisted of three multi-linear regressions (Supplementary Section S3.1) and a non-linear approach based on artificial neural networks (Supplementary Section S3.2), defined as follows.

For details about the models, see Supplementary Section S3.

Our comprehensive analysis of observational data has allowed the identification of many of the physical drivers for the CHL in the six analyzed regions, as well as the teleconnection patterns involved in the large-scale patterns. Figure 9 shows a summary of these main results of our analysis. The analyzed regions are characterized by a vigorous interannual variability; comparing the regional CHL SD before (Figure 1) and after removing the monthly climatology (Figure 9) indicates that the interannual variability has roughly half of the amplitude of the mean seasonal cycle in the tropical regions (1 to 4), but they are comparable in the southernmost regions (5 and 6). This interannual variability is a result from different drivers, herein represented by those oceanographic variables that explain a significant fraction of the observed CHL variability.

In Region 1, SSH and Sssh are variables which can serve as proxies for the regional surface circulations, the North Brazil Current (NBC) (Johns et al., 1998; Vallès-Casanova et al., 2022), the NBC retroflection (Johns et al., 1990; Richardson et al., 1994), and the anticyclonic rings that are shed from it (Richardson et al., 1994; Fratantoni and Glickson, 2002; Bueno et al., 2022; Garzoli et al., 2003) that can interact with other eddies (Subirade et al., 2023). Surface mean circulation is responsible for horizontal export of nutrients (Santos et al., 2008), while eddy features can have different local effects on the CHL (e.g., Pereira et al., 2019). These effects are caused by eddy-wind driven upwelling, vertically transporting nutrients to the surface levels and hence modifying the phytoplankton growth (McGillicuddy et al., 2007; Pereira et al., 2019), or by eddy-interactions which can modify the biogeochemical properties of the water column and shape the CHL distribution (Hernández-Carrasco et al., 2020). Other transient motions like coastal-trapped waves (CTW), both intraseasonal (Imbol Koungue and Brandt, 2021) and interannual (Bachèlery et al., 2020), would affect the local phytoplankton growth. This effect can be by upward nitrate flux induced by nitrocline displacements combined with intense vertical mixing (Körner et al., 2024), or by enhanced CTW-induced along shelf flows (e.g., Rivas, 2017) that might export nutrients and CHL to other zones. Also, the wind fluctuations that generate the CTW might have an enhancing effect of surface seawater turbulence favorable to the phytoplankton (e.g., by turbulent supply of nutrients; Hales et al., 2009). Precipitation rate (TPR) can impact CHL (r = −0.21) through modification of the freshwater input and runoff from land. This runoff carries a large amount of nutrients essential for phytoplankton growth (Santos et al., 2008), and also an important amount of suspended particles that limit the incoming light necessary for the phytoplankton (DeMaster et al., 1986), in addition to the SWR reduced by the prevalent cloudiness that comes with precipitation (e.g., McQuate and Hayden, 1984; Chaves and Cavalcanti, 2001; Yan, 2005).

In Region 2, which is part of the four major Eastern Boundary Upwelling Systems (e.g., Quiñones, 2010), the wind-driven offshore Ekman transport generates a divergence over the shelf (e.g., Trautman and Walter, 2021) and a sea-level decrease near the coast (e.g., Kruse and Huyer, 1983) while cold and saltier water from deeper levels compensates the offshore transport. Dynamical features linked to the upwelling front, like an equatorward alongshore flow (e.g., Allen, 1973) would affect the nutrient distribution, and the mesoscale structures shed from it can increase the CHL (e.g., Hernández-Hernández et al., 2020), especially the cyclonic structures due to a shoaling of the isopycnals and the nitrocline (e.g., Hernández-Hernández et al., 2020), and hence the subsurface chlorophyll-a maximum (e.g., Almazán-Becerril et al., 2012). Rainfall events can interact with upwelling intrusions to produce temporal variations in phytoplankton biomass (Valentin et al., 2021), and cross-shore winds can induce a flux of nutrients and micronutrients from the Saharan desert (Mills et al., 2004; Mahowald et al., 2009; Shi et al., 2022), modulated by the passage of African easterly waves (Bloomfield et al., 2022; Nakamae and Shiotani, 2013).

In Region 3, the atmospheric circulation is influenced by the ITCZ (e.g., De Albuquerque Cavalcanti, 2015) which modulates the wind field, rainfall and solar radiation thus affecting the regional phytoplankton growth. SSH is a proxy for the driving effect by the NBC and the North Brazil Undercurrent (Dossa et al., 2021), together with a weak influence of the Amazon River and the Saõ Francisco River (Rodrigues Holanda et al., 2021).

In Region 4, SST is linked to atmopsheric patterns like the ITCZ which creates strong convective activity, characterized by atmospheric-pressure reductions and wind convergence at surface that generates the thunderstorms that characterize the tropical regions (Huntley, 2019). These storms can lead to a deepening of the mixed-layer depth and a dilution of the deep chlorophyll-a maximum (Barrillon et al., 2023).

In Region 5, cross-shore velocity Uo represents the circulation variability dominated by the Malvinas Current (MC) in the south and the Brazil Current (BC) in the north, which converge near 34^°S (Piola et al., 2018; Artana et al., 2021), and the mesoscale dynamics enhanced as a consequence of the MC-BC interation (+EKE, +CHL) (as in Frey and Kubryakov, 2023). These hydrodynamic features have an influence on the regional transport of nutrients, especially those originated in the La Plata River’s plume and the lagoon systems (Niencheski, 2018). Tx represents the variability of cross-shore winds, responsible for an atmospheric input of limiting micro-nutrients (i.e., iron) from land, mostly from the Patagonian desert (Johnson et al., 2010), and which are driven by the winds associated to high and low pressure systems characteristic of South Brazil and the nearly-weekly southwesterly cold fronts (Rodrigues et al., 2004; Wahrlich et al., 2018).

In Region 6, shortwave-irradiance SWR variability is driven by the characteristic fog and low-level clouds in the Namib area (Andersen et al., 2019, 2020), and this cloudiness would produce some precipitation. Tx represents the cross-shore winds which plays an important role transporting nutrients from desert areas to the ocean surface, in the form of dust plumes linked to salt pans and dry river valleys in the Namib desert (Eckardt and Kuring, 2005). Windblown dust within the river valleys is easily transportable offshore from the Namib over the Benguela Upwelling System (Dansie et al., 2017c), containing great concentrations of bioavailable iron and enriched nitrogen and phosphorus, which may positively impact marine productivity in the southern Atlantic (Dansie et al., 2017a, b).

SST anomalies in the tropical Atlantic (TSA and TNA; Figures 5a–c) have a marked effect on the rainfall over northeastern Brazil (and the Amazon basin) through the modulation of the latitudinal position of the ITCZ (Nobre and Shukla, 1996; Hounsou-Gbo et al., 2019). The El Niño Southern Oscillation (ENSO) also has a significant effect on the rainfall in the western near-equatorial Atlantic; along the ENSO cycle (Figures 6a–c but for SOI < 0), the ascendant branch of the Pacific Walker cell shifts to the east and the Pacific-Atlantic cell shows enhanced subsidence over equatorial South America and the western Atlantic, causing a more stable atmosphere which reduces the moist convection (Rodriguez-Fonseca et al., 2023). Marine resources in the northwest African coast are well affected by the ENSO variability, for example, the El Niño has been used as a predictor of the sardinella population biomass and distribution in the central-southern portion of the Canary Current Upwelling System (López-Parages et al., 2020).

The eastern Atlantic off Mauritania exhibits intense SST interannual varability associated with the Dakar Niño/Niña involving ocean-land-atmosphere coupled processes (Oettli et al., 2016; Koseki et al., 2024), which would have their signature in variables like SSH and others (consistent with the ranked variables in Figure 3). On the eastern tropical Atlantic, positive Dipole Mode Index (DMI) episodes (Figures 5d–f) drive tropospheric circulation changes that lead to an increase in moisture convergence and convection over the Congo basin and an increase in Congo River discharge (Jarugula and McPhaden, 2023), as evidenced in the reduced salinity in that region (Figure 5f). In the western south Atlantic, salinity and CHL anomalies show significant correlation with the Atlantic Multidecadal Oscillation (AMO) (Pacheco et al., 2022), apparently associated with an onshore advection of moist and relatively-warm air.

On the other hand, winds are responsible for upper-ocean nutrient-enrichment processes like coastal upwelling, a determinant factor for phytoplankton growth, and their changes can modify the spatial distribution (and even stability) of the upper-ocean circulation, and hence affect the transports of nutrients and phytoplankton. The wind-field patterns can be also modified by large-scale teleconnection mechanisms. For example, near-equatorial warmings (TSA and TNA) produce a weakening of the trade winds (Nobre and Shukla, 1996; Hounsou-Gbo et al., 2019) (Figure 5b). Also, there is a weak but significant correlation between the ENSO and the Brazil Current (BC) SST (with an 8-month lag when ENSO leads) in austral summer, linked to changes in upwelling-favorable winds and the intensity of the trade winds, associated with changes of the South Atlantic anticyclone intensity and position (Rouault and Tomety, 2022). The ENSO also induces a subsidence over the equatorial western Atlantic, which produces a weakening of the trade winds (Figure 6b for SOI < 0) and a reduction of the upwelling in that area (Rodriguez-Fonseca et al., 2023). These wind-field variations, specifically the wind-stress curl (Fonseca et al., 2004), have a modification effect on the surface equatorial currents (Figures 5a, 6a); the North Brazil Current (NBC) also presents a multidecadal variability (Tretkoff, 2011). Detachment of anticyclonic eddies from the NBC retroflection occurs during periods when this retroflection is weakened (Didden and Schott, 1993), but these eddies’ formation mechanism has been associated with short equatorial Rossby waves, generated by wind variations in the tropical Atlantic (Ma, 1996). The detached eddies can interact with the Amazon River plume, resulting in a CHL increase (Fratantoni and Glickson, 2002), and responsible for the transport of water properties to the northwest (Richardson et al., 1994; Fratantoni and Richardson, 2006).

South of the tropical band, the BC transport is significantly correlated with the sea-level pressure (Schmid and Majumder, 2018) and the wind-stress curl in the western Atlantic (Goes et al., 2019) (Figures 5b, e, 6b–e). The BC transport is also correlated with the ENSO and teleconnection indices not considered in our analysis (namely, the South Atlantic Subtropical Dipole Mode and the Southern Annular Mode) (Schmid and Majumder, 2018). These BC’s variations must modulate the formation of eddies in the Brazil-Malvinas confluence zone (Frey and Kubryakov, 2023), which can be detected by their enhanced CHL levels (Garcia et al., 2004).

On the subtropical eastern coast (around Region 6), Andersen et al. (2020) report that synoptic-scale patterns control the fog and low-level cloud variability off Namibia. They conclude that the fog/clouds are formed by an increased longwave cooling under the dry anomaly close to the coast, and an onshore flow anomaly of marine boundary-layer air masses, modulated by coastal winds and a heat low over southern Africa (generated by greenhouse warming by moist air masses and northerly warm air advection). Our results suggest that the ENSO may modulate the interannual variability of this synoptically-controlled phenomenon. During La Niña conditions (SOI > 0) three dynamical features occur which resemble the conceptual model for April-June proposed by Andersen et al. (2020): a strong synoptic-scale SLP anomaly southwest of Namibia (Figure 6b), a smaller-scale low-pressure system next to the Namibian coast (Figure 6b), and a warming over the Angola-Benguela front (Figure 6a). Thus, this suggests that an El Niño episode would decrease (with a lag of 4–6 months) the occurrence of fog and low-level clouds off Namibia, especially during the austral fall, whereas a La Niña episode would have an opposite effect; this is consistent with the significant correlation between the SOI and this Region’s precipitation TPR (r = 0.18). Fog is an essential component of Namib-region ecosystems (Andersen et al., 2019); as described in Section 3.1.6, the fog and low-level clouds modulate the incoming solar light, limiting the regional phytoplankton growth and the CHL.

A relevant application of the identification of the main driver-linked series for the regional CHL is using these series to reconstruct the CHL series by different model approaches. Each CHL-reconstruction approach adopted in this paper has its merits. The multi-linear regression models are able to reproduce significant fractions of the observed CHL variance (especially lower frequencies), with the climate-index model showing poorer results and the SST-EOF model showing better results, but the oceanographic-variable model is the most physically sound. Although the climate-index model shows a lower performance to reproduce the regional CHL (Table 3, Figure 7), it can be easier to evaluate (if the indices are available) compared to the data processing required to the use of the satellite/reanalysis data.

The SST-EOF model, on the other hand, not only shows the best results among the multi-linear models (Table 3, Figure 7) but it requires the processing of one single variable (SST); however, the validity and physical meaning of the higher EOF modes remain to be clarified. Applying a significance criterion like the “N-test” used in our analysis (Supplementary Section S1.3) is most probably not enough to guarantee that the last selected modes are indeed physically sound, but dynamical criteria are necessary to determine this notion, which is out of the scope of this paper. None of the selected EOF modes was used in all the regional models, only two of them, 1 and 15, were used in five of the six regions (Figure 8). EOF-mode 1 is consistent with the SST-anomaly induced by the Tropical Northern Atlantic Index (TNA) and in lesser extent that by the Tropical Southern Atlantic Index (TSA), with the characteristic warming in the Tropical North Atlantic, but with more clear warmings off Benguela and off La Plata River (not shown). EOF-mode 15 shows a different pattern, with an intense signal along the South American eastern coast, from Region 3 to Region 5 where it is stronger, and weaker signals over the southern portion of Region 2 and over the northern portion of Region 6, as well as a signal of opposite sign off Region 4 (not shown).

The neural-network model shows a high performance to reproduce the CHL variability (Figure 7), even if only one predictor variable (Regions 2 and 4) is used. However, deducing the relative importance of one specific input variable (when more than one predictors are used) is complicated since there is no explicit way to clarify the relation-ship between inputs and outputs in these models, which causes difficulty in interpreting the results (Torres-Faurrieta et al., 2016).

Previous papers have also used non-linear CHL-reconstruction models. For example, Martinez et al. (2020a, 2020b) used a machine learning approach (a non-linear statistical approach based on Support Vector Regression) to reconstruct global spatiotemporal CHL variability from selected surface oceanic and atmospheric physical parameters taken from a numerical model. Using a 13 year-long training period, they skillfully reproduced the CHL variability from a 32-years (1979-2010) global physical-biogeochemical simulation, and their reconstructed CHL more accurately reproduced some aspects (e.g., El Niño signature in the tropical Pacific and Indian Oceans) of the satellite-observed CHL variability and trends compared to the model simulation.

The CHL-reconstruction models present the possibility to reproduce CHL in periods when its observation is unavailable. As described in Section 3.2.2, the models were used to project for two periods beyond the fitting/training (1998-2021) period: 1993–1997 and 2022-2023. The comparisons between the projected CHL-reconstructed series and the satellite-observed CHL, and between each other, show that the multi-linear models have a performance comparable to that for the fitting period, but the index-based model and the EOF-based model show a slightly poorer performance compared to the model based on oceanographic variables (Figure 7). The non-linear model (NARX), on the other hand, shows nearly-perfect correlations in the training period (1998-2021) but shows results slightly poorer than those obtained with the multi-linear models (Figure 7). This issue occurs because the neural-network presents a non-linear response during the unconstrained period which produce significant deviations from the observed CHL, which might make this model inconvenient to estimate long-term CHL series. Then, all these results point to a potential use of a model based on oceanographic variables to reproduce and even project (at least for 2–5 years beyond the fitting period if it is this period is long enough), the regional CHL.