Our approach uses a triply nested configuration of the Coastal and Regional Ocean Community (CROCO; Penven et al., 2006; Debreu et al., 2012) version of the Regional Ocean Modelling System (ROMS; Shchepetkin and McWilliams, 2005) for numerical simulations. The model has 60 terrain-following vertical levels, which were scaled down from 1/4° to 1/12° and finally to a 1/36° grid. Outputs from the 1/36° grid model are used in this study because of the model’s ability to capture and resolve smaller-scale features of the shelf around the South African coastline. The CROCO simulation used for our study is that of Tedesco et al., 2019, who set up the model to investigate key sub-mesoscale processes in the Agulhas Current. The model’s realism has been validated against existing observations in Tedesco et al. (2019) and has subsequently been used in a number of studies (Schubert et al., 2021; Heye et al., 2022; Pfaff et al., 2022; Suthers et al., 2023). CROCO was forced by climatological fields of wind stress, surface temperature and humidity from global reanalysis products and other atmospheric models. Daily and monthly means of the 1/36° grid datasets are available for the period 1993 to 2014. However, for the purposes of this study, we extracted and used daily outputs for a 5-year period spanning from 2010 to 2014. The model outputs are freely available and can be accessed at http://dap.saeon.ac.za/thredds/catalog/SAEON.EGAGASINI/2019.Penven/DAILY_MEANS/1_36_degree/catalog.html.

The 1/12° grid global ocean physics reanalysis (GLORYS12) from the Copernicus Marine Environment Monitoring Service (CMEMS) is used for comparison with the model in this study. The GLORYS12 product is a global eddy-resolving physical ocean and sea ice reanalysis product, with 50 vertical levels and covering the period from 1993 to the present (Jean-Michel et al., 2021). The product is based on the framework of the real-time forecasting CMEMS nested in the Nucleus for European Modelling of the Ocean (NEMO) platform and has been driven at the surface by the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim and ERA5 reanalysis in recent years. In this product, ocean observations are assimilated by means of a reduced-order Kalman filter. GLORYS12 uses the 3D-VAR algorithm to resolve corrections for the slowly evolving large-scale biases in temperature and salinity. The product resolves small-scale variability in surface dynamics and compares well with nonassimilated data. It provides a reliable physical ocean state for climate variability and supports a variety of applications, such as seasonal forecasts and marine biogeochemistry. In this study, we extracted the daily products of GLORYS12 for a 5-year period spanning 2010-2014. The GLORYS12 product can be accessed freely from https://data.marine.copernicus.eu/product/GLOBAL_MULTIYEAR_PHY_001_030/.

The level-4 sea surface temperature (SST) product from the European Space Agency Climate Change Initiative (ESA-CCI) is used in this study to evaluate the performance of the CROCO model in representing the main physical processes around the South African coastline. The level 4 ESA-CCI product provides gap-free maps of daily average SSTs at a depth of 20 cm in the global ocean at a horizontal grid resolution of 5 km. It consists of observations from the Along-Track Scanning Radiometers (A)ATSRs), Sea and Land Surface Temperature Radiometer (SLSTR) and Advanced Very High-Resolution Radiometer (AVHRR) series of sensors (Merchant et al., 2019; Good et al., 2020). The Level-4 product is available from 1 September 1981 to October 2022 and can be accessed from the CMEMS via the following link https://doi.org/10.48670/moi-00169.

The study also uses a global 1km Group for High Resolution Sea Surface Temperature (GHRSST) Level 2P dataset based on multi-channel SST retrievals from the AVHRR on the European Meteorological Operational-A satellite which was first launched on 19 October 2006. This SST product is produced and retrieved in near real time from Metop/AVHRR infrared channels (3.7, 10.8 and 12.0 micrometer) using multispectral algorithm by the Ocean and Sea Ice Satellite Application Facility (OSI SAF). To correct the multispectral algorithm for regional and seasonal biases, atmospheric profiles of water vapor and temperature from numerical weather prediction model and radiative transfer model are applied (OSI SAF, 2014). Data access is provided free of charge via THREDDS services on Catalog https://www.ncei.noaa.gov/thredds-ocean/catalog/ghrsst/L2P/AVHRR_SST_METOP_A/OSISAF/catalog.html.

