Tozuka et al.
Progress in Earth and Planetary Science(2022) 9:70
https://doi.org/10.1186/s40645-022-00532-0

RESEARCH ARTICLE

Progress in Earth andPlanetary Science

Open Access

Simulated decadal variations of surface
and subsurface phytoplankton in the upstream
Kuroshio Extension region
Tomoki Tozuka1,2* , Yoshikazu Sasai3 , Sayaka Yasunaka4 , Hideharu Sasaki2   and Masami Nonaka2   

Abstract
Using outputs from an ecosystem model embedded in an eddy-resolving ocean general circulation model that canrealistically simulate decadal modulations of the Kuroshio Extension (KE) between stable and unstable states, decadal
variations of phytoplankton concentration in the upstream KE region are investigated. During stable states of the KE,
surface phytoplankton concentrations are anomalously suppressed to the south of the KE front, while those to thenorth are anomalously enhanced. Although the surface phytoplankton concentration anomalies are prominent only
during winter to spring, significant subsurface anomalies centered around 60 m depth persist even in summer andautumn. Anomalies persist throughout the year in phytoplankton biomass integrated over the upper 100 m despite
the strong surface anomalies during the spring bloom season. An analysis of nitrogen concentration anomalies suggests that the vertical movement of the isopycnal surfaces, vertical mixing of nutrients, and meridional shifts in the KE
jet contribute to the anomalous phytoplankton biomass.
Keywords: Phytoplankton, Kuroshio Extension, Decadal variation, Ecosystem model, Eddy-resolving ocean generalcirculation model
1 Introduction
The Kuroshio Extension (KE) in the northwestern Pacificundergoes large decadal variations and is known to shift
between stable and unstable states (Qiu and Chen 2005).
It has been shown that wind stress curl anomalies in thecentral North Pacific that excite Rossby waves (Latif andBarnett 1996; Qiu and Chen 2005; Taguchi et al. 2007;
Nonaka et al. 2006; Qiu et al. 2014) are the root cause ofthis decadal variability. When these Rossby waves reachthe KE, they are trapped along the jet and lead to a stable state when downwelling Rossby waves impinge andan unstable state when upwelling Rossby waves arrive(Sasaki and Schneider 2011; Sasaki et al. 2013).

*Correspondence: tozuka@eps.s.u-tokyo.ac.jp
1
Department of Earth and Planetary Science, Graduate School of Science, The
University of Tokyo, Tokyo, Japan
Full list of author information is available at the end of the article

To understand how such decadal variations in the
KE affect marine ecosystem, which are also known toundergo large decadal variations (Noto and Yasuda 1999;
Chiba et al. 2013), it is necessary to understand decadal
variations in the lower trophic level marine ecosystem.
In this regard, Lin et al. (2014) investigated decadal variations in surface chlorophyll associated with decadal
modulations of the KE using satellite data. They revealed
that the surface chlorophyll south of 35°N becomesanomalously low when the KE is in the stable state, while
it increases anomalously during the KE unstable state. In
some extreme years, anomalies can reach up to 35% compared to the climatological value during the spring bloom
season. They suggested that anomalous vertical nutrient
supply induced by westward propagating Rossby waves
related to the North Pacific Gyre Oscillation (Di Lorenzo et al. 2008) and Pacific Decadal Oscillation (Mantua et al. 1997; Newman et al. 2016) contributes to thedecadal surface chlorophyll variations. Lin et al. (2020)© The Author(s) 2022. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

Tozuka et al. Progress in Earth and Planetary Science(2022) 9:70

later successfully reproduced the above decadal variations with their coupled physical–biological model andconducted a mixed layer nitrogen budget analysis. It isshown that the vertical advection mainly caused by nutricline heaving plays an important role. However, theiranalyses mostly focused on the impacts of the decadal
variations of the KE on the surface chlorophyll in theannual mean and the spring bloom. Since chlorophyll a(Chl-a) often retains its maximum in the subsurface layer
in this region (Fujiki et al. 2020; Yasunaka et al. 2022), it
will be illuminating to also examine subsurface variability
for more comprehensive understanding.
This paper is organized as follows. A brief descriptionof the model and datasets used in this study is provided
in the next section. In Sect. 3, simulated seasonal variations in phytoplankton are described first, which is aprerequisite to understand decadal variations. Then, simulated surface and subsurface phytoplankton concentration anomalies during the KE stable state are presented
and their possible mechanisms are discussed. Conclusions and discussions are given in the final section.

2 Methods
2.1 Model

We have used outputs from a simple nitrogen-based
NPZD (nitrogen, phytoplankton, zooplankton, and detritus) pelagic model (Oschlies 2000) embedded in an eddyresolving ocean general circulation model called the
Ocean general circulation model For the Earth Simulator(OFES2) (Masumoto et al. 2004; Sasaki et al. 2020; Sasai
et al. 2006, 2007a, 2007b, 2010, 2022) configured for the
North Pacific. The horizontal resolution is 0.1° × 0.1°, andthere are 105 levels in the vertical. The physical fields are
spun up with the monthly climatology of the JRA55-do(Tsujino et al. 2018) for 50 years and are forced by the
3-hourly data of the JRA55-do from 1958 to 1979. The
last day of 1979 is used in initial physical conditions forthe coupled physical-biological model, and then, it isforced from 1980 to 2019. In this study, we have used
40-year outputs of the coupled physical-biological modelfrom 1980 to 2019. For more details about this ecosystemmodel embedded in the eddy-resolving model including
its more thorough validations in the North Pacific, readers are referred to a paper by Sasai et al. (2022).
2.2 Data

To validate the model simulation, we use daily satellite altimetry data for sea surface height (SSH), which
was downloaded from the Copernicus Climate DataStore. It covers the period from 1993 to 2019, and its
horizontal resolution is 0.25° × 0.25°. Also, we use ahigh-resolution data-assimilation product known as thefour-dimensional variational ocean re-analysis for the

Page 2 of 14

