The Physical-Biogeochemical Responses to a Subsurface Anticyclonic Eddy in the Northwest Pacific

Introduction

Mesoscale eddies are common physical processes in oceans characterized by energetic, swirling, time-dependent circulations that are approximately 100 km wide and include both warm-core and cold-core eddies (Rhines, 2001). They play an important role in biogeochemical cycles. A well-known feature on interannual to decadal timescales is the bimodality of the Kuroshio path to the south of Japan, which flows either in a large meander (LM) or a non-large meander (NLM) path (Yoshida, 1964; Taft, 1972; Kawabe, 1985, 1995). The drastically different changes that occur in the shapes and locations of these two paths significantly impact the climate, fisheries, navigation safety, and marine resources in the northwestern pacific region.

When the Kuroshio takes an LM path, the enhanced cyclonic eddy near the Kii Peninsula moves southward, with the east side of a long-standing anticyclonic eddy accompanying the Kuroshio meander off Shikoku (Figure 1; Sugimoto et al., 1986; Barreto et al., 2021). However, this anticyclonic eddy should be distinguished from the Shikoku Recirculation Gyre (SRG, 25°−35°N, 132°−140°E) that is associated with basin-scale circulation typically on the order of several cm s^–1 (Imawaki et al., 1997). Specifically, the anticyclonic eddy is smaller (∼200 km in diameter at the velocity core) but has a greater velocity (∼20 cm s^–1) than the SRG (Mitsudera et al., 2001; Li et al., 2014). The vertical and horizontal mesoscale dynamics associated with eddies can affect the rates at which nutrients are supplied to the euphotic zone (McGillicuddy et al., 1998; Sasai et al., 2010), further altering phytoplankton productivity and the export of particles to the deep ocean (Bidigare et al., 2003; Xiu et al., 2011; Xing et al., 2020). Previous studies have shown that 10–50% of the global new primary production is caused by eddy-induced nutrient fluxes (Falkowski et al., 1991), so the processes associated with these two strong eddies play crucial roles in modifying the local-scale nutrient distribution and biological production.

During the past few decades, scientists have paid increasing attention to the coupling relations between physical and biological processes (Lewis, 2002). The mechanisms by which mesoscale eddies influence marine biogeochemical processes can be categorized into three primary processes, namely, horizontal advection, vertical flux, and the influence of eddies on stratification and hence on upper ocean mixing. More specifically, mesoscale eddies can influence biogeochemical processes through the following mechanisms: eddy advection of nutrients, phytoplankton, and other biological tracers by the horizontal rotational velocity of eddies (Chelton et al., 2011); eddy pumping presented as an upward or downward nutrient flux induced by isopycnal displacement as an eddy forms (McGillicuddy et al., 1998, 1999; Siegel et al., 1999); the trapping of ecosystems by nonlinear eddies, of which the rotational velocities are faster than the eddy propagation speed, from their region of formation (Early et al., 2011; Gaube et al., 2014); eddy-induced Ekman pumping, which arises from convergence or divergence induced by wind and surface current interactions (Stern, 1965; Martin and Richards, 2001; McGillicuddy et al., 2007; Mahadevan et al., 2008; Morel and Thomas, 2009); eddy strain-induced pumping resulting from a cross-front ageostrophic secondary circulation, which causes upwelling/downwelling along the light/dense side of the front, generating upward nutrient fluxes and increased chlorophyll (CHL) concentrations along the light side of the front (Zhang et al., 2019); and the influence of eddies on mixing through alterations to the mixed-layer depth that impacts nutrient fluxes and, potentially, also modulates light availability (Levy et al., 1998; Gaube et al., 2013).

Chen et al. (2007) found that the cold eddy generated in the Luzon Strait can significantly increase the primary productivity. Lin et al. (2010) studied a phytoplankton bloom produced by an anticyclonic eddy injection in the oligotrophic area of the northern South China Sea. Xiu and Chai (2011), using a coupled three-dimensional physical-biogeochemical model, concluded that the modeled depth-integrated (0–125 m) CHL, zooplankton, new production, and silicate uptake were significantly enhanced in the cyclonic eddies and reduced in the anticyclonic eddies. Kouketsu et al. (2016) reported that high (low) surface area-averaged CHL concentrations were frequently observed in the core of cyclonic (anticyclonic) eddies in the Kuroshio Extension region. Based on historical cruise observations with mesoscale eddies detected by altimeters in the North Pacific, Huang and Xu (2018) showed that the depth-integrated CHL anomalies in the euphotic layer of eddies are less than ± 10%, and the integrated responses of CHL to eddies decline from ocean boundaries to the oligotrophic open ocean in the North Pacific. Using biogeochemical (BGC-Argo) floats, Xiu and Chai (2020) showed that mesoscale eddies can significantly affect the subsurface CHL maximum and subsurface biogenic particles in the central North Pacific Subtropical Gyre.