Hourly ERA5 wind data from the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis for the global climate and weather are used in this study. The hourly ERA5 wind data are available from 1940 onwards, and they are regridded to a regular lat-lon grid of 0.25 degrees. The ERA5 data combines model data with global observations to produce a complete dataset (Hersbach et al., 2023). Conversely, Ekman currents derived from the Copernicus Level 4 multiyear global total velocity product are used to characterize the mean circulation along the west coast of South Africa. The Level 4 multiyear product provides velocity fields at the surface and 15 m depth, with a 3-hour frequency and a 0.25-degree regular grid (product ID: MULTIOBS_GLO_PHY_REP_015_004). The velocity fields are generated by combining the Ekman currents simulated using the approach of Rio et al. (2014) and the CMEMS geostrophic surface currents from satellites. Products are available in hourly, daily and monthly means. In this study, the daily averaged fields from 2010 to 2014 are used. For the purposes of this study, we extracted daily average fields from 2010 to 2014 (Rouault et al., 2010).

To simulate the larval dispersal patterns of the anchovy species Engraulis capensis and sardine Sardinops sagax in South Africa’s coastal marine protected areas, the Ichthyop Lagrangian tool was used (Ichthyop v3.15; Lett et al., 2008). Ichthyop is an open-source java-based tool designed to study physical and biological drivers influencing larval dynamics throughout the water column. Ichthyop utilizes current velocity vector fields, temperature and salinity from hydrodynamic models to simulate advection and connecting pathways of virtual particles (i.e., larvae) in a three-dimensional physical environment. Ichthyop was built with the assumption that ichthyoplankton are advected and diffuse similarly to water particles and can be used in various applications in physical oceanography and marine ecology. The tool follows the particle’s location and records its associated ocean properties, such as temperature and salinity. One of the advantages of the tool is that it provides users with the option to code virtual particles with biological behavior (e.g., diel vertical migration, swimming, etc.) to assess its impact on larval dispersion.

Ichthyop allows users to run simulations offline in either hindcast (backtracking) or forecasting modes. However, in this study, particle trajectories were tracked in forecasting mode at a constant time step using a fourth-order Runge-Kutta model, which enables particle trajectory tracking with adaptive time steps without compromising computational efficiency (Van Sebille et al., 2018). In total, 140, 000 particles were released in the first 30 m of the water column every October (coinciding with the peak sardine and anchovy spawning periods) in 14 marine protected areas located within the 200 m isobath ( Figure 1 ). We released 10, 000 particles in each of the marine protected areas on the first day of October. Simulation experiments were run for 30 days (~ pelagic duration of anchovy and sardine in South Africa) and particles were tracked every half an hour (1800 seconds). The number of virtual particles released was found to be statistically sufficient to reproduce stable connectivity estimates, as similar results were obtained when using twice as many particles. We configured and parameterized both coarse and fine meshes such that particles were prevented from overlapping with land; therefore, particles were considered dead when they reached land or drifted outside the domain. This approach was used to avoid bias in our analysis. Moreover, we released particles on the surface and allowed them to be transported passively by ocean currents derived from both the high-resolution CROCO model and the coarse-resolution GLORYS products. The objective here was to test the sensitivity of particle dispersion to spatial resolution when advected passively without incorporating biological behaviors.

We constructed connectivity matrices by evaluating larval settlement (i.e., potential recruitment) every 15 days along the 30-day particle trajectories. In this study, larval anchovies and sardines were considered to have reached a suitable competency period at the age of 15 days, which is the period when larvae are suitable for settlement (Beckley and Van der Lingen, 1999; Somarakis and Nikolioudakis, 2010; Hampton, 2014). This exercise was performed every half an hour during the competency period. We considered particles to have successfully settled if they were found within a coastal marine protected area at any given time during the competency period. Essentially, connectivity matrices were constructed to summarize the number of virtual larvae exchanged between any two MPAs, as were obtained after 30 days of simulation. Moreover, from any connectivity matrix, we derived the strength of the connection from one MPA to another. The strength of the connection from MPA _i to _j corresponds to the proportion of virtual larvae released on MPA _i that then settled onto MPA _j . Row _i of the matrix corresponds to the source MPA, and column _j corresponds to the settlement MPA. A nonzero value obtained from the connectivity matrix denotes an existing connection between MPA _i and MPA _j (D’Agostini et al., 2015; Andrello et al., 2017; Belharet et al., 2020). Essentially, positive (negative) values of the connectivity matrix denote a strong (weak) connection between MPAs. Self-retention (also known as self-seeding) was also considered in this study and is defined as the number of virtual particles retained at the release MPA (i.e., origin location).