Western North Pacific over 30 years (FORA-WNP30)(Usui et al. 2017). It is based on the MeteorologicalResearch Institute (MRI) Community Ocean Model version 2.4 (Tsujino et al. 2006), and the four-dimensionalvariational analysis scheme version of the MRI Multivariate Ocean Variational Estimation system (Usui et al.
2015) is used. The horizontal resolution is 0.1° × 0.1°, andwe use the product from 1982 to 2015. Subsurface Chl-adata are obtained from the World Ocean Database 2018
(WOD2018; Boyer et al. 2018) and the Global Ocean
Data Analysis Project version 2.2019 Release (GLODAPv2.2019; Olsen et al. 2019). Monthly mean Chl-a
concentrations with a horizontal resolution of 1° × 1°
are calculated at 10 different vertical levels in the upper
200 m (see Yasunaka et al. (2022) for more details). The
World Ocean Atlas 2018 (WOA2018; Garcia et al. 2018)is used for nitrate concentration, while Mixed Layer data
set of Argo, Grid Point Value (MILA GPV; Hosoda et al.
2010) based on the Argo profiles is used for mixed layer
depth (MLD). The horizontal resolution is 1° × 1° for bothdatasets. The data-derived euphotic layer is estimatedfrom the satellite-based photosynthetically available
radiation (PAR) and 1% light level depth from the GlobColour project website (GlobColour_R2018; Frouin et al.
2003; Morel et al. 2007) (see Yasunaka et al. (2022) for
more details).

3 Results
The study area of this paper is schematically shown in
Fig. 1. The Kuroshio that flows along the southern coastof Japan separates from the coast and flows into the interior ocean of the North Pacific as the KE. Horizontally,
nutrients are rich on the northern side of the KE andpoor on the southern side, while in the vertical direction, nutrient-rich water is found at depth, which may
be brought to the surface mixed layer by vertical mixing.
The shortwave radiation from the sun penetrates into theocean, and it is sufficient to support photosynthesis ofphytoplankton in the euphotic layer.
3.1 Seasonal variation

Prior to investigating decadal variations in phytoplankton concentrations, it is necessary to understand its
seasonal characteristics (Tozuka et al. 2017). Figure 2ashows the climatology of phytoplankton concentration
simulated in the upstream KE region averaged between
142°E and 153°E. As expected, the phytoplankton concentration is much higher to the north of the KE jet
throughout the year and reaches its maximum in thespring bloom season in April. Although the maximum
is located at the surface from December to April, themaximum appears in the subsurface centered around
60 m depth from June to October (Fig. 2a). This isTozuka et al. Progress in Earth and Planetary Science(2022) 9:70

Page 3 of 14

Fig. 1 Schematic diagram of the study area showing three-dimensional structure of nitrogen concentration as well as the major circulation pattern
and relevant processes/layers. For nitrogen concentration (shading), nutricline (thick solid line), mixed layer depth (thick dashed line), and the baseof euphotic layer (thin solid line), the climatological mean of October from the OFES are usedbecause nitrogen is depleted after the spring bloom inthe surface mixed layer, defined by the depth at which
the potential density increases by 0.03 kg ­m−3 from its
surface value, but sufficient nitrogen is found in thesubsurface even in summer and autumn (nutricline
defined by the depth of 1 mmol N ­m−3 isoline is deeper
than the mixed layer) and the euphotic layer (i.e., thedepth of the 1 W ­m−2 irradiance; Pastor et al. 2013) isdeep enough. Such a subsurface maximum is seen in
the observed Chl-a concentration in summer/autumn,
although the maximum is located at shallower depths
to the north of the KE jet and a deeper level to the southof the KE jet (Fig. 2b; Itoh et al. 2015; Yasunaka et al.
2022). Also, the vertical migration of the simulated
nutricline is in agreement with the observation (Fig. 2);
more specifically, the nutricline outcrops in February,
but it deepens from April to October and starts to outcrop again to the north of the KE jet in December. Furthermore, the euphotic layer depth in both the modeland observational estimates decreases with latitude
and reaches its seasonal maximum in summer withstronger shortwave radiation, although the depth tends
to be deeper in the model, especially from December toApril. Since subsurface observations of Chl-a concentrations are still limited, Fig. 2b is noisier than Fig. 2aand data are missing in some months. For this reason,
the observational subsurface Chl-a data is used only tovalidate the simulated seasonal variations.

3.2 Decadal variations

To define the stable and unstable states of the KE, wehave first applied an empirical orthogonal function (EOF)
analysis to SSH anomalies in 140°E-170°E, 30°N-40°N following Sasaki et al. (2013). The first EOF mode successfully captures the well-known decadal modulations of the
KE between stable and unstable states; i.e., SSH anomalies with opposite signs are found to the south and northof the KE jet (Fig. 3a). This leading mode explains 14.8%of the total variance. Based on the criterion of Northet al. (1982), it is significantly separated from the second
EOF mode, whose variance contribution is 7.4%. The spatial distribution and variance contribution are consistent
with those obtained from the same analysis with the satellite observation (Fig. 3b) and FORA (Fig. 3c).
We have then computed the annual mean value of thefirst principal component (PC) (Fig. 3d). The PC obtained
from the model is well correlated with those from satellite observations and FORA with their correlation coefficients of 0.70 and 0.62, respectively. The relatively highcorrelations not only suggest that the model has highskill in simulating decadal variability in the KE, but also
indicate that decadal fluctuations in SSH over the KE
have a deterministic component with respect to atmospheric forcing (Qiu et al. 2014). The first PC without
any temporal filter undergoes large decadal variations
with its value above 1 standard deviation in eight years(1990, 2002, 2003, 2004, 2005, 2011, 2015, and 2018) andbelow -1 standard deviations in eight years (1983, 1985,
1986, 1996, 1997, 1998, 2007, and 2008). These years areTozuka et al. Progress in Earth and Planetary Science(2022) 9:70

Page 4 of 14

Fig. 2 Monthly climatology of phytoplankton concentration in the upstream Kuroshio Extension (KE) averaged between 142°E and 153°E for a
OFES2 and b observation. The thick solid lines represent nutricline, the thick dashed lines represent the base of the mixed layer, and the thin solidlines represent the base of euphotic layer

Tozuka et al. Progress in Earth and Planetary Science(2022) 9:70

