Satellite-derived products

A satellite-derived Level-3 data set of Chl a concentration (mg m−3) was obtained from the European Space Agency’s GlobColour project (http://www.globcolour.info). The 8-day composite Chl a concentrations using standard Case 1 water algorithms were used (i.e., OC4Me for Medium-Resolution Imaging Spectrometer, and OC3v5 for Moderate Resolution Imaging Spectroradiometer/Visible Infrared Imaging Radiometer Suite sensors). The altimetry-derived geostrophic velocities (AVISO MADT daily product) were produced by CLS/AVISO (Collecte Localization Satellites/Archiving, Validation, and Interpretation of Satellite Oceanographic data), with the support from the CNES (Centre National d’Etudes Spatiales; http://www.aviso.altimetry.fr/duacs/).

BGC-Argo network

The BGC-Argo dataset used in the present study is publicly available at ftp://ftp.ifremer.fr/ifremer/argo/dac/coriolis/, an Argo Global Data Assembly Center. This dataset represents an international initiative by compiling 132 BGC-Argo floats (a total of 14,415 stations) from AOML-NOAA, BODC, CORIOLIS, CSIRO, and INCOIS. It covers the time period from July 2010 to June 2017 and involves a variety of float platforms (NKE’s PROVOR, Webb’s APEX, and SeaBird’s NAVIS) performing different missions (in dive depths, frequency of profiling, and data acquisition) and equipped with different sensors for measuring Chl a and backscattering at 700 nm (Seabird MCOMS, WETLabs Eco-Triplet, and Eco-FLBB).

A quality control procedure was achieved on the CTD22, Chl a23 and backscattering24 data. Fluorescence data were corrected for non-photochemical quenching on daytime profiles following the method of Xing et al.25 as follows: the maximum Chl a value above the mixed layer depth (MLD), defined as a density difference of 0.03 kg m−3 from a reference value at 10 m, is extrapolated toward the surface. Fluorescence data are then converted to Chl a concentration (mg m−3) by first applying the calibration (dark value and slope) and then multiplying by a SO-specific correction factor.

To decide which calibration factor to apply, we carried out a robust analysis by applying the radiometric method of Xing et al.6 to retrieve F 490 (a refined calibration factor with respect to factory calibration) based on all available BGC-floats equipped with a downward irradiance sensor at 490 nm (OC4 radiometer, Satlantic). We allocated the various BGC-Argo profiles according to SO provinces to detect any potential intra-regional variability in F 490 . The average F 490 within provinces ranges between 0.26 and 0.33 (Supplementary Fig. 1, analysis of 3321 profiles), which translates into an overestimation factor of 3 to 4 with respect to the factory calibration. This value is actually lower than overestimates derived from HPLC Chl a26,27. We note here that these HPLC-based estimates (1) are relevant to a spatio-temporal domain restricted to the float deployment (the estimated correction factor might change as environment and community composition change during the float journey) and (2) present a large (yet unexplained) variability. Here, we use a conservative value of 0.3 for F 490 (corresponding to a factory calibration overestimation of 3.3). This value has the advantage of integrating a broad spatio-temporal domain (e.g., winter conditions) and a large dataset for the estimation of this correction (more than 3300 profiles).

ANDRO dataset

The ANDRO atlas ASCII file (available at http://www.coriolis.eu.org/) contains the float parking pressure (actually a representative parking pressure which is generally an average of the measured pressures during float drift at depth) and temperature, deep and surface displacements, and associated times, deep and surface associated velocities with their estimated errors (see Ollitrault and Rannou28). ANDRO data originate from AOML, Coriolis, JMA, CSIRO, BODC, MEDS, INCOIS, KORDI, KMA, and CSIO and represent a total of 6271 floats contributing to 612,462 displacements.

Each float cycle (deep displacement between two profiles) provides an estimate of the zonal and meridional current velocities at their drifting depth (mostly around 1000 m). From 2002 to 2016, around 21,300 cycles were available close to the SWIR (Supplementary Fig. 2). These velocities were binned into 1° by 1° boxes and then averaged in space and time. Only the box containing more than 5 data points were kept. Standard deviation of the zonal (u’) and meridional (v’) velocity components in each box were used to calculate the mean deep eddy kinetic energy (EKE) as follows:

$$\overline {{\mathrm{EKE}}} = \frac{1}{2}\left( {u\prime ^2 + v\prime ^2} \right)$$ (1)

where the overbar denotes the time average over the whole period (2002–2016). EKE was then interpolated on a finer grid by a Gaussian correlation function, weighted by the local number of data, with a decorrelation radius of 100 km.

Helium dataset

The helium data were extracted from the GLobal Ocean Data Analysis Project (GLODAP) Version 2, a cooperative effort to coordinate global synthesis projects funded through NOAA/DOE and NSF as part of the Joint Global Ocean Flux Study–Synthesis and Modeling Project (JGOFS-SMP)29. The GLODAP Version 2 data product (available at https://www.glodap.info) is composed of data from 724 scientific cruises covering the global ocean. Here, two different expeditions (06AQ19960317; March 23–31, 1996 on the Polarstern) and (35MF19960220; February 28–March 25, 1996 on the Marion Dufresne) were used to generate the two sections of δ3He.

Isopycnals between meridional section

A combination of hydrographic profiles from different sources are used to construct a 3-D climatology of potential density in the region 0–55°E and 55–45°S, with a half degree horizontal resolution, and a 25 m vertical resolution.

We use three distinct sources of observations to maximize the number of profiles. The first set of observations is conductivity-temperature-depth (CTD) data from ship-recorded observations during the period 1906–2016 from the NOAA World Ocean Database (https://www.nodc.noaa.gov/OC5/SELECT/dbsearch/dbsearch.html). We only use profiles that have a quality control flag of 1, containing information on their position, date, temperature, and salinity. The second set of observations we use is float observations from the Argo international program. The Argo float profiles of pressure, salinity, and temperature used in this study were gathered in the period 2002–2016. They provide temperature and salinity between 0 and 2000 m. We only use profiles that have a quality control flag of 1, and contain information on their position, date, temperature, and salinity. As a final data set, we use profiles derived from the animal-borne sensor program MEOP (http://www.meop.net/; Treasure et al.30). Similar to the other datasets, we only use profiles with control flag of 1, and that contains position, date, temperature, and salinity. Altogether, we gathered 33096 profiles in the region 0–55°E–55–45°S.

From this dataset, we computed potential density for each 25 m vertical interval between the sea surface and 2000 m. Then, for each interval, we produced maps of climatological fields of potential density using an Optimal Interpolation procedure. The Optimal Interpolation and gridding method are described in detail in Schmidtko et al.31. As a brief summary, we interpolate onto a 0.5° grid in the region 0–55°E–55–45°S. We used a 550 km isotropic decorrelation scale, incorporating an anisotropic isobath-following component using a “Fast Marching” algorithm, as well as front-sharpening components. In addition, recent data are emphasized in the mapping, which produces a climatology typical of the years 2000–2010 (see Schmidtko et al.31 for more details on the mapping).

Two vertical sections of density are extracted at 28°E and 38°E from the produced climatology: σ 28° (P,lat) and σ 38° (P,lat). Because the ACC is not entirely zonal, comparing the density structure of these two sections as a zonal difference would compare density structure from south and north of given ACC fronts. Instead, we determined how the density structure differed between 28°E to 38°E but alongstream, i.e., following the ACC structures. For that, we use the dynamical height (dh) provided by AVISO for the period 1993–2012 (http://www.aviso.altimetry.fr/) at these two sections to convert σ 28° (P,lat) and σ 38° (P,lat) into σ 28° (P,dh) and σ 38° (P,dh). Because fronts and jets tend to follow individual contours of dynamical height (e.g., Sokolov and Rintoul32; Sallée et al.33), comparing these two sections in dynamic height coordinate ensures an alongstream comparison, i.e., dynamically consistent with regards to the ACC. We therefore produce a difference section: ∆σ(P,dh), and using the mean dynamic height in the sector 28–38°E, we produce a mean relationship between dynamic height and latitude: \(\widetilde {{\mathrm{lat}}} = {\mathrm{f}}({\mathrm{dh}})\), where \(\widetilde {{\mathrm{lat}}}\) is referred to as pseudo-latitude, which we use to produce \(\Delta {\mathrm{\sigma }}({\mathrm{P}},\widetilde {{\mathrm{lat}}})\).

Lagrangian modeling of horizontal iron delivery

An advection scheme based on altimetry was used here to estimate iron delivery due to horizontal stirring from (1) shallow bathymetry (<500 m; as shown in Fig. 1 and in Ardyna et al.8) and from (2) the initial location of the iron pathways in the surface layer along the SWIR (Fig. 4e, f). According to the analysis on vertical divergence (Fig. 4d), the origin of the iron pathways in the surface layer is located in the area of vertical displacement of isopycnals, suggesting the upwelling on the SWIR corresponds to the Andrew Bain fracture zone (see Fig. 4c). This feature is located between (24.5°E–56°S) and (34.2°E–44°S) (http://www.ngdc.noaa.gov/gazetteer/). Thus, the enrichment of the two BGC-Argo floats was identified in the region where the floats crossed this geological structure with a high isopycnal adjustment (Fig. 4d). This area has been represented by two overlapping disks centered in (30°E–50°S and 30.5°E–49.5°S), and with a 1° radius, as shown in Fig. 4e, f. We note that another possible upwelling region may occur north or south to this area, according to the analysis of vertical displacement of isopycnals (Fig. 4d).

The advection scheme seeds each location of the study region with a spatial resolution of 1/4° (Fig. 1) and 1/8° (Fig. 4e–f). Trajectories are derived from surface velocities by applying a Runge–Kutta fourth-order scheme with a time step of 6 h, in which velocity fields have been linearly interpolated in both space and time. The advection scheme then finds the particle’s most recent contact with an iron source (shallow bathymetry and upwelled input along the SWIR) and provides the time at which the contact took place. The iron content of each particle that was in contact with a potential source of iron was estimated with an exponential scavenging relation. This relation reproduces the decreasing concentration of bioavailable iron along the trajectory after the contact with the iron source. This approach was initially developed for predicting the development of the Kerguelen phytoplanktonic plume19 and thereafter extended to the entire Southern Ocean8. The model has been calibrated and validated in the Crozet and Kerguelen regions by combining satellite data (altimetry and ocean color), lithogenic isotopes, iron measurements, and drifters19,20,21 (see d’Ovidio et al.19 for further details). Here the advection scheme is applied to the period 2010–2015 and for the planktonic bloom season, November to March, in order to obtain a mean climatological signal.

Bloom characterization

To determine the bloom magnitude, each float time series was divided into individual annual cycles, starting on 1 July. Cycles that do not cover the theoretical bloom period (from early November to late February the next year) with at least eight float profiles were discarded. For each remaining cycle, float profiles were binned into a 20-day period corresponding to the decorrelation scale of float Chl a records34. The magnitude of the bloom was then computed as the maximum in integrated Chl a biomass from the surface down to either the MLD or the euphotic depth, whichever was deeper. The euphotic depth is the depth at which light is 1% of its surface value, based on the surface Chl a according to Morel et al.35.