Critically, most previous studies have focused on the effect of surface eddies on CHL, or failed to distinguish different eddies, while two types of anticyclonic eddies can be differentiated based on the vertical distribution of their hydrographic signals: regular surface eddies and subsurface-intensified eddies. Regular surface eddies are characterized by expressing the largest isopycnal deformation near the surface, either upward or downward, depending on eddy polarity. Subsurface-intensified anticyclonic eddies, however, appear as positive sea level anomalies with deepening of the main pycnocline and shallowing of the seasonal pycnocline (Sweeney et al., 2003; McGillicuddy et al., 2007; Assassi et al., 2016); similar to surface cyclonic eddies, these eddies can also result in phytoplankton blooms and enhance production via the convex of the isopycnals. For example, in the southern Bay of Biscay (Rodriguez et al., 2003), a twofold enhancement of phytoplankton biomass was detected in the subsurface CHL maximum layers at the core of an anticyclonic eddy, and this enhancement was related to the injection of nutrients into the euphotic layer due to pycnocline doming. There has been a long-standing anticyclonic eddy exists in the region south of Japan; however, due to the limited availability of temporally continuous biogeochemical parameter profile data, the three-dimensional biogeochemical structure of this eddy and its contributions to biogeochemical cycles have never been systematically investigated. A complete depiction of the biogeochemical structure of this subsurface eddy is critically required, especially if eddy-associated biogeochemical cycles are to be understood. In this study, three-dimensional biogeochemical model data were used to evaluate the impacts of this long-standing eddy. We first tracked the subsurface anticyclonic eddy and then analyzed the corresponding physical and biogeochemical responses. The remainder of this paper is organized as follows: Section “DATA AND METHODS” briefly describes the information of the data and methods used in this study. The physical and biogeochemical properties of the anticyclonic eddy are in presented in Section “Physical Properties” and Section “Biogeochemical Properties,” respectively. Section “Chlorophyll Aggregation Below the ZEU in Winter” includes the mechanisms by which the anticyclonic eddy-induced CHL aggregates below the euphotic layer in winter. Summary and discussion are provided in Section “SUMMARY AND DISCUSSION.”

Data and Methods

BGC-Argo

Measurements from two BGC-Argo floats (ID: 2902750 and 2902754) were obtained from the Coriolis Global Data Assembly Centre (GDAC) FTP server¹ to validate the accuracy of the biogeochemistry model data. In addition to temperature and salinity, these two floats simultaneously measured CHL concentrations (units: mg m^–3) and particulate backscattering coefficient at 700 nm [bbp (700); units: m^–1]. Nitrate (NO₃; μmol kg^–1) and dissolved oxygen (DO; units: μmol kg^–1) were also measured by float 2902754. The average sampling rates of the two floats were approximately 2 days and ranged from 1-5 days. The two floats sampled from the surface to ∼1,000 m. The observed trajectories of these two BGC-Argo floats and the simultaneous satellite sea level anomaly (SLA) associated with eddy are shown in Figure 2. The BGC-Argo2902750, which was deployed in the region south of Japan on September 9, 2018, was quickly entrained into the 200−300-km anticyclonic eddy, and we studied and stayed in the eddy until the end of April in 2019; this float was successfully recovered in early June (Figures 2A,C). BGC-Argo 2902754 was released to the Northwest Pacific on August 30, 2018 and traveled for more than 2 years until March 12, 2021. From April to November of 2019, the float was rapidly captured by a cyclonic eddy intensified with the formation of the LM path (Figures 2B,D).

Common and specific quality control tests have been performed on the BGC-Argo-measured variables, and we used only data with “good quality” flags. Additionally, the floats profiled at night (around 22:00–0:00 local time) to avoid phytoplankton fluorescence non-photochemical quenching and CHL profiles were dark corrected based on the on-float measurement before deployment, and further corrected based on match-up satellite CHL data using nonlinear regression function (Xing et al., 2020).

Model Data

Hydrological Data

The three-dimensional daily mean fields of temperature: salinity, zonal, and meridional velocities, as well as two-dimensional daily mean fields of mixed layer depth, were obtained from the Copernicus Marine Environment Monitoring Service (CMEMS) global physical analysis and coupled forecasting product. The coupled modeling system being used to provide the analysis and forecast products is the Met Office Coupled Atmosphere-Land-Ocean-Ice data assimilation system (CPLDA). The scientific and technical implementation of the V3.1 system is described in Lea et al. (2015). These data have a horizontal resolution of 1/4 degree and 43 vertical levels and covered the period from December 2015 to December 2020. The model includes multivariate data assimilation, which consists of a SEEK filter analysis of along-track SLA and sea surface temperature (SST), together with in situ temperature and salinity observation profiles. The altimeter reference period for the assimilated SLA is 20 years (Rio et al., 2014). The assimilated SST is taken from the CMEMS Thematic Assembly Centre (TAC) daily level-4 Operational Sea Surface Temperature and Ice Analysis (OSTIA) composite product (Donlon et al., 2012). The model bathymetry is DRAKKAR v3.3, which is based on the ETOPO1 dataset with additional data in coastal regions from the General Bathymetric Chart of the Oceans (GEBCO).

Compared with observation from a subsurface, mooring system was deployed at 25°N, 146°E from April 2017 to June 2018 (Supplementary Figure 1); the product agrees well with it in terms of temporal variation in temperature [correlation coefficient: 0.76; root mean square difference (RMSD): 0.46°C] and salinity (correlation coefficient: 0.66; RMSD: 0.04 psu). Meanwhile, this product could reproduce the important features of the Kuroshio system, including its spatial pattern of sea surface velocity in NLM and LM path south of Japan (Supplementary Figure 2).

Biogeochemistry Data

