The ocean’s general circulation is simulated using the open-source “Nucleus for European Modelling of the Ocean” (NEMO v.4.2.2; Madec and the NEMO System Team, 2023). This three-dimensional, finite-difference ocean model solves the Navier-Stokes primitive equations and employs a non-linear equation of state that couples active tracers with fluid velocity, following hydrostatic and Boussinesq assumptions.

The physical component of the modeling system is constituted by a global ocean configuration that utilizes a global ORCA tripolar grid (Madec and Imbard, 1996) with a nominal horizontal resolution of 1/4° (~27.75 km) at the equator, refining to approximately 14 km at 60° latitude. A regional domain representing the BUS is nested within the global configuration using the “Adaptive Grid Refinement in Fortran” tool (AGRIF; Debreu et al., 2012). This regional domain achieves a higher horizontal resolution of 1/16° (~7 km) and extends from 15.71°S to 44.89°S and from 7.03°W to 27.28°E ( Figures 1A, B ). A two-way nesting method is adopted, enabling the regional configuration to receive lateral open boundary conditions from the global one, whose physical fields are then improved by means of the regional high-resolution oceanic conditions. Herein, the global configuration is referred to as ORCA025L75 and the regional configuration as BENGUELA16. Both configurations share the same vertical z-coordinate system, consisting of 75 layers ranging from a few meters near the surface to approximately 200 meters in the deep ocean. Data from the General Bathymetric Chart of the Oceans grid (GEBCO; GEBCO Compilation Group, 2021) were used to generate the bathymetry of each configuration.

Vertical mixing is parameterized using the turbulent kinetic energy scheme (Blanke and Delecluse, 1993), while the advection of active tracers (temperature and salinity) is computed using a total variance dissipation scheme, and the lateral mixing is applied through a Laplacian formulation along iso-neutral surfaces. Bottom friction is modeled using a quadratic formulation. In the global configuration, additional parameterizations account for the subgrid-scale vertical mixing and the diffusive bottom boundary layer. Momentum and heat exchange at the air-sea interface are simulated through the ECMWF bulk parameterization scheme (Brodeau et al., 2017), while near-surface atmospheric forcing data were obtained from the ERA5 reanalysis dataset (Hersbach et al., 2020), which has a horizontal resolution of 31 km and an hourly frequency. A spatially explicit, time-invariant geothermal heat flux is prescribed at the seafloor using data from Goutorbe et al. (2011). Riverine freshwater inputs are derived from the Global Flood Awareness System reanalysis (GloFAS v2.1; Alfieri et al., 2013), setting as inputs values zero salinity and the same temperature of the corresponding ocean grid-point.

A spin-up simulation is conducted by repeating the forcing conditions from 1980 over 15 annual cycles. The initial three-dimensional temperature and salinity fields are obtained from the 1981–2010 climatology of the World Ocean Atlas 2018 (Locarnini et al., 2018; Zweng et al., 2019). Temperature and salinity are not restored toward surface reference fields. The time integration of the ocean model is performed with a timestep of 1200 s for the global configuration and 400 s for the regional one. The final physical state at the end of the spin-up period serves as the initial condition for a 41-year hindcast simulation, covering the period from January 1980 to December 2020.

The biogeochemical component is provided by the “Biogeochemical Flux Model” (BFM v.5.3; Vichi et al., 2023). BFM is an open-source, community-developed, generic biomass-based marine biogeochemical model. It represents the system in Eulerian coordinates, simulating the pelagic dynamics of the marine ecosystem through a suite of chemical and biological processes. The pelagic-benthic coupling is implemented via a simple benthic closure, where fixed quotes of the sedimented organic matter are remineralized and returned to the water column as nutrients.