Page 5 of 14

Fig. 3 First empirical orthogonal function (EOF) mode of sea surface height (SSH) anomalies in 140°E-170°E, 30°N-40°N obtained from a OFES2,
b satellite observation, and c FORA. d Normalized annual mean of the first principal component (PC) of SSH anomalies in 140°E-170°E, 30°N-40°N
obtained from the OFES2 (black line), satellite observation (light blue line), and FORA (green line). The horizontal black dotted lines represent ± 1
standard deviation, and the red (blue) dots represent the stable (unstable) states

defined, respectively, as the stable and unstable states.
Since the unstable state is close to a mirror image of thestable state in general, we focus on the stable state in thisstudy.
To investigate decadal variations in phytoplankton
concentrations and their controlling mechanisms, composites of their surface anomalies during stable states
are constructed (Fig. 4). Throughout the year, negative
phytoplankton concentration anomalies are found to thesouth of the KE jet, while positive anomalies are foundto the north, but these anomalies are more conspicuous
during winter and spring. The results are consistent withthose of Lin et al. (2014, 2020).
To examine whether such anomalies are also found in
the subsurface, composites of phytoplankton concentration anomalies along 145°E, where the strongest surface
negative anomalies are found to the south of KE in winter, are prepared (Fig. 5). It is clear that significant phytoplankton concentration anomalies are not confined tothe surface, but extend down to around 100 m depth inTozuka et al. Progress in Earth and Planetary Science(2022) 9:70

Page 6 of 14

Fig. 4 Composites of phytoplankton concentration anomalies at the surface during the stable states. The contour intervals are 0.02 mmol N ­m−3.
Anomalies significant at the 95% confidence level by a two-tailed t test are shaded

winter to the south of the KE jet located around 36°N. It
is intriguing that significant phytoplankton concentration anomalies remain in the subsurface layer even in
August and October despite the disappearance of significant anomalies near the surface during these months.
As a result, significant anomalies exist throughout theyear in the phytoplankton biomass integrated over theupper 100 m even though the strong surface anomalies
are found during the spring bloom season (Fig. 6). Thus,
the decadal variations of the KE may have larger impacts
on marine ecosystems and carbon dioxide uptake thanpreviously thought, as the impacts are not limited to thespring bloom season.
To understand the mechanism of this interestingseasonal evolution of phytoplankton concentration
anomalies, we have constructed composites of nitrogenconcentration anomalies along 145°E. Figure 7a reveals
that significant negative concentration anomalies extend
to the surface south of the KE jet throughout the year.
This may explain why phytoplankton concentrations
decreased anomalously to the south of the KE jet during
stable states. Since the climatological mean nitrogen concentration is already too low to sustain the phytoplankton
activity near the surface in summer and autumn, nitrogen concentration anomalies do not lead to significant
phytoplankton concentration anomalies during these two
seasons.

To check the relative importance of anomalous verticalmovement of the isopycnal surfaces and vertical mixingof nitrogen, composites of nitrogen concentration anomalies in a density coordinate are examined (Fig. 7b). If thevertical process were important as suggested by Lin et al.
(2020) and nitrogen concentration anomalies were solely
due to vertical shifts in the nutricline, we would expect
significant anomalies to disappear in the density coordinate. However, we find significant but weaker nitrogenconcentration anomalies throughout the year. This suggests that adiabatic vertical movements make a substantial contribution, but diapycnal mixing also plays a role in
the negative nitrogen and phytoplankton concentration
anomalies to the south of the KE jet. We note that phytoplankton concentration anomalies may also contribute tonitrogen concentration anomalies through consumption,
but negative phytoplankton concentration anomalies are
expected to lead to less consumption of nitrogen andcannot explain the negative anomalies in the nitrogenconcentration.
We suggest two possible mechanisms explaining how
vertical movements and mixing contribute to the negative nitrogen concentration anomalies to the south of
the KE jet during the KE stable state. First, positive
anomalies in the vertical gradient of nitrogen concentration (i.e., negative nitrogen concentration anomalies
are smaller near the surface) (Fig. 7a) lead to a reduced
supply by the vertical movement and mixing. Second,Tozuka et al. Progress in Earth and Planetary Science(2022) 9:70

Page 7 of 14

Fig. 5 Composites of phytoplankton concentration anomalies (in mmol N ­m−3) along 145°E in the upstream KE during the stable states in a thez-coordinate and b the density coordinate. The black contour intervals are 0.01 mmol N ­m−3. Also shown in (a) with red contours are composites ofpotential density. Anomalies significant at the 95% confidence level by a two-tailed t test are shadedTozuka et al. Progress in Earth and Planetary Science(2022) 9:70

Page 8 of 14

Fig. 6 Composites of phytoplankton concentration anomalies integrated over the upper 100 m. The contour intervals are 2 mmol N ­m−2.
Anomalies significant at the 95% confidence level by a two-tailed t test are shaded

downwelling anomalies to the south of the KE jet also
lead to less upward transport of nitrogen-rich water fromthe subsurface.
Furthermore, when the KE jet latitude coordinate isused, the nitrogen concentration anomalies become
smaller in amplitude (Fig. 8). This suggests that at least
some portion of the negative nitrogen concentration
anomalies to the south of the KE jet may be explained by
the northward migration of the KE jet during the KE stable state.
Although we have focused on the effect of nutrient concentration anomalies, we note that other effects such as
MLD variations, light availability, and temperature variations may also influence phytoplankton concentrations.
However, it turns out that these effects act as a negative
feedback. First, MLD is known to undergo decadal variations associated with the KE state (Sugimoto and Kako
2016) and wintertime MLD anomalies may contributeto nitrogen concentration anomalies because deepening of the mixed layer is associated with entrainment of
nutrient-rich subsurface water. However, their contribution seems to not be dominant, as the MLD becomesanomalously deep to the south of the KE jet during the
KE stable state (figure not shown), which should favor
positive nitrogen concentration anomalies. Also, becauseof the negative phytoplankton concentration anomalies
to the south of the KE jet during stable states, penetrative

