1 Introduction

The North Pacific subtropical mode water (STMW) is a distinct water mass that exists between the upper pycnocline and the lower main pycnocline in the western subtropical North Pacific (e.g., Masuzawa 1969; Suga et al. 1989; Hanawa and Talley 2001; Oka and Qiu 2012; Feucher et al. 2019). It appears as a layer of nearly vertically homogeneous properties and is dynamically characterized by a vertical minimum in potential vorticity (PV; Fig. 1a). The STMW forms in the deep mixed layer south of the Kuroshio and the Kuroshio Extension in winter, which arises from strong ocean surface cooling and wind stirring in conjunction with weak vertical stratification below the mixed layer (Tozuka et al. 2017). The STMW is then subducted and transported to a wide region of the subtropical gyre by the mean gyre circulation (Huang and Qiu 1994; Suga et al. 2004) and mesoscale eddies (Nishikawa et al. 2010; Xu et al. 2014, 2016).

Fig. 1
figure 1

Climatological mean cross sections of (a) potential vorticity (PV) and (b) meridional potential density (\({\sigma }_{\theta }\)) gradients, \(F={\left(\partial {\sigma }_{\theta }/\partial y\right)}_{z}\), in the summer 137°E section from the Japan Meteorological Agency (JMA) repeat hydrographic observations. Contours in both panels are the mean \({\sigma }_{\theta }\)

The STMW has dynamic effects on surface circulation (Kobashi and Kubokawa 2012). Low-PV waters subducted from different locations in the Kuroshio Extension are advected on different isopycnals to the south and overlap vertically because of the beta spiral effect. They become stacked and push the upper pycnocline upward, forming a subtropical front (STF) along the southern flank of the STMW (Kubokawa and Inui 1999; Kubokawa 1999). The STF is manifested as the northward shoaling of the upper pycnocline at subsurface depths of 50–200 m in the latitudinal range of 17°–27°N (Fig. 1b; Uda and Hasunuma 1969). The STF can also be generated by stationary Rossby waves emanated in the presence of the STMW (Kubokawa 1997). The STF creates an eastward shear to the flow in the upper ocean, maintaining the surface current called the subtropical countercurrent (STCC). The STCC causes warm advection and heat flux out of the ocean along the STCC (Xie et al. 2011; Kida et al. 2015). A sea surface temperature (SST) front that forms in the STCC region in winter to spring increases the near-surface baroclinicity in the atmosphere (Kobashi et al. 2008). These atmospheric effects exert a great influence on the precipitation and winds along the STCC (Kobashi et al. 2008; Xie et al. 2011; Zhang et al. 2017). Mode waters also influence large-scale structures of the subtropical gyre (Kimizuka et al. 2015).

The formation of the STMW changes substantially at decadal timescales due to the variability in the wintertime mixed layer depth (MLD) in the STMW formation region. Several mechanisms have been proposed for the decadal MLD variations. One mechanism emphasizes local surface cooling in controlling the wintertime MLD (Bingham 1992; Iwamaru et al. 2010). The other mechanisms are related to the ocean dynamics associated with westward-propagating Rossby waves from the central North Pacific. Basin-scale changes in wind generate pycnocline depth anomalies, which propagate into the Kuroshio Extension region and lead to decadal changes in the Kuroshio Extension jet, its southern recirculation gyre and associated eddies (Qiu 2003; Qiu and Chen 2005, 2010b; Ceballos et al. 2009). The changes in the pycnocline depth and eddy activity modify the upper-ocean stratification and thus regulate the development of the wintertime mixed layer in the STMW formation region (Qiu et al. 2007; Sugimoto and Hanawa 2010). Recently, Sugimoto and Kako (2016) analyzed the historical temperature profiles and revealed the decadal-scale variations in the wintertime MLD in the period of 1968–2014. They showed that the MLD variations were mainly caused by local surface cooling in the late 1970s and 1980s and by ocean dynamics in and after the 1990s. The result was supported by Kim et al. (2020), who analyzed an ocean general circulation model hindcast simulation for 1960–2009.

The decadal variations in the STMW spread downstream away from the formation region (Qiu and Chen 2006). Oka et al. (2019) analyzed more than 50-year-long time series of the Japan Meteorological Agency (JMA) repeat hydrographic observations along the 137°E meridian (Oka et al. 2018) and found a significant correlation at decadal timescales between the winter cross-sectional area of the STMW at 137°E and the winter MLD in the STMW formation region with a lag of 1 year. The decadal signals are also evident in the biogeochemical properties of the STMW observed at 137°E (Oka et al. 2019) and further downstream near the western boundary (Oka et al. 2015).

Mode waters shape the subtropical pycnocline through ventilation, in which PV is conserved while they are advected downstream (Kubokawa 1999; Dewar et al. 2005). Eddy-resolving ocean general circulation models (Yamanaka et al. 2008; Nonaka et al. 2012; Sugimoto et al. 2012) and coupled climate models (Xie et al. 2011; Xu et al. 2012) revealed that the STF is anchored by mode waters in the climatology, in agreement with observations (Aoki et al. 2002; Kobashi et al. 2006). They also showed that the advection of thick, low-PV mode waters causes the upper pycnocline to shoal more than normal, intensifying the STF and STCC on decadal timescales. This mode water-induced change in the STCC is the dominant mechanism for ocean temperature variability in the subtropical gyre (Xie et al. 2011). Xie et al. (2011) proposed a second mode baroclinic adjustment of the pycnocline as the mechanism, in which thick mode waters shoal the upper pycnocline while deepening the lower main pycnocline. The mechanism, however, has not been verified yet. Furthermore, the relationship between the changes in the STF and the mode waters has not been explored from observations. The eddy-resolving models analyzed in the previous studies cannot reproduce the STF in the western North Pacific possibly due to strong dissipation in some models (Nonaka et al. 2012). The mode water simulated by climate models is too strong, and the mode water dynamics may be exaggerated (Xie et al. 2011). The results need to be verified with observations.