The model organizes elemental components—carbon, nitrogen, phosphorus, silicon, iron, and oxygen—into chemical functional families, forming the structural elements of both living functional groups and non-living organic matter. Biogeochemical cycling is solved independently for the following dissolved components: nitrate (NO₃ ^-), ammonium (NH₄ ⁺), phosphate (PO₄ ^3-), silicate (SiO₄ ^4-), iron (Fe), and oxygen (O₂). The ocean carbonate system is modeled with dissolved inorganic carbon (DIC) and total alkalinity as prognostic state variables, following the best practices outlined by Orr et al. (2017). The lower-trophic level ecosystem interactions represented in this study by the BFM can be considered of intermediate complexity (see, e.g., Lovato and Butenschön, 2023). The ecosystem structure accounts for five living functional groups: heterotrophic bacteria, two phytoplankton groups (diatoms and nanoflagellates), and two zooplankton groups (microzooplankton and mesozooplankton). A schematic diagram illustrating the relationships and interactions between the different model components is provided in Figure 2 .

The coupling between the physical and biogeochemical components of the modeling system is achieved through an offline procedure, utilizing high-frequency physical fields previously computed for the BUS regional configuration. These data are processed by the NEMO interface “Tracers in Ocean Paradigm” (TOP; Aumont et al., 2019) to resolve the advective and diffusive transport of biogeochemical tracers, while reaction terms computed by the BFM are integrated using a source-splitting method (see details in Lovato and Butenschön, 2023). The advection of passive tracers is solved using the monotonic upstream scheme for conservative laws (MUSCL), while lateral diffusion follows the same approach used for physical active tracers. A key advantage of the offline coupling procedure is the ability to use a longer integration timestep which was set to 900 s.

Similar to the physical component, a spin-up simulation of marine biogeochemical quantities was conducted for 15 years by repeatedly applying the physical circulation fields of the year 1980. Initial conditions for O₂, SiO₄ ^4-, PO₄ ^3-, and NO₃ ^- were derived from the World Ocean Atlas 2018 (WOA18; Garcia et al., 2019a, 2019b), while total alkalinity and dissolved inorganic carbon concentrations were sourced from the climatological gridded dataset of the Global Ocean Data Analysis Project (GLODAP v2; Lauvset et al., 2022). All other variables were initialized with uniform values in the upper ocean (0–100 m), as described in Orr et al. (2017). Climatological river nutrient loads were prescribed using fields produced by the Global Nutrient Export from Watersheds 2 model (NEWS2; Mayorga et al., 2010). Lateral open boundary conditions for biogeochemical variables were obtained from monthly fields generated in the historical and ssp245 experiments by the CMCC-ESM2 model (Lovato et al., 2023a, 2023b).

After the spin-up simulations, the regional BUS configuration was used to simulate the evolution of marine biogeochemistry over the period 1980–2020.

A comprehensive set of multiplatform observational datasets, including global gridded products and local in-situ measurements, was used to assess the performance of the coupled modeling system.

Global scale gridded datasets employed in the validation of marine physics accounted for three main seawater parameters, namely temperature, salinity, and near-surface current speed. Satellite-derived data products from the European Space Agency-Climate Change Initiative (ESA-CCI) comprised the sea surface temperature (SST, Merchant et al., 2019) at 5 km of spatial resolution over the period 1982–2020 and the sea surface salinity (SSS, Boutin et al., 2021) covering the years range 2010–2019 with a nominal resolution of 25 km. Reference fields of the near-surface (15 meters depth) current velocities were obtained from the Ocean Surface Current Analyses Real-time (OSCAR) from 2012 to 2020 and with a third of degree (~37 Km) horizontal resolution (ESR, 2009). The subsurface ocean validation of three-dimensional temperature and salinity fields was addressed using the multi-observation global ocean ARMOR3D analysis at a resolution of 25 km (CMEMS, 2022), which merges in-situ and satellite observations over the 1993–2020 period using multiple linear regression methods that relate surface and subsurface data (Guinehut et al., 2012).

