Progress in Earth andPlanetary Science

Nagano et al.
Progress in Earth and Planetary Science(2022) 9:57
https://doi.org/10.1186/s40645-022-00513-3

Open Access

RESEARCH ARTICLE

Spatiotemporal vertical velocity variationin the western tropical Pacific and its relation
to decadal ocean variability
Akira Nagano1* , Takuya Hasegawa2,3 and Masahide Wakita4

Abstract: Prior to the cold phase of the Pacific decadal oscillation (PDO), sea surface water in the western tropicalNorth Pacific is heated. To evaluate the impact of the subsurface water on the upper-ocean (< 100 m depth) temperature variation, we examined the vertical velocity variation. We applied altimetry-based gravest empirical mode (AGEM)
method to Argo hydrographic profile (2001–2017) and altimetric sea surface height (1993–2017) data in the western
tropical Pacific and obtained a time-varying three-dimensional (3D) potential density field in the layer from 100 to1000 m depth. From the AGEM-derived potential density, we derived a 3D field of vertical velocity below 100 m depth
by the P-vector inverse method, which is based on potential density and vorticity conservations, and examined thespatiotemporal characteristics, focusing on the seasonal to quasi-decadal (QD) variations. In the mean state, substantial upward current (> 10−6 m s−1) is observed in the region of 10◦–17◦N, 140◦E–180◦ in the North Equatorial Current,southern part of the North Pacific subtropical gyre. Possibly besides the Ekman upwelling, the upward current contributes to the cooling of the upper ocean in the western tropical North Pacific through the upward advection of the
deep cold water. The upward current was found to be intensified (attenuated) by the strengthening (weakening) of
the subtropical gyre interior southward flow on the QD timescale in association with frequent occurrences of La Niña
(El Niño) events. The weakening of the upward current is responsible for the QD increase in upper-ocean temperaturein the western tropical North Pacific, which is equivalent to or larger than 0.2 ◦C, several years after the warm phaseof the PDO. The QD temperature variation may affect the phase reversal of the PDO through the heat transport of thewestern boundary current of the subtropical gyre, i.e., the Kuroshio.
Keywords: Western tropical Pacific, ENSO, PDO, Vertical velocity, P-vector inverse method, AGEM

1 Introduction
The subtropical gyre of the North Pacific is a basin-scale
anticyclonic circulation, consisting of the northward
swift western boundary current, so-called the Kuroshio, and the southward weak flow in the interior region.
The subtropical gyre transports a huge amount of heat
northward (e.g., Bryden and Imawaki 2001; Nagano et al.
2009, 2010). Through the variation in the northward*Correspondence: nagano@jamstec.go.jp
1Center for Coupled Ocean‑Atmosphere Research (CCOAR), Japan Agency
for Marine-Earth Science and Technology (JAMSTEC), 2‑15 Natsushima‑cho,
Yokosuka, Kanagawa 237‑0061, Japan
Full list of author information is available at the end of the article

heat transport, the subtropical gyre of the North Pacificis considered to play a critical role in interannual andlonger timescale midlatitude climate variations such as
the Pacific decadal oscillation (PDO) (e.g., Kawai et al.
2008), the most dominant variation in sea surface temperature (SST) in the North Pacific, (e.g., Mantua et al.
1997; Mantua and Hare 2002). Similar decadal oscillations were produced in ocean and atmosphere coupled
numerical models (Latif and Barnett 1994, 1996; Gu andPhilander 1997).
The spatial pattern of the PDO is defined as the first
empirical orthogonal function (EOF) mode of SST in the
North Pacific. A map of regression slope of SST to the
PDO index (the time coefficient of the first EOF mode of© 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/.

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

the SST) is illustrated in Fig. 1a. In the positive (warm)
phase of the PDO, the pattern exhibits a negative SST
anomaly in the central part of the North Pacific and apositive SST anomaly in the central and eastern equatorial Pacific (Mantua et al. 1997; Mantua and Hare 2002).
Thermal disturbances excited at the northern peripheryof the subtropical gyre due to the strengthening of the
Aleutian Low are observed to subduct in the upper main
thermocline (pycnocline), to propagate southward in theinterior region of the subtropical gyre, and, longer than
10 years after the subduction, reaches a subsurface (upto approximately 100 m depth) layer around a latitudeof 18◦N of the western tropical North Pacific (e.g., Deseret al. 1996; Schneider et al. 1999).
By combined use of altimetric sea surface height (SSH)
data, expendable bathythermograph temperature data,
and Argo salinity data at lines from San Francisco to
Hawaii and from Hawaii to Japan, Nagano et al. (2016a)
studied the quasi-decadal (QD) variation in the subtropical gyre net heat transport due to the variation of
the interior flow with periods of 8–12 years (e.g., Tourre
et al. 2001; White et al. 2003; Hasegawa and Hanawa
2003; Hasegawa et al. 2013). Up to 4 years after the positive (warm) PDO phase, the volume transport-weighted
temperature in the subtropical gyre interior region
increases due to the enhancement in volume transport ina potential density layer between 25.0 and 25.5 σθ (σθ =potential density −1000 kg m−3). As a result, because thetemperature difference of the water masses transportedby the Kuroshio and the gyre interior flow decreases, thenorthward net heat transport decreases. The reductionin the northward heat transport accumulates heat in thewestern tropical North Pacific, where climatological SST
is very high, and enhances the warm pool SST (Nagano
et al. 2016a), as is recognized in the map of the regression slope with a 3-year lag (Fig. 1d). In other words, thewarm pool SST increases prior to the negative phase of
the PDO.
The QD variation in the subsurface volume transport
in the interior region of the subtropical gyre is suggested to contribute to the SST variation in the regionof 10◦–20◦N, 125◦–160◦E in the western North Pacificwarm pool (Nagano et al. 2016a) possibly throughthe enhancement of the upwelling of deep water fromlayers deeper than 100 m depth. If we obtain a threedimensional (3D) field of vertical velocity, we can confirm this idea. The spatiotemporal characteristics of
upward currents in the western tropical North Pacificare, however, unknown because vertical current velocity is too small (< 10−4 m s−1) to be directly measured
by oceanographic instruments such as current meters,
the accuracy of current measurements being of theorder of 10−3 m s−1 at best for now (e.g., Thomson andPage 2 of 26

Emery 2014). Therefore, vertical velocity should be estimated from other physical variables such as 3D water
density field under some hydrodynamic assumptions.
Based on the conservations of potential density andpotential vorticity, Needler (1985) derived a formula of
3D absolute velocity vector ( v ), in the interior region ofsubtropical gyres asρv =gk·(∇ρ ×∇q)(∇ρ ×∇q),∇(f ∂ρ/∂q)·(∇ρ ×∇q)(1)

where ρ is the potential density, f is the Coriolis parameter, g is the gravitational acceleration, q = −f /ρ(∂ρ/∂z) is
Ertel’s potential vorticity neglecting the relative vorticity,and k is the unit vertical vector. In addition to potential
density and potential vorticity, the Bernoulli function of awater parcel advected by the flow is conserved if the flow
is steady and the diffusion can be neglected (e.g., Needler 1985; Pedlosky 1987a, b). For the estimation of the
3D current velocity field from hydrographic data, several
methods have been proposed by past investigators (e.g.,
Stommel and Schott 1977; Olbers et al. 1985; Killworth
1983, 1986; Chu 1995, 2006). These methods are expected
to provide similar 3D current velocity fields because they
are based on the same assumptions, i.e., conservations ofpotential density and potential vorticity.
Chu (1995, 2006) developed a sophisticated method,
called the P-vector inverse method, to calculate 3D
absolute velocity vector from observed hydrographic
data, avoiding cases that isosurfaces of potential density and potential vorticity are parallel (i.e., ∇ρ ×∇q =0), as described in the following section. For example,
he applied the P-vector inverse method to the circulations in the upper layer of the North Atlantic and demonstrated a very strong downward velocity around adepth of the main pycnocline (500 m depth) in the area
south of Iceland, where the subduction zone of the
global ocean conveyor belt is possibly located (Schmitz
1995). Therefore, the P-vector inverse method provides
reasonable 3D distributions of vertical velocity and may
be used to diagnose vertical velocity in the western
tropical and subtropical North Pacific.
By applying the P-vector inverse method to the 3D
potential density field in the tropical and subtropicalNorth Pacific, we expect to obtain a vertical velocity
field in the upper ocean layer and to reveal its seasonal
to decadal variations. In other words, we can identify
the location where the QD variation in the subsurface
volume transport in the interior region of the subtropical gyre reported by Nagano et al. (2016a) emerges in
the western Pacific region and how the QD variation of
the vertical transport to the east of the Philippines isrelated to the variation in the interior region. Further,

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

Page 3 of 26

Fig. 1 Regression slope ( ◦ C) between the Pacific decadal oscillation (PDO) index and sea surface temperature (SST) anomaly with lags of a 0, b 1, c
2, and d 3 years. Dashed square shows the study region of this paper

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