The CHL (units: mg m^––3) concentration data were downloaded by file transfer protocol (FTP) from the CMEMS² global ocean biogeochemistry analysis and forecast product (001-028), which was produced at Mercator-Ocean (Toulouse, France). This product is displayed on a regular grid at 1/4 degree resolution, with 50 vertical levels ranging from 0 to 5,500 m, on the global ocean, and includes daily and monthly mean fields of the 10 following biogeochemical variables over the period from January 2017 to the present: CHL, nitrate, phosphate, silicate, dissolved oxygen, dissolved iron, primary production, phytoplankton, PH, and surface partial pressure of carbon dioxide. This product is based on the pelagic interactions scheme for the carbon and ecosystem studies (PISCES) biogeochemical model [available on the Nucleus for European Modeling of the Ocean (NEMO) modeling platform³ ] and was forced offline at a daily frequency by CMEMS global physical analysis product coarsened at 1/4 degree, with singular evolutive extended Kalman/incremental analysis updating (SEEK/IAU) data assimilation of global ocean CHL data obtained from satellite observations. In addition, the time series is aggregated in time in order to reach a two full years’ time series sliding window. In this study, the CHL fields were log10 transformed, which accounts for the highly skewed distributions of the untransformed data that can occur in many regions of the global ocean, especially in near-coastal regions (Campbell, 1995). The detailed validation of the model data was conducted in Section “Biogeochemical Properties.”

ZEU Data

The ZEU is the depth of the euphotic layer (m), i.e., this subsurface parameter used in this study was obtained from the Glob-Color dataset⁴. The ZEU was computed from the corresponding merged CHL-OC5 products (Morel et al., 2007) using the following empirical equation with y = log10 (CHL-OC5):

Ancillary Datasets

The following ancillary data were used in this paper; the SLA and absolute dynamic topography (ADT) data used in this study were provided by the Archiving, Validation, and Interpretation of Satellite Oceanographic data (AVISO) from January 2017 to December 2019 with spatial and temporal resolutions of 0.25° × 0.25° and 1 day, respectively; temperature and salinity data provided by the World Ocean Atlas 2018 (WOA18, 1° × 1°,⁵) were used to estimate the background profile; data on the daily mainstream location of the Kuroshio were provided by the Hydrographic and Oceanographic Department of Japan⁶; and 12 h, 0.25° × 0.25° air-sea net heat flux (Qnet) and 10-m wind speed were provided by the ERA5 and were sourced from the European Centre for Medium-Range Weather Forecasts (ECMWF).

Eddy-Tracking and Construction Method

The eddy identification followed in this study mainly refers to the method proposed by Nencioli et al. (2010). This method is based on eddy vector geometry and has great advantages over the commonly used Okubo-Weiss method in terms of its correctness and accuracy (Souza et al., 2011). First, four constraints on the geometry of velocity vectors are defined to determine an eddy center, which was defined as the minimum velocity in a given region that extends up to a certain grid point. For more details, readers can refer to Nencioli et al. (2010). Once the center is identified, we then search for the outermost closed streamline around the center, and a tracking algorithm is conducted. If the velocity magnitudes cease to increase in the radial direction, this location can be regarded as the eddy boundary. We applied this method to obtain eddy parameters, including the location, amplitude, radius, intensity, and boundary locations of the eddy.

In order to obtain the three-dimensional physical and biogeochemical structure of the eddy, we searched the profiles located in the range of twofold radius and calculated the relative zonal and meridional distance, ΔX_E and ΔY_E, respectively, of profiles to the eddy center (defined at ΔX_E = Y_E = 0). All of the profiles were then used for composition. To eliminate the influence of radius on the composite result, ΔX_E and ΔY_E are divided by the corresponding radius R so that they are dimensionless. After that, for profile surfacing at the boundary of eddy, it shall meet the condition $\sqrt{\mathrm{△}{X}_E^2+\mathrm{△}{Y}_E^2}=1$. To obtain the eddy-induced three-dimensional structures in different months, all profiles in the respective months on the eddy-coordinate space (ΔX_E, ΔY_E) were mapped onto 0.1 × 0.1 grid using the objective analysis (Barth et al., 2014). Also, at each depth, anomaly data were treated as outliers and removed if they were more than three times away from either the first or the third quartile. Thus, for each month, a unique three-dimensional composite structure was obtained.

Results

Physical Properties

There was a long-standing anticyclonic eddy in the Shikoku Basin, although its strength and location will be significantly influenced by the Kuroshio paths south of Japan; this eddy is closely related to the SRG, which is associated with the particular ocean bottom topography south of Japan.

From January 2016 to December 2019, this long-standing anticyclonic eddy was detected from ADT and velocity images, with an average radius of approximately120 km and amplitude of 20 cm (defined as SLA difference between the eddy center and periphery), locating around 30°−32°N, 134°−136°E off Shikoku (Supplementary Figure 3). The variations in the radius, amplitude, strength, latitude, and longitude of the center of this anticyclonic eddy over its lifetime are shown in Figure 3. If the southernmost position of the Kuroshio axis is located south (north) of 32°N, it is classified as a meander (non-meander) path (Sugimoto and Hanawa, 2014). With the development of the LM path (gray background in Figure 3C), the radius, amplitude, and intensity (the ratio of the amplitude to the radius) of the anticyclonic eddy increase gradually as the position shifts southwest. While, according to the idea of the self-sustained oscillation mechanism (Qiu and Miao, 2000), the strength of the SRG increases continuously when the Kuroshio is in its straight path state, corresponding to the decreasing trend of the spatially averaged relative vorticity. During the transition period from the NLM to the LM path, a sudden release of velocity shear corresponds well to the weakening of the SRG, which plays a key role in modulating the Kuroshio path variations. We tried to explain that the trends observed in the anticyclonic eddy intensity are opposite to those of the SRG from the perspective of energy conversion.