The validation of the biogeochemical component focused on dissolved inorganic nutrients (NO₃ ^- and PO₄ ^3-), dissolved oxygen, and Chl-a from either gridded climatologies or in-situ datasets. Annual climatological data of dissolved oxygen, NO₃ ^-, and PO₄ ^3- were sourced from the CSIRO Atlas of Regional Seas (CARS, 2009), at half a degree (~55 Km) of horizontal resolution and over 79 depth levels, with vertical intervals ranging from 5 m near the surface to 250 m at deep layers. Upper ocean Chl-a concentrations were obtained from satellite observations based on the Copernicus-GlobColour processor, using the Garver-Siegel-Maritorena (GSM) algorithm, covering the period 1998–2020 at a resolution of 4 km (CMEMS, 2023).

In-situ observations were retrieved from different sampling cruises and monitoring platforms, as illustrated in Figures 3A, B . Dissolved inorganic nutrients (NO₃ ^- and PO₄ ^3-) profiles come from two oceanographic cruises. The RV James Cook cruise (JC032) held in April 2009, conducted at a latitudinal transect near 24°S, extending longitudinally from 6°W to 14°E, including 23 vertical profiles, each with approximately 24 depth levels from the surface down to 5500 m depth (Schuster et al., 2013). The RV Africana cruise (Afr258) in December 2009 accounts for 16 profiles along the Benguela coastal zone between 17°S and 33°S (Rixen and Flohr, 2012; Flohr et al., 2014).

Chl-a vertical distribution was assessed using 325 profiles collected by Biogeochemical Argo (BGC-ARGO) float number 6901582 between December 2014 and September 2018 in the open waters of the southern Benguela subregion, extending from the surface to 1000 m of depth (Wong et al., 2020). Daily dissolved oxygen observations were derived from a high-resolution glider deployed from February to June 2018 in the NBUS (Lovecchio et al., 2022). The glider followed a ~12 km triangular path and profiled to 1000 m depth with ~20 cm vertical resolution, surfacing 5–6 times daily and sampling primarily on upward dives.

A brief description of the main physical results from the global ocean configuration is provided here to complement the overall modeling system assessment, with a focus on the reliability of the simulated fields in the South Atlantic region surrounding the Benguela upwelling area.

The global ocean configuration ORCA025L75 demonstrates satisfactory skill in reproducing the annual climatology of SST patterns when compared to ESA-CCI satellite observations ( Figures 4A–C ). The model accurately captures the horizontal distribution of SST, with a spatial mean bias of about 0.02 °C. However, localized positive biases of approximately 3 °C are evident along western boundary current regions, where strong horizontal circulation patterns exhibit substantial variability. Seasonal patterns of SST climatological fields are also fairly well represented by the model ( Supplementary Figure S1 ), with mean seasonal biases ranging from 0.11 °C in austral summer (DJF) to -0.05 °C in austral spring (September to November, SON). A quantitative evaluation of the 3D temperature climatology is conducted using a normalized Taylor diagram ( Figure 5A ), compared against ARMOR3D data. Both annual and seasonal climatological values consistently cluster below 0.40 RMSD and show remarkably high positive correlation coefficients exceeding 0.95 across all depth ranges. The close agreement between ORCA025L75 and ARMOR3D is further supported by the proximity of normalized STD values, indicating the model effectively captures the amplitude of variability in observed data.

The comparison of simulated SSS annual climatology with ESA-CCI satellite observations ( Figures 4D–F ) also reveals a rather good agreement, with a spatial mean bias of approximately -0.03 PSU and seasonal differences ranging from -0.04 to -0.02 PSU in JJA and DJF, respectively ( Supplementary Figure S2 ). Nonetheless, the comparison highlights both positive and negative SSS biases near river mouths. In the vicinity of the BUS, the SSS bias remains minimal, within ±0.50 PSU ( Figures 4D–F ). The normalized Taylor diagram ( Figure 5B ) shows that values are more scattered in near-surface layers than at depth, with RMSD values below 0.50 and generally high positive correlations (> 0.80) in most instances.

