Gravity Core CF1 (water depth 1,180 m; 127.43°E, 28.42°N) (Figure 1), with a length of 351 cm, was retrieved from the western slope of the middle Okinawa Trough during a cruise in 2012. The core was sliced into 1-cm-thick subsamples after the cruise. The lithology of Core CF1 consists of homogeneous gray silty clay. No obvious depositional hiatus or turbidite layers were found within Core CF1. The planktonic foraminiferal species Neogloboquadrina dutertrei from the >150 μm size fraction of eight layers were picked up for accelerator mass spectrometry radiocarbon (AMS ¹⁴C) dating at Beta Analytic Inc.(Florida, United States) (Figure 2).

Considering the top segment of Core CF1 was damaged and lost during the sampling process, the age model of Core CF1 between 8 and 348 cm (Figure 2A) was constructed using the recently published age-depth modeling routine “Undatable” (Lougheed and Obrochta, 2019). The deterministic routines of Undatable, with a positive sedimentation rate prior and bootstrapping, result in median and mean age-depth models with age error estimates of non-normally 68% and 95% percentiles (Lougheed and Obrochta, 2019). Here, the AMS ¹⁴C ages were calibrated to a calendar age before 1950 CE (cal BP) using the Marine20 calibration curve (Heaton et al., 2020) and the ΔR value of 29 ± 18 a (Table 1). Default values for bootstrapping percentage and sedimentation rate uncertainty were set to 20% and 0.1 cm yr⁻¹, respectively. The sedimentation rate of Core CF1 varied from 45 cm/kyr to 190 cm/kyr, with an average value of ∼90 cm/kyr.

Grain-size analysis was conducted using a Mastersizer-2000 laser particle-size analyzer at Qingdao Institute of Marine Geology, China Geological Survey, with a measurement range of 0.02–2,000 μm and a size resolution of 0.01ϕ. The measuring error was within 3%. Before the grain-size analyses, the samples were pretreated with 10% H₂O₂ and 0.5 mol/L HCl for 24 h to remove organic matter and biogenic carbonate, respectively. The high-resolution analysis provided a resolution of ∼12 years on average, with a lower resolution of ∼22 years for the interval between 2,400 and 3,460 years BP.

Varimax-rotated Principal Component Analysis (V-PCA) is a statistical procedure that uses orthogonal transformation to convert a set of possibly correlated variables into a set of linearly uncorrelated variables called principal components. This method allowed us to separate out orthogonal modesand independent grain-size spectra from the grain size matrix that are related to potential input functions and sensitive to specific transport mechanisms (Darby et al., 2009; Hu et al., 2012; Yi et al., 2012; Zheng et al., 2014; Zheng et al., 2016). It has been used to successfully identify the transport mechanisms of the Yellow Sea (Hu et al., 2012) and the Okinawa Trough (Zheng et al., 2014; Zheng et al., 2016). In this study, V-PCA was applied to temporal variation of Core CF1 grain size spectrum with the input grain size matrix ranging from 0.41 to 707 μm. The mode of each extracted grain size component is defined as the grain size sets with largest factor loading, which is most representative for the grain size spectra (Figure 3).