the impact of the vertical transport on the upper ocean
to the east of the Philippines may be estimated.
In general, potential vorticity is not conserved along
streamlines with horizontal and vertical diffusion. To
obtain the absolute velocity vector field by the P-vector
method and other methods based on the potential vorticity conservation, a time-varying 3D gridded water density field with sufficient spatial resolution is required. As
a candidate of such a dataset, we may use gridded hydrographic data such as the grid point value of the Monthly
Objective Analysis using the Argo (MOAA GPV) data
(Hosoda et al. 2008) based on global Argo float data collected since 2000. Using this dataset, we first attempted
to compute the P-vector-derived vertical velocity. However, the vertical velocity computed from the MOAAGPV data applying the P-vector method is unreasonably
large (> 10−4 m s−1) probably because the effective resolution (∼1000 km) is too low to apply the conservationof potential vorticity. Hence, we must use a finer gridded
dataset to obtain reliable vertical velocity estimates by
the P-vector inverse method.
Satellite observations have extensive coverage andsufficient resolutions in space and time to resolve theocean variations on scales ranging from mesoscale toglobal scale. SSH maps with a horizontal resolution of
1/4◦ ×1/4◦ are available to examine current variations(Pujol and Mertz 2020). Particularly, because SSH is aphysical parameter, which corresponds to integrated specific volume (i.e., reciprocal of water density) throughthe water column from the bottom to the sea surface,variations in altimetric SSH should be related to those inwater density through the water column. Thus, combined
use of the Argo hydrographic profile data and altimetric
SSH data might produce a sufficiently fine resolution 3Dhydrographic dataset.
The altimetry-based gravest empirical mode (AGEM)
method has been applied to spatiotemporally sparsehydrographic data to estimate variations of hydrographic
and current fields all over the world from SSH data (e.g.,
Sun and Watts 2001; Swart et al. 2010; Nagano et al.
2016a, b). Because AGEM-derived hydrographic fields
have the same horizontal resolution as that of the altimetric SSH data, they are sufficient to resolve regional
scale circulations and their temporal variations. Nagano
et al. (2016a, b) applied the AGEM method to hydrographic data and altimetric SSH data, and estimated sections of potential temperature and salinity (therefore,
potential density) at a line from San Francisco to near
Japan via Hawaii across the interior region of the North
Pacific subtropical gyre and 3D hydrographic fields in thewestern subarctic North Pacific.
In this study, we compute the 3D monthly potential temperature, salinity, and therefore potential density fields fromPage 4 of 26

the AGEM for areas divided by applying cluster analysis to
Argo hydrographic profile data. We then estimate the vertical velocity field by the P-vector inverse method and examine
the variations on seasonal to decadal timescales. In Sect. 2,
we describe Argo conductivity–temperature–depth (CTD)
profile data, altimetric SSH data, and the AGEM method. In
Sects. 3 and 4, we analyze the AGEM-derived potential density and the P-vector-based vertical velocity fields, focusing
on their variations related to the QD variations in the tropical and subtropical North Pacific and examine the impact ofthe vertical velocity variation on the upper-ocean temperature through the vertical advection of the deep cold water.
Summary and conclusion are provided in Sect. 5.

2 Data and methods
2.1 Argo data

To construct AGEM fields, we used more than 7000
temperature and salinity profiles from the near sea surface (10 dbar) to a pressure level of 1000 dbar (equivalent
to a depth of approximately 1000 m) collected for eachmonth by Argo floats (Advanced automatic QC Argo
Data version 1.2) from January 2001 to May 2017 in theregion of 10◦S–25◦N, 120◦E–180◦. A total of 90192 hydrographic profiles were used in this study. The accuraciesof the Argo temperature, salinity, and pressure data werereported to be ±0.005 ◦ C, ±0.01 (psu), and ±5 dbar by
Argo Science Team, respectively (Hosoda et al. 2008).
Following the Akima (1970) method, temperature andsalinity profiles were interpolated onto the vertical grids
at intervals of 10 dbar from 10 to 1000 dbar.
From the Argo temperature and salinity profile data,
we calculated spiciness parameter by using the polynomial proposed by Flament (2002), which isopleths are
orthogonal to those of potential density on the temperature–salinity diagram. The orthogonal character of the
spiciness parameter permits us to identify water masses
in terms of potential density by the single parameter.
Because of this, spiciness analyses are useful to distinguish water masses, as demonstrated in the Kuroshio andNorth Pacific subtropical gyre interior regions by Nagano
et al. (2009, 2010, 2012). We performed cluster analysis,
applying Ward’s criteria (e.g., Ward 1963; Anderberg
1973) to spiciness values averaged between isopycnal layers at intervals of 0.01 kg m−3 from 24.06 to 26.46 σθ andchanging the number of clusters from 3 to 9. The numberof clusters being set to 7, the classification based on thecluster analysis is found to be roughly consistent with the
distributions of water masses in the western tropical andsubtropical Pacific such as the North/South Pacific Tropical Waters and the North Pacific/Antarctic Intermediate
Waters (e.g., Fine et al. 1994; Kashino et al. 2011; Talley
et al. 2011). Thus, we classified Argo hydrographic profiles into 7 clusters.

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

We divided the study region (10◦S–25◦N, 120◦E–180◦)
into pieces of 5◦ on each side based on this cluster analysis result. Next, the dominant (the most frequent) clusterof hydrographic profiles for each piece was determined.
Based on the dominant clusters in all pieces, the study
region was separated into 7 areas in individual months
except for January (6 areas) and February (5 areas). The
area division based on the cluster analysis is found to beroughly consistent with the regional current and watermass pathways (Fig. 2) (e.g., Reid 1997; Talley et al. 2011).

Page 5 of 26

Using the Argo temperature and salinity profiles, we
calculated the dynamic height with the reference pressure level of 1000 dbar, HCTD, in meters as

HCTD =

1g1000 dbar

δ dp,
10 dbar(2)

where p is pressure and δ = ρ −1 − ρ −1 (0, 35, p) is the
specific volume anomaly relative to the specific volumeof seawater with temperature 0 ◦ C and salinity 35 (e.g.,
Gill 1982) calculated from Argo CTD data. The verticalFig. 2 Divisions of the study region up to 7 areas applying the cluster analysis to Argo CTD profiles collected from a January to l December. Area
numbers from 1 to 7 are indicated in the color bar in the lower right corner of each panel. Numeric values in areas divided by the cluster analysis
indicate correlation coefficients between altimetric SSH and sea surface dynamic height relative to the pressure level of 1000 dbar

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

Page 6 of 26

Fig. 2 continued

integral in Eq. (2) was discretized into the summationof 10-dbar averaged δ multiplied by the vertical intervalof 10 dbar from 10 to 1000 dbar. The error in the Argoderived HCTD was estimated to be smaller than 0.8 cm,
which is small compared with the altimetric SSH data
error (∼3 cm).
2.2 Altimetric SSH data

We used daily altimetric absolute SSH ( HALT) 1/4◦ × 1/4◦
gridded map data in the region of 10◦S–25◦N, 120◦E–180◦
from January 1993 to May 2017. The absolute SSH data

are derived from the delayed-time updated SSH (Level 4)anomaly data (DT_SLA-H) (Pujol and Mertz 2020) added
to the mean dynamic topography data (MDT_CNESCLS13) provided by Rio et al. (2011). Because the SSH
data are erroneous near the coasts, we did not use SSH
data in regions of water depths less than 1000 m to estimate AGEM fields and AGEM-derived 3D potential temperature and salinity fields.
HALT is linearly correlated with the dynamic height
derived from Argo CTD data, i.e., HCTD, with correlation coefficients of 0.60–0.98 (Fig. 2), which are muchhigher than the 1% significance level (i.e., within the 99%

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

confidence limit) based on the Student’s t test. Therefore, we performed linear regressions and based on thoseequations, converted values of HALT to dynamic height,
which is required to construct 3D potential temperatureand salinity fields by the AGEM method described in thenext subsection. Note that the correlation coefficients inthe equatorial and South Pacific are relatively low possibly because of high-frequency (time scales shorter thandays) variations by internal gravity waves or small scale
patchy water masses as observed by Richards et al. (2015)
and others, which cannot be fully resolved by the altimetric SSH observations.
2.3 Area‑wise AGEM method

Vertical structures of temperature and salinity are considered to be greatly different among the areas classified
by the cluster analysis and seasons. Thus, to consider
the spatial and seasonal variations, for each month andeach area, we constructed AGEM fields, i.e., the relations between hydrographic profiles and dynamic height,
HCTD, based on all the Argo CTD data in the area during
the month. The AGEM fields were smoothed for HCTD
by a Gaussian weight function with an e-folding scaleof 2.5 cm to remove noisy features due to the high-frequency variability.
As an example, the AGEM potential temperature andsalinity fields down to a depth of 1000 m for 6 areas inJanuary are shown in Figs. 3a–f and 3g–l, respectively.
The main thermocline illustrated by the sharp potential temperature decrease from 25 ◦ C to 8 ◦ C commonly
deepens as the dynamic height HCTD increases. Alongwith the main thermocline, the salinity maximum layer
deepens with HCTD although the salinity values are different for each area. Salinity in the areas in the Southern
Hemisphere (Fig. 3i, k) is higher than that in other areas(Fig. 3g, h, j, l), as illustrated by Reid (1997) and others.
Corresponding to the regional differences in salinity,
there are also regional variations in potential temperatureand density.
We identify differences in the vertical distributions ofpotential temperature and salinity for the season (notshown). In the upper layer shallower than 100 m, potential temperature, and salinity seasonally vary associatedwith the development of the mixed layer in winter andthe stratification in summer. Because, in this AGEM estimation, we used Argo hydrographic data during approximately 17 years from 2001 to 2017, the constructed
AGEM fields should account the interannual to decadal
variations. Such spatial and seasonal variations of the
AGEM fields beside the SSH variations are responsible
for those of the horizontal circulation and vertical current and their interannual to decadal modulations, as
described in the following sections.

Page 7 of 26