Recently, Cerovečki and Giglio (2016) and Cerovečki et al. (2019) analyzed Argo temperature and salinity profiles for 2005 to 2012 and showed a marked decrease in the STMW volume from 2006 to 2009. The decrease coincided with the shoaling of the bottom of the STMW layer and the deepening of the top of the STMW layer, which apparently agrees with the baroclinic adjustment. Their analysis, however, focused on the entire STMW region, including not only the downstream region but also the formation region. As mentioned above, in the formation region, the change in the main pycnocline causes that in the STMW. It is still unclear how the STMW alters the downstream pycnocline and the STF on decadal timescales.

In this study, we analyze the JMA repeat hydrographic observations at 137°E, the same as that used by Oka et al. (2019). Figure 2 shows the climatology of the STMW thickness and the acceleration potential on the core isopycnal surface of the STMW in summer (June to August) and the MLD in late winter (February to March). The 137°E section crosses the southwest of the major portion of the STMW and is located downstream of the formation region of the deep mixed layer south of the Kuroshio Extension. Because the 137°E observations have continued for more than 50 years, they enable us to study long-term oceanic variations. The present study addresses the following questions: how does the change in the STMW thickness alter the upper and lower main pycnoclines at 137°E on a decadal timescale? Does the STMW affect the variability of the STF? We focus on the summer season when the upper pycnocline develops well above the STMW.

Fig. 2
figure 2

Climatological mean North Pacific subtropical mode water (STMW) thickness (colors) and acceleration potential referenced to 1000 dbar (dashed contours every 0.5 m2s−2) from June to August. The thickness is drawn for the values greater than 20 m. Black thick contours represent 140 and 180 m of the climatological mean mixed layer depth (MLD) from February to March. All the plots are based on the temperature and salinity profiles from the Grid Point Value of the Monthly Objective Analysis using the Argo data (MOAA-GPV). Circles denote grid points of the 137°E repeat hydrographic observations (color figure online)

The rest of the paper is organized as follows. Section 2 describes the data and methods used in this study. In Sect. 3, we first examine the variations in the STMW thickness and their relationships with the mixed layer in the formation region south of the Kuroshio Extension. Then, we explore the variations in the pycnocline at 137°E and make a comparison with those of the STMW thickness. Because a pycnocline fluctuates due to basin-scale wind-driven baroclinic Rossby waves (e.g., Miller et al. 1998; Qiu 2003; Taguchi et al. 2007; Ceballos et al. 2009), we also examine the effect of Rossby waves on the pycnocline variations with a focus on first-mode baroclinic Rossby wave processes that are important for large-scale pycnocline changes (Qiu 2002). We show the importance of the STMW in changing the pycnocline at 137°E. Then, we investigate the variations of the STF and show that they are explained by variations of the STMW on decadal timescales. Section 4 provides a summary and discusses the results.

2 Data and methods

2.1 Observation data and method

The JMA 137°E repeat hydrographic observations have been conducted biannually in winter (January) since 1967 and in summer (mostly July) since 1972. We use optimally interpolated temperature and salinity profiles in summer and winter for the period of 1972 to 2019, which are obtained from the JMA website. They are gridded for each cruise at 1 dbar vertical intervals, every 1° meridional interval south of 30°N, 30′ intervals between 30° and 31°N and every 20′ intervals between 31° and 34°N (Fig. 2). The deepest level is 1250 dbar before 1992 and 2000 dbar afterwards. The interpolation procedure is described by Nakano et al. (2007).

In this study, the temperature and salinity profiles are first smoothed vertically by applying a three-point Hanning filter 10 times to reduce features smaller than 20 dbar in the vertical direction. Then, we calculate the potential temperature (\(\theta\)) and potential density (\({\sigma }_{\theta }\)). PV is defined as \(PV=gf \partial {\sigma }_{\theta }/\partial p\), ignoring relative vorticity, where \(g\) is the gravitational acceleration, \(f\) is the Coriolis parameter and \(p\) is the pressure. The MLD is computed from the shallower depth of either where \(\theta\) or \({\sigma }_{\theta }\) differs by 0.5 °C or 0.125 kg m−3, respectively, from the corresponding values at a depth of 10 m (Hosoda et al. 2010).

The MLD in the STMW formation region and the basin-scale change in the pycnocline are examined using the Grid Point Value of the Monthly Objective Analysis using the Argo data (MOAA GPV) produced by Hosoda et al. (2008). The data consist of monthly temperature and salinity profiles at standard pressure levels on a 1° grid in space from 2001, which is produced from an optimal interpolation of temperature and salinity profiles from real-time quality-controlled Argo data and available ship observations. We use the data from 2005 to 2019, when the Argo observations cover the STMW region well (Toyama et al. 2015).

In this study, the STMW is defined as a layer of PV less than a threshold value, with the vertical PV minimum core in the STMW density range. The method is the same as that used by Oka et al. (2019), except that these authors used the temperature range instead of the density range. The threshold value used is 2.0 \(\times\) 10–10 m−1 s−1. The density range is determined from the histogram of the vertical PV minimum plotted against \({\sigma }_{\theta }\)(not shown) and is between \({\sigma }_{\theta }\) = 24.8 and 25.8 kg m−3. The STMW thickness is defined as the total thickness of the STMW at each observation point in each year. When there is no STMW that satisfies our criterion, we set the STMW thickness to be zero.

