The decadal to multidecadal mixed layer variability is investigated in a region south of the Kuroshio Extension (130°E–180°, 25°–35°N), an area where the North Pacific subtropical mode water forms, during 1948–2012. By analyzing the mixed layer heat budget with different observational and reanalysis data, here we show that the decadal to multidecadal variability of the mixed layer temperature and mixed layer depth is covaried with the Atlantic multidecadal oscillation (AMO), instead of the Pacific decadal oscillation (PDO). The mixed layer temperature has strong decadal to multidecadal variability, being warm before 1970 and after 1990 (AMO positive phase) and cold during 1970–90 (AMO negative phase), and so does the mixed layer depth. The dominant process for the mixed layer temperature decadal to multidecadal variability is the Ekman advection, which is controlled by the zonal wind changes related to the AMO. The net heat flux into the ocean surface Qnet acts as a damping term and it is mainly from the effect of latent heat flux and partially from sensible heat flux. While the wind as well as mixed layer temperature decadal changes related to the PDO are weak in the western Pacific Ocean. Our finding proposes the possible influence of the AMO on the northwestern Pacific Ocean mixed layer variability, and could be a potential predictor for the decadal to multidecadal climate variability in the western Pacific Ocean.
The upper-ocean mixed layer plays an important role in climate variability through connecting the atmosphere and ocean. For example, even the seasonal cycle of mixed layer depth (hereafter hm) can induce “winter-to-winter” memory and thus persistence of sea surface temperature (SST) anomalies via the “re-emergence” mechanism, through the changes in the maximum depth of the mixed layer from one winter to the next (Alexander and Deser 1995; Deser et al. 2003; Sugimoto and Hanawa 2005; Byju et al. 2018). Most of the water masses properties in the global ocean are controlled by the air–sea interaction and turbulent mixing in the surface mixed layer and by subduction process down to the ventilated thermocline layer (Pedlosky 1996). Hence, the variability of the hm in the water-mass formation regions has been studied in association with the analysis of water masses, such as mode water (e.g., Speer and Forget 2013).
In the North Pacific subtropical gyre, the deepest mixed layer is located in a region south of the Kuroshio Extension (KE; black contours in Fig. 1a) where the North Pacific Subtropical Mode Water (STMW) is formed. It is believed that the change of STMW has important impacts on midlatitude climate variability through re-emergence of wintertime SST anomalies (e.g., Alexander et al. 1999; Hanawa and Sugimoto 2004; Sugimoto and Hanawa 2005) and regulates biogeochemical process, for example, via the nutrient cycling in the oligotrophic subtropical gyres (e.g., Oka et al. 2015) and oceanic uptake of atmospheric CO2 (Bates et al. 2002). Numerous authors (e.g., Qiu and Chen 2005, 2006; Qiu et al. 2007; Oka et al. 2011, 2012) have investigated the importance of hm variability on the STMW in the region south of the KE. One of the major factors controlling the hm changes is the intensity of the upper-ocean stratification (e.g., Qiu and Chen 2006), as strong upper-ocean stratification is unfavorable to the mixed layer development. It should be noticed that hm is usually defined by the mixed layer temperature (hereafter as Tm; Fig. 1b) south of the KE (e.g., Qiu and Kelly 1993; Qiu and Chen 2006; Sugimoto and Hanawa 2010), since the temperature is more important than the salinity in the density here (Johnson et al. 2012).
In the North Pacific, the most prominent interannual–decadal variation is known to be the Pacific decadal oscillation (PDO), which is normally characterized by the “PDO index” (Mantua et al. 1997). The PDO index is represented as the intensity change of Aleutian low (AL) and associated with the Pacific–North American teleconnection pattern. The PDO would generate the wind stress curl anomaly in the central North Pacific, and then the ocean response (e.g., SST anomaly and thermocline depth anomaly) would propagate westward at the speed of first-mode baroclinic Rossby waves. It roughly takes 3–5 years to reach the western Pacific and then influence the hm as well as the STMW in the region south of the KE (e.g., Miller et al. 1994; Nakamura et al. 1997; Qiu and Chen 2005; Qiu 2003; Taguchi et al. 2007; Oka and Qiu 2012).
However, a recent study carried by Wu et al. (2020) shows that the decadal to multidecadal variability of the hm and STMW is controlled by the AMO rather than PDO in the northwestern Pacific Ocean. The STMW mean temperature varies with the AMO index. They hypothesize that this decadal to multidecadal variation of the hm and STMW originates from the Tm variability south of the KE, caused by the AMO, and the surface signal would propagate into the subsurface ocean via subduction process. However, the above hypothesis needs to be further tested, in a quantitively way, since there are several processes could affect the change of Tm and then the hm. Another related question here is why the AMO but not the PDO controls the mixed layer decadal to multidecadal variability south of the KE (Fig. 1).
To address the above questions, in the present study, we quantitively estimate the contribution of each process in the mixed layer heat budget analysis and explain how the AMO-related Ekman heat transport controls the Tm decadal to multidecadal variability. The presentation is organized as follows. Descriptions of the mixed layer heat budget governing equation and data are given in section 2, in which seasonal variations of all the terms included are also presented. The contribution of each term on the Tm decadal to multidecadal variation is shown in section 3. In section 4, the relationship between the AMO related zonal wind anomaly and Ekman advection term, which is the major factor controlling the Tm decadal to multidecadal variability, is investigated. We also discuss the relative importance of AMO and PDO on the Tm variability and highlight the dominant role of AMO in this section.
2. Mixed layer heat budget
a. Governing equation
Following Qiu and Kelly (1993), the vertically integrated equation for the mixed layer temperature can be expressed by
where hm denotes the mixed layer depth and Tm the mixed layer temperature, respectively. The other terms are as follows: Um is the horizonal transport in the mixed layer, ΔT is the temperature difference between the mixed layer and the layer just below the mixed layer, we is the entrainment velocity, Ah is the horizontal eddy diffusivity, ∇ = (∂/∂x, ∂/∂y) is the horizontal gradient operator, and ρ0 and c are the reference density and specific heat of seawater, which are taken to be 1026 kg m−3 and 4000 J kg−1 K−1.
The net surface heat flux Qnet includes shortwave and longwave radiation and latent and the sensible heat fluxes (Qnet is positive into the ocean). The downward heat penetration at depth z, q(z) is given by
where q(−hm) is the downward heat penetration at the bottom of the mixed layer; q(0) is the shortwave radiative flux at the sea surface; R is a separation constant; and γ1 and γ2 are the attenuation length scales depending on the water properties. In this study, we used R = 0.77, γ1 = 1.5 m and γ2 = 14 m, which correspond to type II surface water (Paulson and Simpson 1977). In the classification of Jerlov (1968), there are five types of ocean water (I, IA, IB, II, and III) and observational data indicate that the water type in the western Pacific is type II.
where we is computed as the forward difference. Notice that during the detraining period (when the mixed layer is in the shoaling phase; i.e., ΔH/Δt < 0), we = 0.
The horizontal transport Um can be separated into an Ekman component UEK, which is driven by the wind stress τ, and a geostrophic component Ug, which is computed from the thermal wind equations. Therefore, the equation for Um is written as
where f is the Coriolis parameter, ρ is the seawater density, and τx and τy are the surface wind stress in the x and y direction, respectively. Here, we assume a level of no motion at 1500-m depth (H = 1500 m).
b. Data and methods
The south of the KE front (as shown in Fig. 1a with dashed box region) is the analysis region in this study, where the Tm (Fig. 1b) variability along with the front here is important for the hm variability through influencing the upper-ocean stratification intensity. According to Sugimoto and Hanawa (2010), the upper-ocean stratification intensity is defined as the mean value of the vertical temperature gradient (°C m−1) calculated from the surface to 200 m. Here 200 m is chosen to cover the change of seasonal mixed layer adequately, but to exclude the influence from the main thermocline depth change. To compute all the terms of the equations in section 2a, data fields of hm, Tm, Qnet, q(z), and τ are needed. The respective data sources are listed below.
Three observational datasets that cover the 1948–2012 period (longer than other available datasets), referred to as Ishii data, EN4 data, and IAP data, are used in the following calculations. Ishii data are from the research led by Ishii (Ishii et al. 2006); the second set is the EN4 data (version EN.4.2.1-analyses-g10) produced by Met Office Hadley Center (Good et al. 2013); and the last is the product by Institute of Atmospheric Physics, Chinese Academy of Sciences (called IAP data in this study; Cheng et al. 2017). All the datasets are global gridded monthly objective analyses of in situ observations and derived from quality-controlled ocean data that ingest all types of ocean profiling instruments, for example, expendable bathythermographs (XBT), conductivity–temperature–depth (CTD) measurements from research ships and Argo floats. Ishii data cover the global oceans with horizontal resolution of 1° × 1° and 24 levels in the upper 1500 m for the period from 1945 until 2012 (https://rda.ucar.edu/datasets/ds285.3/). Since 1900, EN4 data offer monthly temperature and salinity information in a 1° × 1° horizontal grid with 42 depth levels, with a maximum depth of 5350 m (https://www.metoffice.gov.uk/hadobs/en4/). The IAP data provide monthly three-dimensional subsurface temperature and salinity data with 41 levels up to 2000 m in the same horizontal resolution as Ishii and EN4 datasets (1° × 1°) from 1940 to the present (http://126.96.36.199/cheng/). Detailed information on input data sources used in the compilation and associated quantification of sampling errors and statistical errors is provided in Ishii et al. (2006), Good et al. (2013), and Cheng et al. (2017). To further support our results, the ocean reanalysis product, Simple Ocean Data Assimilation (SODA) version 2.2.4 (Carton and Giese 2008), is used in this study with more AMO cycles. It is an extend ocean reanalysis product using the POP ocean model and the NOAA/National Centers for Environmental Prediction Twentieth Century Reconstruction version 2 winds (Giese and Ray 2011). It provides an estimate of the ocean on a 0.5° × 0.5° grid with 40 vertical levels for the period of 1871–2010 (https://climatedataguide.ucar.edu/climate-data/soda-simple-ocean-data-assimilation).
The hm is derived from the Ishii (or EN4, IAP, SODA) monthly temperature field and defined to be the upper-ocean layer, whose depth-averaged temperature is 1°C higher than the water temperature just below (Td), namely,
The same definition for hm was used by Qiu and Kelly (1993). Besides, we calculated hm by density definition: the base of the mixed layer is the depth at which density is 0.125 kg m3 higher than the surface density (Monterey and Levitus 1997), as shown in Fig. S1 in the online supplemental material. Both the value and pattern of temperature-defined hm are similar to the hm calculated by density. The seasonal cycle of these two hm is also very similar, with the hm defined by temperature being slightly larger than that by density (Fig. S1c). Since the temperature observation is more accurate with larger profiles than the salinity observation, we use the hm defined by temperature hereafter. The temperature averaged within hm is taken as the Tm.
Daily surface fluxes and the surface wind stress during 1948–2012 are from the National Centers for Environmental Prediction (NCEP) and the National Center for Atmospheric Research (NCAR) global reanalysis (simply called the NCEP data in this study; Kalnay et al. 1996; https://www.esrl.noaa.gov/psd/). It covers the global oceans with a T62 Gaussian grid. The net surface heat flux Qnet (positive into the ocean) is the sum of the shortwave radiation, the longwave radiation, the sensible heat flux, and the latent heat flux from the NCEP data. The shortwave radiative flux at the sea surface q(0) from NCEP is given to compute the q(z) in Eq. (2). The wind stress τ (τx, τy) used in this study is directly from NCEP data. The Ssalto/Duacs delayed-time multimission maps of the absolute dynamic topography is used in this study (Fig. 1a) to show the climatological KE position, with 0.25° × 0.25° spatial resolution from 1993 to 2018 (available at http://www.aviso.oceanobs.com).
Two other datasets are used in the following analysis to double-check the accuracy of the flux data. One is the NOAA/Cooperative Institute for Research in Environmental Sciences (CIRES) Twentieth Century Reanalysis version 2 (NOAA 20CRv2; Compo et al. 2011) on a 2° (longitude) × 2° (latitude) grid. NOAA 20CRv2 data cover the period from 1851 to 2014 and are available at https://www.esrl.noaa.gov. The other dataset the is ECMWF twentieth-century reanalysis (ERA-20C), which provides global atmospheric data for the period 1900–2010 with approximately 125 km (spectral truncation T159) horizontal resolution and is available at http://www.ecmwf.int (Poli et al. 2016). The time periods of 1948–2012 and 1948–2010 are selected from NOAA 20CRv2 and ERA-20C, respectively.
The PDO and AMO [sometimes also referred to as the Atlantic multidecadal variability (AMV)] indices represent climate modes in the North Pacific and North Atlantic oceans, respectively. Monthly time series of the PDO index for the 1948–2012 period are obtained from the website of the Joint Institute for the Study of the Atmosphere and Ocean, University of Washington (http://jisao.washington.edu/pdo/PDO.latest). The PDO index is defined as the leading principal component of North Pacific monthly SST variability (poleward of 20°N). The AMO index is based upon the area-weight average anomalies of SST over the North Atlantic basin (0°–65°N, 80°W–0°) with respect to the 1961–90 climatology (https://climatedataguide.ucar.edu/climate-data/atlantic-multi-decadal-oscillation-AMO-index-station-based).
In this study, annual means are constructed from monthly means by averaging the data from January to December, winter means from January to March, and summer means from July to September. The statistical significance of the linear regression coefficient and correlation between two time series is based on a two-tailed Student’s t test in which the degrees of freedom are estimated by considering the lag-1 autocorrelation, calculated following Bretherton et al. (1999).
c. Seasonal variation on the upper-ocean heat balance
To quantitatively calculate the contribution of each term in determining the upper-ocean thermal structures, it is helpful to rewrite Eq. (1) as follows:
Considering that the Ekman and the geostrophic flows (uEK and ug) have disparate effects on the Tm changes, we have separated the advection term into the Ekman and geostrophic components in Eq. (9). Here the Ekman current uEK is obtained by scaling the Ekman transport in Eq. (4) with an effective depth, which is given as 26.5 ± 3 m for the tropical Pacific (Ralph and Niiler 1999). In this study, the effective depth is equivalent to the hm if the latter is not deeper than 30 m; otherwise it is set to 30 m.
The climatological map of each term in the governing equation of the Tm [Eq. (9)], including the meridional and zonal Tm gradient, Ekman velocity, geostrophic velocity, net heat flux (Qnet; positive into the ocean), entrainment velocity, and horizonal eddy diffusion term from 1948 to 2012 is shown in Fig. 2. On the annual-mean basis, the spatial distribution of the meridional temperature gradient (Fig. 2a) is larger than the zonal temperature gradient (Fig. 2d) in the southern KE region. For the Ekman velocity field, the meridional velocity is stronger than the zonal one (Figs. 2b,e), since the zonal wind is much stronger here than the meridional wind (Fig. 3a in Wu et al. 2018). On the contrary, the climatological meridional geostrophic current is weaker than its zonal component (Figs. 2c,f), which is consistent with the large-scale wind and pressure field. The above mean patterns of the temperature gradient and oceanic currents (both Ekman and geostrophic velocity) result in a larger contribution of the Ekman advection term (−uEK ⋅ ∇Tm) than the geostrophic advection term (−ug ⋅ ∇Tm) to the temperature tendency term ∂Tm/∂t. The spatial pattern of the Qnet averaged in 1948–2012 is shown in Fig. 2g and it agrees well with previous study (e.g., Hsiung 1985). As we defined, the entrainment velocity we is always positive (Fig. 2h). South of the KE, we is relative weak and contributes little to the Tm seasonal variability. The diffusive term [Ah(hm∇2Tm − ΔT∇2hm)] shows a large value and a dipole pattern in the KE region, with positive (warming) south of the KE and negative (cooling) north of the KE (Fig. 2i).
Based on Eq. (8), Tm and hm are calculated using the monthly mean 1° × 1° temperature data from different sources (the Ishii, EN4, and IAP datasets). The climatological Tm and hm values averaged in the analysis area are shown in Figs. 3a and 3b (black line for Ishii data, blue line for EN4 data, and yellow line for IAP data). Both time series have well-defined annual cycles. The Tm reaches its maximum, about 25°, in summer and decreases steadily in the subsequent months until the following March. Accompanying the Tm decreasing, there is a gradual increase in the standard deviation of Tm. On the contrary, hm is shallowest, about 20 m, in summer and deepens in the succeeding months until the next March, reaching 200 m. For the spatial patterns of seasonal Tm and hm, we plot their climatological mean in winter and summer as displayed in Figs. 3c–f. It is clear that the spatial pattern of the climatology of Tm (hm) is dominated by the summer (winter) mean and less influenced by the winter (summer) mean. Examining the term balance for the Tm change [Eq. (1), Fig. 3g] reveals that the temporal seasonal change of Tm is more sensitive to terms such as the surface heat flux Qnet and the horizontal advection than to the cold water entrainment through the bottom of the mixed layer.
In Fig. 4, we represent the zonally averaged term balances of Eq. (9) as a function of time and latitude, as in Qiu and Kelly (1993). The seasonal temperature tendency term ∂Tm/∂t displays a pattern very similar to the surface thermal forcing term Qnet − q(−hm)/ρ0chm (Figs. 4a,b). It is not surprising that the surface thermal forcing determines the seasonal cycle of the upper-ocean water temperature, because of the strong surface thermal forcing in the KE region. It is obvious that the annually integrated contribution of the surface thermal forcing to the Tm is positive; that is, the atmosphere warms the mixed layer. This warming effect of surface thermal forcing on the Tm is nearly counterbalanced by the other terms in Eq. (9). We should notice that a perfect annual equilibrium of the Tm would require the integrated ∂Tm/∂t value to be zero, but it cannot be exactly zero because of the interannual changes in Tm associated with the El Niño–Southern Oscillation event as discussed in Qiu and Kelly (1993). There is always a cold effect from the vertical entrainment term −ΔTwe/hm as shown in Fig. 4c, larger in the warm season and smaller in the cold season. Horizontally, the position of annual-mean zero wind stress curl line is located along 30°N (see Fig. 3a in Wu et al. 2018), with westerly wind north of 30°N inducing southward cold Ekman transport and resulting in temperature flux divergence (−uEK ⋅ ∇Tm < 0), and vice versa south of 30°N (Fig. 4d). The temperature advection associated with geostrophic current is not coherent in the meridional direction. In areas north of the KE (latitude > 33°N), −ug ⋅ ∇Tm is generally positive, reflecting the northward heat transfer, and vice versa south of the KE (Fig. 4e). Finally, Fig. 4f indicates that the local fluctuations of the geostrophic advection term are partially balanced by the horizontal eddy diffusion term, with negative north of 33°N and positive south of 33°N.
3. Domination of the Ekman advection on the Tm decadal to multidecadal variation
As we analyzed in the last section about the seasonality of the upper-ocean heat balance, all the time series and patterns are consistent with previous study (e.g., Qiu and Kelly 1993). It is confirmed that the data processing in this study is all correct, and we can further analyze the decadal to multidecadal variability of the mixed layer heat budget in this section.
When we focus on the decadal to multidecadal time scales, the PDO shows no significant impact on the Tm variation (Fig. 5a). But the Tm decadal to multidecadal variability is significantly related with the AMO (Fig. 5b) during 1948–2012 (the correlation coefficient between the 7-yr running-averaged Tm and AMO index is 0.89). Even we use the other data source (from EN4 and IAP in Fig. 5c) and longer ocean reanalysis data (SODA; green line in Fig. 5c) back to 1871, the good correlation between the Tm and the AMO is still robust (r = 0.70 for EN4 data during 1940–2018; r = 0.84 for IAP data during 1940–2018; r = 0.80 for SODA during 1871–2010). It confirms that the basin-averaged Tm well varies with the AMO (Fig. 5), no matter the data source or length.
To better understand the contribution of each term on the Tm anomaly on decadal to multidecadal time scales, we decompose each field into an annual mean (denoted by an overbar) and the departure from the mean (denoted by a prime). If the contributions from the mean values of anomalies are neglected, the equation for the Tm anomalies is written as
Term is the rate of the Tm anomaly change, term is the change of the Tm by the surface thermal forcing anomaly, term represents advection of the mean Tm by anomalous currents, term denotes advection of the Tm anomaly by mean currents, term indicates advection of the Tm anomaly by currents anomaly, term is the entrainment of the Tm through the bottom of the mixed layer, and term is horizontal diffusion. Equation (10) represents a lower-order dynamical description of variations of Tm. The advantage of using this equation is that the key balance terms can all be calculated from existing datasets (see section 2b).
Similar to Fig. 4, the zonally averaged (130°E–180°, as shown in Fig. 1) Tm anomaly (Fig. 6a) and the contribution from each term (Figs. 6b–f) are plotted as a function of time (7-yr running-averaged time series during 1948–2012) and latitude. We can find that the zonally averaged Tm anomaly shows clear decadal to multidecadal variability, undergoing a warm phase before 1970 and after 1990, and a cold phase of 1970–90. This “high–low–high” feature is consistent with the basin-averaged Tm variation as shown in Fig. 5. Not only the horizonal pattern, but also the vertical structure of the water temperature changes represents a clear high–low–high feature, with propagation of the temperature anomaly signal from surface into the subsurface ocean (as shown in Figs. 3c,d in Wu et al. 2020). According to Eq. (10), we calculate the contribution of the Ekman advection term, geostrophic advection term, net surface heat flux term, entrainment velocity term, and horizontal eddy diffusion term on the Tm changes. It is interesting that the Tm anomaly induced by Ekman advection term also has an obvious decadal to multidecadal high–low–high feature (Fig. 6b) as the Tm, whereas the other terms do not.
Instead, the contribution from surface thermal forcing (Qnet term) is nearly opposite to the Tm changes with a low–high–low feature (Fig. 6d). It should be emphasized that the surface thermal forcing is dominant in the seasonal cycle, as we discussed in the last section, whereas it is not true when we focus on the decadal to multidecadal time scales. In fact, the SST variability here is not determined by the air–sea heat flux or by the Tm changes. The air–sea interaction south of the KE reflects ocean driving atmosphere in the decadal to multidecadal time scale as discussed in Wu et al. (2019). For example, an upward heat release increases while Tm anomalies are positive during the warm phase after 1990 (Fig. 6d), and vice versa during the cold phase for 1970–90. Figures 6c and 6f indicate that the contribution from geostrophic advection term and horizontal eddy diffusion term can somewhat balance with each other, which is similar to the seasonal cycle (Figs. 4e,f). The contribution of entrainment velocity to the Tm changes is little, with strong interannual variability but no obvious decadal to multidecadal variability (Fig. 6e). Figure 7 is the domain-averaged heat budget, which further displays the dominant role of Ekman advection in determining the decadal to multidecadal variability of the Tm, while the net heat flux into the ocean surface Qnet acts as a damping term.
From the above analysis, we can conclude that the decadal to multidecadal variability of the Tm is dominated by the Ekman advection. To better understand the physical process, we further examine the decadal changes between two contrasting periods of the Tm anomaly, the 1970–80 and 2000–10 periods (Fig. 7). During these two periods, the AMO is stable for both the negative (1970–80) and positive (2000–10) phase, while the PDO just changes its phase and would have very little impact (the gray shadings in Fig. 5), so that we can remove the PDO effect on the Tm decadal changes in the following analysis. First, we compare the decadal change of Tm (1970–80 and 2000–10, latter minus former) between the observation and integrated term in Eq. (10), as shown in Figs. 8a and 8b. They show exactly the same warm pattern, which confirms that the calculation processing is correct. Then we compute the contribution from each term on the Tm decadal change as shown in Figs. 8c–f, including the total advection term, Qnet term, entrainment velocity term, and diffusion term. Not surprisingly, the total advection term is responsible for the major warm anomalies of the Tm (Fig. 8c), although the locations of the warming do not match well in details. The different warming patterns between Figs. 8a and 8c may be first caused by the inconsistency of the datasets during the calculation (i.e., temperature from the Ishii data but wind stress from the NCEP data). Second, the contributions from other terms (i.e., the Qnet term, entrainment velocity term, and diffusion term as shown in Figs. 8d–f) also account for the decadal change of Tm, although not as much as advection term. On the contrary, both the surface net heat flux term and entrainment velocity term dampen the warm anomaly with a cooling effect (Figs. 8d,e). The diffusion term does not have any spatial pattern with patches of warm and cold anomalies (Fig. 8f), and the sum of it in the whole region is also negligible. Here, the aforementioned damping effect of the Qnet term as exhibited in Figs. 6–8 is double-checked with other flux datasets, including NOAA 20CRv2 (Compo et al. 2011) and ERA-20C (Poli et al. 2016). It is clear that the results from NCEP (Figs. 6–8) are in agreement with NOAA 20CRv2 and ERA-20C as shown in Figs. S2 and S3.
As in Fig. 6, we separate the total advection term into an Ekman advection term and a geostrophic advection term (Figs. 9a,b). The Ekman advection term contributes a lot to the anomalous warm pattern of the Tm and then dominates the decadal change of the Tm. Based on Eq. (11), the Ekman advection term can be decomposed into the , , and terms, which indicate the Tm decadal change caused by the anomalous Ekman velocity, the anomalous temperature gradient, and the anomalous Ekman velocity and temperature gradient, respectively (Figs. 9c–e). It is clear that the term mostly controls the Ekman advection term as shown in Fig. 9c, with a warm pattern similar to Fig. 9a, and is much more important than the other two terms ( and in Figs. 9d,e). The heat transport by the mean Ekman velocity ( term; Fig. 9d) is patchy and also smaller than the anomalous Ekman velocity term (Fig. 9c). The advection by two anomalous terms is very small and can be neglected (Fig. 9e). It should be noted that the Ekman velocity anomaly in the term is closely connected with the wind stress changes as shown in Eqs. (5) and (6). We check the wind stress decadal change between 1970–80 and 2000–10 (latter minus former). As expected, it is the easterly wind anomaly that causes the northward Ekman heat transport, warming the south of the KE and controlling the Tm decadal to multidecadal variation (Fig. 9f).
As analyzed above, the zonal wind anomaly is the key point for the Ekman advection term controlling the Tm heat budget. The question remaining here is: What is the relationship and mechanism between the zonal wind anomaly and the AMO? Wu et al. (2019) pointed out the similar connection between wind field in the Oyashio Extension region and the AMO. In the past 37 years (1982–2018), it is the AMO-induced zonal wind anomaly influencing the meridional shift of the Oyashio Extension fronts. In this study, we will further check the possible link between the wind field changes and the AMO south of the KE. Figure 10a shows the zonal wind speed anomaly (black line) averaged in the study area, which is well correlated with the AMO index (red line; r = −0.82). During the AMO positive phase (before 1970 or after 1990), there is an easterly wind anomaly, and vice versa for the AMO negative phase. The AMO impact on the zonal wind anomalies are better shown in Figs. 10b and 10c. In the AMO negative phase (Fig. 10b), there is a broad westerly wind anomaly, from the subtropical to subpolar North Pacific Ocean, and enhanced in the western Pacific Ocean. It is similar for the AMO positive phase with a broad easterly wind anomaly (Fig. 10c). Due to the zonal wind anomalies, the surface Ekman transport will carry cold or warm water meridionally, changing the Tm and controlling the mixed layer decadal to multidecadal variability south of the KE region.
In this section, the decadal to multidecadal variation of the Tm is analyzed by detailed heat budget. It is confirmed that the Ekman advection, which is controlled by the wind stress changes associated with the AMO (Fig. 10), dominates the Tm variability.
4. Discussion and conclusions
Through the above detailed heat budget analysis, we found that the time-integrated Ekman advective flux convergence contributes the most low-frequency variability of the Tm in the past 60 years, while the Qnet forcing acts as a damping term (e.g., Fig. 7). Here, we will further decompose Qnet into all its components, net shortwave radiation, sensible and latent heat fluxes, and upward and downward longwave radiation, to evaluate their contributions to the decadal to multidecadal variability of the Tm anomalies. First, we check the climatology mean of each term as shown in Fig. S4. Notice that the values of upward (Fig. S4d) and downward longwave radiation flux (Fig. S4e) are roughly 3 times larger than the others, while the sum of the upward and downward longwave radiation flux (i.e., net longwave radiation flux; Fig. S4f) is small. It means that the upward and downward longwave radiation almost compensate each other, which is consistent with previous studies (e.g., Yamanouchi and Kawaguchi 1984; Kotsuki et al. 2015). Second, we decompose the Qnet term based on Eq. (10) as follows:
where QULWRF denotes the upward longwave radiation flux, QDLWRF the downward longwave radiation flux, QSHF the sensible heat flux, QLHF the latent heat flux, and QNSWRF the net shortwave radiation flux.
It is found that the damping effect of the Qnet term on the Tm anomalies on decadal to multidecadal time scales is mainly from the QLHF term (red line in Fig. 11) and partially from QSHF term (green line in Fig. 11). The sum of the upward and downward longwave radiation fluxes (i.e., net longwave radiation flux QNLWRF; orange line in Fig. S5), exhibits clear warming effect on the Tm without multidecadal variation, which may be related to the global warming. The QULWRF, QDLWRF, QSHF, and QNSWRF terms all counteract each other. The turbulent exchange of heat between the ocean and atmosphere is due to sensible heat (temperature) and latent heat (moisture) fluxes (Large and Pond 1982). We found that not only the QLHF term but also the QSHF term (green line in Fig. 11) have similar variations to the Qnet term on decadal to multidecadal time scales: negative before 1970 and after 1980, and positive during 1970–80. The longwave radiation is not only controlled by the SST (more precisely, skin temperature) but also is sensitive to the variation of other meteorological parameters, such as clouds, water vapor amount, and so on (Yamanouchi and Kawaguchi 1984). Therefore, it is understandable that the QULWRF term does not function as a pure damping term on the decadal to multidecadal variability of the Tm.
The Tm is considered as the key index to show the changes of the mixed layer in the past 60–70 years. Actually, we found that not only Tm but also hm vary with the AMO. Figure 12 gives the hm (black line) variability during 1948–2012. It is coherent with the AMO index (red line), with deep hm during negative AMO phase and shallow hm during AMO positive phase. It should be noticed that the hm in this study is defined by the Tm according to Eq. (8) and hence the change of hm is strongly related to the change of Tm, namely, the upper-ocean stratification. Furthermore, we examined the variability of the salinity including the time-varying and spatial structure, and its impact on the mixed layer density (or hm), on the decadal to multidecadal time scale, comparing to temperature as shown in Fig. S6. On the one hand, it can be readily seen that the time-varying mixed layer density (black line in Fig. S6a) fluctuates coherently with Tm (red line in Fig. S6a) very well over the entire period (r = −0.96). The mixed layer salinity (yellow line in Fig. S6a) exhibits an obvious decreasing trend from 1948 to 2012 without any decadal to multidecadal feature. On the other hand, the spatial structure of mixed layer density decadal change (Fig. S6b) is also consistent with the Tm decadal change (Fig. 8a). The 2D correlation between the mixed layer density and Tm decadal change pattern is high (−0.87), while there is no significant correlation with the mixed layer salinity (Fig. S6c). The above results indicate that the mixed layer density (or hm) is controlled by Tm, from both the time-varying variation and spatial structure perspective, on the decadal to multidecadal time scales.
As we discussed above, the Tm anomaly induced by Ekman heat transport will transfer into the subsurface ocean, change the upper-ocean stratification, and further affect the hm variability. It is shown in Fig. 12 that the upper-ocean stratification (blue line) also varies with the AMO. When the AMO is in its negative phase, the upper-ocean stratification is weak and hm deepens, and vice versa for the AMO positive phase. In addition, the 200-yr outputs from a PI-Control simulation [see Wu et al. (2020) for details], with more AMO cycles, also support our findings. The evolution of the 200-yr hm south of the KE varies with the modeled AMO significantly (r = −0.63; Fig. S7), which further confirms that the good correlation between hm and the AMO is not a coincidence.
Considering the well-known seasonality of the Tm and hm, we investigate the responses of seasonal Tm and hm to the AMO and PDO forcing. Figure 13 demonstrates the 7-yr low-pass filtered Tm and hm regressed upon the similarly filtered normalized AMO and PDO in winter and summer during 1948–2012. Consistent with the previous study (Figs. 2a,b in Wu et al. 2020), south of the KE (130°E–180°, 25°–35°N), the decadal to multidecadal Tm and hm pattern associated with the PDO shows nearly no signal (Figs. 13e–h), while displaying the largest anomalies associated with the AMO (Figs. 13a–d), which indicates that a positive AMO is the potential cause of a warmer Tm and shallower hm in both winter and summer. On the other hand, the PDO-induced Tm and hm anomalies are concentrated in central eastern Pacific, resembling the classic PDO pattern (Mantua et al. 1997). Specifically, the variance of Tm (hm) anomalies in summer (winter) is larger than those in winter (summer), and the patterns correspond to their seasonal climatology mean (Figs. 3a–f). For example, the largest hm anomalies associated with AMO in winter are located in the area of 140°–160°E, 30°–35°N (Fig. 13b), similar to the deepest mean hm south of the KE (Fig. 3d). Our above results suggest that the effect of AMO on the Tm, which is caused by local wind forcing via Ekman heat transport south of the KE, is important, especially in summer. It is worth emphasizing that the winter wind and its variability are usually larger than those in summer (Fig. 4 in Wu et al. 2019). The Ekman heat transport is more important in summer than that in winter due to the seasonality of the hm (Fig. 3). Because the Tm change by Ekman heat transport is inversely proportional to hm [Eq. (6) in Wu et al. 2019], the thin hm in summer makes the Ekman heat transport much more predominant as compared to the winter season.
It is interesting that the AMO rather than the PDO dominates the decadal to multidecadal variability of the mixed layer in the south of the KE. Figure 14 shows the SST and wind stress regressed upon the AMO and PDO index, respectively. South of the KE, there is strong easterly wind anomaly induced by the AMO with obvious SST anomaly (Fig. 14a). But the PDO-induced wind and SST changes are mainly confined in the central and eastern Pacific Ocean, with very weak loading in the western Pacific Ocean (Fig. 14b). It indicates that the decadal to multidecadal variability of zonal wind anomaly south of the KE (Fig. 10a) is controlled by the AMO rather than the PDO, and hence will result in the mixed layer variation. It is worth noting that, on an interannual scale, the correlation between Tm (or zonal wind anomaly) and the AMO is not so good (thin lines in Figs. 5 and 10a), which means that on the different time scales, the physical mechanism is different. As shown by previous studies, the PDO dominates the state of the KE on interannual to decadal time scales (e.g., Qiu and Chen 2005, 2010), via oceanic Rossby wave propagation and eddy–mean flow interaction. In this study, when we focus on the decadal to multidecadal variability, large-scale wind forcing and related Ekman advection become important. The reason why the dominant mechanism is different on different time scales is an interesting topic and needs to be further addressed.
The whole physical process can be described as in Fig. 15. When the AMO in its positive phase, warmer SSTs in the whole North Atlantic Ocean (the red circle over North Atlantic in Fig. 15a) may influence the western North Pacific through the atmospheric teleconnections, inducing the easterly wind anomaly south of the KE. This easterly wind anomaly would further generate the northward Ekman heat transport and increases the Tm. The increased Tm could even transfer into the subsurface ocean, enhancing the upper-ocean stratification, and finally shallow the hm, and vice versa for the AMO negative phase (Figs. 15a and 12). But for the PDO, the wind anomaly is weak and mainly in the meridional direction to the south of the KE (Fig. 15b). It should be emphasized that the zonal Ekman transport, caused by meridional wind, has little impact on the Tm variability, since the SST front is nearly zonal here (Fig. 1). So, the hm does not change under the PDO forcing.
It should be noted that the aforementioned “atmospheric teleconnections” from the Atlantic to the subtropical North Pacific (the dashed arrows in Fig. 15a) have not been well or fully explained. Sun et al. (2017) demonstrated that the SST variability in the tropical western Pacific is largely explained by the AMO. They suggest that the Atlantic Ocean acts as a key pacemaker for the western Pacific decadal climate variability through the cross-basin teleconnection. Further, more and more researchers have noticed that the Atlantic Ocean has important effects on the tropical Pacific Ocean (e.g., Li et al. 2013, 2016; Chikamoto et al. 2016; Levine et al. 2017; Zhang et al. 2017), subtropical Pacific Ocean (Zhang and Delworth 2007; Wu et al. 2020), and subpolar Pacific Ocean (Gong et al. 2020). To some extent, the above studies provide the possible links between the Atlantic and Pacific through so-called atmospheric teleconnections, although they focus on the regional Pacific Ocean. To our knowledge, no literature presents a well-explained atmospheric mechanism for the teleconnection from the Atlantic to the entire western Pacific yet; this is important work for the next step.
In conclusion, we report the findings that the decadal to multidecadal variation of the Tm (and the hm) south of the KE is dominated by the Ekman advection (anomalous Ekman velocity, the term) through the detailed heat budget, and that it is controlled by the AMO-induced zonal wind changes. The net heat flux into the ocean surface Qnet acts as a damping term and results mainly from the effect of latent heat flux and partially from sensible heat flux. The Tm has strong decadal to multidecadal variability, being warm before 1970 and after 1990, and cold during 1970–90. Unlike in previous studies, the changes of Tm are not due to the PDO, but vary with the AMO on decadal to multidecadal time scales. The salient contribution of this study is that we propose the possible influence of the AMO on the mixed layer variability in the region south of the KE. Due to the close connection between the mixed layer variation and the STMW formation, the AMO could be an important driver for the decadal climate and could provide memory for the prediction in the Pacific Ocean (Wu et al. 2020). In addition, it is interesting to note that the fish catch variations in the Kuroshio region, dropping before the 1970s and after the 1990s (Fig. 2e in Tian et al. 2014), are coincident with the Tm variation and the AMO index, but not with the PDO. From the above considerations, we believe that our findings will be a step toward the possible predictor for both the climate and ecosystem decadal to multidecadal variability in the western Pacific Ocean.
This study benefited from discussion with Steven J. Lentz. Three anonymous reviewers provided useful and constructive comments in revising the manuscript. Xiaopei Lin is supported by the China’s national key research and development projects (2016YFA0601803) and the National Natural Science Foundation of China (41925025 and U1606402). Baolan Wu is supported by the China Scholarship Council (201806330010). Lisan Yu thanks NOAA for support for her study on climate change and variability.