By converting monthly altimetric SSH, HALT, to HCTD
and referring to the AGEM fields, we estimated the temporal changes of the 3D fields of potential temperatureand salinity, and that of potential density on the horizontal grids of 1/4◦ at the vertical grids of 10 m from10 to 1000 m. Ranges of converted SSH, i.e., dynamic
height, variations in individual areas were within thoseof the dynamic height of the AGEM fields (abscissa in
Fig. 3). This suggests that the AGEM-derived 3D potential temperature and salinity fields involve the seasonal todecadal variations that are consistent with the SSH variations. To eliminate noisy features and to suppress artificial gaps between the areas in AGEM-derived potential
temperature and salinity fields for reasonable estimationby the P-vector inverse method, spatial smoothing was
performed using a Gaussian filter with a damping scale of
150 km. For the following analysis, we used the smoothed
AGEM-derived potential temperature and salinity fields.
In Fig. 4, we show the root-mean-square (RMS) potential density differences between Argo and AGEM-derived
data at depths of 100, 200, 400, and 800 m; the RMS
potential density difference at each grid was computed byusing data collected in the month by Argo floats withina radius of 1000 km. The RMS potential density differences are reduced with the depth from approximately
1.2 kg m−3 at 100 m to 0.1 kg m−3 at 800 m due to theweakening of the variations on shorter timescales thanapproximately 1 month. In addition, note that there are
no artificial spatial gaps in the RMS potential density
differences due to the area division based on the cluster
analysis.
Figure 5a shows the horizontal distribution of potential
density at a depth of 200 m averaged in all months during the study period. In addition, according to Rintoul
et al. (2002) and Nagano et al. (2016b), we calculated thestreamfunction defined as
1 1000 dbar
χ1000 =p δ dp.
(3)
g 10 dbar
The streamfunction in the mean state integrated fromthe sea surface to a pressure level of 1000 dbar (approximately 1000 m depth) is illustrated by contours in Fig. 5a.
The estimated mean current circulations are similar to
those observed by past investigators (e.g., Nitani 1972;
Reid 1997; Qu et al. 1998; Qu and Lukas 2003). The North
Pacific subtropical and tropical gyres to the north andsouth of approximately 13◦N are characterized by highand low potential density. Associated with the northward
sharp decrease in potential density at a latitudinal band
between 10◦N and 15◦N, there exists a strong current
toward the western boundary of the Pacific basin, i.e., theNorth Equatorial Current (NEC), which consists of theNagano et al. Progress in Earth and Planetary Science(2022) 9:57

Page 8 of 26

Fig. 3 AGEM (a–f) potential temperature ( ◦ C) and (g–h) salinity fields for January in terms of dynamic height relative to a pressure level of1000 dbar, HCTD, (m). AGEM potential temperature (salinity) fields from area 1 to area 6 are displayed in panels from a (g) to f (l). Ordinate is depth(m) and abscissa is HCTD. Dots at the top of each panel denote values of the dynamic height based on Argo-observed CTD data

westward current of the southern part of the subtropicalgyre and that of the northern part of the tropical gyre.
The Mindanao Dome (MD) is produced in the region
off the eastern coast of the Mindanao Island ( 4◦–10◦N)and the eastern edge of the circulation extends to the eastof 140◦E in boreal spring (Fig. 5c) and to the east of 150◦E
in boreal autumn (Fig. 5e). The remarkable seasonalchange of the circulation observed by Qu et al. (1998),
Qu and Lukas (2003), and others, and produced in anumerical model by Tozuka et al. (2002) are reproducedin this AGEM estimation.

2.4 P‑vector inverse method

By using a unit vector P , called the perfect vector (P-vector); and the magnitude of the velocity vector γ , called
the speed parameter, at location x = (x, y, z), Chu (1995,
2006) expressed Eq. (1) as v(x) = γ P . Here, using potential density ρ and potential vorticity q = −f /ρ(∂ρ/∂z), P
is defined as

P = (Px , Py , Pz ) = (∇ρ ×∇q)/|∇ρ ×∇q|.

(4)

Using the thermal wind equation, γ at a depth level of zk(γk ) is calculated as

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

Page 9 of 26

Fig. 4 Root-mean-square difference between AGEM-derived and Argo-observed potential density at depths of a 100, b 200, c 400, and d 800 m

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

Page 10 of 26

Fig. 5 Horizontal distributions of AGEM-derived potential density σθ at a depth of 200 m (kg m−3, color shading) and streamfunction χ1000(M kg s−1, contours) averaged in a all months, b January–March, c April–June, d July–September, and e October–December during 1993–2016.
Contour intervals of streamfunction are 0.1 M kg s−1. Letters H and L indicate maxima and minima of streamfunction, respectively


 


(r)  (k) (r)  �uk Px   Px Px 
γk =  ,
 �vk Py(r)   Py(k) Py(r) (5)where (uk , vk ) is the vertical shear vector of horizontal velocity relative to the reference depth level, zr . In thisstudy, we set zr to 1000 m.
Using the monthly 3D potential density fields estimated by the area-wise AGEM method, we calculatedNagano et al. Progress in Earth and Planetary Science(2022) 9:57

the 3D velocity vector below a depth of 100 m, where
the AGEM-based potential density is reliable. Note that
we did not conduct the calculation in the equatorial
region between 3◦S and 3◦N because |f| is close to zero
approaching the equator. There are two necessary conditions for the solution to exist. Isosurfaces of potential
density and those of potential vorticity are not parallel(∇ρ ×∇q = 0), and the denominator of Eq. (5) does notvanish (for the existence of the P-spiral). By examining
the two necessary conditions for the existence of theabsolute velocity vector, we can remove inappropriate
solutions. In this study, we did not calculate the P-vector and the speed parameter in the locations wherethese necessary conditions are not satisfied. The necessary conditions are satisfied in areas exceeding 86% of
the calculation region every month. Values of P and γ
in areas with no data were interpolated with a spatial
Gaussian weight function with a horizontal e-folding
scale of 150 km, which is much smaller than spatial
scales of the phenomena of our interests (> 500 km).
The estimation of the 3D absolute current velocity
field was conducted under some assumptions. In particular, the neglect of relative vorticity may produce
significant error in the vertical velocity. Taking thecurl of the momentum equation and using the continuity equation for incompressible fluid, we obtain∂w/∂z = (1/f )(∂ζ /∂t), where ζ = ∂v/∂x − ∂u/∂y isrelative vorticity (e.g., Gill 1982). The horizontal velocity is no more than 50 cm s−1 in the NEC and the horizontal scale of the variation is larger than 1000 km, sothat the magnitude of ζ is order of 10−7 s−1 except in
the western boundary region. Taking into account thevertical and time scales of the w variation, 100 m and1 month, respectively, the contribution of relative vorticity except in the western boundary region is lessthan 10−6 m s−1, which is smaller than the w variation.
In the western boundary region, the error larger than
10−6 m s−1 is considered to be included in the estimatedvertical velocity data.
By using the thermal wind equation, the vorv = (f /β)∂w/∂z ,tex-stretchingequation,givesρ∂ 2 w/∂z 2 = −(βg/f 2 )∂ρ/∂x , where β is the latitudinal variation of f (Killworth 1983). On the basis of thisequation, the error of the estimated vertical velocity
attributed to the potential density error can be evaluated. The RMS difference between Argo and AGEMderived potential density data in layers deeper thanthe main thermocline, where variations possibly due
to high-frequency internal gravity waves are relatively
small, may be regarded as the AGEM potential density error. If we set the error of the potential density to
0.1 kg m−3 , we evaluate the error of the w variation witha horizontal scale of 1000 km and a vertical scale ofPage 11 of 26

100 m to be 10−7 m s−1. The error is one order of magnitude smaller than the variation in w on the seasonal
and longer time scales, as will be discussed in the following sections.
To examine the characteristics of interannual andlonger timescale variations in vertical velocity, wesmoothed the data by a 15-month running mean filter.
We decomposed the interannual vertical velocity variation by EOF modes to reveal regional-scale variation
patterns, as the EOF mode decomposition is helpful to
examine the spatiotemporal characteristics of the variations (e.g., Thomson and Emery 2014). Based on the EOF
mode patterns, we determine the upwelling region where
we compute volume and temperature transports.
2.5 SST data

To examine the impact of vertical velocity on upperocean (< 100 m depth) temperature in the region eastof the Philippines, we used monthly US National Oceanic and Atmospheric Administration (NOAA) Optimum Interpolation Sea Surface Temperature (SST) V2(NOAA_OI_SST_V2) data (Reynolds et al. 2002) fromJanuary 1993 to May 2017 in the region of 10◦S–25◦N,
120◦E–180◦. We smoothed the data by a 15-month running mean filter to eliminate the seasonal variation andevaluated the temporal tendency of the temperature bytaking the smoothed data to the water temperature averaged in the top 100 m.
2.6 Climate indicesWe used the Niño 3.4 index from January 1993 to
December 2016, which is defined as the monthly mean
SST anomaly from the seasonal variation in the central
and eastern equatorial region of 5◦N–5◦S, 170◦W–120◦W
as an El Niño–Southern Oscillation (ENSO) index by
Trenberth (1997). Anomalously high and low SST eventsin the central and eastern equatorial Pacific, i.e., El Niño
and La Niña events, are indicated by the Niño 3.4 index.
In addition, the PDO index from January 1993 to December 2016, defined as the time coefficient of the first EOF
mode of SST in the region poleward of 20◦N for theperiod from 1900 to 1993 (Mantua et al. 1997), was used
to calculate the regression slope between the PDO index
and SST in the Pacific (Fig. 1) and to conduct comparison
with the time series of the vertical current transports inthe western tropical North Pacific.

3 Results
3.1 Mean circulation and vertical velocity fields

The area-wise AGEM estimation produced the horizontal upper-ocean circulations in the western parts of
the North Pacific subtropical gyre and the North/SouthNagano et al. Progress in Earth and Planetary Science(2022) 9:57