The STF is detected based on the strength of a density front (\(F\)), which is defined as \(F={\left(\partial {\sigma }_{\theta }/\partial y\right)}_{z}\) on the depth coordinate \(z\). It is estimated by a least square fitting using five grid-point values including two adjacent grid points north and south of each grid. This method has an advantage in that it is not as sensitive to small-scale features compared to the difference method between the two grid points; thus, this method can reduce the contamination of vigorous eddies in the STF region (Qiu 1999; Kobashi and Kawamura 2001). Figure 1b shows the mean section of \(F\) in summer, where the broad STF band appears along the southern flank of the low-PV STMW (Fig. 1a) with two subsurface narrow bands at 19°N and 25°N. These two bands are notable characteristics of the STF (Aoki et al. 2002; Kobashi et al. 2006), successfully captured by the present method. The subsurface STF between below the MLD and above the top of the STMW is hereafter referred to as the STF.

In the present study, we adopt a three-year running mean as the low-pass filter to highlight decadal signals. We conduct cross correlation analysis and assess the statistical significance of the results at the 95% confidence limit using the degrees of freedom estimated from the length of the time series on the assumption of the dominant timescale of a decade.

2.2 Linear baroclinic Rossby wave model

Wind-driven baroclinic Rossby waves are examined by adopting a 1.5-layer reduced gravity model in which the main pycnocline is represented by the interface between the upper layer and the infinitely deep bottom layer. The model captures first-mode baroclinic Rossby wave processes. Under the longwave approximation, the linear vorticity equation is expressed as follows:

$$\frac{\partial h}{{\partial t}} - C_{R} \frac{\partial h}{{\partial x}} = - w_{e} - \varepsilon h,$$
(1)

where \(h\) is the upper layer thickness anomaly,\({C}_{R}\) is the zonal phase speed of the long baroclinic Rossby wave,\({w}_{e}=curl\left({\varvec{\tau}}/{\rho }_{0}f\right)\) is the Ekman pumping velocity, \({\varvec{\tau}}\) is the wind stress vector, \({\rho }_{0}\) is the mean density of sea water, \(f\) is the Coriolis parameter, and \(\varepsilon\) is the Newtonian dissipation rate (Meyers 1979; Qiu 2002; Qiu and Chen 2010b). We regard \(h\) as the main pycnocline depth. Given the wind stress vector, \(h\) is solved by integrating Eq. (1) from the eastern boundary (\(x={x}_{e}\)) along the baroclinic Rossby wave characteristic:

$$h\left( {x,y,t} \right) = \int\limits_{{x_{e} }}^{x} {\frac{1}{{C_{R} }}} w_{e} \left( {x^{\prime},y,t + \frac{{x - x^{\prime}}}{{C_{R} }}} \right)\exp \left[ {\frac{ \varepsilon }{{C_{R} }}\left( {x - x^{\prime}} \right)} \right]dx^{\prime}.$$
(2)

We ignore the upper layer thickness anomalies along the eastern boundary, \(h\left({x}_{e},y,t\right)=0\), because its influence is limited to the area near the boundary (Fu and Qiu 2002). The Rossby wave speed is estimated using the relation, \({C}_{R}\left(x,y\right)=\beta {g}^{^{\prime}}{h}_{0}/{f}^{2}\), where \(\beta\) is the meridional derivative of \(f\),\({g}^{^{\prime}}\) is the reduced gravity and \({h}_{0}\) is the mean thickness of the upper layer. We set \({g}^{^{\prime}}=0.03\) m s−2 following Qiu (2002). \({h}_{0}\left(x,y\right)\) is calculated as the depth of the isopycnal surface using the annual climatology of temperature and salinity from the World Ocean Atlas (WOA) 2013 (Locarnini et al. 2013; Zweng et al. 2013) . The isopycnal surface is decided so that the \({C}_{R}\) values become close to those of Killworth et al. (1997), who showed Rossby wave speed in the presence of a mean flow. As a result, we choose an isopycnal depth of \({\sigma }_{\theta }\)= 26.5 kg m−3 for \({h}_{0}\).\({C}_{R}\) is a key parameter for the estimation of \(h\). We carried out the same analysis by changing the \({C}_{R}\) values by 5% and 10% and found that the results do not change substantially in a qualitative sense.

The main pycnocline depth anomalies are calculated using the monthly surface wind stress curl computed from JRA–55 monthly surface wind stress vector on a 1.25° grid from 1958 to 2019 (Kobayashi et al. 2015). The dissipation rate is set to \(\varepsilon =1/6\) yr−1 following Qiu (2002). To compare with the observation at 137°E, the pycnocline depth anomalies on the JRA–55 grid points are spatially interpolated onto the 137°E observation points for every month using a Gaussian filter with an e-folding scale of 150 km.

3 Results

3.1 Variations in the STMW thickness

Figure 3 shows the time–latitude plot of the STMW thickness anomalies in the summer 137°E section together with the average and standard deviation. The anomaly is the deviation from the average over the whole period from 1972 to 2019. The mean thickness and the standard deviation are large between 26° and 32°N and decrease to the south (Fig. 3b, c). The thickness variations seem to be similar in the region between 20° and 32°N (Fig. 3a). Figure 4a shows the time series of the thickness anomalies averaged between 26° and 29°N, displaying decadal-scale variations. The autocorrelation function was calculated at each latitude between 20° and 32°N, showing a dominant timescale of approximately 9–15 years (not shown). The increases from the early 1990s, early 2000s and early 2010s and their subsequent decreases are in agreement with those identified from the Argo float and other hydrographic observations by Qiu and Chen (2006); Oka (2009); Qiu and Chen (2013); Rainville et al. (2014); Oka et al. (2015) and Cerovečki and Giglio (2016).

Fig. 3
figure 3

a Latitude-time plot of low-pass filtered STMW thickness anomalies in the summer 137°E section, along with b the mean thickness and c the standard deviation. The standard deviation is calculated from the low-pass filtered time series. Horizontal bars at the top in (a) denote the large-meander periods of the Kuroshio

Fig. 4
figure 4