The ensemble empirical mode decomposition (EEMD) method is based on the noise-assisted data analysis, which takes the average value of multiple measurements and obtains the ensemble means of corresponding intrinsic model functions (IMFs) as the results (Wu and Huang, 2009). EEMD is implemented through a sifting process that uses only local extrema. A complete sifting process stops when the residue becomes a monotonic function from which no more IMFs can be extracted. The total number of IMFs of a data set is lose to log₂N (N is the number of total data points). Detail process and explanation can be found in the MATLAB EEMD code program instructions (Wu and Huang, 2009). Before performing the EEMD, the PC1 loading of Core CF1 is linearly interpolated at a 10-years interval. Eight IMFs are generated from our data by the EEMD method and spectral analyses is applied to determine the periodicities and periodic stabilities using the PAST software (details can be found in https://palaeo-electronica.org) (Figure 4). We aim to quantify the relative contributions of both the long-term trend and centennial-scale oscillations. Considering the sample resolution used in this study, the IMF1 and IMF2 components with periodicity less than 100 years are thus considered as decadal-scale noise (Figure 5A). At the same time, the IMF3, IMF4, and IMF5 components exhibit centennial-scale cyclic oscillations with a total variance of 39% (Figure 5B), while the IMF6, IMF7, and IMF8 components are regarded as millennial variability or long-term trend, contributing altogether 29% to the total variance (Figure 5C).

Previous studies suggested the Okinawa Trough has multiple potential terrigenous sediment sources, mainly including the large rivers in China (Yellow and Yangtze Rivers) and small mountain rivers in Taiwan and Kyushu Islands, all of which display large temporal-spatial variations influenced by sea-level fluctuations, oceanic circulations, and East Asian monsoons over the last glacial-interglacial cycles (Dou et al., 2010; Dou et al., 2012; Wang et al., 2015; Dou et al., 2016; Beny et al., 2018; Li et al., 2019; Xu et al., 2019; Hu et al., 2020). Based on the results of previous studies, the sediment source-to-sink process in the ECS-Okinawa Trough since the LGM can be broadly summarized as follows: During the glacial periods marked by the low sea level (e.g., about 135 m below the present sea level during the LGM), the paleo-channels of Yellow and Yangtze Rivers were located on the YS and ECS shelves and may have directly entered the northern and middle regions of the Okinawa Trough, respectively. In contrast, Taiwan-derived sediments are confined to the southern region of the Okinawa Trough due to the weak or even absent KC during the glacial periods. Subsequently, sea level gradually rose from -135 m to -40 m during the deglaciation period (22–10 ka), resulting in the YS and ECS shelves being flooded by seawater along with the paleo-mouths of Yellow and Yangtze River which retreated quickly towards the YS shelf. During the process of sea level rise, terrigenous sediment supply became trapped on the shelf, although more intense continental erosion occurred resulting from the strengthened East Asian Summer Monsoon (EASM). Meanwhile, deglaciation transgression caused some reworked materials to be released from the ECS shelf into the Okinawa Trough, forming a broad tidal sand ridge covering the ECS middle-outer shelf. Conversely, the mainstream of the KC gradually deflected to the west of the Ryukyu Islands due to the rapid sea-level rise and then delivered the Taiwan Rivers sediment northward during this period. Provenance proxies indicated a significant change from a dominance of the paleo-Yangtze or paleo-Yellow River to an increased influence from Taiwanese rivers after 9–10 ka or ∼7 ka, associated with the establishment of a “water barrier” of strong KC intensity and a modern circulation pattern since then. Combining the results of previous studies, it can be found that the sea level has stabilized or changed very little and the ocean circulation system of ECS has been formed since the Late Holocene. At that time, the sediment transport dynamics affecting the middle Okinawa Trough are mainly controlled by changes in the strength of the KC.

For Core CF1, results of the V-PCA analyses on the grain-size matrix are shown in Figure 2A. Three principal grain-size components (PC1: 69%, PC2: 16%, and PC3: 10%) were identified for Core CF1, which totally account for 95% of the variance. The first mode (PC1) has a broad negative peak at 4–16 μm and a positive peak at 32–148 μm (Figure 2A). The second mode (PC2) has a positive peak at 0.4–3.2 μm and a trough around 18–26 μm (Figure 2A). The third mode (PC3) has a strongly positive peak at 297–707 μm (Figure 2A). Two nearby cores A7 and OKI02 gave the same modes with a consistent grain size structure, and such consistency strongly indicates that there must be identifiable mechanisms to account for these specific modes (Zheng et al., 2014; Zheng et al., 2016). As suggested by Zheng et al. (2014, 2016), three PCs of the Okinawa Trough cores (A7 and OKI02) were interpreted as KC silt (PC1), bottom nepheloid clay (PC2), and turbidity sand (PC3), respectively. Sr-Nd isotopic compositions of 5–18 μm fractions from Core AF2/OF3 were further compared with the potential endmembers (the Yangtze River, Yellow River, and rivers of Taiwan Island) to assess the provenance, which clearly indicated that the 5–18 μm fractions originated from rivers of Taiwan Island and are transported by, and sensitive to, the strength of the KC (Zheng et al., 2016). Specially, sea level has stabilized or changed very little and the ocean circulation system of ECS has been formed since the Middle Holocene, and at that time, the sediment transport dynamics affecting the middle Okinawa Trough are mainly controlled by changes in the strength of the KC (Zheng et al., 2014; Zheng et al., 2016). These lines of evidence support the interpretation of PC1 loading of Core CF1 as a dynamic proxy for the KC intensity.

To quantify the relative contributions of long-term trend and centennial-scale oscillations of the KC, we performed a noise-assisted data analysis EEMD on the PC1 loading of Core CF1 after linearly interpolated at a 10-years interval. Based on the results of the EEMD, we combined the IMF6, IMF7 and IMF8 components (IMF6+7 + 8) to reveal the millennial variability or long-term trend of KC intensity (Figure 6A). Our reconstruction of the KC intensity contradicts recent inferences based on low-resolution magnetic parameters and paired organic paleothermometers (Li et al., 2019; Li et al., 2020), but is consistent with several other lines of evidence. For example, abundance of P. obliquiloculatata (Jian et al., 2000; Lin et al., 2006; Xiang et al., 2007), sediment mercury (Hg) enrichment factor (Lim et al., 2017), and sediment provenance changes (Jiang et al., 2019; Xu et al., 2019), as well as SST variations in the South Yellow Sea (Wang et al., 2011; Jia et al., 2019) (Figure 6B). Within the age uncertainty estimates, our IMF6+7 + 8 record of Core CF1 is broadly consistent with the abovementioned proxy evidence and suggests a slow-down of the KC during the period of 4.6–2.0 ka, followed by quickly enhanced KC and attaining the highest level at approximately 0.5 ka (Figure 6A). Additionally, climate records from the Taiwan Island (Selvaraj et al., 2007; Wang et al., 2015; Ding et al., 2020), Cheju Island (Park et al., 2016; Park, 2017) and Northeastern Asia (Hong et al., 2005; Zhao et al., 2021) reflect gradual precession-driven cooling/drying since the Early-Middle Holocene with rapid climate amelioration since ∼2.0 ka. Their similarity to the KC intensity evolution suggests the KC transports huge amounts of latent heat and water vapor to the atmosphere along its path and plays an important role in modulating surrounding climate change.

What, then, explains Late Holocene KC evolution? Firstly, although a few uncertainties remain with respect to the Holocene ENSO activities, most of the records support a strong relationship between ENSO variability and the Equatorial Pacific mean state and suggests more frequent El Niño events during the Late Holocene (Figures 6C,D) (Moy et al., 2002; Conroy et al., 2008; Koutavas and Joanides, 2012; Sadekov et al., 2013; Gill et al., 2016; Barr et al., 2019). Enhanced El Niño activities during the Late Holocene, characterized by weakened Walker circulation and northernmost position of NEC bifurcation latitude, is not conducive to KC reinforcement at that time.

Secondly, the variability of the KC along the ECS margin at interdecadal scales is more closely related to PDO than to ENSO, and KC generally decelerated during the negative PDO phase, and vice versa (Andres et al., 2009; Wu, 2013; Wu et al., 2019). Recent study shows that PDO-like pattern shifts from a positive to a negative phase at ∼2 ka based on ΔSST_E-W between Cores MD01-2,412 and ODP1019C (Figure 6E) (Chen et al., 2021). Similar results are also found in the PDO record from the Santa Barbara Basin over the last 2,700 years (Beaufort and Grelaud, 2017). Commonly, negative PDO phase since ∼2 ka, characterized by decreased negative wind stress curl (WSC), should decelerate the southward Sverdrup flow in the subtropical gyre and thus weaken return flow (i.e., KC), which conflicts with the KC records from the Okinawa Trough.

Finally, the East Asian monsoon systems, including the southeasterly EASM and the northwesterly East Asian Winter Monsoon (EAWM), are another important factor influencing the variation of KC intensity (Jian et al., 2000; Zheng et al., 2016). High-resolution reconstruction of EASM and EAWM from three loess sections in the Loess Plateau display continuously weakened EASM (Figure 6F) and persistently strengthened EAWM during the Late Holocene (Figure 6G) (Kang et al., 2020). These situations would trigger less negative or even positive WSC over the North Pacific subtropics, suppressing the KC intensity, which also contradicts the enhanced KC since ∼2.0 ka.

Overall, low-latitude processes (e.g., East Asian monsoons and/or changes in ENSO/PDO-like phases) should be a secondary factor or negligible for the long-term trend of the KC intensity on millennial timescale during the Late Holocene.

As discussed earlier, forcings other than low-latitude processes are required to explain the long-term trend of the KC intensity during the Late Holocene. Indeed, the KC is driven by a basin-scale negative WSC over the North Pacific, depending on the combined effects of subtropical easterly and mid-latitude westerly winds (Seager and Simpson, 2016). Enhanced negative WSC over the North Pacific forces stronger southward Sverdrup flow in the inner ocean, which is compensated by strengthened northward flow of the KC at the western boundary (Seager and Simpson, 2016). Accordingly, we hypothesize that the mid-latitude Westerly Jet (WJ) and its associated WSC in the North Pacific play a dominant role in the long-term trend of KC intensity since the Late Holocene. Recent studies have highlighted teleconnections between North Pacific atmospheric and oceanic circulations over the glacial-interglacial cycles, with equatorward (poleward) and intensified (weaker) WJ during the glacial (interglacial) periods, coupling with the expansion and strengthening of circulation gyres during cold periods and vice versa (Gray et al., 2018; Gray et al., 2020; Abell et al., 2021).

Moreover, dust records from the Tarim Basin (Han et al., 2019), the Cheju Island (Lim and Matsumoto, 2006, 2008), and the Japan Sea (Nagashima et al., 2013) suggested strengthening and southward shift of WJ occurred during the Late Holocene (Figure 6H), which is further confirmed by Holocene spatiotemporal precipitation patterns (Figure 6I) and climate model simulations (Chen et al., 2019; Herzschuh et al., 2019; Routson et al., 2019; Zhou et al., 2020). The forcing mechanism may be that seasonal difference between the winter and summer insolation across the latitudes results in gradually increased pole-to-equator thermal gradient during the Late Holocene (Chen et al., 2019; Routson et al., 2019). Thus, we provide direct evidence for the predominance of accelerating and southward shift of mid-latitude WJ, rather than low-latitude processes, causing more negative WSC over the North Pacific and thus stronger KC since ∼2 ka.

Superimposed on the long-term trend of the KC is significant centennial-scale variability, which displays a larger amplitude and a longer anomaly duration after ∼2.5 ka (Figure 5B). We found a periodicity of ∼1,800-years (at the 95% significance level), ∼700-years (at the 90% significance level) and ∼500-years (at the 95% significance level) of the original PC1 record from Core CF1 by the PAST software (Figure 7). However, it seems inappropriate to determine the 1,800-years cycle as one of the credible climate oscillations due to their insufficient temporal length of the original PC1 records. We therefore consider only the 500–700-year cycles because of underlying solar forcing. Moreover, we also performed wavelet power spectrum analyses on the IMF3+4 + 5 record and reflected statistically significant centennial periodicities on approximately 500-years and 700-years (>90% confidence level) (Figure 8). Similar periodicities are also found in the solar activity records and exhibit strong power after ∼2.5 ka (Wanner et al., 2008; Steinhilber et al., 2012). This further suggests a link between the KC variability and solar activity on a centennial-scale, with reduced solar activity corresponding with weak KC intensity (Figure 9). Not all KC weakening events covary with sunspot number (SN) on a centennial-scale for the Late Holocene, which may stem from age model uncertainties and sampling resolution.

Previous studies highlight that solar variability (500–700-year cycles), amplified by low-latitude oceanic-atmospheric interactions (e.g., ENSO), could have played an important role in regulating centennial-scales climate variations (Liu et al., 2014; Zhu et al., 2017; Xu et al., 2019; Xu et al., 2020). Model simulations suggest that the ENSO system acts as a mediator of solar activities on the climate system’s low-latitude heat engine (Emile-Geay et al., 2007). This has been further confirmed by paleoclimate reconstructions, which show that the tropical Pacific mean states shift in response to the solar activities of the last 1,500 years (Mann et al., 2009), for the early-middle Holocene (Marchitto et al., 2010), as well as the entire Holocene (Ersek et al., 2012). In addition, solar irradiance variations can also trigger changes in the strength and direction of WJ across the North Atlantic via a “top-down” mechanism (Ineson et al., 2011), with the southward migration and enhanced intensity of WJ responding to lower solar irradiance (Olsen et al., 2012; Wirth et al., 2013; Lan et al., 2020). Nonetheless, if some direct (though amplified) solar forcing of the WJ was the dominant control on the KC intensity, we would expect enhancement of KC intensity correlated with reduced solar activity, conflicting with the fact observed (Figure 3). Taken together, we can conclude that, since the Late Holocene, weak solar irradiation makes the tropical Pacific to be an El Niño-like state, associated with weakened easterly trade winds and attenuated EASM, which overwhelms the enhanced intensity of mid-latitude WJ and contributes to the decreased KC intensity on centennial-scale.