Previous studies have indicated that baroclinic instability is a necessary component on the formation of the LM path (Yoon and Yasuda, 1987; Tsujino et al., 2006). Eddy energetics analysis is a useful tool for evaluating flow instabilities (Lorenz, 1955, 1960). However, it is difficult to define the background field due to the great variability that exists in the path of the main Kuroshio axis south of Japan. Tsujino et al. (2006) suggested decomposing an instantaneous field into the field at the previous time step as a mean field and the deviation from this previous field as an eddy component. Here, we adopted the same method for the monthly mean velocity and analyzed the local energetics at each time step. The barotropic (BT) conversion rate from the mean kinetic energy (MKE) to eddy kinetic energy (EKE), which indicates the barotropic instability, is defined as follows:

where ρ₀ denotes the density, $\overline{u}$($\overline{v}$) is the mean zonal (meridional) velocity, and u′(v′) is the perturbed velocity in the eddy field. The baroclinic (BC) conversion rate from mean potential energy (MPE) to eddy potential energy (EPE), which is associated with the baroclinic instability, is calculated as follows:

where g is the gravity acceleration, N² is the Brunt-Vaisala frequency squared, and ρ with the overbar and prime denote the mean and eddy components, respectively. More details of the above eddy energy conversion rate equations and the specific calculation process can refer to Tsujino et al. (2006).

The spatial distributions of the two energy conversion rates described above were calculated by Equations (2) and (3), and the energy conversion distribution areas were located near the main axis of the Kuroshio, which may have been due to the relatively large velocity in this region (Li et al., 2014). In addition, the volume-integrated (upper 1,000 m) energy conversion rates were calculated over the area (30°–35°N, 132°–140°E) (Figure 4). We found that the energy conversion rate of BC over the volume integral was larger than that of BT, indicating that the BC instability plays a major role in the formation of the LM path. Sharp increases in these two conversion rates can be seen before the onset of the LM path, demonstrating the transfer of energy from the mean field to the eddy field. It is possible that the enhanced strength of the anticyclonic eddy resulting from the eddy kinetic energy was provided by the dramatic increase observed in the positive BC conversion rate.

After having analyzed the variations in the anticyclonic eddy intensity when Kuroshio takes different paths, we subsequently tried to investigate differences in temperature and salinity structure of it. The composited monthly temperature and salinity anomalies vertical section structures of the studied anticyclonic eddy in winter (February), spring (May), summer (August), and autumn (November) during the LM path period (2018–2019) are displayed in Figure 5, respectively. As the figure shows, this is a typical subsurface anticyclonic eddy (SSAE), characterized by maximum temperature and salinity anomalies core as well as maximum velocity (Supplementary Figure 4) in the subsurface. And the SSAE is fairly strong, with magnitude of maximum temperature and salinity anomalies exceeding 4°C and 0.3, respectively. The comparison results showed that the positive temperature and salinity anomalies above 800 m caused by the SSAE during the NLM path years (2016–2017, Figure 6) had smaller magnitudes and wider horizontal scales of the core; robust negative temperature and salinity anomalies resulting from the influence of the Kuroshio are clearly seen on the right side of the figure in the LM path period (Figure 5, February, May, and August). Notably, significant seasonal variation exists in the central potential density structure in the LM path period; the SSAE shows a normal concave isopycnal structure from January to April and a convex lens structure from May to December but is less pronounced in NLM path periods.

The lens structure of isopycnals is typical for subsurface eddies, and the existence of the low-potential vorticity (PV) inside the eddy center was thought to be the major reason for this structure. Relative vorticity is usually much smaller than planetary vorticity in most of the interior of the ocean, if ignoring the effects of the relative vorticity; the PV has the following expression (Pedlosky, 1996; Stewart, 2008; Nan et al., 2017):

where ρ denotes the potential density. The relative vorticity (defined as ζ$=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$) of the SSAE is negative (clockwise horizontal velocity, Supplementary Figure 4); thus, it always has the opposite sign to the planetary vorticity, which means the PV anomaly caused by relative vorticity always has the same sign with the PV anomaly caused by the lens-shaped density field at the eddy center (Zhang et al., 2017). Therefore, failing to consider the influence of the relative vorticity causes only the amplitude of the absolute PV anomaly to be underestimated, indicating that only some weak subsurface eddies may be missed; at the same time, the stronger SSAE detected by Equation (4) is reliable. From the composite vertical PV sections shown in the bottom row of Figure 5, the seasonal variation in the central potential density structure of SSAE may have been related to the subduction of low-PV North Pacific Subtropical Mode Water (STMW).

Generally, anticyclones deepen the mixed layer depth (MLD), whereas cyclones thin it, with the magnitude of these eddy-induced MLD anomalies being largest in winter (Gaube et al., 2019). Furthermore, the extent to which eddies modulate the MLD is linearly related to the sea surface height amplitude. Eddy-centric composite averages from Figures 5, 6 also reveal that the largest MLD anomalies occurred at the eddy center and decreased with increasing distance from the center. As a result, due to the unique seasonal variations in the central potential density of SSAE, this fairly strong eddy could deepen the mixed layer by an average of more than 100 m from January to April; in other seasons, however, it has the contrary effects.