Time series of a the STMW thickness anomalies averaged between 26° and 29°N in the summer 137°E section. b February–March mean MLD in the region south of the Kuroshio Extension (29.5°–33.5°N, 141.5°–155.5°E), based on the MOAA-GPV data. Gray lines denote the yearly values, and black lines denote the low-pass filtered values (color figure online)

The composite PV sections are made with respect to the STMW thickness anomalies averaged in 26°–29°N (Fig. 4a): one is made by averaging the data when the positive thickness anomalies exceed one standard deviation (thick STMW years; 1981–1983, 1995–1997 and 2013–2015), and the other is made by averaging the data when the negative anomalies exceed one standard deviation (thin STMW years; 1977–1978, 2000–2001 and 2008–2010). In the thick STMW years, the low-PV STMW appears with the large area north of 26°N and extends to the south around 22°N (Fig. 5a). In contrast, the STMW is almost absent along the 137°E section in the thin STMW years (Fig. 5b).

Fig. 5
figure 5

Composite PV cross sections in a thick STMW years (1981–1983, 1995–1997 and 2013–2015) and b thin STMW years (1977–1978, 2000–2001 and 2008–2010) in the summer 137°E section. These years are chosen based on the time series of the STMW thickness anomalies between 26° and 29°N in Fig. 4a

The STMW in the 137°E section originates in the winter mixed layer south of the Kuroshio Extension (Suga et al. 1989; Suga and Hanawa 1995). Indeed, the comparison between the STMW thickness at 137°E and the February-to-March mean MLD south of the Kuroshio Extension, which is calculated using the MOAA-GPV data, shows the close relationship on a decadal timescale (Fig. 4b), though the time series is short because of the limitation of the available observations. As demonstrated by Oka et al. (2019), the change in the MLD heads slightly that in the STMW thickness and their low-pass filtered time series are significantly correlated (a coefficient of 0.89), with a lag of 1 year. This suggests that the change in the STMW thickness in the summer 137°E section originates in the wintertime mixed layer south of the Kuroshio Extension in the preceding year.

In the region north of 30°N, the STMW thickness changes in relation to the Kuroshio path. The STMW became thin from 1975 to 1979, from 1981 to 1983 and after 2018 during the Kuroshio large meander periods, as indicated by the horizontal bars in Fig. 3a. During the Kuroshio large meander, the Kuroshio takes an offshore detouring path between 135° and 139°E (Kawabe 1995) and separates the recirculation gyre south of the Kuroshio, hindering the advection of the STMW from the east north of 30°N (Suga and Hanawa 1995). The reduction in the STMW thickness is not as obvious during the other large meander periods, which are relatively short (less than 1.5 years) compared to the three meanders mentioned above. This is due to the low-pass filter used in Fig. 3a. Without applying the low-pass filter, the short-term Kuroshio meander is accompanied by the reduction in the STMW thickness (figure not shown).

3.2 Variations in the pycnocline

In this section, we first describe the variations in the pycnocline depth and then examine their relationships with the STMW thickness and wind-driven baroclinic Rossby waves.

The variations in the pycnocline are investigated by analyzing the isopycnal depths on an isopycnal surface. The depth is linearly interpolated with a density increment of 0.1 \({\sigma }_{\theta }\). Figure 6 shows the mean isopycnal depth and its standard deviation in the summer section. The maximum wintertime outcrop density is obtained from the winter 137°E section data and superimposed in Fig. 6. In the STMW region between 20° and 32°N, the standard deviation is relatively large over almost the entire layer below the wintertime outcrop density, indicating that the pycnocline changes considerably in the STMW region.

Fig. 6
figure 6

Climatology (contours) and standard deviation (color) of the depth of isopycnal surfaces in the summer 137°E section, plotted as a function of \({\sigma }_{\theta }\) and latitude. The standard deviation is calculated using the detrended, low-pass filtered time series at the grid points with a time series of more than 30 years. The thick line denotes the wintertime maximum outcrop density among the winter 137°E sections (color figure online)

The correlation between the STMW thickness and the isopycnal depth is calculated for every isopycnal surface at each latitude grid point in the STMW region between 20° and 32°N (Fig. 7a). In the region between 26° and 32°N where the thick STMW is present in the climatology (Fig. 3b), the correlation is generally negative and positive above and below the STMW layer at approximately \({\sigma }_{\theta }\)= 25.0–25.5 kg m−3, respectively, which signifies that the presence of the thick STMW coincides with the shoaling of the upper pycnocline and the deepening of the lower main pycnocline. The negative correlation in the upper pycnocline is significant at the 95% confidence limit and extends to the south around 20°N over the STMW region, with a slight drop in the correlation in the two STF bands at 21°N and 24°–26°N (Fig. 1b). In contrast, the positive correlation in the main pycnocline is low and confined to the north of 25°N with limited significance, except in the region north of 30°N. The pattern of the negative and positive correlation may be suggestive of the second mode baroclinic adjustment proposed by Xie et al. (2011), though it appears only north of 25°N in the northern part of the STMW region.

Fig. 7
figure 7

a Correlation coefficients (color) and regression coefficients (contours) of the isopycnal depth on the STMW thickness at each latitude in the summer 137°E section, plotted in the STMW region between 20° and 32°N. b The same as in (a) but those of the isopycnal depth on the main pycnocline depth calculated from the wind-driven Rossby wave model. The correlation and regression are calculated using the detrended, low-pass filtered time series. Crosses indicate a significant correlation at the 95% confidence limit. The thick line denotes the wintertime maximum outcrop density among the winter 137°E sections (color figure online)

The linear regression of the isopycnal depth is calculated using the STMW thickness as an explanatory variable (see contours in Fig. 7a). The regression coefficient is approximately 0.1–0.2 in the upper pycnocline above the STMW layer, indicating that the upper pycnocline changes in depth by approximately 10–20% of the change in the STMW thickness.

