Episodic subduction patches in the western North Pacific identified from BGC-Argo float data

Subduction associated with mesoscale eddies is an important but difficult-to-observe process that can efficiently export carbon and oxygen to the mesopelagic zone (100–1000 dbar). Using a novel BGC-Argo dataset covering the western North Pacific (20–50^∘ N, 120–180^∘ E), we identified imprints of episodic subduction using anomalies in dissolved oxygen and spicity, a water mass marker. These subduction patches were present in 4.0 % (288) of the total profiles (7120) between 2008 and 2019, situated mainly in the Kuroshio Extension region between March and August (70.6 %). Roughly 31 % and 42 % of the subduction patches were identified below the annual permanent pycnocline depth (300 m vs. 450 m) in the subpolar and subtropical regions, respectively. Around half (52 %) of these episodic events injected oxygen-enriched waters below the maximum annual permanent thermocline depth (450 dbar), with >20 % occurring deeper than 600 dbar. Subduction patches were detected during winter and spring when mixed layers are deep. The oxygen inventory within these subductions is estimated to be on the order of 64 to 152 g O₂/m². These mesoscale events would markedly increase oxygen ventilation as well as carbon removal in the region, both processes helping to support the nutritional and metabolic demands of mesopelagic organisms. Climate-driven patterns of increasing eddy kinetic energies in this region imply that the magnitude of these processes will grow in the future, meaning that these unexpectedly effective small-scale subduction processes need to be better constrained in global climate and biogeochemical models.

1 Introduction

Ocean subduction is the process of transporting water from the wind-mixed surface layer into or below the permanent thermocline, resulting in the efficient injection of heat, carbon, and oxygen to the ocean interior (Fig. 1). Subduction therefore plays an important role in regulating global climate and carbon cycles (Sabine et al., 2004; Qu and Chen, 2009; Stukel et al., 2017, 2018; Boyd et al., 2019; Martin et al., 2020). Many studies focus on the subduction of mode waters driven by large-scale circulation and the seasonal cycle of the mixed-layer dynamics (Williams, 2001; Qu et al., 2002; Qiu et al., 2007; Koch-Larrouy et al., 2010; Kawakami et al., 2015; Nie et al., 2016). But recent advances have highlighted the importance of small-scale (1–100 km) dynamical processes for vertical transport and biogeochemistry in the upper ocean, driven by mesoscale eddies and sub-mesoscale processes (Lévy et al., 2001; Xu et al., 2014; Omand et al., 2015; McGillicuddy, 2016; Llort et al., 2018; Resplandy et al., 2019). Ocean general circulation models typically resolve the large-scale subduction of mode waters (Koch-Larrouy et al., 2010) but cannot accurately capture small-scale, short-term subduction processes because of their episodic characteristics (Xu et al., 2014; Llort et al., 2018).

Figure 1 An illustration of the Kuroshio and Oyashio Extension region depicting the different modes of carbon export below the maximum annual mixed-layer depth; the biological gravitational pump (sinking export, zooplankton migration) and subduction in the region of the Kuroshio and its extension (yellow line) and Oyashio and its extension (grey line). The subducted surface waters, apparently driven by mesoscale eddy processes, travel along isopycnal surfaces transporting water containing high dissolved oxygen (DO), dissolved organic carbon (DOC), and slowly sinking particulate organic carbon (POC) into the mesopelagic zone (low DO, DOC, and POC). The green layer represents the euphotic zone, and the blue layer below is the mesopelagic zone.

Subduction associated with mesoscale and sub-mesoscale dynamics has been observed at higher latitudes in the North Atlantic (Omand et al., 2015) and Southern oceans (Llort et al., 2018), and similar processes have been shown to occur in the Kuroshio Extension (KE) region in the western subtropical Pacific. Shipboard sampling techniques have been used there to identify small water parcels within the main thermocline having low potential vorticity, elevated dissolved oxygen (DO), and anomalous salinity, which are signals indicative of small-scale subduction (Yasuda et al., 1996; Okuda et al., 2001; Oka et al., 2009). Analogous phenomena have been observed in mooring data from the region (Nagano et al., 2016; Inoue et al., 2016a; Kouketsu et al., 2016; Zhu et al., 2021), and more focused sampling of anticyclonic eddies with Argo floats (Zhang et al., 2015; Inoue et al., 2016b) and Seagliders (Hosoda et al., 2021) confirm the existence of discrete subsurface water mass exchanges. These episodic features will contribute to both ventilation of the mesopelagic zone and export of dissolved inorganic and organic carbon from surface waters (i.e., the solubility pump; Sarmiento and Gruber, 2006), but their frequency, spatial extent, and lifetimes remain unknown (Hosoda et al., 2021).

Eddy-associated processes that generate vertical transport of productive and detrital planktonic biomass into the mesopelagic zone affect not only carbon export but also carbon sequestration timescales (i.e., time that carbon remains within the ocean interior). In general, sequestration timescales are proportional to the depth of injection, but the more important factor is whether these injections extend below the annual maximum mixed-layer depth (MLD), or permanent pycnocline, which hinders its return to the atmosphere (Boyd et al., 2019). Although eddy subduction has the potential to contribute significantly to global carbon export, evidence of the subsurface fate of injected carbon has been indirect and patchy (Estapa et al., 2019), highlighting the challenge of detecting and quantifying carbon export associated with mesoscale and sub-mesoscale processes.