Pacific tropical gyres, as illustrated by the streamfunction in Fig. 5a. Feeding the subtropical gyre interior
southward flow, the westward flowing NEC is strengthened toward the western boundary from 148◦E to 135◦E(Fig. 5a). The cyclonic circulation of the MD centeredaround 8◦N, 130◦E is separated from the eastern tropical cyclonic circulation east of 165◦E. Associated with thewestward shrinkage of the MD in boreal winter (Fig. 5b)and spring (Fig. 5c), the westward current of the southernpart of the NEC around a latitude of 8◦N is attenuated inlongitudes between 145◦E and 168◦E as illustrated by thewidening and undulation of the χ1000 contours. When the
MD is weak, the NEC approaches the North Equatorial
Countercurrent (NECC), and it generates wavy patterns
between the currents. The wavy patterns may be attributed to fluctuations due to the barotropic instability (e.g.,
Pedlosky 1987a), which can occur in the enhanced horizontal shear between the opposing currents, as observedin the central and eastern equatorial regions (e.g., Tanaka
and Hibiya 2019). In boreal summer (Fig. 5d) and autumn(Fig. 5e), the MD extends eastward and strengthens theNEC.
Upward (positive) velocities of approximately
1×10−6 m s−1 are found in a part of the NEC, 10–16◦N,
140◦E–180◦, in layers of 100 m (Fig. 6a) and 200 m(Fig. 6b) depths. This upwelling region is located in theboundary between the subtropical and tropical gyres,where the meridional gradients of the potential density
and streamfunction are remarkable (Fig. 5a). Meanwhile,at a depth of 400 m, there is no systematic pattern of vertical velocity. As indicated by the vertical velocity maps at
100 m depth in Fig. 7, the upward velocity is extensively
observed in the region north of approximately 10◦N,
mostly in the subtropical gyre region, in boreal winter(Fig. 7a) and spring (Fig. 7b) but is limited in the latitudes
south of 15◦N in summer (Fig. 7c) and autumn (Fig. 7d).
The magnitudes of the vertical velocity tend to be amplified in low latitudes near the western boundary regions
because of γ ∼f −1 as derived from Eq. (5). Although values of vertical velocity at individual grids are similar toor a little larger than the error (10−6 m s−1) estimated in
Sect. 2.4, the error of mean value across such a regional
scale (> 2000 km) is evaluated to be 10−8 m s−1. Therefore, we focus on the variations with spatial scales equivalent to or larger than the upwelling region in the NEC.

Page 12 of 26

In Sect. 3.3, we will discuss the vertical volume and temperature transports across the upwelling region.
The seasonal changes in vertical velocity are also identifiable in the zonal and meridional sections; for example,
we show zonal sections at the latitude of 12.6◦N (Fig. 8)
and meridional sections at the longitude of 145.4◦E(Fig. 9). The upward and downward velocities were produced in a layer down to a depth of approximately 350 m
in association with the undulation of the isopycnal surfaces of 25.0–26.5 σθ in the longitudinal band between
140◦E and 170◦E. The intensity of the vertical velocity is greatest in boreal winter (Figs. 8b, 9b) and spring(Figs. 8c, 9c). In the region from 140◦E to 160◦E, theupward velocity appears to be stronger than the downward velocity; particularly, in the layer shallower than
150 m, upward velocity is predominant (Fig. 8).
The
considerably
large
verticalvelocity(> 5×10−6 m s−1) is generated in winter and springaround a latitude of 12◦N, where the main pycnocline
sharply shoals toward the equator along with the strongly
baroclinic NEC, as illustrated in Fig. 9b, c. The strong vertical velocity in boreal winter to spring is considered tobe driven by the winter intensification of the NEC, which
may be partly attributed to the variation of the subtropical gyre interior flow volume transport in a potential
density layer of 25.0–25.5 σθ reported by Nagano et al.
(2012, 2016a). Despite the great seasonal change in alayer between 150 and 350 m depths, the seasonal change
in the shallower layer around a latitude of 10◦N between
140◦E and 180◦ is remarkably small. Beneath a depth ofapproximately 350 m, where the vertical potential density
gradients are less than those in the overlying layer, themagnitudes of the vertical velocity are unsubstantial.
3.2 Interannual and longer timescale variations in verticalvelocity

The vertical velocity field varies on various time scales
longer than 1 year. As described in Sect. 2.4, we obtained
the interannual variation in vertical velocity by a
15-month running mean filter and calculated the EOF
modes. The percentages of variance explained by the firstand second EOF modes are 66% and 11%, respectively,
and the sum of these modes explains 77% of the total variance. Each higher EOF mode accounts for no more than
3% of the total variance. Therefore, the leading two EOF(See figure on next page.)
Fig. 6 Horizontal distribution of vertical velocity w at depths of a 100, b 200, and c 300 m. Vertical velocity was estimated by the P-vector inverse
method from the AGEM-derived potential density field. Positive values show upward velocity. Contour intervals are 1×10−6 m s−1. Gray shadings
indicate regions of no data because of small values of the Coriolis parameter and shallow water depth. Zonal and meridional thin solid lines denotelocations of the w sections shown in Figs. 8a and 9a. Dotted square delineates the region where Nagano et al. (2016a) studied the QD variation in
SST. In the solid square, vertical volume and temperature transports and SST tendency due to the vertical current variation were calculated as inFig. 14

Nagano et al. Progress in Earth and Planetary ScienceFig. 6 (See legend on previous page.)(2022) 9:57

Page 13 of 26

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

Page 14 of 26

Fig. 7 Horizontal distribution of vertical velocity w at a depth of 100 m averaged in a January–March, b April–June, c July–September, and d
October–December during 1993–2016. Positive values show upward velocity. Contour intervals are 1×10−6 m s−1. Gray shadings indicate regionsof no data because of small values of the Coriolis parameter and shallow water depth. Zonal and meridional thin solid lines denote locations of the
w sections shown in Figs. 8b–e and 9b–e. Dotted square delineates the region where Nagano et al. (2016a) studied the QD variation in SST. In thesolid square, vertical volume and temperature transports and SST tendency due to the vertical current variation were calculated as in Fig. 14

modes can explain most of the interannual to decadalvertical velocity variations.
The horizontal patterns of the first EOF modeat depths of 100 and 200 m are shown in Fig. 10a, b,
respectively. At 100 m depth, apart from the western boundary region, the first EOF mode varies witha zonal scale exceeding 2000 km such as in the negative loading region of 10◦–15◦N, 150◦E–180◦ (Fig. 10a),where the mean vertical velocity is largely upward, containing smaller scale (∼1000 km) variations (Fig. 6a).
The positive and negative alternating large loadings arepresent to the east of Mindanao Island (3◦–10◦N, 125◦–140◦E), but, because relative vorticity was not taken
into account, the vertical velocity estimated in thewestern boundary current region may be unreliable.
At a depth of 200 m (Fig. 10b), remarkably large
positive and negative loadings are aligned alternatively along a latitude of 12◦N with a zonal scale ofapproximately 500 km. The zonal scale of the verticalvelocity variation at this depth is smaller than that at
100 m depth. The area of either sign of the loading at
200 m depth appears not to be dominant in the region
10◦–15◦N, 150◦E–180◦. In both zonal (Fig. 10c) andmeridional (Fig. 10d) sections of the EOF mode, theremarkable variation is observed down to a depth ofapproximately 350 m, below which vertical gradient ofpotential density is small and the vertical velocity variation is weak. Therefore, the first EOF mode variationof the vertical velocity is responsible for the intensitychange of the upward current above approximately
350 m depth in the region of 10◦–17◦N, 140◦E–180◦
shown in Figs. 6, 8a, and 9a.
The second EOF mode of the vertical velocity at a depthof 100 m has a large negative loading in a region similar
to the first EOF mode but is limited within the region of
10◦–15◦N, 155◦E–180◦ (Fig. 11a). Regions equatorward of
10◦N and 10◦S show mostly positive values. At a depthof 200 m (Fig. 11b), positive and negative loadings with

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

Page 15 of 26

Fig. 8 Zonal sections of vertical velocity w (m s−1, color shades) and potential density σθ (kg m−3, contours) at the latitude of 12.6◦N averaged in a
all months, b January–March, c April–June, d July–September, and e October–December during 1993–2016. Positive values show upward velocity.
Contour intervals of σθ are 0.25 kg m−3

spatial scales smaller than approximately 500 km are sporadically distributed in the region south of 15◦N. The spatial patterns of the vertical velocity are identifiable in thezonal and meridional sections down to a depth of 290 m,
but in layers deeper than 300 m, we find noisy featureswith vertical scales of 10–20 m (Fig. 11c, d).
The time coefficient of the first EOF mode (blacksolid line in Fig. 12a) has peaks several months after
peaks of the Niño 3.4 index (gray solid line), i.e., El
Niño events, as in 1997–1998, 2002–2003, 2004–2005,

2009–2010, and 2015–2016. Similarly, troughs of thetime coefficients are found after La Niña events such as
in 1995–1996, 1998–2001, 2007–2008, and 2010–2011.
The time-lagged correlation coefficient between thefirst EOF mode and the Niño 3.4 index shows the maximum of 0.52 with a lag of 8 months (black solid line in
Fig. 13, Table 1). The maximum correlation coefficientis higher than the 5% significance level (i.e., within the
95% confidence limit) based on the Student’s t test.

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

Page 16 of 26

Fig. 9 Meridional sections of vertical velocity w (m s−1, color shades) and potential density σθ (kg m−3, contours) at the longitude of 145.4◦E
averaged in a all months, b January–March, c April–June, d July–September, and e October–December during 1993–2016. Positive values show
upward velocity. Contour intervals of σθ are 0.25 kg m−3

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

Page 17 of 26

Fig. 10 First EOF mode of vertical velocity w in the region of 10◦–25◦N, 120◦E–180◦ and the layer from 100 to 1000 m depths. The w data weresmoothed by a 15-month running mean filter. Loadings of the EOF mode at depths of 100 m and 200 m are shown in panels a and b, respectively.
Dotted square delineate the region where Nagano et al. (2016a) studied the QD variation and solid square the region (Region A) where we
calculated vertical volume and temperature transports and SST tendency due to vertical current variation as in Fig. 14. Sections of the EOF modeloading at zonal 12.6◦N and meridional 145.4◦E lines (thin solid lines) in a and b are indicated in panels c and d, respectively. Contour intervals are
0.001