The relationship between the pycnocline and the STMW is recognizable in the time-vertical density section shown in Fig. 8. When the thick, low-PV STMW appears at depths of 200–350 m, the isopycnal surfaces move upward above the STMW with positive density anomalies near the surface. On the other hand, the deepening of the main pycnocline is not so obvious, but it may be identifiable below the low-PV STMW in the periods from 1980 to 1982 and from 1990 to 1997. In 1975 and from 2010 to 2015, the main pycnocline rises despite the presence of the low-PV STMW. The rise occurs over the entire depth above 600 m.

Fig. 8
figure 8

Yearly time series of \({\sigma }_{\theta }\)(contours), detrended \({\sigma }_{\theta }\) anomalies (color) and PV (blue contours of 1.5 and 1.8 \(\times\) 10–10 m−1 s−1) averaged between 26° and 29°N in the summer 137°E section. The low PV contours represent a thick, low-PV STMW. The thick black line denotes the wintertime MLD computed from the winter 137°E sections. All the plots are smoothed by the low-pass filter (color figure online)

The heaving of isopycnal surfaces is noticeable even at the depth where the mixed layer forms in the preceding winter as shown in Fig. 8. This feature is also confirmed by the negative correlation above the maximum winter outcrop density shown in Fig. 7a. In spring, the sea surface heat flux changes from cooling to heating, and the seasonal pycnocline, which consists of the upper pycnocline, forms, resulting in a rapid shoaling of the mixed layer (Qiu and Kelly 1993). The seasonal pycnocline grows and extends downward from spring to summer, steadily eroding the upper boundary of the STMW mainly through eddy-related processes such as diapycnal eddy diffusion (Qiu et al. 2006; Sukigara et al. 2011) . Because a thick STMW is transported continuously from the formation region on decadal timescales, the seasonal pycnocline could be prevented from developing more deeply during the thick STMW period.

The strong positive correlation below the STMW layer between 30° and 32°N might represent not only the adjustment to the STMW but also the influence of the Kuroshio on the STMW. As mentioned in Sect. 3.1, during the periods of the Kuroshio large meander, the Kuroshio shifts to the south around 30°N at 137°E, which shoals the main pycnocline and reduces the STMW thickness there.

Figure 7b shows the correlation between the observed change in the isopycnal depth and the pycnocline depth anomalies at 137°E computed from the wind-driven Rossby wave model (Sect. 2.2). The correlation is significantly positive over the entire pycnocline in the North Equatorial Current region south of 15°N, where wind forcing plays a key role in interannual-to-decadal changes in the upper ocean through Rossby wave dynamics (Qiu and Chen 2010b). In contrast, in the STMW region, the correlation is overall positive over the entire pycnocline but is not significant between 20° and 29°N. The depth variations of the upper pycnocline are related more to the STMW thickness (Fig. 7a) than to the main pycnocline depth predicted by the Rossby wave model. As for the main pycnocline layer, the correlation is overall not significant with either the STMW thickness (Fig. 7a) or the Rossby waves (Fig. 7b). The mechanisms of the variations in the main pycnocline need further investigation.

3.3 Variations in the STF and their cause

In this section, we investigate the relationship between variations in the STF and the STMW on isopycnal surfaces. Figure 9 shows the mean and standard deviation of the meridional density gradient, F on an isopycnal surface in the summer section. F is linearly interpolated with a density increment of 0.1 \({\sigma }_{\theta }\) for each year and is then averaged for the entire analysis period. As can be seen in Fig. 1b, the STF appears with the two bands of the positive \(F\) at 19°N and 25°N above the STMW layer around \({\sigma }_{\theta }\)= 25.0–25.5 kg m−3 (Fig. 9a). The standard deviation is large in the regions north of each STF band (Fig. 9b). We focus on these two large-variability regions and estimate the changes in the STF strength by averaging \(F\) vertically between the isopycnals at the bottom of the mixed layer and the top of the STMW (\({\sigma }_{\theta }\)= 24.8 kg m−3) and then by averaging it between 20° and 22°N and between 26° and 28°N (Fig. 10). The STF shows decadal-scale variability in the two regions. The autocorrelation function estimates a dominant period of approximately 8–15 years, which agrees well with that of the STMW thickness (Sect. 3.1).

Fig. 9
figure 9

a Climatology of the meridional potential density \({\sigma }_{\theta }\) gradients, \(F={\left(\partial {\sigma }_{\theta }/\partial y\right)}_{z}\), in the summer 137°E section, plotted as a function of \({\sigma }_{\theta }\) and latitude (color). b Standard deviation of the detrended, low-pass filtered time series of \(F\)(color). The contours in both panels are the climatology of \(F\). The thick line denotes the wintertime maximum outcrop density among the winter 137°E sections (color figure online)

Fig. 10
figure 10

Time series of a the STF strength anomalies, \({F}^{^{\prime}}\), in the summer 137°E section, which are obtained by averaging \(F\) in the vertical direction between the isopycnals at the bottom of the mixed layer and the top of the STMW (\({\sigma }_{\theta }\)= 24.8 kg m−3) and across latitudes between 20° and 22°N and then by subtracting the time mean value. Gray lines denote the yearly values, and thick lines denote the low-pass filtered values with a three-year running mean. Dashed lines are the same as the thick lines but are calculated using the time mean squared Brunt-Väisälä frequency and the time-varying isopycnal slope in Eq. (5) (see the text for details). b The same as in (a) but between 26° and 28°N. The correlation at the top denotes the correlation coefficient between the thick and dashed lines (color figure online)