Annual climatological fields of near-surface (15 m depth) horizontal current speeds, as simulated by ORCA025L75, are compared with OSCAR analysis data ( Figures 4G–I ). Mean biases remain well below -0.01 m/s, and the model effectively reproduces the main patterns of oceanic circulation. However, notable deviations are observed primarily in the equatorial band and along western boundary currents, consistent with the SST bias patterns.

Overall, the mean oceanic physical conditions simulated by the global component of the modeling system are characterized by relatively modest biases in the South Atlantic region surrounding the Benguela upwelling area, helping to provide improved boundary conditions.

The long-term trends of sea surface temperature and upper ocean Chl-a concentrations ( Figure 13 ) have been estimated from the simulated data over the BENGUELA16 domain and the satellite-based datasets ESA-CCI and Copernicus-GlobColour (see Section 2.3 for details). Data were averaged longitudinally over a 1° band from the coastline to retain the upwelling activity zone and scaled to express decadal rates of change. We used linear least-squares regression to obtain the trends at each grid point over the latitudinal range. For the trend analysis, only statistically significant trends with p-values < 0.05 are considered.

In general, the trends derived from the BENGUELA16 SST data show strong coherence with those from the ESA-CCI dataset ( Figure 13A ). A cooling trend of -0.12 °C per decade is identified between 20°S and 26°S, with the model indicating a broader seasonal span from August to October. This negative trend extends into the central Benguela subregion. However, in this area, the two datasets show contrasting signs in the trends during the austral summer. Both datasets indicate cooling across the southern subregion during the austral summer, though ESA-CCI reveals a more intense and seasonally extensive cooling, with trends reaching up to -0.64 °C per decade, compared to a maximum of -0.24 °C per decade in the model.

The BENGUELA16 simulation also captures three statistically significant warming patterns: i) in the Angola-Benguela Frontal Zone (ABFZ, 15.7°S–18.7°S) during late austral spring and summer (October to January), with values up to 0.40 °C per decade; ii) in central Benguela (24.5°S–27.5°S), a weaker warming of about 0.25 °C per decade is observed in January; iii) in southern Benguela (30.7°S–34.7°S), a more modest warming trend of 0.15 °C per decade extends from austral autumn through spring (May to November).

For upper ocean Chl-a trends, the comparison between BENGUELA16 and CopernicusGlobColour reveals several similarities, particularly in regions with positive trends ( Figure 13B ). Both datasets show statistically significant increases in Chl-a exceeding 0.80 mg/m³ per decade during late austral winter and spring (August to October), especially over southern Benguela (27°S–29°S). BENGUELA16 also shows a positive trend at 29°S during austral autumn, although with lower intensity (0.20–0.40 mg/m³) compared to satellite estimates (0.60–0.80 mg/m³). Differences between the datasets become more pronounced in the central and southern Benguela during austral summer (February), while the positive trend emerging in April in satellite data also appears in the model, though delayed by one month. Additionally, Copernicus-GlobColour shows a small but statistically significant negative trend over the ABFZ during austral summer, which is not captured by the model.

Overall, the BENGUELA16 regional model shows good skill in capturing the spatio-temporal patterns and intensity of SST trends, and to a lower extent, those of Chl-a concentrations.