Therefore, the first EOF mode variation is significantly
related to the ENSO.
Peaks of the time coefficient of the second EOF mode
are found almost simultaneously with those of the
Niño 3.4 index (Fig. 12b). With a lag of 1 month, thecorrelation coefficient reached the maximum of 0.74,
which is higher than the 5% significance level and istherefore statistically significant (black dashed line in
Fig. 13, Table 1). The second EOF mode variation is
excited by the ENSO before the first EOF mode variation and is more highly correlated with the ENSO than
the first EOF mode variation. The first EOF mode lags
9 months behind the second EOF mode with the maximum correlation coefficient of 0.71, which is muchhigher than the 5% significance level and is statisticallysignificant (pink solid line in Fig. 13).
3.3 Impacts of vertical velocity variation on SST

In this subsection, we examined the impact of the interannual and longer timescale variations in the verticalvelocity on the upper-ocean temperature in the region
10–17◦N, 140◦E–180◦, referred to as Region A in what
follows, where the vertical velocity at 100 m depth isprominently upward on average (shown by a solid squarein Fig. 6a) and quite variable (Fig. 10a).
We calculated vertical volume transport (Vz ) and temperature transport ( Fz ) at a depth of 100 m in Region A,
which are defined as

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

Page 18 of 26

Fig. 11 Same as Fig. 10, but for the second EOF mode of vertical velocity w

Vz =



w dxdy

(6)



(7)

A

and

Fz = ρCpA

θw dxdy,

where Cp is the specific heat capacity of seawater at constant pressure and θ is the potential temperature. ρCp is
set to 3.99×106 J K−1 m−3.
Time series of the vertical volume and temperaturetransports are shown by thick black and green lines
in Fig. 14a, respectively. Both modes of transports in
Region A are always positive, i.e., upward, and fluctuates inversely related to the PDO and Niño 3.4 indiceswith correlation coefficients of −0.80 and −0.79, respectively (Table 1), being consistent with the EOF analysis
in Sect. 3.2. In addition, the transport is relatively low(high) due to frequent occurrences of El Niño (La Niña)

events during periods such as 1994–1998 and 2002–2005(1999–2001 and 2007–2014). This relationship is clearlyshown by the QD component of the Niño 3.4 indexsmoothed by a 37-month running mean filter (Hasegawa
et al. 2013) (thin magenta line in Fig. 14a), which isinversely related to the vertical volume transport (thinblack line in Fig. 14a).
When the volume transport reached the maximum
over 1.5 Sv in summer (July) of the 2008 La Niña
year (downward arrow in Fig. 14a), the upward current velocity exceeded 1×10−6 m s−1 in the region
from 150◦E to 180◦ around 12◦N (Fig. 15b). The strong
upwelling was driven principally by the southwardcurrent greater than 1 cm s−1 at 100 m depth (contours in Fig. 15b) existing above the isopycnal zone of
25.0–25.5 σθ at 200 m depth (vertical line mesh), which
impinged on the slope of the isopycnals shoaling southward. In this period, the southward current in the isopycnal layer of 25.0–25.5 σθ was mostly distributed in theNagano et al. Progress in Earth and Planetary Science(2022) 9:57

Fig. 12 Time coefficients of the a first and b second EOF modes ofvertical velocity w (black lines) and Niño 3.4 index (gray lines). For
comparison, Niño 3.4 index was smoothed by a 15-month runningmean filter. All the time series were normalized by their standard
deviations

subtropical gyre interior region (10◦–25◦N, 145◦E–180◦)(Fig. 16b). The westward velocity of the NEC in thisisopycnal layer was higher than 15 cm s−1 to the east of160◦E.
Meanwhile, in the summer (July) of the 1997 El Niño
year (upward arrow in Fig. 14a), the vertical velocity
in Region A was weakly upward or downward, and thecorresponding southward current was observed only ina limited region centered around 14◦N, 175◦E (Fig. 15a).
The area of the subtropical gyre interior southwardcurrent in the isopycnal layer of 25.0–25.5 σθ and theintensity of the NEC are both smaller than those in the2008 summer (Fig. 16a). The QD variations in the vertical velocity just beneath a depth of 100 m are expected
to contribute to the QD variation in SST in the western
tropical North Pacific related to the PDO as illustrated
in Fig. 1.
Because the water temperature and vertical velocity
at 100 m depth are not horizontally uniform, the temperature of upwelling water is not necessarily identicalto the area-averaged temperature within Region A. The
water temperature of water mass transported by thevertical current can be examined by the volume transport-weighted temperature () at a depth of 100 m
defined asPage 19 of 26

Fig. 13 Time-lagged correlation coefficients of the first (black solidline) and second (black dashed line) EOF modes of vertical velocity w
with respect to Niño 3.4 index and that between the first and second
EOF modes (pink solid line). The time series were smoothed by a
15-month running mean filter. The positive time lag with respect toNiño 3.4 index indicates that the EOF modes of vertical velocity lag
the Niño 3.4 index. The positive time lag between the EOF modes
denotes that the first mode lags the second mode. The vertical solidline denotes the zero time lag. The horizontal dotted and solid lines
show the 95% confidence interval (±0.46) and zero correlation,
respectively

Table 1 Correlation coefficients of the first and second EOF
mode time coefficients of vertical velocity w, vertical volumetransport (Vz ), and vertical temperature transport (Fz ) at 100 m
depth with respect to Niño 3.4 index; and those of Vz , Fz , thetemperature tendency (∂Tw /∂t ) and temperature variation (Tw′ ) in
the top 100 m layer in Region A with respect to PDO index
Climate index
Niño 3.4

PDO

w-related variable

Correlation
coefficientTime lag(months)

EOF1_w

0.52

8

EOF2_w

0.74

1

Vz ,

−0.79

0

Fz

−0.79

0

Vz

−0.80

0

Fz

−0.80

0

∂Tw /∂t

0.63

0

Tw′

−0.65

−34

Positive time lag indicates that the w-related variables lag behind the climate
indexNagano et al. Progress in Earth and Planetary Science(2022) 9:57

Page 20 of 26

Fig. 14 a Volume transport Vz (Sv, thick black line) and temperature transport Fz (TW, green line) at 100 m depth in Region A shown by the solid
square in Figs. 6, 7, 10, and 11. Vz , Niño 3.4, and PDO indices smoothed by a 37-month running mean filter are shown by thin black, magenta, andcyan lines, respectively. Upward and downward arrows indicate July 1997 and July 2008, respectively. b Volume transport-weighted temperature ( ◦ C, gray solid line) and area-averaged temperature T100 ( ◦ C, gray dashed line) at 100 m depth and Niño 3.4 (magenta line) and PDO (cyan line)
indices. c Upper-ocean (< 100 m depth) temperature tendency, ∂Tw /∂t , averaged in Region A (◦ C month−1, red line) and temperature variation,
Tw′ , (◦ C, blue line) in the top 100 m layer. The horizontal black thin line indicates the mean value of the upper-ocean temperature tendency(−0.031◦ C month−1). The climate indices in panels a and b were normalized by the standard deviations

�=

Fz
.
ρCp Vz(8)

The volume transport-weighted temperature (gray solidline in Fig. 14b) varies on the QD timescale in opposite

phase with the Niño 3.4 index (magenta line in Fig. 14b)
because of the anomalously warm (cold) water in theupper layer of the western Pacific warm pool in La Niña
(El Niño). Therefore, the increase in  just below a depthof 100 m in La Niña years acts to reduce the cooling effectNagano et al. Progress in Earth and Planetary Science(2022) 9:57

Fig. 15 Horizontal distributions of smoothed upward velocity w(m s−1, color shading) and northward velocity v (cm s−1, contours)
at a depth of 100 m in a July 1997 and b July 2008, i.e., periods ofthe weak and strong upward transports in Region A (solid square),which are denoted by upward and downward arrows in Fig. 14a,
respectively. Data were smoothed by a 15-month running mean filter.
Contour intervals of v are 1 cm s−1 and cross line meshing indicatesnorthward velocity. Vertical line meshing illustrates the isopycnal
zone of at a depth of 200 mof the enhancement of the upward velocity. Because of
the nonuniformity of w,  is approximately 0.1 ◦ C higher
than the area-averaged temperature, T100 (gray dashedline in Fig. 14b). However, note that, as will be shown
below, the volume transport-weighted temperature variation is much less effective in the upper-ocean temperature variation than the upward current variation.
The upward volume and temperature transports in
Region A decrease with the depth because, in layers
deeper than 150 m, upward currents are compensatedby similar magnitude downward currents (Fig. 8). For
example, the volume and temperature transports at
200 m depth are mostly positive (upward) but smaller
than 0.8 Sv and 50 TW, respectively. Their correlation coefficients with respect to the Niño 3.4 index are−0.44 , absolute values of which are smaller than the 5%
significance level.
To examine the impact of the vertical velocity andtemperature variations at 100 m depth on the temperature in the overlying layer, we estimate temperature

Page 21 of 26

Fig. 16 Distributions of eastward velocity u (cm s−1, contours)
and northward velocity v (cm s−1, color shading) averaged within
isopycnal layers of 25.0 and 25.5 σθ in a July 1997 and b July 2008,
i.e., periods of the weak and strong upward transports in Region A(solid square), which are denoted by upward and downward arrows
in Fig. 14a, respectively. Data were smoothed by a 15-month runningmean filter. Contour intervals of u are 5 cm s−1 and dashed contours
indicate negative values (westward velocity)tendency (∂Tw /∂t ) averaged in the top 100 m due to thevertical velocity variation, w, i.e., vertical advection oftemperature (w∂T /∂z ), by the following formula:w(Tu − Td )
∂Tw
=−
,
∂t
hu

(9)