shortwave radiation anomalies are positive and shortwave radiation can reach deeper (figure not shown). This
would favor an increase in phytoplankton concentrationsat depth, but the fact that phytoplankton concentration
anomalies are negative suggests that the penetrative
shortwave radiation anomalies act as a negative feedback, but their impact is smaller compared with the effectof anomalous nutrient concentrations. Furthermore,because of the anomalously deep thermocline, positive
temperature anomalies are found to the south of the KE
jet (figure not shown). Again, this favors enhanced concentrations of phytoplankton at depth and opposes thenegative phytoplankton concentration anomalies associated with reduced nutrient concentration. For moredetails on their effects, readers may also refer to a paper
by Sasai et al. (2022).
Finally, although our analyses have mostly focused on
the 145°E section, we present composites of nitrogen andphytoplankton concentration anomalies in February andAugust that compare cross sections at 142.5°E, 145°E,
147.5°E, and 150°E (Figs. 9 and 10). Negative nitrogenconcentration anomalies are found to the south of the KE
jet axis from 142.5°E to 150°E with the strongest anomalies at 145°E in both February and August (Fig. 9). In
February, strong negative phytoplankton concentration
anomalies extending from the surface to about 150 m
depth at 142.5°E are further strengthened at 145°E, butTozuka et al. Progress in Earth and Planetary Science(2022) 9:70

Fig. 7 As in Fig. 5, but for nitrogen concentration anomalies (in mmol N ­m−3). The black contour intervals are 1 mmol N ­m−3

Page 9 of 14

Tozuka et al. Progress in Earth and Planetary Science(2022) 9:70

Page 10 of 14

Fig. 8 As in Fig. 7a, but in the KE jet latitude coordinate

weaken gradually from 147.5°E to 150°E (Fig. 10). Since
the euphotic layer depth is about 80 m, negative phytoplankton concentration anomalies below this levelmay be caused by intense wintertime mixing within themixed layer that extends to about 150 m. In August, thepeak negative phytoplankton concentration anomalies found between 24.0 and 24.5 potential density surfaces strengthen from 142.5°E to 145°E, but disappear at
147.5°E and reappear at 150°E (Fig. 10). The disappearance at 147.5°E suggests that the negative phytoplankton
concentration anomalies are not just advected zonally by
the KE.

4 Conclusions and discussions
In this study, we have investigated decadal variations ofphytoplankton concentrations in the upstream KE region
using outputs from a simple nitrogen-based NPZD
pelagic model embedded in an eddy-resolving ocean
general circulation model, which can represent the KE
jet and associated frontal structures. Observed seasonalvariability in the subsurface structure of Chl-a concentrations and vertical migration of nutricline is also well simulated in the model. In agreement with past studies basedon satellite observations and model simulations (Lin et al.
2014, 2020), the surface phytoplankton concentrations
are anomalously suppressed to the south of the KE front,
while those to the north are anomalously enhanced during winter and spring of the stable state of the KE. However, we have shown for the first time that significant
subsurface anomalies centered around 60 m depth persist
even in summer and autumn. Although the surface phytoplankton concentration anomalies are prominent only
during winter to spring, significant subsurface anomalies centered around 60 m depth persist even in summer and autumn. As a result, significant anomalies exist
throughout the year in the phytoplankton biomass integrated over the upper 100 m despite the strong surface
anomalies during the spring bloom season. A comparison of nitrogen concentration anomalies in two different
coordinate systems suggests that both anomalous verticalmovement and mixing contribute to decadal variationsin the nitrogen and phytoplankton concentrations. Also,
negative nitrogen concentration anomalies to the southof the KE jet associated with stable states of the KE may
partially be explained by the northward migration of theTozuka et al. Progress in Earth and Planetary Science(2022) 9:70

Page 11 of 14

Fig. 9 Composites of nitrogen concentration anomalies (in mmol N ­m−3) along (left to right) 142.5°E, 145°E, 147.5°E, and 150°E during the stablestates in (upper row) February and (bottom row) August. The contour intervals are 1 mmol N ­m−3. Also shown with red contours are composites ofpotential density. Anomalies significant at the 95% confidence level by a two-tailed t test are shaded

Fig. 10 As in Fig. 9, but for phytoplankton concentration anomalies (in mmol N ­m−3). The contour intervals are 0.01 mmol N ­m−3

Tozuka et al. Progress in Earth and Planetary Science(2022) 9:70

KE jet, which results in a northward expansion of nutrient poor waters from the subtropical gyre.
Although three possible mechanisms have been proposed, it is difficult to quantify their relative contributions. For example, anomalies in vertical velocity and thelatitude of the KE jet axis are negatively correlated witheach other. Thus, if we were to quantify their respective
contributions with a multiple regression analysis, collinearity becomes an issue. In other words, a northward shift
in the KE jet not only leads to a northward extension ofnutrient poor waters from the subtropical gyre, but also
results in a northward extension of the subtropical gyreassociated with downwelling. Another reason is that we
do not have the vertical diffusion coefficient as a model
output and it is not possible to quantify vertical mixing of
nutrients. Conducting an online nitrogen budget analysis(i.e., saving all terms in the nitrogen budget equation ateach time step and accumulating them, as was done by
Nakazato et al. (2021) for the heat budget equation) may
provide a solution.
The KE region is known for its enhanced turbulent
diapycnal mixing, and it is attributed to the strong confluence of the Kuroshio and Oyashio currents, and anageostrophic secondary circulation (D’Asaro et al. 2011;
Nagai et al. 2012). The importance of diapycnal mixing
suggested by this study is in agreement with Kaneko et al.
(2013, 2021), who provided observational evidence thatvertical mixing plays an important role in nutrient supply
to the euphotic layer.
Although intrinsic variability greatly reduces potential predictability of the KE (Nonaka et al. 2016), westward propagating Rossby waves have been suggested toprovide a source of memory for decadal prediction (Qiu
et al. 2014; Joh et al. 2022). The relatively high correlations in simulated decadal SSH variations with satellite
observations and an ocean re-analysis product (Fig. 3)
suggest that decadal variations of phytoplankton concentration in the KE that are closely related to the KE state
are potentially predictable.
The region south of the KE is also known as the formation area of the North Pacific subtropical mode water(Masuzawa 1969). Since it is circulated in the western
North Pacific (Suga and Hanawa 1995; Suga et al. 2008),
decadal variability of phytoplankton concentrations in
the KE may have profound impacts on the water properties in vast areas of the western North Pacific (e.g., Oka
et al. 2015, 2018) and the uptake of carbon dioxide. Oka
et al. (2015) showed some hint of decadal downstream
signatures in biogeochemical properties along isopycnal
surfaces that correspond to the subtropical mode water,
but denser observations and more sophisticated analyses such as water mass transformation diagnostics (e.g.,
Iudicone et al. 2011) are required to more clearly detect