The uncertainty about the contribution of eddy subduction to carbon and oxygen transport into the mesopelagic and deeper ocean interior has ramifications for both biogeochemical and ecological processes (Fig. 1). The transport of freshly produced particulate and dissolved organic carbon, along with oxygen, from surface waters to the mesopelagic zone is critical for balancing upper-ocean carbon budgets (Emerson, 2014) and supporting the nutritional demands of mesopelagic organisms (Dall'Olmo et al., 2016). The knowledge gap in these episodic processes is particularly evident in the mid-latitude western North Pacific, where mesoscale eddies, recirculation gyres, fronts, and jets are amplified under the influence of the Kuroshio and Oyashio currents and their extensions (Nishikawa et al., 2010). Shoaling of the maximum annual MLD in this region relative to higher latitudes (Cronin et al., 2013; Palevsky and Doney, 2018) has the potential to increase carbon sequestration efficiency and oxygenation of the deep mesopelagic zone (Bushinsky and Emerson, 2018).

Here we investigate small-scale subduction events in the western North Pacific region over the past decade (2008–2019). These events were identified with a new algorithm utilizing anomalies of apparent oxygen utilization (AOU; a proxy for dissolved and particulate organic matter degradation) and potential spicity (π; a characteristic water mass marker) obtained from multiple biogeochemical Argo (BGC-Argo) datasets (Claustre et al., 2020; Chai et al., 2020). These findings show the spatial and temporal distributions of subduction patches reflecting episodic injection processes that contribute to the missing fraction of carbon and oxygen export into the deep twilight zone (Emerson, 2014; Martin et al., 2020) but also have the potential to become increasingly significant under future climate scenarios.

2 Data and methods

2.1 Data

After standard data quality control, 7120 profiles from 43 BGC-Argo floats in the western North Pacific (20–50^∘ N, 120–180^∘ E) between 2008 and 2019 were selected (Fig. 2). All of these profiles contained measurements of temperature, salinity, pressure, and dissolved oxygen (DO, µmol/kg). The upper 1000 dbar of the ocean was sampled in each profile, and the typical profiling interval was between 5–10 d, with the floats parking at 1000 dbar depth in between. The typical vertical sampling frequency was every 5, 10, and 50 dbar for depth intervals of 0–100, 100–500, and 500–1000 dbar, respectively. Some floats were set with daily profiling and higher vertical frequency (e.g., every 2 dbar) for specific purposes.

Figure 2 Horizontal distribution of the quality-controlled BGC-Argo data profiles between 2008 and 2019 in the western North Pacific. The Argo profiling tracks are color-coded by the Julian day (a) and data density (number of available profiles) for each grid ($\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$) (b). The location of Station no. 234 from float MD5904034 is denoted by the black dot in (a) (see Fig. 3); the white curve inside the dashed box represents the trajectory of float MR2901556 between 28 July and 18 August 2014, and the black star indicates the beginning of the float trajectory during this period (see Fig. 4). The white box in (b) denotes the region with strong energetic ocean processes (i.e., Kuroshio–Oyashio Extension, associated with eddy activities).

All BGC-Argo variables were vertically smoothed with a three-bin running average to remove sharp noises or spikes (Llort et al., 2018). Two key variables, apparent oxygen utilization (AOU) and potential spicity (π), were derived from the direct measurements. Specifically, AOU is defined as the difference between saturated oxygen concentration (O_sat) and DO, and O_sat is estimated from temperature and salinity (Garcia and Gordon, 1992). AOU is a proxy for water mass age which reflects the microbial respiration of dissolved and particulate organic matter (Sarmiento and Gruber, 2006). Potential spicity referenced to the surface pressure is calculated from pressure, temperature, and salinity following Huang et al. (2018). Seawater is a two-component system. The water mass anomaly is commonly analyzed in terms of the (potential) temperature and salinity anomaly, and isopycnal analysis is also widely used. By definition, the temperature and salinity anomaly on an isopycnal surface is density compensated; thus, the water mass anomaly on an isopycnal surface is commonly described in terms of another thermodynamic variable, which is called spice, spiciness, or spicity. Over the past decades, there have been different definitions of such a thermodynamic variable; however, a most desirable property of such a thermodynamic function is that it is orthogonal to the density. Recently, Huang et al. (2018) proposed a potential spicity function (π) by the least-squares method, which is practically orthogonal to the potential density, with the root mean square of the angle deviation from orthogonality at the value of 0.0001^∘. Therefore, combining density and spicity gives rise to an orthogonal coordinate system. It is the thermodynamic variable we used in this study, which allows differentiating water masses with distinct thermohaline properties but similar density. In addition, potential density (σ) referenced to the surface pressure was derived from pressure, temperature, and salinity based on the thermodynamic equation (TEOS-10; McDougall and Barker, 2011); and MLD was estimated based on a threshold (0.05 kg/m³) of the difference in density from a near-surface value (i.e., at 10 dbar) (Brainerd and Gregg, 1995). All these derived variables were calculated for each of the 7120 profiles.

In addition to the BGC-Argo float data, satellite data of daily sea level anomalies (SLAs) and daily geostrophic velocity anomalies (u^′ and v^′) between 1993 and 2018 were also processed. The geostrophic velocity anomalies were used to calculate the eddy kinetic energy (EKE) as $\mathrm{EKE}=\frac{\mathrm{1}}{\mathrm{2}}\sqrt{{u^{\prime }}^{\mathrm{2}}+{v^{\prime }}^{\mathrm{2}}}$. These data were used to identify the spatial relationship between surface mesoscale circulation and the float profiles.

2.2 Methods

2.2.1 Subduction detection

When a BGC-Argo float passes through a parcel of water injected from the mixed layer, it captures coherent anomalous features in AOU and π distinct from the surrounding waters (Fig. 1). These anomalies can be used to identify subduction patches that are indicators of subduction events occurring in the vicinity (Omand et al., 2015; Llort et al., 2018). Quantifying anomalies in AOU and π (denoted as Δ_AOU and Δ_π) requires defining the reference values of AOU and π at the mean state of the profile without subduction. Llort et al. (2018) used the 20-bin running averages of the profiles as the references; however, we found that this approach could dampen the subduction signal and thus miss subduction patches as well as misidentify other signals as subduction (see Fig. S1). To avoid misreporting these anomalies, a revised detection method was developed by trial and error, as shown in example profiles of AOU, π, DO, and σ for Station no. 234 of float MD5904034 (Fig. 3; see Fig. 2a for its sampling location). Two subduction patches are visually apparent at ∼ 230 and ∼ 300 dbar (yellow shading in Fig. 3a, b). The identification of the lower subduction patch at ∼ 300 dbar from the spicity profile is briefly described below and is illustrated in Fig. 3c:

1. Calculate the slopes (i.e., first-order derivative) for profiles of AOU and π against depth.

2. Locate the peaks in AOU and π profiles (e.g., the blue star in Fig. 3c) based on their slopes. Specifically, if at one sampling point the slope changes from positive to negative when moving downwards, it is called a negative peak and vice versa. Only the negative or positive peaks in π associated with a negative peak in AOU are considered, as only negative AOU anomalies indicate potential water transport from the surface mixed layer (Llort et al., 2018).

3. Locate the coherent peaks in both AOU and π, and mark their depths as the targeted locations (represented by pressure, p) for potential subduction patches.

4. Calculate the peak Δ_π at each targeted pressure. For the case of a negative (positive) peak, identify the maximum (minimum) values of π within the depth ranges of [$p-\mathrm{\Delta }p,p$] and [$p,p+\mathrm{\Delta }p$], respectively (green triangles in Fig. 3c), and the depth interval Δp=100 dbar is chosen, considering the general vertical scale (i.e., a few tens of meters) of the eddy-induced subduction features (Zhang et al., 2015; Hosoda et al., 2021); the reference profile is defined by the straight line in between. The anomaly Δ_π (red bracket in Fig. 3c) is defined as the difference between the reference profile and the original profile of π at pressure p (green star in Fig. 3c).

5. Calculate Δ_AOU using the same method, independent of Δ_π.

6. The thresholds used to determine whether the signals meet the criteria of a subduction patch or not were set to −10 µmol/kg for Δ_AOU and ±0.05 kg/m³ for Δ_π following Llort et al. (2018).

Figure 3 Vertical property distributions of profile no. 234 (on 24 June 2016) of float MD5904034 (the black dot in Fig. 2a) demonstrate subduction patches observed by the BGC-Argo floats. (a) The profiles of potential density (σ, blue line) and potential spicity (π, dotted red line), (b) the profiles of DO (blue line) and AOU (red dotted line); (c) the same spicity profile as in (a), which is used to demonstrate the steps to detect subduction signals described in Methods. Note that the red dots in each panel represent the raw field observations; the overlaid red curves are the three-bin running averages to remove sharp noises or spikes, which are used to calculate the anomalies in AOU and π; and the black line represents the MLD. The yellow shading in (a) and (b) highlight the subduction features identified using the detection method in (c).

The refined algorithm presented here had improved performance for detecting subduction patches in these BGC-Argo profile data compared to that used in previous studies (Llort et al., 2018) (see Fig. S1). The main difference in our approach is in selecting the frame of reference for identifying AOU and π anomalies from irregular features in “typical” vertical profiles.

The sensitivity of the method to the interval of Δp (in step 4) was investigated by varying Δp between 70 and 130 dbar (see Table S1). For Δp of 100±3 dbar (i.e., 97, 98, 99, 101, 102, and 103 dbar), fewer than seven (≤2 %) subduction patches were missed, and the resulting Δ_AOU and Δ_π show a RMSD of ≤ 3.8 µmol/kg (≤8.3 %) and ≤0.03 kg/m³ (≤ 9.2 %). More details are provided in Text S1. The sensitivity analysis suggests the validity and robustness in the choice of Δp of 100 dbar. After verifying that our approach better captured subduction indicators in a subset of BGC-Argo data from this region, the algorithm was applied to all profiles to identify the locations, depths, timing, and strengths (i.e., Δ_AOU, Δ_DO, and Δ_π) of the subduction patches.

2.2.2 Quantification of oxygen export

For all the subduction patches identified using the method developed above, we obtain a first-order estimate of oxygen export based on the DO anomalies (ΔDO) with the assumptions that (1) the surface processes initiating these subduction events generated similar levels of DO (i.e., surface phytoplankton production) and (2) the water parcels containing this DO are subducted into the ocean's interior.

We estimated the average oxygen inventories within the water column based on the BGC-Argo profiles. We calculated DO inventories (per m²) through these features in two ways: by integration of the anomaly above the estimated baseline (Eq. 1) and by using the anomaly peak height (Eq. 2) (see Fig. 3c).

The equation for the integrated estimates for each profile is

$$\begin{array}{}\text{(1)}& {\text{Oxygen Inventory}}_{\text{IA}}(\mathrm{g}{\mathrm{O}}_{\mathrm{2}}/{\mathrm{m}}^{\mathrm{2}})={\sum }_{z=p\mathrm{1}}^{z=p\mathrm{2}}\mathrm{\Delta }{\text{DO}}_z,\end{array}$$

where ΔDO_z is the DO anomaly at depth z within the water column of the subduction patch, and the integrated areas (IA) of DO anomalies are converted from µmol/kg to mg/m² based on seawater density.

The inventory calculated using the peak height (PH) approach is

$$\begin{array}{}\text{(2)}& {\text{Oxygen Inventory}}_{\text{PH}}(\mathrm{g}{\mathrm{O}}_{\mathrm{2}}/{\text{m}}^{\mathrm{2}})={\mathrm{\Delta }}_{\text{DO\_peak}}\times H,\end{array}$$

where H is the thickness (i.e., vertical height between the green triangles in Fig. 3c, in units of m) of the subduction patch and the Δ_{DO_peak} is the maximum anomalous value of DO converted to mg/m² as above. The oxygen inventory using the peak height method represents the maximum potential of the anomalous DO inventory within the subduction patch.

3 Results and discussion

3.1 Case study – detecting subduction in BGC-Argo datasets

Subduction associated with eddy pumping is a recognized important contributor to the transfer of carbon and other materials from the surface euphotic layer to the ocean interior (McGillicuddy, 2016; Bord et al., 2019), but investigating the spatial distributions, physical dynamics, and biogeochemical consequences of these episodic small-scale processes is difficult. The BGC-Argo program provides an exceptional data resource for this purpose (Claustre et al., 2020; Chai et al., 2020), but detecting subduction signals where differences among water masses are small is challenging.

Subduction patches below the seasonal and permanent pycnoclines can be identified in vertical profiles by anomaly matrices of temperature, salinity, and dissolved oxygen (DO). Examples of these events are illustrated in time series from the BGC-Argo profiling float (MR2901556), between 28 July and 18 August 2014 (Fig. 4). Positioned on the southern perimeter of the Kuroshio Extension region, the surrounding ocean conditions are less energetic, with fewer eddy activities and small sea level anomalies. During the observation period, the float was trapped in a warm-core eddy (Fig. S2) and moved from the margin of the eddy to the eddy core, as indicated by the gradual increase in depth of the anomalous patch (right panels of Fig. 4). Due to the vertical cruising of the Argo profiler in an environment with velocity shear, it may sample different parts of the same subduction patch, as indicated by slightly different depths, and anomalous spicity and oxygen concentration. Here, intermittent patches of elevated spicity (π), lower AOU, and greater dissolved oxygen are visible in the upper 600 dbar (boxes 1–3, Fig. 4). Potential spicity (π), a parameter dependent on pressure, temperature, and salinity (Huang et al., 2018), is a sensitive indicator of water mass differences. AOU is the difference between the measured dissolved oxygen concentration and its equilibrium saturation concentration in water with the same physical and chemical properties. It reflects the degree of progressive microbial decomposition of organic matter since the water was last at the surface in contact with the atmosphere (Garcia and Gordon, 1992; Sarmiento and Gruber, 2006). Despite this oxygen consumption, these injected waters retain excess net oxygen concentrations relative to the surrounding mesopelagic zone (Fig. 4d).

Figure 4 Trajectory of float MR2901556 between 28 July 2014 (Station no. 123) and 18 August 2014 (Station no. 144) (a) and its time series of π (b) AOU (c), and DO (d). Vertical lines in (b), (c), and (d) represent the BGC-Argo profiles, and the section distance along the x axis is the path distance from Station no. 123 (the red star in a). The three boxes (Box1, Box2, and Box3) in (b, c, d) outline the coherent anomalous features in π, AOU, and DO, which were identified as subduction patches following the detection procedure in Sect. 2.2.1. The red lines in (b, c, d) indicate the MLD, and the horizontal black lines are the isopycnals. Anomalies of magnitude less than −10 µmol/kg for Δ_AOU and ±0.05 kg/m³ for Δ_π (e.g., at section distances of ∼ 25, 125, 175, 275, and 475 km) were below our conservative thresholds for identifying intrusions (−10 µmol/kg for Δ_AOU and ±0.05 kg/m³ for Δ_π).

Llort et al. (2018) successfully identified eddy subduction in BGC-Argo data from the Southern Ocean using anomalies in spiciness (Flament, 2002; Huang, 2011; McDougall and Krzysik, 2015), a parameter derived from a different function of pressure, temperature, and salinity than potential spicity (Huang et al., 2018). However, we found that spiciness frequently missed signs of subduction while misidentifying other signals as subduction, and the 20-bin method used by Llort et al. (2018) significantly dampened the subduction signals in our data. Potential spicity (π) (Huang et al., 2018), on the other hand, greatly improves the ability to distinguish among similar water masses due to its orthogonal coordination with density, a feature that spiciness lacks. This added sensitivity revealed reliable signals of subduction in these BGC-Argo data. The algorithm based on peak detection here shows better capabilities in capturing and quantifying the subduction signals (see Methods, Fig. S1).

For the same subduction event, continuous subduction patches are expected to be identified from the Argo profiles. The discrete anomalous π and AOU signals, highlighted in boxes 1–3 in the example time series (Fig. 4a–d; 31 July, 10 and 12 to 15 August), indicate that they stemmed from distinct subduction events, opportunistically captured by this BGC-Argo float. The first two anomalies (July and early August) each appeared in only a single profile, perhaps indicating a limited spatial scale of these subduction events. In contrast, the mid-August anomaly persisted over four consecutive profiles. We further examined the corresponding time series of temperature, salinity, and potential density and found salinity also showed a similar anomalous signal. As such, we suspect the consecutive subduction patches were most likely from a more sustained or a larger spatial subduction event. It should be noted that the high detection chance of subduction within 21 d was likely given that the float was trapped in a warm-core eddy (Fig. S2), which does not indicate this area is a place where it is easy to observe subduction events (see Sect. 3.2).

3.2 Spatial and temporal distributions of subduction

We used our peak detection algorithm with the π and AOU data and applied it to all 7120 BGC-Argo profiles (2008–2019) in the western North Pacific (Fig. 5). Our algorithm resolved 335 subduction patches, spread over an unexpectedly large area in the western North Pacific. Overall, subduction patches were identified in 288 profiles (4.0 %) (some profiles have multiple patches at different depths), with approximately 83 % of these being concentrated in the Kuroshio–Oyashio Extension region (Fig. 5a). High (≥6 cm) climatologic sea level anomalies (SLAs) and the corresponding distribution of eddy kinetic energy (EKE) are evidence of the strong energetic ocean processes in this region (Fig. 5a and b). By contrast, far fewer subduction patches were identified in the less energetic region to the south of 29^∘ N despite a higher BGC Argo sampling density (Fig. 2b), consistent with eddy-related processes being important for driving these subduction events. Even so, the true frequency of these events across the entire region is certain to have been under-sampled given their small scales relative to the dispersed BGC-Argo float positions.

Figure 5 Horizontal distribution of the BGC-Argo data profiles associated with subduction patches (a, b) between 2008 and 2019 in the western North Pacific. The profiles with detected subduction patches are color-coded by different intervals of depths of the subduction patches (a), AOU anomalies (b), and seasons (c), with percentages of detected patches in each interval annotated. The purple background data in (a) represent all the analyzed profiles as shown in Fig. 2a. The greyscale background map in (a) is the annual mean EKE climatology, with EKE contour lines of 0.3, 0.2, and 0.1 m²/s² shown in magenta, cyan, and black, respectively, and the greyscale background map in (b) is the annual mean SLA climatology, with SLA contour lines of ≥0.06, 0.04, and 0.02 m shown in magenta, cyan, and black, respectively. The seasons in (c) are divided into the spring of March–May, summer of June–August, fall of September–November, and winter of December–February.

Discrete signals of subduction were detected throughout the mesopelagic depth range (∼ 100–1000 dbar), with the majority detected below 300 dbar (green and yellow dots in Fig. 5a). The deepest penetrations (≥450 dbar) occurred largely in areas experiencing the highest EKE while the shallowest (100–300 dbar) were largely restricted to areas with lower EKE (Fig. 5a). Based on 16 years of Argo float data (N=1 226 177) in the global ocean, Feucher et al. (2019) found that the depth of the permanent pycnocline differs between the subtropical (i.e., <35^∘ N) and subpolar (>35^∘ N) regions, with the depth of the permanent pycnocline being 300 and 450 m in the subpolar and subtropical sections of the western North Pacific. Similarly, using the limited BGC-Argo dataset used in this study (Fig. 2), we also found shallower annual maximum MLD in the subpolar section than that in the subtropical section (see Fig. S3). As a result, 56 (16.7 %) and 104 (31.0 %) subduction patches were found to be above and below the depth of the permanent pycnocline (i.e., 450 m) in the subtropical section; and in the subpolar section, 34 (10.1 %) and 141 (42.1 %) subduction patches were above and below the permanent pycnocline (i.e., 300 m). Overall, roughly half (52 %) of the detected subduction signals were below 450 dbar in this region of the western North Pacific, while 22 % penetrated far deeper (up to 800 dbar; Table S2).

There is a distinct seasonality in subduction, with most (∼ 70 %) signals being observed between March (the maximum) and August (Figs. 6 and S4), with the peaks appearing in early March. Since these anomalous patches can be observed after their formation, there is a delayed period between the peak of formation and the peak of observation.

Ocean eddies are sustained by many processes including small-scale mixed-layer instabilities set up by large-scale atmospheric forcing in winter (Sasaki1 et al., 2014); in particular, the extreme mixed-layer deepening events are attributable to the accumulation of excessive surface cooling driven by synoptic storms in connection with cold-air outbreaks (Yu et al., 2020). Each cooling episode brought by the synoptic storm may lead to the formation of subduction patches due to strong mixed-layer deepening. Thus, subduction patches are primarily formed during winter and spring when deep mixing processes occur. Although only 8.3 % of the total profiles were obtained in March, they accounted for 17.3 % of all observed subduction patches (Fig. S4a); correspondingly, the monthly subduction detection rate (i.e., the number of profiles with identified subduction patches divided by the total number of profiles available) was the highest in March, at ∼ 10 % (Fig. S5). In a pioneering work, Stommel (1979) argued that a demon works in the ocean by selecting the later winter (typically later March in the Northern Hemisphere) water mass properties and injecting them into the subsurface ocean. This mechanism is now called the Stommel Demon in dynamical oceanography (Huang, 2010). The high detection rate of episodic subduction patches in March was consistent with observations of large-scale subduction in this region during late winter because mesoscale and sub-mesoscale eddy activities are prevalent when large-scale subduction occurs (Qu et al., 2002; Qiu et al., 2007; Nishikawa et al., 2010; Liu and Huang, 2012; Zhang et al., 2014; Xu et al., 2014). The March-to-August time frame also coincides with the onset and establishment of warming-induced shoaling of the mixed-layer depth, when winter-subducted waters are less likely to be re-entrained into surface waters by winds (Dall'Olmo et al., 2016; Palevsky and Doney, 2018). Indeed, based on the entire BGC-Argo dataset in Fig. 2, we found that the monthly MLD reached a maximum in February and March and then decreased until August (Fig. S3). It should be noted that, despite the number of subduction patches identified in the time frame of April–August being slightly larger than those in September–December (Fig. S4a), the detection rates did not vary much between these time frames (Fig. S5). In contrast, comparatively few (3.0 %) of the subduction patches were detected in January and February, in which time the detection rates were also low (< 2 %, Fig. S5). Although specific timelines between the observed subduction patches and their formation could not be determined, it is reasonable to anticipate that more energetic winds and the accumulated strong heat loss during mid-winter contributed to the peak in subduction signatures observed in March. However, there were no spatial patterns of the subduction patches detected in each season (Fig. 5c). The current BGC-Argo profiling asset is not sufficient for studying how those subduction patches change on interannual scales.

Figure 6 Temporal distribution of the number of patches (a), integrated AOU anomaly (b), integrated π anomaly (c), and integrated DO anomaly (d), by the Julian day based on 7-point smoothing. Spicity in subducted patches can be lower or higher than the surrounding waters, resulting in negative Δπ (red lines) or positive Δπ (blue lines) anomalies, respectively (see text in Sect. 3.4). The integrated anomalies indicate the significance and prevalence of the episodic subduction events over time.

The Kuroshio–Oyashio Extension zone lies between the subtropical and subpolar gyres in the North Pacific, and it is a recognized hot spot for water mass exchange via eddy transport (Yasuda et al., 1996; Talley, 1997; Joyce et al., 2001; Zhang et al., 2014; Xu et al., 2016) and substantial ocean-to-atmosphere heat flux (Jing et al., 2020). It is not surprising then that the majority of subduction signals were observed in this region in spite of less float coverage (Fig. 5). In the southern area of the Kuroshio Extension, where the subtropical mode water is formed, the MLD could reach 300–400 m between January and March (Fig. S3). Subduction patches are formed after the water parcels are detached from the base of the mixed layer, and they could be identified below the base of the winter mixing layer (Fig. 5c). Large-scale circulation and seasonal variability in the mixed-layer depth here typically result in late-winter subduction of subtropical mode waters (Qiu et al., 2007; Oka et al., 2009; Oka and Qiu, 2012; Xu et al., 2014, 2016), and sharp horizontal density gradients can enhance strong vertical exchanges (Marshall et al., 1993; Hurlburt et al., 1996; Liu et al., 2012; Ma et al., 2017). Rapid heat loss to the wintertime cool, dry continental air masses flowing across the Kuroshio–Oyashio Extension erodes the seasonal thermocline to its maximum depth in February–March (Cronin et al., 2013), the latter portion in which the subduction patches were most frequently observed (Fig. 6).

Ascertaining the frequency and spatial extent of these lower-latitude episodic events will be important for establishing their overall contribution to the transport of surface waters into the mesopelagic zone, but this goal is challenged by the presently limited distribution of BGC-Argo floats. It may be possible to obtain a first-order estimate of their frequency by linking the subduction signals here to surface-expressed indicators of mesoscale circulation processes. Moreover, our findings suggest that spicity should be adopted more generally in probing BGC-Argo datasets to improve our understanding of the spatial and temporal distribution of subduction processes.

3.3 Properties of subduction patch

Beyond being a water mass indicator, AOU is a proxy for cumulative net community respiration and a sensitive indicator of carbon export in the upper mesopelagic zone (Emerson et al., 2001; Pan et al., 2014; Catala et al., 2018; Bushinsky and Emerson, 2018). This export comprises remineralized carbon as well as dissolved and slowly sinking particulate organic matter carried by the subducting waters (Stukel et al., 2017). The magnitude of AOU may be used as an indicator of the time since subduction, with the first-order assumption being that the larger-scale processes initiating these subduction events generated similar surface production. Values of Δ_AOU at the anomalous peak depth ranged between −10 (the minimal threshold used) and −81 µmol/kg (Fig. 7a). This proxy was highly variable over the space–time domain, similarly to the variations in Δ_π (Fig. 7b). In general, 61.7 % of the subduction patches had Δ_AOU values in the range of −30 to −10 µmol/kg with the remainder having greater oxygen depletions (i.e., $\le -\mathrm{30}$ µmol/kg) (Fig. 5b). Water masses subducted below 450 dbar (i.e., the permanent pycnocline) had an average AOU anomaly of −25.7 ± 15.3 µmol/kg.

Figure 7 Vertical spatial distribution of the detected subduction patches in the western North Pacific, color-coded by the magnitudes of the subduction strengths in terms of the AOU anomaly (a) and π anomaly (b).

There was no clear relationship between the depth of subduction and Δ_AOU (Fig. 7a), suggesting either surface conditions (e.g., water temperature, primary productivity) were substantially different when the seawater parcels were subducted or these signatures stem from non-systematic differences in the time since subducted waters were last at the surface. On the other hand, it is noted in Fig. 7 that the depth positions of the subduction patches appear to somewhat extend from northeast to southwest and are deeper along the isopycnal surface as illustrated in Fig. 1. This phenomenon is clearly shown when averaging the depth of subduction patches both latitudinally and longitudinally (Fig. S6). Along the latitude, despite a few deep subduction patches identified at 42–43^∘ N (at around 550 m), the mean depths of the subduction patches show a clear increasing pattern from latitude 37–42^∘ N to a latitude of 32–37^∘ N, i.e., 300 m vs. 500 m. However, the depth positions tend to become shallower and shallower south of 32^∘ N. Along the longitude, the depth positions generally appear to be deeper from east to west. As such, it is most likely that the subduction that occurred in the northern KE (37–42^∘ N) could travel southwestward from shallow to deep depth, and these waters could reach 32^∘ N. The increasing depth positions of subduction patches from 26 to 32^∘ N tend to suggest the gradually downward movements of the subducted water masses carried by the general trend of the anticyclonic gyre-scale circulation, yet further investigation is needed.

In the subpolar region, for the subduction patches identified above and below the depth of the permanent pycnocline (i.e., 300 m), respectively, the averaged Δ_AOU values are −32.9 and −25.8 µmol/kg, averaged Δ_DO values are 42.5 and 32.5 µmol/kg, and averaged thicknesses (i.e., vertical extension of the subduction patch) are 127.5 and 126.6 m (Table 1). In the subtropical region, the depth of the permanent pycnocline was deeper (i.e., 450 m) and the subduction patches above and below this layer were associated with a mean Δ_AOU of −27.2 and −28.5 µmol/kg, mean Δ_DO of 31.2 and 36.4 µmol/kg, and mean thickness of 128.7 and 128.1 m (Table 1). In general, the vertical extension (i.e., thickness) of the subduction patches identified in each layer and in each region did not vary much between 126.6 and 128.7 m. The mean Δ_AOU and Δ_DO were stronger above the depth of the permanent pycnocline than those below the depth of the permanent pycnocline in the subpolar region, yet the opposite case is shown for the subtropical region, where the mean Δ_AOU and Δ_DO were weaker above the depth of the permanent pycnocline than those below the depth of the permanent pycnocline. Interestingly, it is noted that the mean Δ_AOU and Δ_DO values in the subtropical region below 450 m were also weaker than those in the subpolar region above 300 m, which further supports the potential northeast-to-southwest pathway of subducted waters shown in Fig. 7.

Table 1 Statistics of the subduction patches and the associated oxygen exports into the ocean's interior. See Sect. 2.2.2 for details on the calculation of the DO inventory. Note that these statistics are based on the subduction patches identified, without considering their episodic characteristics and spatial and temporal inhomogeneity. Note that the DO inventories are the average values of the amount of oxygen injected into the ocean interior by an episodic subduction event.

Most subduction patches with strong AOU anomalies were observed between March and August (particularly March; see Fig. S4), after the seasonal mixed layer began to shoal, consistent with expected higher levels of phytoplankton production, which results in a greater degree of respiration in the subducted waters. More respiration means a great degree of oxygen consumption and thus a more negative offset from the surface-saturated concentrations before subduction. Only 0.6 % of the total subduction patches had Δ_AOU values of $\le -\mathrm{30}$ µmol/kg in January and February (Fig. S4b). It should be noted that Δ_AOU also strongly depends on the surface water temperature (which determines the solubility of oxygen) when it is subducted. The π anomalies show similar variation patterns with months (peaked in March), with stronger Δ_π coupled with stronger Δ_AOU (Fig. S4c).

Table 2 Summary of the subduction patches associated with positive and negative π anomalies; bold numbers indicate statistics of the sum and mean based on absolute values of π anomalies.

3.4 Oxygen injections into the twilight zone

Global ocean inventories of oxygen have been decreasing, and current climate models predict this trend is likely to accelerate over the next century (Oschlies et al., 2018). However, these models suffer from considerable gaps in understanding, one of which is the absence of small-scale transport processes such as the events captured here (Oschlies et al., 2018). The average residual DO enrichment in the subduction patches, defined as the difference in DO concentrations within and adjacent to the subducted waters, was 34.5±19.8 µmol O₂/kg, with levels as high as ∼ 88 µmol O₂/kg below 450 dbar during March. These differences reflected ∼ 20 % higher oxygen concentrations than in the surrounding mesopelagic waters. The integrated DO enrichment reached a maximum in March (Fig. 6d). Based on these residual excess oxygen concentrations, the oxygen inventory within these features was estimated to be on the order of 64 to 152 g O₂/m² (Eqs. 1 and 2). Specifically, the DO inventories below the permanent pycnocline in the subtropical and subpolar regions were on the order of 64.3–161.5 and 61.2–142.1 g O₂/m², respectively (Table 1). Note that the DO inventories here are the average values of the amount of oxygen injected into the ocean interior by an episodic subduction event. These oxygen inventories may represent a significant source of ventilation to our study region.

Co-injection of oxygen below the permanent pycnocline by eddy pumping has not been given close consideration in previous studies, largely because it is less relevant for high-latitude, oxygen-rich waters. However, weak ocean ventilation in the tropical and subtropical mesopelagic zone is leading to declining oxygen concentrations (Karstensen et al., 2008; Oschlies et al., 2018; Robinson, 2019) and expansion of oxygen minimum zones in many regions of the oceans (Stramma et al., 2008; Breitburg et al., 2018). These episodic, dispersed subduction events likely represent a significant source of ventilation to help offset the de-oxygenation phenomenon and to support the expected climate-driven effects of increasing temperature on the metabolic oxygen demand of mesopelagic organisms (Wohlers et al., 2009). Enriched oxygen supplies into the mesopelagic zone will also influence remineralization rates of sinking particulate organic carbon in the ocean's twilight zone (Buesseler et al., 2007; Steinberg et al., 2008), affecting carbon sequestration timescales. Current global-scale biogeochemical models are too coarse to capture the effect that these sub-mesoscale processes may have on mesoscale oxygen variability (Takano et al., 2018) or to account for this additional oxygen supply. Overall, the intensity of these export events below the permanent pycnocline is remarkable, and they should be adequately considered in biogeochemical models.

Eddy-associated pumping is also one of several processes contributing to net global ocean carbon export (McGillicuddy, 2016; Boyd et al., 2019), but its importance is generally thought to be comparatively small because the relatively shallow penetration leads to shorter carbon sequestration times (Lévy et al., 2001; Karleskind et al., 2011a, b; Omand et al., 2015; Nagai et al., 2015; Boyd et al., 2019). That is, much of the carbon “exported” to the upper mesopelagic zone over spring and summer is returned to the atmosphere by deep winter mixing. At higher latitudes, where eddy pumping has been most studied, subduction must extend up to > ∼ 1000 m to reach below the permanent pycnocline (Palevsky and Doney, 2018; Boyd et al., 2019). However, the permanent pycnocline in the western North Pacific is much shallower – on the order of ∼ 300–450 dbar (Qiu and Huang, 1995; Feucher et al., 2019) – and most of the observed subduction signals here extended far below this depth (Tables 1 and S2). Thus although the subduction depths shown here are similar to those observed at higher latitudes, they represent much longer carbon sequestration timescales than those previously associated with eddy pumping (Boyd et al., 2019). As such, in addition to oxygen exports, the observed subduction patches seem to also transport large amounts of carbon into the ocean interior particularly below the permanent pycnocline. However, the lack of carbon measurements on the BGC-Argo floats used in this study impeded us from quantifying the carbon inventory within the subduction patches.

Because the BGC-Argo profiler only captures snapshots of subduction events, it is impossible to quantify the vertical transporting rate, which is needed to quantify export fluxes, of subduction from the BGC-Argo float data alone. Alternatively, the lifetime of subduction patches could be used to infer subduction rates, yet due to the dynamics and episodic characteristics of eddy subduction, there are currently no estimates of how much time these water masses maintain differentiated properties in the mesopelagic zone for, and there are numerous physical and biogeochemical processes influencing them.

3.5 Surface forcing of subduction

The AOU, DO, and π anomalies were integrated within the study domain over the year to assess the extent of subduction in the western North Pacific (Fig. 6, Table 2). π anomalies were divided into negative or positive Δ_π – i.e., π being greater or less than that in surrounding waters – which can suggest their modes of formation. Negative Δ_π would correspond with the subduction of colder and less saline waters, such as along the edges of cyclonic eddies, while positive Δ_π would be associated with the eddy pumping of warmer-core, anticyclonic eddies. The subduction patches were clearly dominated by negative Δ_π, and more negative Δ_π values corresponded with much larger Δ_AOU and Δ_DO values (Fig. 6, Table 2), suggesting they were associated with cyclonic, cold-core, upwelling-dominated eddies that have higher oxygen solubilities, higher nutrient flux to the surface, and thus higher plankton production. Conversely, the association of lower Δ_AOU and Δ_DO values with positive Δ_π values would align with the lower oxygen solubility, nutrient flux, and plankton production expected for warmer-core, downwelling, anticyclonic eddies. Moreover, the majority of deep intrusions had negative Δ_π values (Fig. 6, Table 2) consistent with colder waters following deeper isoclines. In contrast, anticyclonic eddies would push warm, lower-oxygen, and less-biomass-containing waters to shallower depths. These findings suggest that tracking the activity of cyclonic eddies in regions with shoaling permanent pycnoclines (Chelton et al., 2011; McGillicuddy, 2016) may be particularly important for quantifying these deeper subduction processes.

The findings here indicate that eddy-associated subduction is an important mechanism driving oxygen enrichment below the permanent pycnocline across the western subtropical Pacific region, particularly near the Kuroshio Extension (KE). Moreover, the abundance of these discrete, small-scale subduction events is almost certainly under-sampled in the BGC-Argo dataset. The frequency of this subduction is expected to vary as the KE oscillates between two dynamic states – quasi-stable and unstable – linked to the Pacific Decadal Oscillation (PDO) or North Pacific Gyre Oscillation (NPGO) (Di Lorenzo et al., 2008). When quasi-stable, the KE jet shifts north and generates less eddy activity than the unstable, highly meandering, southward KE jet, which reduces eastward transport and sharply increases eddy kinetic energy (Qiu and Chen, 2010; Lin et al., 2014). Superimposed on these KE oscillations has been an increase in the ratio of cyclonic to anticyclonic eddies associated with a climate-driven intensification of tropical storms in the western Pacific and the multidecadal trend of acceleration in Kuroshio flow (Zhang et al., 2020), suggesting that the importance of eddy-associated subduction processes in this region has been increasing and may continue to increase in the future. This linkage needs to be considered in designing future ocean observation programs and the modeling of global biogeochemical cycles to adequately capture the damping effects that eddy-associated subduction may exert on increasing atmospheric CO₂ and de-oxygenation in the tropical and subtropical ocean.

4 Conclusions

Biogeochemical measurements obtained from the BGC-Argo float data provide new insights into the small-scale vertical water mass exchange in the ocean. In particular, spicity and AOU are key parameters in capturing the episodic subduction events and their significance. Although these floats cannot capture the full pathways of subduction, they provide first-hand data on the locations, depths, timing, and strengths of episodic subduction patches. Here we analyze float data in the western North Pacific and show significant subduction export of dissolved oxygen to the mesopelagic zone particularly below the permanent pycnocline; thus, the BGC-Argo data available over the global oceans can be used to extend the current study to other oceanic regions. Carbon measurements are needed to quantify the carbon export associated with the subduction patches. These two factors – carbon export and re-oxygenation – would help to offset the apparent budget imbalance between the biological gravitational pump and mesopelagic carbon demand and would support the increasing metabolic oxygen demand of mesopelagic organisms as ocean warming continues.