where hu is the depth of the uppermost level of the estimation of w (i.e., 100 m), Td is the potential temperature just below 100 m depth, and Tu is the upper-ocean(< 100 m depth) temperature. Td and Tu are set to the
AGEM-derived potential temperature at the 110 m depth
and monthly gridded value of SST prepared by Reynolds
et al. (2002), respectively. We assigned to w the timevarying vertical velocity at a depth of 100 m derived fromthe P-vector inverse method. When the vertical velocity
is negative (or downward), w is set to zero, i.e., w = 0.
The temperature tendency, ∂Tw /∂t , based on Eq. (9)and averaged in Region A is shown in Fig. 14c. The QD
variability (red line in Fig. 14c) is remarkable and is almostin phase with the PDO index (cyan line in Fig. 14b) with

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

a significant correlation coefficient of 0.63 (Table 1).
Associated with the enhancement of the vertical volumetransport (Fig. 14a), the temperature tendency has substantial troughs in 2000, 2009, and 2014, which occur
6–8 years after the peaks of the volume transport in the
isopycnal layer between 25.0 and 25.5 σθ in the subtropical gyre interior region in 1994, 2001, and 2006 estimatedby Nagano et al. (2016a). The minimum value of the temperature tendency (−0.04 ◦ C month−1) in the spring of2000 is approximately 30% of the annual mean coolingrates of the mixed layer due to the vertical entrainment
and the net surface heat flux, which were individuallyestimated to be approximately −0.13 ◦ C month−1 by Qu(2003); thus, the variation in the upward advection of the
deep cold water due to the vertical velocity variation has
a critical influence on the QD upper-ocean temperature
variation.

4 Discussion
4.1 Relation of the upward current cooling to the PDO

The temperature variation obtained by temporallyinte
grating the anomaly of ∂Tw /∂t , i.e., Tw′ = ∂Tw′ /∂t dt ,exhibits a smooth QD variation with peaks in 1999 and2006 (blue line in Fig. 14c). These peaks largely coincidewith those of the QD SST variation in the North Pacificwarm pool region calculated by Nagano et al. (2016a) in
the region marked by the dotted square (Figs. 6, 7, 10a,11a) and the troughs of the QD variation in the net heattransport examined by Nagano et al. (2016a). The peakvalues of Tw′ reach or exceed 0.2 ◦ C, which is equivalent
to the regression slope in the western tropical North
Pacific with a lag of 3 years as illustrated in Fig. 1d. Tw′
inversely varies leading 34 months to the PDO index witha significant correlation coefficient of −0.65 (Table 1). The
leading time (34 months) is a little longer than a leading
time (2 years) of the variation in coral nitrogen isotope
concentrations in Tatsukushi Bay, Japan, a possible proxy
for the surface Kuroshio water temperature, to the PDO
index reported based on 36-year time series by Yamazaki
et al. (2016). The discrepancy (8 months) between thetime leads of this study and that of Yamazaki et al. (2016)
can be explained by the advection time of approximately
9 months evaluated by Nagano et al. (2014, 2017) fromthe northward propagation of sea surface salinity variation along the western subtropical North Pacific.
Therefore, 6–8 years after the variation in the volumetransport in the isopycnal layer of 25.0–25.5 σθ in thesubtropical gyre (Nagano et al. 2016a), the vertical volume transport in the region east of the Philippines (10◦–17◦N, 140◦E–180◦) changes and it modifies the upperocean temperature through the vertical advection of thecold deep water to the upper ocean. Considering thePage 22 of 26

horizontal advection time scale of water masses in thesubtropical gyre interior region (∼10 years) (e.g., Sasaki
et al. 2010), the time lag between the variations in thevertical volume transport to the east of the Philippines
and the volume transport in an isopycnal layer of 25.0–
25.5 σθ (6–8 years) appears to be reasonable.
The relations of the upwelling strength in the NEC tothe positive and negative phases of the PDO are summarized in Fig. 17a, b, respectively. In the positive phase of
the PDO, the subtropical gyre interior flow in the isopycnal layer between 25.0 and 25.5 σθ is strengthened(Fig. 17a). The strengthened interior flow disturbanceproceeds southwestward and, 6–8 years after, reaches theeast of the Philippines, and cools the upper ocean in thewestern tropical North Pacific (Fig. 17b), leading to thenegative phase of the PDO.
The excessive accumulated heat in the western tropicalNorth Pacific during approximately 4 years by the attenuation of the vertical velocity is considered to be advected
northward by the Kuroshio (Fig. 17a) as inferred fromthe extension of the positive regression slope (high SST
after the positive PDO phase) in the western subtropical
Pacific (Fig. 1d). Note that the high SST anomaly is attributable not to the positive PDO phase 3–4 years before
but to the negative PDO phase approximately 10 years
before. The anomalously warm water is transported
to the midlatitude region of the North Pacific by the
Kuroshio and the Kuroshio Extension, yielding the next
negative phase of the PDO. The strengthening (weakening) of the upward transport in Region A in the negative(positive) phase of the PDO is critically related to thereversal of the PDO phase. Note that, in the region eastof the Philippines, there are mixed layer depth and barrier layer thickness variations (Sato et al. 2004; Nagano
et al. 2018), which may change responses of SST to thesurface heat flux and vertical current variations but cannot be accounted in this study. In particular, it is unclear
how the barrier layer thickness varies and affects the SST
variation. Future comprehensive studies are necessary tobetter understand the SST variation in the western North
Pacific related to the reversal of the PDO phase.
4.2 Effect of vertical eddy diffusion

The upward transport in Region A can be attributed toboth the adiabatic and diabatic processes. The adiabatic
process contributes to the upward transport through the
upward displacement of the main pycnocline (around
10◦N latitude) associated with the intensification of thenorthern part of the NEC. On the other hand, the isopycnal layer of 25.0–25.5 σθ , where Nagano et al. (2016a)
found the QD volume transport variation in the subtropical gyre interior region, does not reach the mixed layer

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

Page 23 of 26

Fig. 17 Schematic representation of the upwelling strength in the NEC region and temperature anomalies in the a positive and b negative PDO
phases. SST anomalies are denoted by solid shades and subsurface temperature anomalies dotted shades. High and low temperature anomalies arecolored red and blue, respectively. The sea surface Kuroshio and Mindanao currents are illustrated by solid arrows and the subsurface subtropicalgyre interior currents dashed arrows

base (< 100 m depth) (de Boyer Montégut et al. 2004)
in Region A, although the upward velocity in this region
extends from 250 to 100 m in all seasons (Figs. 8, 9). The
mechanism of the water mass transport across the isopycnal surfaces, i.e., the diabatic process, is considered
to be of importance in the process that the subsurface
water-mass volume fluctuations are entrained into thesea surface mixed layer.
In this subsection, we examine whether or not theeddy diffusive mechanism is important, by estimating

the order of the magnitude of diapycnal eddy diffusivity
from the AGEM potential density and vertical velocity.
Assuming that the vertical advection is balanced with thediapycnal diffusion:

w

∂ 2ρ
∂ρ
= Kρ 2 ,∂z∂zwhere Kρ is diapycnal eddy diffusivity (e.g., Pedlosky
1987b), we can estimate Kρ . On the basis of the AGEMbased potential density, we set ∂ρ/∂z and ∂ 2 ρ/∂z 2 toNagano et al. Progress in Earth and Planetary Science(2022) 9:57

1.0 × 10−2 kg m−4 and 1.0 × 10−4 kg m−5, respectively,
and for w = 1.0 × 10−6 m s−1, Kρ is estimated to be
1.0 × 10−4 m2 s−1. Considering that the isopycnals canmove adiabatically in the real ocean, this value is considered to be the upper limit of the diffusivity that should beobserved.
Whalen et al. (2012) estimated global distribution of
diapycnal eddy diffusivity based on Argo hydrographic
data obtained in the period from 2006 to 2011, during
which La Niña events frequently occurred. The estimatedvalue of Kρ in the present study is of the same order ofmagnitude as globally averaged diffusivity at a depthof 200 m (Fig. S.5a of Whalen et al. 2012) calculated by
Whalen et al. (2012). Note that the vertical transports are
relatively large (black and green lines in Fig. 14a) in Whalen et al.’s study period (2011–2016). There exists noticeably enhanced diffusivity in a layer from 250 to 500 m
depths over the rough ocean bottom in the western tropical North Pacific (Figs. 2a and S.1 of Whalen et al. 2012),
which covers Region A. Even though the diapycnal process is not accounted for in the P-vector inverse method,the enhancement of the eddy diffusivity in La Niña condition by the rough bottom topography is consistent with
the La Niña-related strengthening of the vertical velocity
estimated in this study. Future observations of verticaleddy diffusion by microstructure profilers in addition to
CTD data in both phases of the ENSO may be useful toclarify the relationship between the vertical velocity in
Region A and the bottom topography.
4.3 Vertical currents by the Ekman current divergencenear sea surface

Note that the vertical currents near the sea surface cannot be studied by the P-vector method because they are
considered to be strongly affected by the Ekman currents.
Near the sea surface in the tropical (south of 13◦N) andsubtropical (north of 13◦N) western North Pacific, thewind-driven Ekman currents are divergent and convergent, causing upwelling and downwelling, respectively(e.g., Fig. 7 of Tozuka et al. 2002). The upward velocity
induced by the Ekman currents is particularly strong in
winter and reaches the order of 10−6 m s−1. Although the
thickness of the Ekman layer is approximately 20 m evenin the tropical North Pacific (Fig. 4 of Ralph and Niiler
1999), which is remarkably less than 100 m depth, theupward current at 100 m depth in the western tropicalregion may be strengthened by the Ekman upwelling.

5 Conclusions
Prior to the negative (positive) phase of the PDO, the SST
in the western tropical North Pacific is heated (cooled).
The SST variation is suggested to be related to the QDtimescale (8–12 years) variation in the volume transport

Page 24 of 26