Because the STF is associated with the northward shoaling of the upper pycnocline (e.g., Kobashi and Kubokawa 2012), we rewrite \(F\) in terms of the meridional slope of the isopycnal surfaces. Considering a small change in density \(\delta \rho\) from the reference density \({\rho }_{0}\) in the y, z plane, we transform the meridional density gradient from z to isopycnal coordinates as follows:

$$\left[ {\frac{{\left( {\rho_{0} + \delta \rho } \right) - \rho_{0} }}{\delta y}} \right]_{z} = \left[ {\frac{{\left( {\rho_{0} + \delta \rho } \right) - \rho_{0} }}{\delta z}} \right]_{y} \left( {\frac{\delta z}{{\delta y}}} \right)_{\rho } .$$
(3)

This transformation is described in some textbooks (e.g., Holton 2004). Taking the limits \(\delta y, \delta z\to 0\) yields.

$$\left( {\frac{\partial \rho }{{\partial y}}} \right)_{z} = - \left( {\frac{\partial \rho }{{\partial z}}} \right)_{y} \left( {\frac{\partial z}{{\partial y}}} \right)_{\rho } ,$$
(4)

\({\mathrm{where} \left(\partial z/\partial y\right)}_{\rho }\) on the right-hand side of the equation represents the slope of the isopycnal surfaces. Using \({\sigma }_{\theta }\) for \(\rho\), the depth of the isopycnal surfaces \(Z\left({\sigma }_{\theta }\right)\) and the Brunt–Väisälä frequency \(N={\left\{-\left(g/{\rho }_{0}\right)\left(\partial {\sigma }_{\theta }/\partial z\right)\right\}}^{1/2}\), Eq. (4) can be rewritten as.

$$F = \left( {\frac{{\partial \sigma_{\theta } }}{\partial y}} \right)_{z} = - \frac{{\rho_{0} }}{g}N^{2} \left( {\frac{{\partial Z\left( {\sigma_{\theta } } \right)}}{\partial y}} \right)_{{\sigma_{\theta } }} ,$$
(5)

where \(F\) is proportional to the product of the squared Brunt–Väisälä frequency and the isopycnal slope. We replace \({N}^{2}\) in Eq. (5) with the time mean value,\(\stackrel{-}{{N}^{2}}\), and calculate the STF strength in the same manner (Fig. 10). The time series are almost identical to the original one, which indicates that the change in the STF strength is closely associated with that in the isopycnal slope.

We estimate the change in the STF strength due to the STMW, \({F}_{S}^{^{\prime}}\), by using the \(\stackrel{-}{{N}^{2}}\) and the anomalies of isopycnal depth associated with the STMW thickness,\({Z}_{S}^{^{\prime}}\left({\sigma }_{\theta }\right)\), at each latitude point as follows:

$$F_{s}^{^{\prime}} = - \frac{{\rho_{0} }}{g}\overline{{N^{2} }} \left( {\frac{{\partial Z_{S}^{^{\prime}} \left( {\sigma_{\theta } } \right)}}{\partial y}} \right)_{{\sigma_{\theta } }} ,$$
(6)

where the prime denotes the anomaly from the mean. \({Z}_{S}^{^{\prime}}\left({\sigma }_{\theta }\right)\) is calculated from the linear regression of the isopycnal depths on the STMW thickness. For comparison, we also estimate the change in the STF strength due to wind-driven Rossby waves, \({F}_{w}^{^{\prime}}\), in the same manner, but by using the isopycnal depth anomalies regressed on the main pycnocline depth calculated from the Rossby wave model,\({Z}_{w}^{^{\prime}}\left({\sigma }_{\theta }\right)\) as follows:

$$F_{w}^{^{\prime}} = - \frac{{\rho_{0} }}{g}\overline{{N^{2} }} \left( {\frac{{\partial Z_{w}^{^{\prime}} \left( {\sigma_{\theta } } \right)}}{\partial y}} \right)_{{\sigma_{\theta } }} .$$
(7)

Figure 11 shows the correlation coefficients between \({F}^{^{\prime}}\) and \({F}_{S}^{^{\prime}}\) and between \({F}^{^{\prime}}\) and \({F}_{w}^{^{\prime}}\) on an isopycnal surface. The correlation between \({F}^{^{\prime}}\) and \({F}_{S}^{^{\prime}}\) is generally positive and significant in the upper pycnocline above the STMW layer around \({\sigma }_{\theta }\)= 25.0–25.5 kg m−3 north of 20°N, with the two high-correlation bands around 21°N and 28°N. These bands correspond roughly to those of the large STF variability shown in Fig. 9b. In contrast, the correlation between \({F}^{^{\prime}}\) and \({F}_{w}^{^{\prime}}\) is weak and not significant in the STF region. These results indicate the importance of the STMW in the STF variations.

Fig. 11
figure 11

Correlation coefficients between the STF strength anomalies of a \({F}^{^{\prime}}\) and \({F}_{S}^{^{\prime}}\) and b \({F}^{^{\prime}}\) and \({F}_{w}^{^{\prime}}\) in the summer 137°E section. a is plotted only in the STMW region between 20° and 32°N. The correlation is calculated using the detrended, low-pass filtered time series. Crosses denote a significant correlation at the 95% confidence limit. The thick line is the wintertime maximum outcrop density among the winter 137°E sections

The time series of \({F}^{^{\prime}}\) in the two high-correlation bands are plotted in Fig. 12, together with the contributions of the STMW thickness \({F}_{S}^{^{\prime}}\) and the Rossby wave \({F}_{w}^{^{\prime}}\). The contributions of \({F}_{S}^{^{\prime}}\) and \({F}_{w}^{^{\prime}}\) are calculated by linearly regressing \({F}^{^{\prime}}\) on \({F}_{S}^{^{\prime}}\) and \({F}^{^{\prime}}\) on \({F}_{w}^{^{\prime}}\), respectively. The variations in \({F}^{^{\prime}}\) agree well with those in \({F}_{S}^{^{\prime}}\) in each region, with significant positive correlation coefficients of 0.73 and 0.69 for the northern and southern regions, respectively, while the contribution of \({F}_{w}^{^{\prime}}\) to \({F}^{^{\prime}}\) is small, with a small amplitude in both regions.