Biogeochemical Properties

A comparison with existing observations was conducted to demonstrate the accuracy of the model data before using it to study the biogeochemical structure of SSAE. The model could reasonably reproduce the general spatial pattern of CHL, including the relatively high CHL concentrations in the coastal region, the low values in the deep basin, and the westward pointing “arrowhead” shape of the CHL tongue in the tropical band of the Pacific Ocean (Supplementary Figures 5, 6). The modeled vertical distributions trend of CHL, NO_3, and DO agrees well with the compared BGC-Argo, but was overall underestimated or overestimated when compared to BGC-Argo (Figure 7). The temporal variabilities, especially on the seasonal cycle at mid-latitudes and high-latitudes, of model CHL and two BGC-Argos were similar (r = 0.5, p < 0.01; correlation coefficient even reaching.68 above 100 m): spring or winter bloom. The model also showed fairly well temporal correlation with BGC-Argo DO (r = 0.86, p < 0.01) and NO₃ (r = 0.85, p < 0.01). As the Argo observations are discrete in both temporal and spatial, the results of the vertical profiles comparison for this long time-series are valuable, although they have some discrepancies. And, from April to November of 2019, the Argo2902754 was captured by a cyclonic eddy (Figure 2D) with significantly elevated isopycnals (Figures 7E,F), and the model also responded well to this mesoscale eddy. Therefore, it is plausible to investigate the ecological effects of mesoscale eddies using model data.

In addition to the above qualitative comparison of the biogeochemistry parameters, we should also quantify the bias between the model and the observed data. To achieve meaningful quantitative results of variables, we combined measurements of other 9 BGC-Argo floats in the region of Kuroshio Extension (120°−180°E and 25°−45°N, Supplementary Table 1). Figure 8 shows the comparison between the model and profiles of BGC-Argo CHL, NO_3, and DO concentrations. To estimate errors more accurately, the BGC-Argo CHL profiles were firstly corrected based on match-up satellite CHL data using nonlinear regression function (Xing et al., 2020; Xiu and Chai, 2020; Figures 8A,B), and temporally averaged CHL profiles were further compared with the model (Figures 8C,D). Overall, the model and float CHL vertical concentrations distributions compare reasonably well with a correlation coefficient of 0.71, a negative bias of 0.02 mg m^–3 and an RMSE of 0.06 mg m^–3, which means the model tends to underestimate the CHL concentrations; for NO₃ concentrations, the model is in good agreement with a correlation coefficient of.89, a negative bias of 1.89 μmol kg^–1, and an RMSE of 4.99 μmol kg^–1 (Figure 8E). The model tends to overestimate NO₃ for NO₃ < ∼20 μmol kg^–1 and to underestimate high NO₃ concentrations; and for DO concentrations, the model and float are in favorable agreement with a correlation coefficient of 0.92, a positive bias of 15.56 μmol kg^–1, and an RMSE of 18.74 μmol kg^–1 (Figure 8F). The model tends to underestimate DO concentrations at lower levels.

In the region south of Japan, there was a significant seasonal change in CHL, which increases every winter and early spring and decreased in summer (Supplementary Figure 6). Usually, persistent thermal stratification of the upper ocean means that the surface-mixed layer rarely reaches the depth of the nutricline. Therefore, the vertical distribution of CHL is not typically homogeneous and displays a pronounced maximum close to the base of the euphotic zone, which is called the subsurface chlorophyll maximum (SCM) (Figures 7A–D; Anderson, 1969; Furuya, 1990; Dore et al., 2008; Mignot et al., 2011; Xiu and Chai, 2020). Usually, in the open ocean, the depth of SCM varies over a range of 80−130 m, becoming shallower with increasing latitude; while, in near-shore regions, it ranges between 5 and 50 m (Gong et al., 2012). It should be stressed that CHL concentrations are greater in region south of the Kuroshio main axis than near-shore regions at depths below 50 m (Supplementary Figure 6). The SCM exists throughout the year in tropical and subtropical regions. The seasonal variation of the SCM depth is significant; it is shallow in winter and deep in summer, which is proportional to the light intensity (Radenac and Rodier, 1996; Hense and Beckmann, 2008). Other than that, the depth of the SCM is negatively correlated with the depth of the mixed layer (Figures 7A,C; Qiu et al., 2018). Despite being influenced by many other factors, the SCM often results from photo-acclimation of the phytoplankton organisms in stratified waters, which induces an increase in the intracellular CHL in response to low-light conditions (Kiefer et al., 1976; Winn et al., 1995; Fennel and Boss, 2003; Dubinsky and Stambler, 2009).

