Efficient biological carbon export to the mesopelagic ocean induced by submesoscale fronts

Oceanic submesoscale processes are ubiquitous in the North Pacific Subtropical Gyre (NPSG), where the biological carbon pump is generally ineffective. Due to difficulties in collecting continuous observations, however, it remains uncertain whether episodic submesoscale processes can drive significant changes in particulate organic carbon (POC) export into the mesopelagic ocean. Here we present observations from high-frequency Biogeochemical-Argo floats in the NPSG, which captured the enhanced POC export fluxes during the intensifying stages of a submesoscale front and a cyclonic eddy compared to their other life stages. A higher percentage of POC export flux was found to be transferred to the base of mesopelagic layer at the front compared to that at the intensifying eddy and the mean of previous studies (37% vs. ~10%), suggesting that the POC export efficiency was significantly strengthened by submesoscale dynamics. Such findings highlight the importance of submesoscale fronts for carbon export and sequestration in subtropical gyres.


Text S2. Eddy tracking results
Based on SLAs, eddies were tracked, and their trajectory products were downloaded from AVISO+ with the algorithm of Pegliasco et al. 1 for the delayed-time product.Eddy information includes the location of the eddy center, the outline of the eddy edge, the radius R and the amplitude.Using this eddy dataset, the parameters of the eddies sampled by the BGC-Argo floats were computed (Table S2).According to the satellite-derived eddy identification and tracking dataset, the characteristics of eddies CE 1 and CE 2, including the eddy radius, amplitude and propagation speed, were averaged over the period during which they were sampled by the floats, as well as over the whole lifespan of the eddy.

Text S3. Lagrangian particle tracking experiment at the front
To investigate whether the float drifted along the front, a numerical particle tracking experiment was conducted using the TRACMASS Lagrangian trajectory code 2,3 .The particle was released at the same location as the float on 21 April 2019 and transported purely by surface current.Here, we run the code with the current velocity derived from AVISO.The Lagrangian particle trajectory derived from the particle tracking model demonstrated that the virtual particle starting at the location of the float moved along the front (Fig. S15).The float moved slower than the particle but mainly followed its path, suggesting that the float drifted in a quasi-Lagrangian manner during the analysis period.The difference in moving velocity is likely that the particle was only transported by the surface current, whereas the float was moved by both the surface and deep currents.