Fig. 12
figure 12

a Low-pass filtered time series of the STF strength anomalies, \({F}^{^{\prime}}\)(thick lines), and the contribution of the STMW thickness \({F}_{S}^{^{\prime}}\)(dashed lines) and the wind-driven Rossby waves \({F}_{w}^{^{\prime}}\)(thin lines) in the summer 137°E section. These are the averages in the vertical direction between the isopycnals at the bottom of the mixed layer and the top of the STMW (\({\sigma }_{\theta }\)= 24.8 kg m−3) and across latitudes between 20° and 22°N.\({F}^{^{\prime}}\) is the same as the thick lines in Fig. 10. b The same as in (a) but between 26° and 28°N. The correlation at the top denotes the correlation coefficient between the thick and dashed lines (color figure online)

How does the STMW change the STF? The isopycnal depth of the STF,\({Z\left(\rho \right)}_{STF}\), is expressed as.

$$Z\left( \rho \right)_{STF} = Z\left( {\rho_{b} } \right) - Q,$$
(8)

where Z \(\left({\rho }_{b}\right)\) is the main pycnocline depth of the isopycnal surface \({\rho }_{b}\), and \(Q\) is the thickness of the layer between \({Z\left(\rho \right)}_{STF}\) and Z \(\left({\rho }_{b}\right)\). Subtracting the time mean and taking the meridional derivative yields.

$$\frac{{\partial Z^{\prime}\left( \rho \right)_{STF} }}{\partial y} = \frac{{\partial Z^{\prime}\left( {\rho_{b} } \right)}}{\partial y} - \frac{{\partial Q^{\prime}}}{\partial y}.$$
(9)

The left-hand term is the slope of the isopycnal surfaces of the STF. The first term on the right-hand side is the slope of the main pycnocline, and the second term is the meridional gradient of the thickness of the layer. The right-hand terms can be decomposed into the components associated with the STMW thickness and the others as follows:

$$\frac{{\partial Z^{\prime}\left( \rho \right)_{STF} }}{\partial y} = \frac{{\partial Z_{s}^{^{\prime}} \left( {\rho_{b} } \right)}}{\partial y} - \frac{{\partial Q_{s}^{^{\prime}} }}{\partial y} + R,$$
(10)

where \({Q}_{s}^{^{\prime}}\) is the STMW thickness anomaly,\({Z}_{s}^{^{\prime}}\left({\rho }_{b}\right)\) is the change in the main pycnocline depth associated with the STMW thickness and \(R\) is the other factors that include wind-driven Rossby waves. The first two terms on the right-hand side are associated with the STMW. The equation states that a thicker STMW can increase the meridional gradient of the STMW layer and intensify the slope of the STF to the south. A similar diagnostic equation was derived in terms of PV by Kobashi et al. (2006). We calculate the first two terms on the right-hand side of the equation at each latitude using \({\sigma }_{\theta }\)= 26.0 kg m−3 for \({\rho }_{b}\) and the STMW thickness. Figure 13 shows the correlation and regression between the left-hand term and the sum of the first two right-hand terms for every isopycnal surface above the top of the STMW (\({\sigma }_{\theta }\)= 24.8 kg m−3). The correlation is significantly positive, with the maxima around 20°–22°N and 27°–28°N, consistent with the results in Fig. 11a. This demonstrates that the STMW changes the STF strength.

Fig. 13
figure 13

Correlation coefficients (color) and regression coefficients (contours) of the meridional slope of the isopycnal surfaces on the sum of the meridional gradients of the main pycnocline associated with the STMW thickness and the meridional gradient of the STMW thickness (those of the left-hand term on the sum of the first two terms on the right-hand side in Eq. (10)). The correlation and regression are calculated using the detrended, low-pass filtered time series. Crosses denote a significant correlation at the 95% confidence limit. The thick line is the wintertime maximum outcrop density among the winter 137°E sections (color figure online)

The regression coefficients are smaller than unity in Fig. 13, which means that the change in the isopycnal slope of the STF is smaller than the expectation from that in the STMW thickness. In addition, because the correlation coefficients are approximately 0.6–0.7 in the STF region in Fig. 13, they are equivalent to the explained variance of approximately 40–50%. More than half of the variance remains to be understood. The correlation and regression coefficients decrease to the surface above the wintertime maximum outcrop density (Fig. 13). The seasonal upper pycnocline is affected by various processes such as heat and freshwater fluxes at the sea surface and horizontal advection of heat and freshwater by the mean circulation. These processes could alter the STMW-induced slope of the upper pycnocline and thus the STF.

4 Summary and discussion

Decadal-scale variations in the STMW thickness and their influence on the pycnocline are examined by analyzing the JMA 137°E repeat hydrographic observations for the period from 1972 to 2019, with a particular focus on the summer when the seasonal upper pycnocline develops above the STMW. The STMW appears between 20° and 32°N at 137°E, with the thickness varying on decadal timescales. The dominant period is approximately 9–15 years. The analysis of the Argo float observations after 2005 suggests that the change in the STMW thickness at 137°E originates in the wintertime mixed layer south of the Kuroshio Extension in the preceding year.

The pycnocline in the STMW region at 137°E is influenced by the STMW variations. Figure 14 is a schematic summary of the results. The presence of thick STMW shoals the upper pycnocline and is occasionally concurrent with the deepening of the lower main pycnocline. The change is robust in the upper pycnocline where the heaving of isopycnal surfaces is accompanied by density anomalies up near the surface. Because the heaving is obvious even at the depth where the mixed layer forms in the preceding winter, the continuous advection of the STMW from the upstream could change the local development of the seasonal upper pycnocline.