It is demonstrated via our previous analysis that there is a significant seasonal variation in the dynamic structure of the SSAE; this seasonal variation determines the distribution of biochemical parameters. Therefore, we intended to investigate the seasonal variation of biochemical parameters corresponding to the dynamic structure of SSAE. Only the biogeochemical structures of SSAE during the LM path period (2018–2019) were composited here due to the lack of biogeochemical data before 2017. In the same way as before, the composite three-dimensional CHL, NO_3, and DO anomaly structures of the SSAE from the surface to 150 m in February, May, August, and November during the LM path period (2018–2019) are shown in Figure 9. In February, DO within the SSAE showed positive anomalies due to the concave of isopycnals, while the CHL and NO₃ concentrations also showed positive anomalies rather than their typically observed negative anomalies. The consistently higher CHL levels throughout the top and bottom of the ZEU in the center of the SSAE resulted from enhanced winter mixing; the deepened MLD was favorable to the entrance of nutrients into the ZEU from the bottom of the mixed layer. A convex lens isopycnal structure was observed from May to December in the SSAE, with convex isopycnals occurring above 200 m opposite to February. During this time period, the composite CHL, NO_3, and DO structures within the SSAE showed positive anomalies around the depth of SCM (50−80 m), and elevated CHL and DO levels and reduced NO₃ levels along the periphery of the eddy below the depth of SCM. According to Xu et al. (2019), they proposed two potential generation mechanisms for these CHL rings: horizontal advection and wind-current interaction. The role of advection is excluded at this depth; therefore, these CHL, DO, and NO₃ rings of SSAE may have been related to wind-current interaction and sub-mesoscale process-associated vortex Ross by waves (Zhang and Qiu, 2020). Twofold to threefold increase of CHL was detected in the subsurface CHL maximum layers at the core of the SSAE, and this enhancement was related to the injection of nutrients into the ZEU due to winter mixing and the convex of isopycnals.

Chlorophyll Aggregation Below the ZEU in Winter

The near-surface mixed layer is the conduit by which the atmosphere influences the ocean interior, and, conversely, the ocean modulates fluxes into the atmosphere. In addition, primary production in the ocean is modulated by nutrient and phytoplankton fluxes through the base of the mixed layer and the availability of light (Dawson et al., 2018; Frenger et al., 2018; Song et al., 2018). Although there is no doubt that the mixed-layer processes are important for biological processes, the relationship between mixed-layer properties and biogeochemical parameters inside eddies remains unclear. Above, we mentioned that the depth of the SCM is negatively correlated with the MLD, and is shallow in winter and deep in summer, ranging from 50 to 200 m. Thus, it is amazing that high CHL and total phytoplankton concentrations were found at depths far below the euphotic layer (Figure 10). As we know, there is low CHL and productivity below 200 m (where CHL concentration is typically smaller than 0.1 mg m^–3); thus, it is intriguing to explore this unusual event.

As illustrated in Figure 10, the CHL concentrations around the SSAE showed obvious downward propagation with depth from 1 January to 24 March 2017 (during the NLM period). Due to the CHL, concentrations in this area are relatively high in winter; northern coastal nutrients and plankton are entrained by the horizontal advection of the SSAE, resulting in elevated CHL along the periphery of the eddy in the initial stage; then, stirring caused a gradually uniform CHL distribution to be observed in the eddy. As the eddy intensity increased, due to the combined effect of the central downward displaced isopycnals in SSAE and winter mixing, high CHL concentrations were progressively found at depths of 222, 266, and even 318 m. Based on thickness and duration of high-CHL-level occurrences, we roughly estimated that the vertical downward speed caused by the SSAE was 2 m d^–1; this speed is consistent with the magnitude of strong eddy-induced Ekman pumping (Gaube et al., 2014). Additionally, a similar phenomenon was captured from January to March 2018 (during the LM period, Supplementary Figure 7). However, the high CHL level associated with the SSAE occurred at shallower depths during this year compared to 2017. Based on our previous analysis in Section “Physical Properties,” the strength of the SSAE increases resulting from the dramatic increased positive BC conversion rate during the transition period from the NLM to the LM path. It seems to indicate that the weaker SSAE corresponds to more pronounced CHL downward displacement below the ZEU.

In order to better investigate the differences between these two processes, time series of regionally averaged CHL, total phytoplankton (concentration of phytoplankton expressed as carbon in sea water, Phyc), and MLD during January to April 2017 and 2018 are given in Figure 11. The averaged CHL, Phyc, and MLD values were derived by obtaining the daily fields located within 1.5 times the length of the SSAE radius (1.5 R). The more pronounced CHL downward displacement below the ZEU (at depths of ∼90 m) caused by the SSAE reached 320 m in 2017, resulting in the CHL levels being elevated by approximately threefold at that depth, whereas this downward propagation reached only 260 m and a nearly 2-fold-elevated CHL level in 2018. The CHL concentration was used as the best available proxy for the phytoplankton biomass; it means that this downward displacement process enables more phytoplankton at that depth. And the ultimate downward propagation depth closely fitted the MLD in the SSAE, implying that the vertical displacement of CHL may be mainly explained by the variations in the MLD (Figures 11A,B,E,F). Near-surface CHL is known to be strongly dependent on processes that control the MLD, with wind and heat fluxes having significant roles (Chiswell et al., 2015). Generally, the MLD is directly proportional to the sea surface wind stress and heat flux. To demonstrate two factors causing MLD variations, regionally averaged (over the area from 132°−140°E and 28°−35°N) 10-m wind stress and air-sea net heat flux (Qnet) were calculated in 2017 and 2018, and the results revealed that enhanced ocean heat losses and wind stresses in January and February of 2017 allowed a deeper MLD (Figures 11G–J).

In addition to the two factors causing MLD variations mentioned above, physically, stratification in the summer/fall months is important because it preconditions the development of the convective mixed layer during the subsequent cooling season (Qiu and Chen, 2006; Iwamaru et al., 2010). The upper ocean stratification during the preceding warm season as weaker (stronger) stratification results in a deeper (shallower) MLD in the subsequent winter; here, the July–September value is chosen to represent the condition before the onset of the cooling season. The result also suggested that the upper ocean stratification during the preceding warm season in 2017 was weaker compared to 2018 (Supplementary Figure 8), which was favorable to the development of a deeper MLD in the following winter.