The simulation outputs of the Ichthyop were analyzed using open source Python software (v.3.11; https://www.anaconda.com/) implemented in the Windows Anaconda platform for data science as well as using R Statistical Software (v.4.0.4; R Core Team, 2021). We calculated dispersion statistics for all 5 years combined (i.e., 2010–2014) while distinguishing between the results obtained from coarse-resolution GLORYS and the high-resolution CROCO model. We calculated the following from the 5 years of simulations: proportions of virtual particles 1) retained at the release MPA, 2) those settling at adjacent MPAs other than the release site, or 3) those not retained in any MPAs. We computed trajectory density maps, connectivity matrices and residence times for particles retained in any marine protected area.

The performance of the 3-km grid-resolution CROCO model in accurately reproducing physical processes around the South African coastline is assessed by comparing it to the 5-km grid-resolution SST data derived from the European Space Agency Climate Change Initiative (ESA-CCI) product. The ESA-CCI SST product is based only on satellite observations and objective mapping, thus providing us with independent observations to which we can compare the CROCO model SSTs. The SST comparison exercise serves to ascertain the model’s realism.

There is a demonstrable consistency between the results obtained from the model, GLORYS and the ESA-CCI SST ( Figure 3A ). The results show good agreement between the model and the three altimetry products (i.e., ESA-CCI,GLORYS and GHRSST) in terms of the representation of the mean SST and the location of the Agulhas Current ( Supplementary Figure S1 ). The results also show a relatively warm bias of +1°C observed off the Cape Peninsula region in the CROCO model, which is more strongly expressed along the path of the northwards-flowing, shelf-edge Benguela jet current. The location of the jet appears to be slightly more west in the model as compared to the GLORYS product ( Figure 3B ).

Density maps of virtual particles released using a 3-km grid model (CROCO) and 8-km global reanalysis (GLORYS) are used to illustrate the dominant connecting pathways of particles that eventually reach a recruitment zone ( Figure 4 ). The main connecting pathways observed in this study are consistent with those of general circulation along coastlines. The fast, south-westwards-flowing Agulhas Current drives variability off the east coast and virtual particles in that region tend to align with the flow of the Agulhas Current and are transported southwestward.

In the northern Agulhas Current region, a high number of particles are observed in the Natal Bight and along the path of the Agulhas Current in both models. South of 32 ⁰S however, the model shows that particles released on the east coast of South Africa tend to move inshore of the Agulhas Current, whereas those in the global reanalysis (GLORYS) are either entrained by the current or drift to the outer regions of the current. Relatively large alongshore loss of particles is observed in the GLORYS as compared to the model. Relatively higher levels of energy are also observed in the Natal Bight. The core of the Agulhas Current is observed to be more inshore in the CROCO than in GLORYS. The large demonstrable differences in mean kinetic energy (MKE) and eddy kinetic energy (EKE) are observed south of the African continent.

The GLORYS show relatively sporadic (erratic) dispersal patterns to the south of the continent. However, both the model and GLORYS show a coherent high-density area of particles released over the Agulhas Bank. Interannual variability significantly influences dispersal patterns of anchovy and sardines ( Supplementary Figure S3 . The west coast is characterized by a relatively large magnitude of wind-modulated Ekman transport which appears to be influencing dispersion in the region ( Figure 5 ). High variability of particle dispersals observed off the west coast. The model shows a stronger jet at the inshore edge of the Benguela Current accompanied by rapid advection of particles to the north. As a result, the particles follow the mean jet path more leading to less retention in the model. The GLORYS show strong inshore movement of particles leading to a more pronounced local retention off the west coast.

The schematic of the alongshore connectivity matrices illustrates a generally high level of alongshore connectivity in the region ( Figure 6 ). The strong connectivity observed in Figure 6 is consistent with the general hypothesis that connectivity increases from east to west along the South African coastline. The hypothesis is premised largely on the general coastal circulation around southern Africa driven mainly by the 1) south-westwards fast-flowing Agulhas Current on the east, 2) sluggish flow on the south coast driven by the combination of winds and the mesoscale activities of the Agulhas Current, and 3) northwesterly flowing frontal Benguela jet Current off the west coast. While the general pattern of alongshore connectivity was comparable between the model and GLORYS, the strength of the connections between release sites varied significantly. For instance, the model showed a high degree of variability in terms of connectivity between the release sites on the east and south coasts. Relatively high local recruitment and self-retention were observed on the south and west coasts in both the model and GLORYS. The GLORYS show a stronger connection between release sites located on the south and west coasts. Both the model and GLORYS show a coherent upstream connectivity along the west coast with stronger connections in the area between the Cape Peninsula and Cape Columbine (i.e., northwardly from Robben Island to the Childs Bank).