in the subsurface layer of 25.0–25.5 σθ (Nagano et al.
2016a). To reveal where the QD variation in the volumetransport in the subsurface layer emerges in the western
tropical North Pacific and how great the impact of thevariation on the upper-ocean temperature is, we examined the spatiotemporal characteristics of the 3D fieldof vertical velocity. By the area-wise AGEM method,
we estimated 3D potential temperature and salinity
fields in the region of 10◦S–25◦N, 120◦E–180◦ during theperiod from January 2001 to May 2017 using altimetric
SSH data and CTD profile data collected by Argo floats.
The 3D potential density field was computed from the
AGEM-derived potential temperature and salinity. The
streamfunction in the top 1000 dbar calculated from the
AGEM-derived potential density is consistent with previous studies.
Applying the P-vector inverse method to the 3D
AGEM-derived potential density field, we obtained thetime-varying 3D vertical velocity field. Vertical velocity
is observed in a layer down to a depth of approximately
350 m. Associated with the seasonal westward shrinkage of the MD in boreal winter and spring, the westward
flowing NEC approaches to the NECC and wavy patternsof the streamfunction and corresponding vertical velocity
fluctuations are generated between the currents possibly
through the barotropic instability due to the enhanced
horizontal shear. Furthermore, the regional-scale upward
current was found in the NEC, 10◦–17◦N, 140◦E–180◦,
called Region A in this paper. In this region, the southward shoaling of the isopycnal surfaces between 25.0 and26.5 σθ produces a mean substantial upward velocity of
approximately 1×10−6 m s−1 at a depth of 100 m depth.
Possibly in addition to the Ekman upwelling in the western tropical North Pacific, the upward current driven by
the southward current impinged on the isopycnal slope
advects the deep cold water to the upper ocean.
To examine the spatiotemporal characteristics, theinterannual to decadal variations in the vertical velocity were decomposed by the EOFs. The first and second
EOF modes represent the ENSO-timescale enhancement (attenuation) of the upward velocity at 100 m depth
in the region of 10◦–17◦N, 140◦E–180◦ (Region A) withpeaks approximately 8 and 1 months after La Niña (ElNiño) events, explaining 66% and 11% of the total variance, respectively. Furthermore, due to frequent occurrences of La Niña (El Niño) events, the upward velocity
was observed to vary on the longer timescale, i.e., the QD
timescale. The volume and temperature transports of thevertical velocity in the region of 10◦–17◦N, 140◦E–180◦
vary on the QD timescale. The maxima of the volume
and temperature transports are estimated to be approximately 1.5 Sv and 150 TW in the La Niña year of 2008,
respectively. At this time, the upward current reached theNagano et al. Progress in Earth and Planetary Science(2022) 9:57

maximum exceeding 1×10−6 m s−1 in Region A, being
forced by the southward current greater than 1 cm s−1
impinged on the sloping isopycnal surfaces.
The cooling effect of the upper ocean (< 100 m depth)
due to the upward current advection was examined. The
attenuation (intensification) of the upward transport
in the positive (negative) phase of the PDO keeps theupper-ocean temperature tendency, i.e., temporal change
in the temperature, anomalously high (low) approximately 6–8 years after the decrease (increase) in the volume transport in an isopycnal layer of 25.0–25.5 σθ in thesubtropical gyre interior region, as reported by Nagano
et al. (2016a), and increase (decrease) SST in the western tropical North Pacific, as schematically illustrated,
respectively, for the positive and negative PDO phases
in Fig. 17. The decreases in the vertical volume transportcause the upper-ocean temperature increase by 0.2 ◦ C or
greater, which is equivalent to the SST change associatedwith the PDO. The QD variation in the vertical velocity
in Region A can be critically related to the reversal of the
PDO phase through the heat transport of the Kuroshio.
This result, as well as atmospheric teleconnections fromthe tropical region (Alexander et al. 2002), suggests that
there is a significant oceanic link connecting the midand high-latitude regions to the tropical region.
Abbreviations
3D: Three-dimensional; AGEM: Altimetry-based gravest empirical mode; CTD:
Conductivity–temperature–depth; ENSO: El Niño-Southern Oscillation; EOF:
Empirical orthogonal function; MOAA GPV: Grid point value of the Monthly
Objective Analysis using the Argo; MD: Mindanao Dome; NEC: North Equatorial Current; NECC: North Equatorial Countercurrent; NOAA: US National
Oceanic and Atmospheric Administration; PDO: Pacific decadal oscillation;
QD: Quasi-decadal; RMS: Root mean square; SSH: Sea surface height; SST: Sea
surface temperature.
Acknowledgements
The 3D hydrographic field estimation by the AGEM method applying toArgo profile data owes a lot to discussion with Prof. Kimio Hanawa and Prof.
Shoichi Kizu (Tohoku University). The authors are grateful to the editors,
Dr. Michio Kawamiya (JAMSTEC) and Prof. Eiji Ohtani (Tohoku University), andanonymous reviewers for constructive review comments.
Author contributions
TH currently belongs to Japan Meteorological Agency. He contributed to thisstudy when he belonged to Tohoku University and Hokkaido University. AN
proposed the topic, conceived, designed the study, analyzed the data, andconstructed the manuscript. TH analyzed the data and collaborated with thecorresponding author (AN) in the construction of manuscript. MW discussed
the results in cooperation with AN and TH. All authors read and approved thefinal manuscript.
Funding
This work was partly supported by the Japan Society for the Promotion of Science, Grant-in-Aid for Scientific Research (Grant Nos. JP15H04228, JP17K05660,
JP20K04072, JP20H02236, JP20H04349, JP20KK0097).
Availability of data and materials
Quality-controlled Argo profile data (Advanced automatic QC Argo Data
version 1.2) were provided in the Japan Argo database (ftp://​ftp2.​jamst​ec.​
go.​jp/​pub/​argo/). Daily altimetric SSH data and monthly NOAA Optimum

Page 25 of 26