Page 12 of 14

such signatures. Thus, future studies examining the consequence of the decadal variability revealed in this study
may further elucidate decadal variability not only in themarine ecosystem in the North Pacific, but also in theoceanic storage of anthropogenic carbon dioxide (Bateset al. 2002).
Mesoscale eddy activity in the KE region is also known
to undergo decadal variations with enhanced eddy activity during unstable states and reduced activity during
stable states (Qiu and Chen, 2005). Cyclonic eddies tothe south of the KE detached from the KE jet have beenshown to play an important role in supplying nutrientsand enhancing near-surface productivity (Sasai et al.
2010, 2019; Kouketsu et al. 2016; Uchiyama et al. 2017;
Honda et al. 2018; Lin et al. 2020). It will be illuminating to investigate how decadal variations in eddy activity
associated with the KE state contribute to those in phytoplankton concentrations in the KE region. A study in thisdirection is underway.
Although the present results are mostly based on modeling, this study suggests the importance of monitoring
subsurface phytoplankton variabilities in addition to satellite observations that can measure surface chlorophylla concentration. This may be accomplished by expanded
deployment of Biogeochemical Argo floats (Bittig et al.
2019).
Abbreviations
Chl-a: Chlorophyll a; EOF: Empirical orthogonal function; FORA-WNP30: Fourdimensional variational ocean re-analysis for the Western North Pacific over
30 years; KE: Kuroshio Extension; GLODAPv2.2019: Global ocean data analysis
project version 2.2019 release; MILA GPV: Mixed layer data set of Argo, grid
point value; MLD: Mixed layer depth; MRI: Meteorological Research Institute;
MRI.COM: Meteorological Research Institute community ocean model; NPZD:
Nitrogen, phytoplankton, zooplankton, and detritus; OFES: Ocean generalcirculation model for the Earth Simulator; PC: Principle component; SSH: Sea
surface height; WOA2018: World Ocean Atlas 2018; WOD2018: World Ocean
Database 2018.
Acknowledgements
We thank two anonymous reviewers for their constructive comments, whichwere very helpful in improving our manuscript. The OFES simulation was
conducted on the Earth Simulator under the support of JAMSTEC.
Author contributions
TT and YS conceived and designed the study. YS and HS carried out the model
simulations. TT and YS analyzed the model outputs, and SY analyzed theobservational data. TT wrote the first draft of this manuscript. All authors readand approved the final manuscript.
Funding
This study is supported by the Japan Society for Promotion of Science throughGrant-in-Aid for Scientific Research on Innovative Area (MEXT KAKENHI Grant
number 19H05701).
Availability of data and material
The outputs from the ecosystem model are available from Yoshikazu Sasai on
reasonable request. The satellite altimetry data were obtained from https://​cds.​clima​te.​coper​nicus.​eu/​cdsapp#​!/​datas​et/​satel​lite-​sea-​level-​global, climatological mixed layer depth data were obtained from MILA GPV (http://​www.​jamst​ec.​go.​jp/​ARGO/​argo_​web/​argo/?​page_​id=​223an​dlang=​en), and nitrateTozuka et al. Progress in Earth and Planetary Science(2022) 9:70