As illustrated in Figure 4 , the global domain physical component exhibits systematic biases for both SST and SSS in correspondence with the western boundary currents, namely the Gulf Stream and the Kuroshio Current—similar to findings from previous modeling studies (Bryan et al., 2007; He et al., 2014; Renault et al., 2019; Liu et al., 2023). As shown by Iovino et al. (2023), the model’s horizontal resolution poses certain limitations on its ability to reproduce the spatial variability of these highly dynamic systems. Along with grid resolution, certain model parameterizations also contribute to the emergence of such biases. Kamenkovich et al. (2000) demonstrated that horizontal viscosity plays a critical role in shaping the dynamics of western boundary currents. Furthermore, Renault et al. (2019) showed that incorporating air-sea feedbacks, particularly through the eddy killing mechanism, can help mitigate persistent biases in simulations of western boundary currents. Although the model is forced with high-frequency GloFAS runoff data, local positive and negative surface salinity anomalies remain evident along coastal areas. This is likely related to the use of zero-salinity inputs, which may underrepresent the real runoff characteristics and would benefit from the imposition of ad-hoc salinity values at river outlets.

The physical component of the BENGUELA16 configuration is characterized by positive salinity biases along its southern boundary ( Figures 6D–F ; Supplementary Figure S4 ), indicating a deficiency in reproducing freshwater masses entering the BUS from the Southern Ocean. These waters typically feed into the South Atlantic gyre and introduce fresher water along the South African coast (Beal et al., 2011). This bias may also be linked to limitations in representing the northward displacement of the South Atlantic Subtropical Front, which influences both the influx of low-salinity water and the zonal extent of the Agulhas Retroflection.

Marine biogeochemical variables show varying levels of agreement with the selected observational datasets (see Section 3.2.2). The comparison of NO₃ ^- and PO₄ ^3- with both gridded and in-situ data indicates that biases mainly affect the surface and subsurface layers. This discrepancy can be attributed to an intense biological consumption in the model, resulting in lower nutrient concentrations. Additionally, the offshore advection of newly upwelled waters may rapidly transport nutrients away from the coast, reducing their residence time and limiting phytoplankton growth. This mechanism is also linked to the slight overestimation of near-surface current speeds compared to OSCAR data, especially over the narrow coastal shelf ( Figures 6G–I ).

Dissolved oxygen biases along the western and southern boundaries of the regional domain, particularly between 500 and 1000 m depth ( Supplementary Figure S5 ), are primarily attributed to lower accuracy in the boundary condition data used to prescribe oxygen influxes from the Southern Ocean. The comparison with glider measurements in the northern Benguela subregion ( Figure 10 ) indicates that simulated oxygen profiles generally match observations in the upper and deeper layers, but the model underestimates concentrations by about 30 mmol O₂/m³ at intermediate depths around 200 m. This difference likely results from a combination of unrealistic oxygen levels in the Southern Atlantic Central Water entering the region (see Mohrholz et al., 2008; Schmidt and Eggert, 2016) and enhanced oxygen consumption due to heterotrophic bacterial remineralization of organic matter.

Coastal Chl-a values within the first 100 km of the upwelling zone are underestimated, while concentrations in the open ocean—especially in the central and southern regions, are slightly overestimated. This is more evident in the punctual comparison of simulated and satellite-derived Chl-a fields ( Figure 11B ). In southern Benguela, the pattern of over- and underestimated values near the coast deviates notably from the ideal fit, while open ocean data are in better agreement. Nevertheless, the spatially averaged annual cycle ( Figure 11C ) shows that the model adequately captures the seasonal variability of algal blooms, especially in the NBUS, where Chl-a concentrations peak in August, December, and April, as previously reported by Louw et al. (2016).

However, comparing model outputs with satellite data in optically complex coastal waters remains challenging. In fact, satellite measurements of spectral reflectance can be influenced by factors unrelated to phytoplankton abundance, such as the presence of colored dissolved organic matter (Siegel et al., 2005; Aurin and Dierssen, 2012) and suspended particulates that alter light attenuation (Chauhan et al., 2002; Gohin, 2011; Theenathayalan et al., 2022). As a result, satellite retrieval algorithms are subject to greater uncertainty in these regions and may sometimes overestimate Chl-a concentrations (Coelho et al., 2017; Poddar et al., 2019).