Interpolation SST V2 data were produced by the EU Copernicus Marine Environment Monitoring Service (CMEMS) (ftp://​my.​cmems-​du.​eu/) and the US
NOAA/OAR/ESRL (https://​psl.​noaa.​gov/​data/​gridd​ed/​data.​noaa.​oisst.​v2.​html).

Declarations
Competing interests
The authors declare that they have no competing interest.
Author details
1
Center for Coupled Ocean‑Atmosphere Research (CCOAR), Japan Agency
for Marine-Earth Science and Technology (JAMSTEC), 2‑15 Natsushima‑cho,
Yokosuka, Kanagawa 237‑0061, Japan. 2 Faculty of Environmental Earth Science, Hokkaido University, N10‑W5 Kita‑ku, Sapporo, Hokkaido 060‑0810,
Japan. 3 Graduate School of Science, Tohoku University, 6‑3 Aoba Aramaki‑aza
Aoba‑ku, Sendai, Miyagi 980‑8578, Japan. 4 Mutsu Institute for Oceanography(MIO), Japan Agency for Marine-Earth Science and Technology (JAMSTEC), 690
Kitasekine Sekine, Mutsu, Aomori 035‑0022, Japan.
Received: 24 June 2022 Accepted: 29 September 2022

References
Akima H (1970) A new method of interpolation and smooth curve fitting
based on local procedures. J Assoc Comput Mach 17:589–602
Alexander MA, Bladé I, Newman M, Lanzante JR, Lau N-C, Scott JD (2002) The
atmospheric bridge: the influence of ENSO teleconnections on air-sea
interaction over the global oceans. J Clim 15:2205–2231
Anderberg MR (1973) Cluster analysis for applications: probability and mathematical statistics. Academic Press, New York
Bryden HL, Imawaki S (2001) Ocean heat transport. In: Siedler H, Church J,
Gould J (eds) Ocean circulation and climate: observing and modeling the
global ocean. Academic Press, London, pp 455–474
Chu PC (1995) P-vector method for determining absolute velocity from hydrographic data. Mar Technol Soc J 29(3):3–14
Chu PC (2006) P-vector inverse method. Springer, Berlin. https://​doi.​org/​10.​
1007/​978-3-​540-​33386-9
de Boyer Montégut C, Madec G, Fischer AS, Lazar A, Iudicone D (2004) Mixedlayer depth over the global ocean: an examination of profile data and a
profile-based climatology. J Geophys Res. https://​doi.​org/​10.​1029/​2004J​
C0023​78
Deser C, Alexander MA, Timlin MS (1996) Upper-ocean thermal variations in
the North Pacific during 1970–1991. J Clim 9:1840–1855
Fine RA, Lukas R, Bingham FM, Warner MJ, Gammon RH (1994) The westernequatorial Pacific: a water mass crossroads. J Geophys Res 99(C12):25063–
25080. https://​doi.​org/​10.​1029/​94JC0​2277
Flament P (2002) A state variable for characterizing water masses and theirdiffusive stability: spiciness. Prog Oceanogr 54:493–501
Gill AE (1982) Atmosphere-ocean dynamics. Academic Press, London
Gu D, Philander SGH (1997) Interdecadal climate fluctuations that depend on
exchange between the tropics and extratropics. Science 275:805–807.
https://​doi.​org/​10.​1126/​scien​ce.​275.​5301.​805
Hasegawa T, Hanawa K (2003) Decadal-scale variability of upper ocean heat
content in the tropical Pacific. Geophys Res Lett. https://​doi.​org/​10.​1029/​
2002G​L0168​43
Hasegawa T, Ando K, Ueki I, Mizuno K, Hosoda S (2013) Upper-ocean salinity
variability in the tropical Pacific: case study for quasi-decadal shift during
the 2000s using TRITON buoys and Argo floats. J Clim 26:8126–8138.
https://​doi.​org/​10.​1175/​JCLI-D-​12-​00187.1
Hosoda S, Ohira T, Nakamura T (2008) A monthly mean dataset of global
oceanic temperature and salinity derived from Argo float observations.
JAMSTEC Rep Res Dev 8:47–59
Kashino Y, Ishida A, Hosoda S (2011) Observed ocean variability in the Mindanao Dome region. J Phys Oceanogr 41:287–302. https://​doi.​org/​10.​
1175/​2010J​PO4329.1
Kawai Y, Doi T, Tomita H, Sasaki H (2008) Decadal-scale changes in meridional
heat transport across 24°N in the Pacific Ocean. J Geophys Res. https://​doi.​org/​10.​1029/​2007J​C0045​25

Nagano et al. Progress in Earth and Planetary Science(2022) 9:57

Killworth PD (1983) Absolute velocity calculations from single hydrographic
sections. Deep-Sea Res 30(5A):513–542
Killworth PD (1986) A Bernoulli inverse method for determining the ocean
circulation. J Phys Oceanogr 16:2031–2051
Latif M, Barnett TP (1994) Causes of decadal climate variability over the North
Pacific and North America. Science 266:634–637
Latif M, Barnett TP (1996) Decadal climate variability over the North Pacific andNorth America: dynamics and predictability. J Clim 9:2407–2423
Mantua NJ, Hare SR (2002) The Pacific decadal oscillation. J Oceanogr 58:35–44
Mantua NJ, Hare SR, Zhang Y, Wallace JM, Francis RC (1997) A Pacific interdecadal climate oscillation with impacts on salmon production. Bull Am
Meteorol Soc 78(78):1069–1079
Nagano A, Ichikawa K, Ichikawa H, Konda M, Murakami K (2009) Synoptic flow
structures in the confluence region of the Kuroshio and the Ryukyu Current. J Geophys Res. https://​doi.​org/​10.​1029/​2008J​C0052​13
Nagano A, Ichikawa K, Ichikawa H, Tomita H, Tokinaga H, Konda M (2010)
Stable volume and heat transports of the North Pacific subtropical gyrerevealed by identifying the Kuroshio in synoptic hydrography south ofJapan. J Geophys Res. https://​doi.​org/​10.​1029/​2009J​C0057​47
Nagano A, Ichikawa H, Yoshikawa Y, Kizu S, Hanawa K (2012) Variation of thesouthward interior flow of the North Pacific subtropical gyre, as revealed
by a repeat hydrographic survey. J Oceanogr 68(2):361–368. https://​doi.​
org/​10.​1007/​s10872-​012-​0102-3
Nagano A, Uehara K, Suga T, Kawai Y, Ichikawa H, Cronin MF (2014) Origin of near-surface high-salinity water observed in the Kuroshio
Extension region. J Oceanogr 70:389–403. https://​doi.​org/​10.​1007/​
s10872-​014-​0237-5
Nagano A, Kizu S, Hanawa K, Roemmich D (2016a) Heat transport variation
due to change of North Pacific subtropical gyre interior flow during
1993–2012. Ocean Dyn 66(12):1637–1649. https://​doi.​org/​10.​1007/​
s10236-​016-​1007-2
Nagano A, Wakita M, Watanabe S (2016b) Dichothermal layer deepening in
relation with halocline depth change associated with northward shrinkage of North Pacific western subarctic gyre in early 2000s. Ocean Dyn
66(2):163–172. https://​doi.​org/​10.​1007/​s10236-​015-​0917-8
Nagano A, Hasegawa T, Ueki I, Ando K (2017) El Niño-Southern Oscillation-timescale covariation of sea surface salinity and freshwater flux in the western
tropical and northern subtropical Pacific. Geophys Res Lett 44:6895–6903.
https://​doi.​org/​10.​1002/​2017G​L0735​73
Nagano A, Ueki I, Hasegawa T, Ando K (2018) Ocean-atmosphere observations
in Philippine Sea by moored buoy. In: OCEANS - MTS/IEEE Kobe TechnoOcean (OTO), Kobe, Japan, pp 1–6. https://​doi.​org/​10.​1109/​OCEAN​SKOBE.​
2018.​85588​86
Needler GT (1985) The absolute velocity as a function of conserved measurable quantities. Prog Oceanogr 14:421–429
Nitani H (1972) Beginning of the Kuroshio. In: Stommel H, Yoshida K (eds) Kuroshio-its physical aspects. University of Tokyo Press, Tokyo, pp 129–163
Olbers DJ, Wenzel M, Willebrand J (1985) The inference of North Atlantic
circulation patterns from climatological hydrographic data. Rev Geophys
23(4):313–356
Pedlosky J (1987a) Geophysical fluid dynamics, 2nd edn. Springer, New York.
https://​doi.​org/​10.​1007/​978-1-​4612-​4650-3
Pedlosky J (1987b) Thermocline theories. In: Abarbanal HDI, Young WR (eds)
General circulation of the ocean. Springer, Berlin, pp 55–88
Pujol M-I, Mertz F (2020) Product user manual for sea level SLA products.
CMEMS, CMEMSQu T (2003) Mixed layer heat balance in the western North Pacific. J Geophys
Res. https://​doi.​org/​10.​1029/​2002J​C0015​36
Qu T, Lukas R (2003) The bifurcation of the North Equatorial Current in the
Pacific. J Phys Oceanogr 33:5–18
Qu T, Mitsudera H, Yamagata T (1998) On the western boundary currents in thePhilippine Sea. J Geophys Res 103(C4):7537–7548
Ralph EA, Niiler PP (1999) Wind-driven currents in the tropical Pacific. J Phys
Oceanogr 29(9):2121–2129
Reid JL (1997) On the total geostrophic circulation of the Pacific Ocean: flow
patterns, tracers, and transports. Prog Oceanogr 39:263–352
Reynolds RW, Rayner NA, Smith TM, Stokes DC, Wang W (2002) An improvedin situ and satellite SST analysis for climate. J Clim 15:1609–1625
Richards KJ, Natarov A, Firing E, Kashino Y, Soares SM, Ishizu M, Carter GS, Lee
JH, Chang KI (2015) Shear-generated turbulence in the equatorial PacificPage 26 of 26

produced by small vertical scale flow features. J Geophys Res 120:3777–
3791. https://​doi.​org/​10.​1002/​2014J​C0106​73
Rintoul SR, Sokolov S, Church J (2002) A 6 year record of baroclinic transportvariability of the Antarctic Circumpolar Current at 140°E derived fromexpendable bathythermograph and altimeter measurements. J Geophys
Res. https://​doi.​org/​10.​1029/​2001J​C0007​87
Rio MH, Guinehut S, Larnicol G (2011) New CNES-CLS09 global mean dynamic
topography computed from the combination of GRACE data, altimetry,
and in situ measurements. J Geophys Res. https://​doi.​org/​10.​1029/​2010J​
C0065​05
Sasaki YN, Schneider N, Maximenko N, Lebedev K (2010) Observational
evidence for propagation of decadal spiciness anomalies in the North
Pacific. Geophys Res Lett. https://​doi.​org/​10.​1029/​2010G​L0427​16
Sato K, Suga T, Hanawa K (2004) Barrier layer in the North Pacific subtropicalgyre. Geophys Res Lett. https://​doi.​org/​10.​1029/​2003G​L0185​90
Schmitz WJ Jr (1995) On the interbasin-scale thermohaline circulation. Rev
Geophys 33(2):151–173. https://​doi.​org/​10.​1029/​95RG0​0879
Schneider N, Miller AJ, Alexander MA, Deser C (1999) Subduction of decadal
North Pacific temperature anomalies: observations and dynamics. J Phys
Oceanogr 29:1056–1070
Stommel H, Schott F (1977) The beta spiral and the determination of theabsolute velocity field from hydrographic station data. Deep-Sea Res
24:325–329
Sun C, Watts DR (2001) A circumpolar gravest empirical mode for the Southern
Ocean hydrography. J Geophys Res 106(C2):2833–2855
Swart S, Speich S, Ansorge IJ, Lutjeharms JRE (2010) An altimetry-based gravest empirical mode south of Africa: 1. Development and validation. J
Geophys Res. https://​doi.​org/​10.​1029/​2009J​C0052​99
Talley LD, Pickard GL, Emery WJ, Swift JH (2011) Descriptive physical oceanography: an introduction, 6th edn. Academic Press, Boston. https://​doi.​org/​
10.​1016/​C2009-0-​24322-4
Tanaka Y, Hibiya T (2019) Generation mechanism of tropical instability waves inthe equatorial Pacific Ocean. J Phys Oceanogr 49(11):2901–2915. https://​doi.​org/​10.​1175/​JPO-D-​19-​0094.1
Thomson RE, Emery WJ (2014) Data analysis methods in physical oceanography, 3rd edn. Elsevier, Amsterdam. https://​doi.​org/​10.​1016/​
C2010-0-​66362-0
Tourre YM, Rajagopalan B, Kushnir Y, Barlow M, White WB (2001) Patterns ofcoherent decadal and interdecadal climate signals in the Pacific Basin
during the 20th century. Geophys Res Lett 28(10):2069–2072. https://​doi.​
org/​10.​1029/​2000G​L0127​80
Tozuka T, Kagimoto T, Masumoto Y, Yamagata T (2002) Simulated multiscalevariations in the western tropical Pacific: the Mindanao Dome revisited. J
Phys Oceanogr 32:1338–1359
Trenberth KE (1997) The definition of El Niño. Bull Am Meteor Soc
78:2771–2777
Ward J (1963) Hierarchical grouping to optimize an objective function. J Am
Stat Assoc 58:236–244. https://​doi.​org/​10.​2307/​22829​67
Whalen CB, Talley LD, MacKinnon JA (2012) Spatial and temporal variabilityof global ocean mixing inferred from Argo profiles. Geophys Res Lett.
https://​doi.​org/​10.​1029/​2012G​L0531​96
White WB, Tourre YN, Barlow M, Dettinger M (2003) A delayed action oscillator
shared by biennial, interannual, and decadal signals in the Pacific Basin. J
Geophys Res. https://​doi.​org/​10.​1029/​2002J​C0014​90
Yamazaki A, Watanabe T, Tsunogai U, Iwase F, Yamano H (2016) A 150-year
variation of the Kuroshio transport inferred from coral nitrogen isotope
signature. Paleoceanography 31:838–846. https://​doi.​org/​10.​1002/​2015P​
A0028​80

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