The GLORYS show that particles take longer in the northern Agulhas region to move through the system ( Figure 7 ). Similarly, longer residence times (> 10 days) were observed south of the continent in the GLORYS. This is especially true over the Agulhas Bank. Particles released in the model spend on average half of the time spent by those in the GLORYS on the south coast. Large variability of residence times observed off the west coast in the GLORYS with particles spending on average 8 days within the release sites.

Figure 3A shows a remarkable agreement between the results obtained from the model and altimetry products (i.e., ESA-CCI and GLORYS). The results are further supported by quantification of SST error in terms of bias, RMSE and correlation coefficient ( Supplementary Figure S1 ). The results also show a relatively warm bias of +1°C off the Cape Peninsula region in the CROCO model. This warm bias in sea surface temperature (SST) is more strongly expressed along the path of the northwards-flowing, shelf-edge Benguela jet current along the west coast of South Africa. The warm bias may be attributed to the location of the jet current which appears to be slightly far west in the CROCO model than in ESA-CCI product. While the model has a demonstrable ability to resolve fine-scale processes, it seems to position the mean location of the shelf-edge jet slightly far west. The jet appears to be too strong on the leading edge of the Benguela Current off the Peninsula. It follows that the observed warm bias in SST associated with eastern boundary currents is a fairly common phenomenon in regional coupled climate models (Dufois et al., 2012; Mao et al., 2019).

This warm bias typically manifests as higher than actual SSTs in subtropical eastern boundary regions of the ocean, particularly affecting upwelling systems resulting in discrepancies in modeled SST which can influence climate predictions and oceanic processes. Dufois et al. (2012) used the Pathfinder monthly SST climatologies and found out that this warm bias can be significant, with summertime differences up to 3–5°C in regions like the California, Humboldt, Canary, and Benguela upwelling systems. Ma et al. (2019) found that insufficient ocean dynamic upwelling is the dominant cause for SST warm bias in subtropical eastern boundary regions. Weakened eastern boundary currents result in additional warming that surpasses the cooling effects typically provided by upwelling in these regions. The origins of this warm bias, particularly during upwelling seasons, remain a subject of ongoing research.

Figure 4 shows a high number of particles observed in the Natal Bight and along the path of the Agulhas Current in both the model and the GLORYS. The Natal Bight is recognized as a significant area for larval retention, which is vital for various ecological processes, including the successful settlement of fish larvae in the northern Agulhas Current region (Heye et al., 2022). For example, Heye et al. (2022) found that the Natal Bight Coastal Counter increases connectivity between MPAs within the Natal Bight, where the current increases water retention entrapping larvae, thus contributing to increasing the suitability of the habitat for larval development and settlement. However, south of the 31⁰S, the model shows that particles released on the east coast of South Africa tend to move inshore of the Agulhas Current increasing residence times ( Figure 7 ). Conversely, the GLORYS show that particles are either entrained by the current or drift to the outer regions of the current. The study of Tedesco et al. (2019) demonstrated the ability of the 3 km grid CROCO model to represent frontal features of the northern edge of the Agulhas Current such as shear driven meso-and sub-mesoscale cyclonic eddies. Studies by Krug et al. (2017) and Schubert et al. (2021) showed that these shear instabilities extract their energy from the mean Agulhas Current jet and increase the residence times of particles inshore of the current area.