The comparison of long-term mean state physical fields in the global and regional domains of the BUS modeling system shows that biases are generally modest, with a good reproduction of climatological physical patterns. In the validation of the ORCA025L75 global domain ( Figure 4 ), temperature and salinity showed realistic spatial patterns, and a rather low bias was obtained for near-surface current speed fields. The global configuration also provided a reliable representation of the seawater temperature and salinity fields in the subsurface ocean at most of the selected depth levels ( Figure 5 ).

The reproduction of physical fields in the BENGUELA16 regional domain was remarkable, as in the comparison with satellite observations ( Figure 6 ) and ARMOR3D analysis data at different depths ( Figure 7 ), where the obtained biases were comparable and in some cases lower than previous similar modeling applications (Veitch et al., 2009; Gutknecht et al., 2013). The regional configuration effectively captures the horizontal and vertical spatial patterns of seawater physical parameters as determined by the Benguela coastal upwelling dynamics. The low SST values obtained along the Namibian shelf ( Figure 6A ) indicate the ascent of cold waters from deeper layers, with an extension of the coastal upwelling plume comparable to that observed by satellite surveys. Along the vertical, the simulated upward tilting of isolines with the consequent outcropping of deeper water compares well with ARMOR3D analysis fields ( Figure 8 ), particularly over the narrow coastal shelf despite some minor underestimations. Near-surface current speeds produced by the Benguela regional configuration are also largely coherent with those available from the OSCAR Analysis ( Figures 6G–I ) and with the circulation patterns described in previous oceanographic studies (Stramma and Peterson, 1989; Lutjeharms, 2007; Beal et al., 2011; Tim et al., 2018; Majumder and Schmid, 2018). The model clearly shows the northwestward flow of the Benguela Current around 27°S that is fuelled by the Southern Atlantic and Indian Oceans via the South Atlantic Current. Furthermore, the BENGUELA16 physical component consistently simulates the Agulhas retroflection at about 16°E, which generates the Agulhas Return Current directed toward the Indian Ocean subtropical gyre and the Agulhas leakage’s mesoscale eddies propagating into the South Atlantic gyre (Reason et al., 2003).

The long-term mean fields of NO₃ ^- and PO₄ ^3- exhibited reasonably good skills at the different vertical layers ( Figure 7 ), as the BENGUELA16 data don’t significantly deviate from the mean state representation of the CARS climatology. Conversely, the evaluation against in-situ data from the JC032 and Afr258 cruises showed a better agreement ( Figure 9 ). In particular, the comparison of longitudinal sections with JC032 data shows a good representation of nutrient concentration isoline deformation caused by the coastal upwelling ( Figures 9C, D ). Dissolved oxygen long-term mean fields of BENGUELA16 showed to be consistent with CARS climatology, particularly within the biologically active upper layers (0–200 m), while oxygen levels were underestimated at depths over the central and northern Benguela shelf ( Supplementary Figure S5 ). However, the model successfully captured key DO dynamics, including the distribution of hypoxic waters (< 60 mmol O₂/m³), particularly in the northern Benguela, where the Southern Atlantic Central Water transports oxygen-depleted waters via the poleward undercurrent (Monteiro and van der Plas, 2006; Mohrholz et al., 2008; Lass and Mohrholz, 2008). The validation of simulated upper-ocean Chl-a concentrations against satellite data shows how the model is capable of capturing the upwelling-driven response of phytoplankton, as the onshore-offshore Chl-a gradient observed in Copernicus-GlobColour data is well represented ( Figure 11A ), with additional skills in the reproduction of the whole seasonal cycle ( Figure 11C ). Moreover, the comparison with Chl-a vertical profiles collected by the ARGO-BGC float indicated the satisfactory skills of the model in replicating the depth and structure of the vertical Chl-a distribution over space and time.