Then, we have concluded that more pronounced CHL downward displacement below the ZEU associated with the SSAE mainly results from the deeper MLD in winter 2017. To further explore the relationship between the MLD and biochemical parameters inside eddies, we plot in Figure 12 the time series of NO₃ and DO averaged in 1.5 times radius of the SSAE in the winter of 2017. The nutricline represents the concentration of one or several nutrients rapidly increasing or decreasing, which depends on the degree of water column stratification. When physical mixing increased in the winter of 2017, the upper mixed layer penetrated the nutricline, thereby providing a source of nutrients to the euphotic zone (Figure 12A); under these conditions, the phytoplankton became elevated throughout the upper water column. Consequently, a CHL bloom was generated in the SSAE in April 2017 (Figure 12C). Conversely, when thermal stratification increases, the upper water column is deprived of nutrients, leading to the progressive deepening of the nutricline closely paralleling the ZEU. Recent studies have suggested that CHL enrichment in anticyclonic eddies in southeastern Indian Ocean is dominated by winter mixing, which enhances the entrainment of nutrients into the euphotic zone from the bottom of the deepened mixed layer (Dufois et al., 2014, 2016). In this eddy case, both the MLDs in and out of the SSAE deepen in winter. However, the MLD in the eddy is deeper (on average, 100 m) than that out of the eddy. The deeper MLD in the eddy, closer to the deep-nutrient-rich water, may possibly lead to more nutrients being entrained into the euphotic zone and, therefore, higher CHL concentrations.

As we discussed in Section “Physical Properties,” the path of the Kuroshio south of Japan had a significant impact on the SSAE. It is worth noting that these two downward propagation processes associated with the SSAE in winter 2017 and 2018 were at different periods of the Kuroshio path; the former occurred during the NLM period on the eve of the path transition, while the latter was during the stable LM period. According to previous studies, the transition process from a straight to a meandering path is modified by stratification and bottom topography. The former enhances a small meander much more than in the barotropic case through shoaling of the interface depth west of Kyushu; this leads to the small meander, formed during the velocity increase, developing enough to trigger the transition to a meandering path. On the other hand, strengthened stratification has the opposite effect on suppressing the formation of a meandering path (Masuda and Akitomo, 2000; Kurogi and Akitomo, 2006; Nagano et al., 2021). To investigate whether the formation of the Kuroshio LM paths south of Japan in 2017 is influenced by stratification, we plot in Figure 13 the time series of the log10 transformed Brunt-Vaisala frequency squared from January to April. The result shows that the Kuroshio south of Japan gradually shifted from the NLM path to the LM path from the beginning of 2017 (Figure 13A). As the stratification was weaker in early 2017 (Figure 13C), the conditions were favorable for the formation of the LM path; after the formation of the LM path in July 2017, upper ocean stratification during the July–September as stronger results in a shallower MLD in the subsequent winter 2018.

To summarize the two processes associated with the SSAE described above, we have attempted to briefly analyze the mechanism by which wintertime CHL aggregates at depths far below the ZEU. The central downward displaced isopycnals caused by the intensified SSAE, and winter mixing deepened to the nutricline due to the thicker MLD and weaker stratification that occurred in winter 2017 (the NLM path period) (Figure 14A). The combined result of these two processes may lead to numerous nutrients and CHL enrichments in the mixed layer, thus generating a CHL bloom in the following April. However, although the SSAE intensified, the shallower MLD and stronger stratification in winter 2018 (the LM path period) limited the depth of downward displacement of isopycnals.

Summary and Discussion

A long-standing anticyclonic eddy was detected from SLA images with an average radius of approximately 130 km and an average amplitude of 25 cm; this eddy was located at around 30°−32°N, 134°−136°E off Shikoku. With the development of the LM path, the radius, amplitude, and intensity of the anticyclonic eddy gradually increased as its position shifted to the southwest. Different from the idea of the self-sustained oscillation mechanism, the enhanced strength of the anticyclonic eddy that resulted from the eddy kinetic energy was provided by a dramatic increased positive BC conversion rate. The composite anomaly structure suggested that this eddy was a subsurface anticyclonic eddy (SSAE), characterized by having the maximum temperature and salinity anomalies core as well as maximum velocity in the subsurface. In addition, there is significant seasonal variation in the central potential density structure of SSAE during the LM period, with a normal concave isopycnal structure observed from January to April and a convex lens structure observed from May to December; this variation may have been related to the subduction of low-PV STMW.

After that, we further investigated the biogeochemical characteristics of the SSAE; the complete knowledge of which is unclear due to the limited availability of temporally continuous biogeochemical parameter profile data. Based on numerical modeling data, the three-dimensional biogeochemical structure of the SSAE and its contributions to the biogeochemical cycles were investigated. Seasonal variation in eddy dynamic structures decisively influenced its biogeochemical structure. In February, DO within the SSAE showed positive anomalies due to the concave of isopycnals, while CHL and nitrate showed positive anomalies rather than the negative anomalies we usually recognized. From May to December, there was a convex lens structure, and the composite structures of CHL, NO_3, and DO within the SSAE showed positive anomalies around the depth of SCM (50−80 m), while elevated CHL and DO concentrations and reduced NO₃ concentrations were observed along the periphery of the eddy below the depth of SCM, which may be related to wind-current interaction and sub-mesoscale process-associated vortex Rossby waves. Twofold to threefold increase of CHL was detected in the subsurface CHL maximum layers at the core of the SSAE, and this enhancement was related to the injection of nutrients into the ZEU due to winter mixing and the convex of isopycnals.