Another major difference between larval dispersals obtained from the model and GLORYS is that there is a significantly large offshore loss of particles in the GLORYS between the eastern and southern coasts of South Africa ( Figure 4 ). Demonstrable differences in the MKE and EKE are observed south of the African Continent. The differences in MKE can be attributed to the fact that the current is stronger and more inshore in the CROCO model, as opposed to the GLORYS where the current is slightly far offshore. The CROCO model’s realism has been validated by Tedesco et al. (2019) against in-situ Agulhas System Climate Array (ASCA) observations and altimetry from AVISO. The difference in EKE can be attributed to the net source of mesoscale eddy energy in the Agulhas Current region (Lutjeharms, 1989; Backeberg et al., 2008; Pauline et al., 2022; Tedesco et al., 2022). The net source of mesoscale eddy energy is caused by a strong generation of eddies (Tedesco et al., 2022) resulting in increased eddy energy in the Southern Agulhas Current. In the model, the local generation of mesoscale eddies dominates the net EKE budget leading to the observed EKE difference.

Circulation off the southern coast is relatively sluggish and is driven by winds and the Agulhas Current (Schumann et al., 1982; Schumann et al., 1995 and Lutjeharms and De Ruijter, 1996; Schumann, 1999; Lutjeharms, 2006). The region experiences a more variable wind regime especially during October, shifting from southwesterly in early October to southeasterly in late October (Lutjeharms, 2006; Rouault et al., 2010; Schumann, 2010). The GLORYS shows relatively sporadic dispersal patterns to the south of the continent as opposed to the model. However, both the GLORYS and model show a coherent high-density of particles over the Agulhas Bank ( Figure 4 ). The winds influence the cross-shore transport leading to inshore movement of larvae in the region (Boyd and Shillington, 1994; Malan et al., 2019). The Agulhas Current influences the inshore movement of particles through intermittent mesoscale Natal Pulses meanders (Krug and Penven, 2011; Krug et al., 2014; Malan et al., 2018) that bring surface water closer to the shore. Krug et al. (2014) show that the Natal Pulses impact the eastern Agulhas Bank for 110 days per year on average, and through their intrusion, they impact drive the particles to the shelf. Furthermore, Krug et al. (2014) show that these Natal Pulses are responsible for the largest SST and current velocity anomalies leading to increased particles on the west Agulhas Bank.

The circulation along the highly dynamic west coast is largely driven by winds and the Benguela Jet (Boyd and Nelson, 1998; Shannon, 2001). The position and intensity of the shelf-edge Benguela Jet are largely driven by seasonal wind forcings (Strub et al., 1998; Penven et al., 2000; Fennel et al., 2012; Veitch et al., 2018). These dynamics influence the high variability of dispersal patterns observed on the west coast ( Figures 4 , 5 ). The wind-modulated Ekman transport drives the onshore movement of particles in the region ( Figure 5 ; Desbiolles et al., 2015; Ragoasha et al., 2019; Ragoasha et al., 2022; Bordbar et al., 2021). We have anticipated the CROCO model to be able to resolve small-scale processes with more frontal variability at the inshore edge of Benguela Jet coupled with larval retention but our results show that the Jet is stronger there resulting in rapid advection of particles (and less retention) to the north and particles follow the mean jet path more. The model places the jet slightly more west than the GLORYS with implications for dispersals in that region. Conversely, the GLORYS show strong inshore movement leading to overestimation of local larval retention especially near the recognized recruitment area off the west coast.

Both the model and GLORYS show coherent upstream connectivity coupled with relatively high local recruitment and self-retention rates on the southern and western coasts ( Figure 6 ) - which are regions characterized by relatively sluggish circulation driven by winds (Boyd and Nelson, 1998; Schumann et al., 1982; Lutjeharms and De Ruijter, 1996; Shannon, 2001; Lutjeharms, 2006) and mesoscale meanders on the Agulhas Bank (Krug et al., 2014; Malan et al., 2018; Tedesco et al., 2022; Tedesco et al., 2024). The general upstream connectivity in these regions is associated with large residence times ( Figure 7 ). The difference between the model and GLORYS is found in the highly dynamic west coast region where GLORYS show larger connections, especially along the Cape Peninsula to the Cape Columbine (i.e., northwardly from Robben Island to the Childs Bank) stretch ( Figure 5 ).

Stronger connections and large residence times are observed in the MPAs located in the northern Agulhas Current ( Figures 6 , 7 ) in both the model and GLORYS. These results highlight the influence of the frontal activity of the Agulhas Current especially near the Natal Bight region. The results show that there is significant larval retention on the western Agulhas Bank. These findings are consistent with literature which demonstrates that the western Agulhas Bank is a significant retention area for many coastal fish species.