Fig. 14
figure 14

Schematic meridional density section illustrating the baroclinic adjustment of the pycnocline to the change in the STMW thickness and resultant change in the STF on decadal timescales

The thick STMW enhances the northward shoaling of the upper pycnocline and intensifies the STF. The STF exhibits large variability in the regions north of the two STF bands. The dominant timescales are consistent with those of the STMW thickness. The STF variations are explained by the STMW-induced change in the upper pycnocline slope on decadal timescales, consistent with previous findings from numerical models. The STMW accounts for approximately 40–50% of the STF variance. More than half of the variance remains to be understood. Further analysis is necessary for a full understanding of the mechanisms of the STF.

We take a brief look at the spatial distribution of the change in the pycnocline associated with the STMW using MOAA-GPV monthly temperature and salinity profiles from 2005 to 2019. During this period, the STMW exhibits a notable decadal change (Fig. 3). Figure 15 indicates the correlation coefficients between the yearly time series of the June–August mean STMW thickness and depths of the isopycnal surfaces of \({\sigma }_{\theta }\)= 23.3 and 26.0 kg m−3, which represent the upper and lower pycnocline depths, respectively. These two isopycnal surfaces correspond to the peaks in the histogram of the vertical maximum of the squared Brunt–Väisälä frequency plotted against \({\sigma }_{\theta }\) in the STMW region of 24.5°–32.5°N, 134.5°–170.5°E in June to August (not shown).

Fig. 15
figure 15

a Correlation coefficients between yearly time series of the June-to-August mean STMW thickness and depth of the upper pycnocline with the isopycnal surface \({\sigma }_{\theta }\)= 23.3 kg m−3, plotted in the STMW region with a climatological-mean thickness of more than 20 m. Thick contours denote 20 m and 120 m of the climatological mean STMW thickness from Fig. 2. The white area within the thick contours corresponds to the outcrop of the upper pycnocline. b The same as in (a) but the depth of the main pycnocline with the isopycnal surface \({\sigma }_{\theta }\)= 26.0 kg m−3. Dots denote a significant correlation at the 95% confidence limit. The correlation is calculated using the detrended, low-pass filtered time series. All the plots are based on the MOAA-GPV temperature and salinity profiles from 2005 to 2019

The correlation in the upper pycnocline is overall negative, with a significant correlation along the southern and eastern flank of the thick STMW (delineated by 120-m contours in Fig. 15a). On the other hand, the main pycnocline exhibits a significant positive correlation to the east of 140°E, except in some regions along approximately 30°N east of 160°E and in the southwest region south of 28°N and west of 150°E (Fig. 15b). The pattern of the negative correlation in the upper pycnocline and the positive one in the main pycnocline may be similar to that observed at the JMA 137°E section (Fig. 7a). Interestingly, the locations of the significant correlation are hardly overlapped between the upper and main pycnoclines. The feature could be different from the second mode baroclinic adjustment.

The negative correlation in the upper pycnocline is confined to the north of 26°N in the MOAA GPV (Fig. 15a), while it extends significantly to the south at 20°N in the JMA 137°E section (Fig. 7a). Compared to the STMW in the JMA section (Fig. 3b), the STMW is thin and absent to the south of 24°N at 137°E in the MOAA GPV (Fig. 2). This is probably due to the strong spatial smoothing adopted in the MOAA GPV (Hosoda et al. 2008). The STMW may not be well resolved in the MOAA GPV.

Regarding the decrease in the STMW from 2006 to 2009 (Fig. 4a), several studies analyzing Argo float observations emphasize basin-scale wind forcing (Qiu and Chen 2013; Cerovečki and Giglio 2016; Cerovečki et al. 2019). The positive wind curl forcing shoals the main pycnocline and results in vertical shrinking of all isopycnals above the main pycnocline and thus a decrease in the STMW thickness. In that case, the change in the STMW results from that in the pycnocline. Figure 16 shows the time-longitude diagram of monthly wind-driven pycnocline depth anomalies from the Rossby wave model, together with the STMW thickness anomalies observed in the summer 137°E section. The Rossby wave model expects the arrival of the shoaling pycnocline anomalies in 2009 (Fig. 16b), concurrent with the decrease in the STMW (Fig. 16a). This result might be consistent with those of previous studies. However, such a relationship is not present in the other period, and there is no correlation between the Rossby waves and the STMW thickness at 137°E in the period from 1972 to 2019 (coefficients less than 0.2).

Fig. 16
figure 16

a Time series of the STMW thickness anomalies averaged between 26° and 29°N in the summer 137°E section (the same as in Fig. 4a). b Time-longitude plot of monthly main pycnocline depth anomalies averaged for a zonal band of 26°–29°N, derived from the wind-driven Rossby wave model. Positive values denote deep anomalies. All the plots are smoothed by the low-pass filter

The STMW region at 137°E corresponds to the so-called “reemergence” area of winter SST anomalies (Hanawa and Sugimoto 2004). Part of the STMW is re-entrained into the deepening surface mixed layer in the autumn and can alter the SST (Alexander and Deser 1995; Hanawa and Sugimoto 2004). The reemergence process determines oceanic thermal inertia and affects climate variability on interannual to decadal timescales (Newman et al. 2016). The heaving of the upper pycnocline might influence the reemergence process through entrainment at the base of the mixed layer. These analyses are to be performed in the future.

The water properties of the STMW change on decadal and longer timescales (e.g., Yasuda and Hanawa 1997; Sugimoto et al. 2017; Oka et al. 2017; Ogata and Nonaka 2020). They are also expected to alter subtropical pycnoclines in the downstream region. Comprehensive analyses are needed to understand the role of the STMW in changing the pycnocline.