concentration data are obtained from the WOA18 (https://​www.​nodc.​noaa.​
gov/​OC5/​woa18/). Subsurface chlorophyll a (Chl-a) data were obtained fromthe WOD2018 (https://​www.​nodc.​noaa.​gov/​OC5/​WOD/​pr_​wod.​html) and the
GLODAPv2.2019 (https://​www.​glodap.​info/). Satellite-based PAR and 1% lightlevel depth data were obtained from the GlobColour project website (http://​hermes.​acri.​fr/​index.​php). This study used the FORA-WNP30, which was produced by Japan Agency for Marine-Science and Technology (JAMSTEC) andMeteorological Research Institute of Japan Meteorological Agency (JMA/MRI)and is available from http://​synth​esis.​jamst​ec.​go.​jp/​FORA/e/​index.​html.

Declarations
Competing interests
The authors declare that they have no competing interest.
Author details
1
Department of Earth and Planetary Science, Graduate School of Science,
The University of Tokyo, Tokyo, Japan. 2 Application Laboratory (APL), ResearchInstitute for Value‐Added‐Information Generation (VAiG), Japan Agency
for Marine‐Earth Science and Technology (JAMSTEC), Yokohama, Japan.
3
Earth Surface System Research Center, Research Institute for Global Change,
JAMSTEC, Yokosuka, Japan. 4 Graduate School of Science, Tohoku University,
Sendai, Japan.
Received: 29 June 2022 Accepted: 5 December 2022

References
Bates NR, Pequignet AC, Johnson RJ, Gruber N (2002) A variable sink foratmospheric ­CO2 in subtropical mode water of the North Atlantic Ocean.
Nature 420:489–493. https://​doi.​org/​10.​1038/​natur​e01253
Bittig HC, Mauer TL, Plant JN, Schmechtig C, Wong APS, Claustre H, Trull TW,
Udaya Bhaskar TVS, Boss E, Dall’Olmo G, Organelli E, Poteau A, Johnson
KS, Hanstein C, Leymarie E, Le Reste S, Riser SC, Rupan AR, Taillandier V,
Thierry V, Xing X (2019) A BGC-Argo guide: planning, deployment, data
handling and usage. Front Mar Sci 6:502. https://​doi.​org/​10.​3389/​fmars.​2019.​00502
Boyer TP, Baranova OK, Coleman C, Garcia HE, Grodsky A, Locarnini RA, Mishonov AV, Paver CR, Reagan JR, Seidov D, Smolyar IV, Weathers K, Zweng
MM (2018) World Ocean Database 2018, Mishonov AV, Technical Ed.;NOAA Atlas NESDIS 87, 2018
Chiba S, Di Lorenzo E, Davis A, Keister JE, Taguchi B, Sasai Y, Sugisaki H (2013)
Large-scale climate control of zooplankton transport and biogeography
in the Kuroshio-Oyashio Extension region. Geophys Res Lett 40:5182–
5187. https://​doi.​org/​10.​1002/​grl.​50999
D’Asaro E, Lee C, Rainville L, Harcourt R, Thomas L (2011) Enhanced turbulenceand energy dissipation at ocean fronts. Science 332:318–322. https://​doi.​
org/​10.​1126/​scien​ce.​12015​15
Di Lorenzo E, Schneider N, Cobb KM, Franks PJS, Chhak K, Miller AJ, McWilliams
JC, Bograd SJ, Arango H, Curchitser E, Powell TM, Rivière P (2008) North
Pacific Gyre Oscillation links ocean climate and ecosystem change. Geophys Res Lett 35:L08607. https://​doi.​org/​10.​1029/​2007G​L0328​38
Frouin R, Franz B, Werdell P (2003) The SeaWiFS PAR product. Algorithm
updates for the fourth SeaWIFS data reprocessing, 22, 46–50
Fujiki T, Inoue R, Honda MC, Wakita M, Mino Y, Sukigara C, Abe O (2020) Timeseries observations of photosynthetic oxygen production in the subtropical western North Pacific by an underwater profiling buoy system. Limnol
Oceanogr 65:1072–1084. https://​doi.​org/​10.​1002/​lno.​11372
Garcia HE, Weathers K, Paver CR, Smolyar I, Boyer TP, Locarnini RA, Zweng MM,
Mishonov AV, Baranova OK, Seidov D, Reagan JR (2018) World Ocean
Atlas 2018, Volume 4: Dissolved Inorganic Nutrients (phosphate, nitrateand nitrate+nitrite, silicate), Mishonov A Technical Ed.; NOAA Atlas
NESDIS 84, 35pp
Honda MC, Sasai Y, Siswanto E, Kuwano-Yoshida A, Aiki H, Cronin MF (2018)
Impact of cyclonic eddies and typhoons on biogeochemistry in theoligotrophic ocean based on biogeochemical/physical/meteorologicaltime-series at station KEO. Prog Earth Planet Sci 5:42. https://​doi.​org/​10.​
1186/​s40645-​018-​0196-3

Page 13 of 14

Hosoda S, Ohira T, Sato K, Suga T (2010) Improved description of global mixedlayer depth using Argo profiling floats. J Oceanogr 66:773–787. https://​doi.​org/​10.​1007/​s10872-​010-​0063-3
Itoh S, Yasuda I, Saito H, Tsuda A, Komatsu K (2015) Mixed layer depth andchlorophyll a: Profiling float observations in the Kuroshio-Oyashio Extension region. J Mar Sys 151:1–14. https://​doi.​org/​10.​1016/j.​jmars​ys.​2015.​
06.​004
Iudicone D, Rogers KB, Stendardo I, Aumont O, Madec G, Bopp L, Mangoni O,
d’Alcala MR (2011) Water masses as a unifying framework for understanding the Southern Ocean carbon cycle. Biogeosci 8:1031–1052. https://​doi.​
org/​10.​5194/​bg-8-​1031-​2011
Joh Y, Delworth TL, Wittenberg AT, Cooke WF, Yang X, Zeng F, Jia L, Lu F,
Johnson N, Kapnick SB, Rosati A, Zhang L, McHugh C (2022) Seasonal-todecadal variability and prediction of the Kuroshio Extension in the GFDL
coupled ensemble reanalysis and forecasting system. J Clim 35:3515–
3535. https://​doi.​org/​10.​1175/​JCLI-D-​21-​0471.1
Kaneko H, Yasuda I, Komatsu K, Itoh S (2013) Observations of vertical turbulent
nitrate flux across the Kuroshio. Geophys Res Lett 40:3123–3127. https://​doi.​org/​10.​1002/​grl.​50613
Kaneko H, Yasuda I, Itoh S, Ito S (2021) Vertical turbulent nitrate flux fromdirect measurements in the western subarctic and subtropical gyresof the North Pacific. J Oceanogr 77:29–44. https://​doi.​org/​10.​1007/​
s10872-​020-​00576-0
Kouketsu S, Kaneko H, Okunishi T, Sasaoka K, Itoh S, Inoue R, Ueno H (2016)
Mesoscale eddy effects on temporal variability of surface chlorophyll a in
the Kuroshio Extension. J Oceanogr 72:439–451. https://​doi.​org/​10.​1007/​
s10872-​015-​0286-4
Latif M, Barnett TP (1996) Decadal climate variability over the North Pacific andNorth America: dynamics and predictability. J Clim 9:2407–2423. https://​doi.​org/​10.​1175/​1520-​0442(1996)​009%​3c2407:​DCVOTN%​3e2.0.​CO;2
Lin P, Chai F, Xue H, Xiu P (2014) Modulation of decadal oscillation on surface
chlorophyll in the Kuroshio Extension. J Geophys Res Oceans 119:187–
199. https://​doi.​org/​10.​1002/​2013J​C0093​59
Lin P, Ma J, Chai F, Xiu P, Liu H (2020) Decadal variability of nutrients andbiomass in the southern region of Kuroshio Extension. Prog Oceanogr
188:102441. https://​doi.​org/​10.​1016/j.​pocean.​2020.​102441
Mantua NJ, Hare SR, Zhang Y, Wallace JM, Francis RC (1997) A Pacific decadal
climate oscillation with impacts on salmon. Bull Amer Meteorol Soc
78:1069–1079. https://​doi.​org/​10.​1175/​1520-​0477(1997)​078%​3c1069:​
APICOW%​3e2.0.​CO;2
Masuzawa J (1969) Subtropical mode water. Deep-Sea Res 16:463–472.
https://​doi.​org/​10.​1016/​0011-​7471(69)​90034-5
Masumoto, Y, Sasaki H, Kagimoto T, Komori N, Ishida A, Sasai Y, Miyama T, Motoi
T, Mitsudera H, Takahashi K, Sakuma H, Yamagata T (2004) A fifty-year
eddy-resolving simulation of the world ocean—Preliminary outcomes of
OFES (OGCM for the Earth Simulator). J Earth Simulator 1:35–56
Morel A, Huot Y, Gentili B, Werdell PJ, Hooker SB, Franz BA (2007) Examining
the consistency of products derived from various ocean color sensors
in open ocean (Case 1) waters in the perspective of a multi-sensor
approach. Remote Sens Environ 111:69–88. https://​doi.​org/​10.​1016/j.​rse.​
2007.​03.​012
Nagai T, Tandon A, Yamazaki H, Doubell MJ, Gallager S (2012) Direct observations of microscale turbulence and thermohaline structure in the
Kuroshio Front. J Geophys Res 117:C08013. https://​doi.​org/​10.​1029/​2011J​
C0072​28
Nakazato M, Kido S, Tozuka T (2021) Mechanisms of asymmetry in sea surface
temperature anomalies associated with the Indian Ocean Dipole
revealed by closed heat budget. Sci Rep 11(1):22546. https://​doi.​org/​10.​1038/​s41598-​021-​01619-2
Newman M, Alexander MA, Ault TR, Cobb KM, Deser C, Di Lorenzo E, Mantua
NJ, Miller AJ, Minobe S, Nakamura H, Schneider N, Vimont DJ, Phillips AS,
Scott JD, Smith CA (2016) The Pacific Decadal Oscillation, revisited. J Clim
29:4399–4427. https://​doi.​org/​10.​1175/​JCLI-D-​15-​0508.1
Nonaka M, Nakamura H, Tanimoto Y, Kagimoto T, Sasaki H (2006) Decadal variability in the Kuroshio-Oyashio Extension simulated in an eddy-resolving
OGCM. J Clim 19:1970–1989. https://​doi.​org/​10.​1175/​JCLI3​793.1
Nonaka M, Sasai Y, Sasaki H, Taguchi B, Nakamura H (2016) How potentiallypredictable are midlatitude ocean currents? Sci Rep 6:20153. https://​doi.​
org/​10.​1038/​srep2​0153

Tozuka et al. Progress in Earth and Planetary Science(2022) 9:70

North GR, Bell TL, Cahalan RF, Moeng FJ (1982) Sampling errors in the estimation of empirical orthogonal functions. Mon Wea Rev 110:699–706.
https://​doi.​org/​10.​1175/​1520-​0493(1982)​110%​3c0699:​SEITEO%​3e2.0.​CO;2
Noto M, Yasuda I (1999) Population decline of the Japanese sardine, Sardinops
melanostictus, in relation to sea surface temperature in the Kuroshio
Extension. Canad J Fish Aquat Sci 56:973–983. https://​doi.​org/​10.​1139/​
f99-​028
Oka E, Qiu B, Takatani Y, Enyo K, Sasano D, Kosugi N, Ishii M, Nakano T, Suga T(2015) Decadal variability of subtropical mode water subduction and its
impact on biogeochemistry. J Oceanogr 71:389–400. https://​doi.​org/​10.​
1007/​s10872-​015-​0300-x
Oka E, Yamada K, Sasano D, Enyo K, Nakano T, Ishii M (2018) Remotely forced
decadal physical and biogeochemical variability of North Pacific subtropical mode water over the last 40 years. Geophys Res Lett 46:1555–1561.
https://​doi.​org/​10.​1029/​2018G​L0813​30
Olsen A, Lange N, Key RM, Tanhua T, Álvarez M, Becker S, Bittig HC, Carter
BR, Cotrim da Cunha L, Feely RA, van Heuven S, Hoppema M, Ishii M,
Jeansson E, Jones SD, Jutterström S, Karlsen MK, Kozyr A, Lauvset SK, Lo
Monaco C, Murata A, Pérez FF, Pfeil B, Schirnick C, Steinfeldt R, Suzuki T,
Telszewski M, Tilbrook B, Velo A, Wanninkhof R (2019) GLODAPv2.2019—
An update of GLODAPv2. Earth Syst Sci Data 11:1437–1461. https://​doi.​
org/​10.​5194/​essd-​11-​1437-​2019
Oschlies A, Koeve W, Garçon V (2000) An eddy-permitting coupled physicalbiological model of the North Atlantic: 2. Ecosystem dynamics and comparison with satellite and JGOFS local station data. Global Biogeochem
Cycle 14:499–523. https://​doi.​org/​10.​1029/​1999G​B9000​80
Pastor MV, Palter JB, Pelegrí JL, Dunne J (2013) Physical drivers of interannual
chlorophyll variability in the eastern subtropical North Atlantic. J Geophys
Res Oceans 118:3871–3886. https://​doi.​org/​10.​1002/​jgrc.​20254
Qiu B, Chen S (2005) Variability of the Kuroshio Extension jet, recirculation
gyre, and mesoscale eddies on decadal time scales. J Phys Oceanogr
35:2090–2103. https://​doi.​org/​10.​1175/​JPO28​07.1
Qiu B, Chen S, Schneider N, Taguchi B (2014) A coupled decadal predictionof the dynamic state of the Kuroshio Extension system. J Clim 27:1751–
1764. https://​doi.​org/​10.​1175/​JCLI-D-​13-​00318.1
Sasai Y, Ishida A, Sasaki H, Kawahara S, Uehara H, Yamanaka Y (2006) A
global eddy-resolving coupled physical and biological model: Physicalinfluences on a marine ecosystem in the North Pacific. SIMULATION
82:467–474. https://​doi.​org/​10.​1177/​00375​49706​068943
Sasai Y, Sasaoka K, Sasaki H, Ishida A (2007a) Seasonal and intraseasonal variability of chlorophyll-a in the North Pacific: Model and satellite data. J
Earth Simulator 8:3–11
Sasai Y, Sasaki H, Sasaoka K, Ishida A, Yamanaka Y (2007b) Marine ecosystemsimulation in the eastern tropical Pacific with a global eddy resolving
coupled physical-biological model. Geophys Res Lett 34:L23601. https://​doi.​org/​10.​1029/​2007G​L0315​07
Sasai Y, Richards KJ, Ishida A, Sasaki H (2010) Effects of cyclonic mesoscaleeddies on the marine ecosystem in the Kuroshio Extension region
using an eddy-resolving coupled physical-biological model. Ocean Dyn
60:693–704. https://​doi.​org/​10.​1007/​s10236-​010-​0264-8
Sasai Y, Smith SL, Siswanto E, Sasaki H, Nonaka M (2022) Physiological flexibilityof phytoplankton impacts modeled biomass and primary production
across the North Pacific Ocean. Biogeosci 19:4865–4882. https://​doi.​org/​
10.​5194/​bg-​19-​4865-​2022
Sasai Y, Honda MC, Siswanto E, Kato S, Uehara K, Sasaki H, Nonaka M (2019)
Impact of ocean physics on marine ecosystems in the Kuroshio and Kuroshio Extension regions: A high-resolution coupled physical-biologicalmodel study. In: Nagai T, Saito H, Suzuki K, Takahashi M. (eds) Kuroshio
Current: Physical, Biogeochemical and Ecosystem Dynamics. AGU Geophysical Monograph Series. AGU-Wiley, pp 175–188. https://​doi.​org/​10.​
1002/​97811​19428​428.​ch11
Sasaki YN, Schneider N (2011) Decadal shifts of the Kuroshio Extension jet:application of thin-jet theory. J Phys Oceanogr 41:979–993. https://​doi.​
org/​10.​1175/​2011J​PO4550.1
Sasaki YN, Minobe S, Schneider N (2013) Decadal response of the Kuroshio
Extension jet to Rossby waves: observation and thin-jet theory. J Phys
Oceanogr 43:442–456. https://​doi.​org/​10.​1175/​JPO-D-​12-​096.1
Sasaki H, Kida S, Furue R, Aiki H, Komori N, Masumoto Y, Miyama T, Nonaka M,
Sasai Y, Taguchi B (2020) A global eddying hindcast ocean simulation with
OFES2. Geosci Model Develop 13:3319–3336. https://​doi.​org/​10.​5194/​
gmd-​13-​3319-​2020

Page 14 of 14

Suga T, Hanawa K (1995) The subtropical mode water circulation in the North
Pacific. J Phys Oceanogr 25:958–970. https://​doi.​org/​10.​1175/​1520-​
0485(1995)​025%​3c0958:​TSMWCI%​3e2.0.​CO;2
Suga T, Aoki Y, Saito H, Hanawa K (2008) Ventilation of the North Pacificsubtropical pycnocline and mode water formation. Prog Oceanogr
77:285–297. https://​doi.​org/​10.​1016/j.​pocean.​2006.​12.​005
Sugimoto S, Kako S (2016) Decadal variation in winter mixed layer depth southof the Kuroshio Extension and its influence on winter mixed layer temperature. J Clim 29:1237–1252. https://​doi.​org/​10.​1175/​JCLI-D-​15-​0206.1
Taguchi B, Xie SP, Schneider N, Nonaka M, Sasaki H, Sasai Y (2007) Decadal
variability of the Kuroshio Extension: observations and an eddy-resolving
model hindcast. J Clim 20:2357–2377. https://​doi.​org/​10.​1175/​JCLI4​142.1
Tozuka T, Cronin MF, Tomita H (2017) Surface frontogenesis by surface heat
fluxes in the upstream Kuroshio Extension region. Sci Rep 7:10258.
https://​doi.​org/​10.​1038/​s41598-​017-​10268-3
Tsujino H, Usui N, Nakano H (2006) Dynamics of Kuroshio path variations in ahigh-resolution GCM. J Geophys Res Oceans 111:C11001. https://​doi.​org/​
10.​1029/​2005J​C0031​18
Tsujino H, Urakawa S, Nakano H, Small RJ, Kim WM, Yeager SG, Danabasoglu
G, Suzuki T, Bamber JL, Bentsen M, Böning CW, Bozec A, Chassignet EP,
Curchitser E, Dias FB, Durack PJ, Griffies SM, Harada Y, Ilicak M, Josey SA,
Kobayashi C, Kobayashi S, Komuro Y, Large WG, Le Sommer J, Marsland
MJ, Masina S, Scheinert M, Tomita H, Valdivieso M, Yamazaki D (2018) JRA55 based surface dataset for driving ocean-sea-ice models (JRA55-do).
Ocean Model 130:79–139. https://​doi.​org/​10.​1016/j.​ocemod.​2018.​07.​002
Uchiyama Y, Suzue Y, Yamazaki H (2017) Eddy-driven nutrient transport andassociated upper-ocean primary production along the Kuroshio. J Geophys Res Oceans 122:5046–5062. https://​doi.​org/​10.​1002/​2017J​C0128​47
Usui N, Fujii Y, Sakamoto K, Kamachi M (2015) Development of a fourdimensional variational assimilation system for coastal data assimilation
around Japan. Mon Wea Rev 143:3874–3892. https://​doi.​org/​10.​1175/​
MWR-D-​14-​00326.1
Usui N, Wakamatsu T, Tanaka Y, Hirose N, Toyoda T, Nishikawa S, Fujii Y, Takatsuki
Y, Igarashi H, Nishikawa H, Ishikawa Y, Kuragano T, Kamachi M (2017) Fourdimensional variational ocean reanalysis: a 30-year high-resolution dataset in the western North Pacific (FORA-WNP30). J Oceanogr 73:205–233.
https://​doi.​org/​10.​1007/​s10872-​016-​0398-5
Yasunaka S, Ono T, Sasaoka K, Sato K (2022) Global distribution and variabilityof subsurface chlorophyll a concentration. Ocean Sci 18:255–268. https://​doi.​org/​10.​5194/​os-​18-​255-​2022

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.