Text S4. POC transfer efficiency affected by fronts in subtropical gyres
To illustrate the potential implications of our findings, we combined previous studies that have estimated POC attenuation coefficients in the North Pacific and North Atlantic gyres (Table S4, Fig. S12).In much of the world's ocean, temperature variations play a leading role in controlling the spatial near-surface density gradients, while salinity plays a secondary role 4 .Sea-surface temperature gradients derived from satellites at a 4-km resolution have been used to study fronts and filaments at scales of O(1 to 10) km 5,6 .A strong correlation between the b value and GM was found (bGM=-3.74•GM+1.43;R 2 =0.52, p<0.01,Fig. 5).This statistical relationship implies that a strong front is linked with high POC transfer efficiency to the mesopelagic layer, which could be used to estimate the impact of fronts on carbon export in subtropical gyres.
At the global scale, water temperature has been proposed as the key factor controlling the variation in b from low to high latitudes through its influence on microbial remineralization 7,8 .According to Marsay et al. 7 , mesopelagic POC flux attenuation is correlated with water column temperature (bT=0.062•T+0.303,where T is the median temperature of the upper 500 m).For the NPSG, there is substantial submesoscale variability in bGM compared to bT or constant b values in the gyre, and the bGM generally deceased from low to high latitudes in a way that was similar to bT (Fig. S16a and c).
The mesopelagic POC transfer efficiency (Teff), which is the POC flux at 1000 m relative to the flux at the base of the euphotic zone, was estimated using bGM and bT (Text S6).Compared to the Teff estimated from bT (Teff(bT)), the spatial variability in Teff estimated from bGM (Teff(bGM)) was large (Fig. S16b and d), which reflected the influence of frontal dynamics on enhancing local carbon transfer efficiency 9 .Overall, bGM was smaller than bT in the NPSG (on average 1.16±0.04versus 1.64±0.02)during 2003-2019 (Fig. S16e).The mean bGM value was more consistent with the mean observational b of 1.19 (Table S4).Both bGM and bT were lower in winter than in summer, which is in agreement with the seasonal b value in the NPSG estimated using POC flux profiles recorded at ALOHA (b=0.81±0.19 in winter versus b=1.13±0.39 in summer; Fig. S17).The lower b value in winter than in summer is likely attributable to stronger submesoscale activities or lower remineralization rates in winter than in summer 7,10 .As a result, the Teff in winter was higher than that in summer, as estimated from both bGM and bT.Due to the lower value of bGM, the Teff(bGM) at 1000 m was much higher than the Teff(bT) (on average 12.4±1.2%versus 4.7±0.1%)during 2003-2019 (Fig. S16f).The elevated Teff(bGM) is consistent with the BGC-Argo observations in this study and the observations of Guidi et al. 11 .For the entire NPSG, estimations of POC export flux at 1000 m by using bT or a constant b (1.19) may underestimate the flux estimated with bGM by ~58% and ~20% on average, respectively.
When the float was deployed, seawater samples were collected at 10 locations, 6 of which (marked as black dots in Fig. 1a; Table S5) were used for nitrate determination, and the others (marked as yellow pentagrams in Fig. 1a) were used for POC determination.Nitrate from discrete samples were analyzed onboard with a Fourchannel Continuous Flow Technicon AA3 Auto-Analyzer (Bran-Lube, GmbH) at depths of 0-200 m with an interval of 10 m and 200-500 m with an interval of 100 m.
The POC concentrations were determined by analyzing 8 L of seawater collected at intervals of 30 m between depths ranging from 30 to 180 m, as well as at a depth of 200 m using an Elemental vario EL cube CHNS Analyzer 12 .The uncertainty for POC data is < 10%.
The shipboard nitrate concentration was used to calibrate the nitrate estimations derived from BGC-Argo float data.By matching each shipboard sampling station with the nearest float sampling location from the BGC-Argo trajectories, the vertical distributions of nitrate were compared between the shipboard and BGC-Argo measurements (Fig. S13a and b).The nitrate concentrations in the upper layer were below the limit of quantification (BLQ), which is 0.1 μmol L -1 .From 100 m to 300 m depth, the ship-based nitrate concentration increased from 0.1 μmol L -1 to ~9 μmol L -1 .
The measurements from BGC-Argo floats captured the vertical structure of nitrate distributions in the NPSG well, despite a slight mismatch of the values below 300 m for both floats (Fig. S13a and b).The X-Y scatter plots for both BGC-Argo floats demonstrated a high correlation between the shipboard and BGC-Argo derived nitrate concentrations (R 2 =0.98 for float 2902753 and R 2 =0.96 for float 2902756; Fig. S13c     and d).The linear regression relationships for both BGC-Argo floats (Y=0.92•X+0.01 for float 2902753 and Y=1.13•X-1.1 for float 2902756) suggested an overestimation for float 2902753 and an underestimation for float 2902756, which were corrected for all the float profiles.
The POC concentration was estimated by converting the BGC-Argo derived bbp(700) into the POC concentration based on empirical equations [13][14][15][16] .In subtropical gyres, the relationship between POC and bbp(700) is usually described as follows 16 : This relationship was shown to overestimate the POC in our study region compared to in situ measurements (Fig. S14).We thus derived another relationship between POC and bbp(700) by applying a linear fitting method to the in situ measured POC and the BGC-Argo derived bbp(700): Eq. (S2) was shown to better estimate the POC distribution than Eq.(S1) in the NPSG, as the mean square error (MSE) was much smaller than that from Eq. (S2).

Text S6. Estimating the transfer efficiency of POC flux in the NPSG
The relationships between the b value and the SST gradient (bGM) and the temperature (bT) were applied to the NPSG to estimate the POC transfer efficiency.The relationship between bGM and GM was shown in Fig. 5: The bGM accounts for the influence of front on POC transfer efficiency.On the other hand, Marsay et al. 7 suggested that the attenuation coefficient can be estimated based on the relationship between the b value and upper water column temperature: where T is the median temperature of the upper 500 m.
The daily SST values were obtained from the CMEMS at a spatial resolution of ~4.5 km to calculate the spatial distributions of GM following the method provided by Belkin and O'Reilly 17 .The bGM from 2003 to 2019 was then obtained using Eq.(S3).
The bT in the NPSG was obtained using the upper water column (500 m) median temperature provided by the World Ocean Atlas climatological data.
The transfer efficiency of POC flux to the mesopelagic layer was estimated from bGM ( eff ( GM )) and from bT ( eff ( T )) as

Figure S1 .
Figure S1.Scatterplots of ZEu, MLD and ZSCM versus r/R in a CE 1 and versus SLA in b CE 2. Source data are provided as a Source Data file.