The unusual event that high CHL concentration and total phytoplankton were found at depths far below the euphotic layer in the winter of 2017 and 2018 was investigated. More pronounced CHL downward displacement below the ZEU (∼90 m) caused by the SSAE reached 320 m in 2017, whereas it reached only 260 m in 2018. We found that the ultimate downward propagation depth closely fitted the MLD in the SSAE, implying that the vertical displacement of CHL may be mainly explained by the variations in the MLD. In 2017, the enhanced ocean heat losses and wind stresses in January and February, and the weaker upper ocean stratification during the preceding warm season allowed a deeper MLD.

It is worth noting that these two downward propagation processes associated with the SSAE in winter 2017 and 2018 were at different periods of the Kuroshio path; the former occurred during the NLM period on the eve of the path transition, while the latter was during the stable LM period. As the stratification was weaker in early 2017, the conditions were favorable for the formation of the LM path; after the formation of the LM path in July 2017, upper ocean stratification during the July–September as stronger results in a shallower MLD in the subsequent winter 2018. In summary, the combined result of these two processes: (1) the central downward displaced isopycnals caused by intensified SSAE, and (2) winter mixing deepened to the nutricline due to the thicker MLD and weaker stratification in winter 2017 (during the NLM path period) may have led to numerous nutrients and CHL enrichments in the mixed layer, thus generating a CHL bloom in the following April. Whereas, although the SSAE intensified, the shallower MLD and stronger stratification in winter 2018 (during the LM path period) limited the depth of downward displacement of isopycnal.

Chang et al. (2019) suggested that the Kuroshio LM was speculated to be connected to the low eel recruitment in 2018 because of its occurrence when there were lower overall catches of glass eels, while there were more glass eel catches in 2017. As we discussed above, the mixed layer in the SSAE deepened to the nutricline in the winter of 2017, leading to numerous nutrients and CHL enrichments throughout the mixed layer and thus generating a CHL bloom in the following April. This process may have been connected to the increased glass eel catches recorded in 2017. Thus, the physical processes associated with the SSAE during the winter preceding the onset of the LM path may have largely affected local phytoplankton blooms and fish catches.

Oxygen is essential for the survival of individual animals, and it also regulates the global cycles of major nutrients as well as carbon. The oxygen content of the ocean and coastal waters has been declining for at least the past half-century, primarily due to human activities that have increased global temperatures and nutrients discharged to coastal waters. These changes have accelerated the consumption of oxygen by microbial respiration, reduced the solubility of oxygen in water, and decreased the oxygen resupply rate from the atmosphere to the ocean interior; these effects have broad biological and ecological consequences (Breitburg and Levin, 2018). But due to the downward displaced isopycnals caused by intensified SSAE and deepened winter mixing, the dissolved oxygen concentration was elevated throughout the upper water column (Figure 12B), causing local dissolved oxygen levels to increase, and it is important for the survival of individual animals and cycles of major nutrients as well as carbon. Consequently, the physical and biogeochemical process of the studied SSAE and its contributions to the biogeochemical cycles should be taken into account in the future research; there is more waiting for us to explore.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: ftp.ifremer.fr/ifremer/argo/dac/, https://hermes.acri.fr/index.php, https://www1.kaiho.mlit.go.jp/KANKYO/KAIYO/, and http://marine.copernicus.eu.

Author Contributions

YD conducted the analyses and wrote the manuscript. FY read through the manuscript and provided revisions. QR, FN, and RW participated in the observation data collection and formal analysis. YSL and YT participated in the discussions and improvement of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Global Climate Changes and Air-sea Interaction Program (Grant No. GASI-02-PACSTMSspr), NSFC-Shandong Joint Fund for Marine Science Research Centers (Grant No. U1406401), and the National Natural Science Foundation of Shandong Province (Grant No. ZR2020MD057).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to thank Yann Drillet, Marc Tressol, Miloud Rached, Matthieu Clavier, David, and Elena from Copernicus Marine Service for assistance in processing the Global Biogeochemical Analysis and Forecast product, which was obtained from the Copernicus Marine Environment Monitoring Service (CMEMS; http://marine.copernicus.eu). We would like to thank China BGC-Argo Project led by Fei Chai at Second Institute of Oceanography and JAMSTEC (Japan Agency for Marine-Earth Science and Technology) project for the BGC-Argo dataset, which were downloaded from the Coriolis GDAC FTP server (ftp://ftp.ifremer.fr/ifremer/argo). The altimeter products are produced by the Ssalto/Duacs and distributed by the AVISO, with support from the CNES. We would also like to thank the two reviewers FZ and XC for their careful reading and for providing very constructive comments that improved the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.766544/full#supplementary-material

Footnotes

1. ^ ftp.ifremer.fr/ifremer/argo/dac/
2. ^ http://marine.copernicus.eu/
3. ^ https://www.nemo-ocean.eu/
4. ^ https://hermes.acri.fr/index.php
5. ^ https://www.nodc.noaa.gov/OC5/woa18/
6. ^ https://www1.kaiho.mlit.go.jp/KANKYO/KAIYO/

References