Overall, the coupled modeling system effectively captures the key physical and biogeochemical dynamics of the Benguela Upwelling System, and it may be further exploited to support a range of applications aimed at monitoring system dynamics. From a physical perspective, different indices can be implemented to assess the variability of upwelling intensity by exploiting either SST changes along the Namibian shelf (Chen et al., 2012; Santos et al., 2012) or the simulated current velocities to estimate offshore Ekman transport across the different subregions (Veitch et al., 2010). Moreover, the combined use of physical and biogeochemical quantities enables the application of the Biologically Effective Upwelling Transport Index (BEUTI; Jacox et al., 2018), which evaluates nitrate flux into the photic zone. The good performance of the modeling system in representing dissolved oxygen across various depths also makes it well-suited for studying long-term trends in minimum oxygen zones, particularly under intensified upwelling conditions (Bakun et al., 2015). Finally, model outcomes could be further utilized to monitor the evolution of the Benguela marine ecosystem and assess potential impacts on productivity, habitat quality, and fisheries potential, helping to identify oceanographic shifts that support informed, climate-resilient management of coastal areas (Chang et al., 2023; Lamont et al., 2019).

Sea surface temperature is widely recognized as a reliable proxy for assessing both temporal and spatial variability in upwelling intensity (Gutierrez-Guerra et al., 2024), due to its theoretical connection to the ascent of deep, cold water to the surface layer. The model outcomes reveal a significant cooling trend in both the northern and southern Benguela subregions, during their respective upwelling seasons. In the northern Benguela, a cooling trend is evident during late austral winter and austral spring ( Figure 13A ), which is the primary upwelling season for this subregion (Shannon, 1985; Hagen et al., 2001; Tim et al., 2015). Similarly, in southern Benguela, the strongest cooling trend occurs during its upwelling season, namely from December to May (Rixen et al., 2021; Rouault and Tomety, 2022).

When comparing the simulated SST seasonal trends with those from ESA-CCI satellite data, the model demonstrates good skill in detecting cooling patterns in northern and southern subregions, particularly during upwelling seasons. However, some discrepancies in the spatial and temporal distribution and intensity of these trends emerged, particularly in the ABFZ and within southern Benguela. These discrepancies may be explained by the model having produced a weaker water inflow from the Angola Current resulting in a lower extent of the warming trend in the NBUS, while the simulated upwelling signal in the SBUS is weakened by a stronger Agulhas warm-water intrusion during the austral spring and autumn.

Nevertheless, the model’s outcomes align with previous observational studies. For example, Tomety et al. (2024) identified a significant cooling trend in satellite-derived SST data in the northern Benguela subregion (23°S–25°S) during austral spring between 1982 and 2015, with a cooling rate of -0.10 to -0.20 °C per decade, similar to our results ( Figure 13A ).

A warming trend in the ABFZ latitudes (16°S) is consistent with findings from previous studies. Tomety et al. (2024) reported a significant warming trend during austral summer (DJF), in addition, Monteiro and van der Plas (2006) documented a warming trend in SST over the ABFZ, with values exceeding 0.60 °C per decade from 1982 to 2005. This warming trend may be linked to latitudinal shifts in its position, likely driven by the southward intrusion of the warmer Angola Current into the cooler Benguela Upwelling Zone. Vizy et al. (2018) reported a southward latitudinal shift of 0.05° – 0.55° per decade between 1980 and 2014. The ABFZ follows a pronounced semi-annual cycle, peaking between November and December (Koseki et al., 2019), which aligns with our detected warming trends in the region during these months. Variability in its position and intensity can significantly impact the coastal ecosystem, as warm and nutrient-poor intrusions may alter biological communities (Huenerlage and Buchholz, 2013).