Figure S2 .
Figure S2.Comparison of sea surface Chl between satellite and BGC-Argo observations in both spatial and temporal variability.a Radial distribution of sea surface Chl at the float locations derived from both the BGC-Argo and satellite observations.The BGC-Argo results were directly obtained from the surface Chl measured by the float.The satellite-derived Chl values were calculated by averaging the Chl from OC-CCI over the float sampling period and then interpolating to the relative locations of the float in CE 1.Because the eddy propagated over time, the location of the float in the eddy was constructed by calculating the azimuthal and radial positions of the float relative to the eddy center.b Satellite-derived Chl averaged inside CE 1 and sea surface Chl derived from BGC-Argo observations.The BGC-Argo results were directly obtained from the surface Chl measured by the float.The satellite-derived Chl values were calculated by averaging the satellite-derived Chl values inside the eddy at the time when the float was at the position of r/R.The color in b indicates date.

Figure S3 .
Figure S3.Radial distributions of a temperature, b NO3 -, c Chl, and d POC in CE 1 at the depth of ZSCM.The solid lines indicate the linear regression lines.Source data are provided as a Source Data file.

Figure S4 .
Figure S4.Time series of sea level anomaly (SLA) along (orange) the float trajectory and (black) averaged inside CE 2. Source data are provided as a Source Data file.

Figure S5 .
Figure S5.Scatterplots of a temperature, b NO3 -, c Chl, and d POC at the depth of ZSCM versus SLA during (blue) the eddy intensification stage and during (orange) the other observational period in CE 2. The solid lines represent the linear regression lines.Source data are provided as a Source Data file.

Figure S6 .
Figure S6.The SST gradient magnitude (GM), lateral buoyancy frequency ( 2 ) and normalized strain rate (S/f) at the front.a Weekly averaged sea-surface temperature (background color) and the locations of the BGC-Argo floats in the eddy during 21-27 April 2019 (purple points).The purple line is the float trajectory from 21 April to 3 May 2019.b The same as a but during 28 April-3 May 2019.c The GM distribution on 25 April 2019.The purple line is the float trajectory during 21-27 April 2019.d The same as c but on 29 April 2019.The purple line is the float trajectory during 28 April-3 May 2019.e The GM at the BGC-Argo float location during 21 April-3 May 2019.

Figure S7 .
Figure S7.Physical characteristics of the float measurements.During 21 April-3 May 2019, the float was on one side of the front.Time series of a temperature, b salinity, c temperature anomaly and d salinity anomaly measured by float 2902753 from 1 April to 3 May 2019.The black, red and cyan lines indicate the mixed layer depth (MLD), the depth of subsurface chlorophyll maximum (ZSCM) and the euphotic zone (ZEu), respectively.The white lines are the isopycnals of 23, 24 and 25 kg m -3 , respectively.The temperature and salinity anomalies were defined as the instantaneous values subtracting the full time mean.e Temperature-salinity diagram with data points from float 2902753 from 1 April to 3 May 2019.The color indicates the sampling date.

Figure S8 .
Figure S8.The lateral buoyancy frequency ( 2 ), Richardson number (Ri) and finitesize Lyapunov exponent (FSLE) at the front.a BGC-Argo-derived lateral buoyancy frequency during 21 April-3 May 2019. 2 = −   ⁄ , where s is the along-track direction.The black line represents the mixed layer depth.b BGC-Argo-derived Richardson number during 21 April-3 May 2019.  =  2    ⁄  2 , where f and b are the Coriolis parameter and buoyancy, respectively.The black line represents the mixed layer depth.c Surface FSLE in the frontal region (orange line) and the lateral buoyancy frequency (blue dashed line) at 100 m during 21 April-3 May 2019.Blue solid lines are the buoyancy frequency envelope.Purple dash lines indicate the shift from the intensifying to decaying stage.

Figure S9 .
Figure S9.The POC and spike of bbp(700) at the front.a BGC-Argo derived POC distribution during 21 April-3 May 2019.The black, red and cyan lines are the MLD, ZSCM and ZEu, respectively.b Depth-integrated POC in the euphotic zone during 21 April-3 May 2019.Gray and black lines in b represent the time series of the integrated POC without and with a 5-point moving filter, respectively.c Spikes of bbp(700) during 21 April-3 May 2019.The cyan line indicates ZEu.

Figure S10 .
Figure S10.Vertical distribution of spiciness and POC.Time series of a spiciness and c POC distributions measured by float 2902753 from 21 April to 3 May 2019.Time series of b spiciness anomaly and d POC anomaly measured by float 2902753 from 21 April to 3 May 2019.The cyan line indicates the euphotic zone depth.The black lines are the isopycnals of 24 and 25 kg m -3 .The spiciness and POC anomalies were defined as the instantaneous values subtracting the annual mean at the same isopycnal layer.

Figure S11 .
Figure S11.POC flux profiles.a The POC flux profile with the subduction pump (ESP) between ZEu and 250 m accounted for when estimating the b value.b The POC flux profile without the ESP when estimating the b value.The black dashed lines and error bars correspond to the euphotic zone depth and one standard deviation, respectively.The red areas represent the euphotic zone.Source data are provided as a Source Data file.

Figure S12 .
Figure S12.Locations of the POC flux data from previous studies (Fig. 5, TableS4) in subtropical gyres.The background color is the mean chlorophyll concentration from satellite data, and the black contour is the chlorophyll concentration of 0.1 mg m -3 .

Figure S13 .
Figure S13.Calibration of the BGC-Argo derived nitrate concentrations.a-b Vertical profiles and c-d scatterplot showing the comparison between the in situ and BGC-Argo derived nitrate concentrations.a and c are results for float 2902753, and b and d are results for float 2902756.Solid lines represent average results, and shades represent their standard deviations.The number of in situ sampling stations N is indicated in a and b.The orange lines in c and d are the linear regression lines.

Figure S14 .
Figure S14.Relationship between BGC-Argo derived bbp(700) and in situ measured POC concentrations.Scatterplot of BGC-Argo derived bbp versus in situ POC.Eq. (S1) (gray dashed line) was from Stramski et al. 16, and Eq.(S2) (black solid line) was derived from the linear regression of BGC-Argo derived bbp versus in situ measured POC.Orange and blue colors indicate data points from floats 2902753 and 2902756, respectively.The mean square error (MSE) for Eq.(S1) and Eq.(S2) are 3.7 and 0.1 mmol 2 C m -6 , respectively.

Figure S15 .
Figure S15.Trajectories of the (triangle) BGC-Argo and (dots) the virtual particle over the period when the float was at the front (21 April-3 May 2019).The particle was released at the same position as the float on 21 April 2019 and transported by the surface current velocity derived from AVISO.The particle trajectory was tracked using the TRACMASS code.Color indicates the date.Blue and orange lines are the contour lines of the sea level anomaly (SLA).Arrows represent the sea surface current velocity anomaly.Source data are provided as a Source Data file.

Figure S16 .
Figure S16.Transfer efficiency of POC flux at 1000 m estimated based on attenuation coefficient b. a-b Spatial distributions of attenuation coefficient b derived from GM (bGM) on 25 April 2019 (a) and its transfer efficiency of POC flux at 1000 m (Teff(bGM)) (b).c-d Spatial distributions of the attenuation coefficient b derived from the median temperature of the upper 500 m (bT) on 25 April 2019 (c) and its transfer efficiency of POC flux at 1000 m (Teff(bT)) (d).e-f Bar plots quantifying the b value (e) and the transfer efficiency of the POC flux Teff at 1000 m (f) during 2003-2019 in the NPSG.The black contour in a-d is the Chl of 0.1 mg m -3 , which represents the outline of the NPSG.Bars are plotted as the mean ± 1 standard deviation.Source data are provided as a Source Data file.

Figure S17 .
Figure S17.The seasonal b value at the ALOHA station.a POC flux profiles and the corresponding Martin curves in summer.Black dots are observational profiles from sediment traps.Gray lines are the fitted Martin curves.The blue line is the Martin curve derived from the averaged b value and averaged POC flux at 150 m. b is the same as a but in winter.c The b value and R 2 for each observational POC flux profile in both summer and winter.We picked up all observational POC flux profiles with at least three vertical sampling points for Martin curve fitting.The b values in summer are significantly different from those in winter, as h=1 and p=0.036 in the two-sample t test result.The b values in a-b are shown as the mean ± 1 standard deviation.Source data are provided as a Source Data file.

Table S1 .
The operation details of sampling when BGC-Argo encountered

Table S3 .
Correlation coefficients (CR) of temperature, NO3 -, DO, Chl and POC with r/R at both the sea surface and ZSCM for CE 1 and with SLA during the intensifying stage (Intensification) and the whole observational period for CE 2 (ALL).

Table S4 .
Compilation of POC export fluxes in the subtropical North Pacific and North Atlantic gyres a F (1000 m) is derived from the fitted Martin curve.

Table S7 .
Linear least squared estimations for the rate of change of AOU (unit: mmol O2 m -3 d -1 )