Several studies converged on the coastal SST cooling along the BUS at the annual scale, attributing this primarily to upwelling intensification (Santos et al., 2012; Polonsky and Serebrennikov, 2021; Gutierrez-Guerra et al., 2024). Polonsky and Serebrennikov (2021) noted a cooling tendency in the Benguela coastaloceanic thermal gradient, accompanied by an intensification of the Ekman transport. Gutierrez-Guerra et al. (2024) assessed changes across all Eastern Boundary Upwelling Systems, identifying a negative trend (-0.23°C per decade) in the BUS coastal waters. The cooling trends simulated by the BENGUELA16 physical component, which moderately coincide with satellite observations during upwelling seasons, are much more likely pointing to an upwelling intensification.

The simulated surface Chl-a trends ( Figure 13B ) show statistically significant positive trends over a few locations, indicating an overall increase in upper Chl-a concentrations over the studied period (1998–2020). In particular, the comparison with Copernicus-GlobColour data showed that the model is able to detect the main positive trend in southern Benguela. However, discrepancies exist, such as the simulated positive trend over the northern Benguela during the austral spring, which is not evident in satellite observations. These may be related to inherent limitations of satellite-derived Chl-a data (see Section 4.1) or potential shortcomings in the model, indicating further required improvement in the simulation of Chl-a levels.

The model-derived positive Chl-a trends, which are potentially linked to upwelling intensification, are also consistent with previous studies. Gregg et al. (2005) noted a positive Chl-a trend along the Namibian-Angolan coast between 1998 and 2003. Verheye (2000) documented an increase in surface Chl-a concentrations between March and June in St. Helena Bay (32°S). Bordbar et al. (2021) observed a significant correlation between the spatial patterns of Chl-a and upwelling intensity from 2008 to 2018 in the whole BUS region. Similarly, Chen et al. (2012) reported a strong positive correlation between upwelling intensity and Chl-a across latitudes from 23°S to 34°S during the 2003–2008 period.

The simulated trends reveal an inverse proportionality between Chl-a and SST trends. The positive Chl-a trends in northern Benguela correspond to the cooling trend of SST found at the same latitudes and seasons. In northern Benguela, the positive Chl-a trend slightly aligns with the detected cooling, though with a one-month temporal lag ( Figures 13A,B ). This relationship supports the linkage between upwelling intensification, surface water cooling, and increasing activity of primary producers (García-Reyes et al., 2015). However, an increase in primary producers’ activity does not always directly correlate with stronger upwelling. Enhanced upwelling typically draws nutrient-rich water into the mixed layer, but our trend comparisons suggest this conceptual framework may not apply across all the Benguela subregions.

A stronger cooling trend may occur without a corresponding positive Chl-a trend, as observed in southern Benguela during austral summer and autumn. In the proposed modeling framework, this could be due to intensified upwelling activity enhancing dilution, further limiting Chl-a accumulation, or to an already saturated phytoplankton growth regulated by top-down control. While upwelling generally stimulates phytoplankton growth by transporting nutrient-rich water to the surface (Bakun et al., 2015), other factors—including the chemical composition of the upwelled water, can significantly influence primary production (Rykaczewski and Dunne, 2010; Bograd et al., 2015). Jacox et al. (2018) noted that subsurface nutrient concentrations may vary independently of local upwelling vertical transport. Flynn et al. (2020) reported that in southern Benguela, primary production is not exclusively tied to upwelling signals, as high nutrient fluxes are regenerated over the continental shelf (i.e., nutrient trapping), providing variability in primary producers independent of upwelling-derived nutrients.

Finally, intensified upwelling in the BUS may boost primary producers and consumers in surface waters but also increase organic export, enhancing microbial respiration that depletes O₂ and raises CO₂ (Gruber, 2011). This can induce hypoxic conditions over the continental shelf, threatening marine life (Grantham et al., 2004; Chan et al., 2008). The exposure of the upper ocean ecosystem to more acidic conditions during the upwelling events may have detrimental effects on marine biota (Lluch-Cota et al., 2014; Bakun et al., 2015). Moreover, changes in upwelling timing can disrupt predator-prey dynamics and affect productivity in upwelling regions (Schwing et al., 2006).