Muted orbital-scale monsoon variability over the Korean Peninsula

A recent study identified an east-west dipole pattern in Asian Summer Monsoon precipitation and in its oxygen-isotope composition variability in response to orbital forcing (Wen et al., 2024). The Korean Peninsula (KP) is situated near the nodal line of this dipole and should therefore exhibit muted precessional variability in both. So far, this conjecture has not been tested using paleoclimatic data. Here, we present speleothem oxygen isotope (δ¹⁸O) records from the KP, which support the notion of suppressed orbital-scale hydroclimate variability. Conducting a transient model simulation with the isotope-enabled Community Earth System Model (iCESM1.2) covering the past 130 000 years, along with tagging experiments for low and high insolation conditions, we show that, on precessional timescales, isotopic contributions from oceanic and continental moisture sources compensate each other over the KP, resulting in only a weak regional signal in precipitation δ¹⁸O. Based on iCESM1.2 simulations, we further demonstrate that deuterium excess (d-excess) variability over the KP would still capture the moisture source region's zonal seesaw response to precessional forcing. We propose that, given the region's strong sensitivity to changes in moisture sources, paleo-water d-excess reconstructions from East Asia speleothems could provide valuable additional constraints on the drivers of regional monsoon systems. This study offers new insights into the spatiotemporal variability of Pan-Asian hydroclimates and its links to changes in moisture source.

1 Introduction

Speleothem (secondary carbonate deposits in caves) records, leveraged by radiometrically constrained chronologies, have been a valuable resource for reconstructing paleo-hydroclimate variability in tropical and subtropical regions (Braun et al., 2019; Cai et al., 2015; Cheng et al., 2016; Gasse et al., 2011; Jo et al., 2017; Maher and Thompson, 2012). Such work is predicated on stable oxygen isotope (δ¹⁸O) measurements of speleothem calcite (δ¹⁸O_c), which ultimately derives from atmospheric moisture. When speleothem CaCO₃ precipitates under near-equilibrium conditions, δ¹⁸O_c encodes cave air temperature and drip water δ¹⁸O. Drip water typically reflects the δ¹⁸O of precipitation (δ¹⁸O_p), integrating processes such as Rayleigh fractionation along water-vapor pathways and, in some regions, the amount of regional rainfall (Lachniet, 2009). In regions where large changes in δ¹⁸O_p are expected with respect to the effect of temperature, δ¹⁸O_c is often considered as an indicator of atmospheric variability. As such, δ¹⁸O_c-records from Chinese caves (e.g., Cheng et al., 2009; Wang et al., 2001) have provided key insights into seasonal and centennial-to-orbital-scale variability in the Asian monsoon system. Speleothems that grew within the two sub-monsoon branches of this system, the Indian Summer Monsoon (ISM) and the East Asian Summer Monsoon (EASM), generally show a close correspondence in their δ¹⁸O_cvariability (Cai et al., 2015; Cheng et al., 2019; Kathayat et al., 2016), at least on orbital timescales. The relationship was initially interpreted as an illustration of large-scale in-phase variations in atmospheric processes across both monsoons (Cheng et al., 2009, 2016; Kathayat et al., 2016; Wang et al., 2001). However, unlike the ISM, the EASM is characterized by more complex spatial patterns in precipitation and in δ¹⁸O_p, even under present-day climate conditions (Chiang et al., 2017; Wang et al., 2008; Zhao et al., 2012).

Our current understanding is that the EASM precipitation pattern is mainly driven by the position of the subtropical front, which responds to Northern Hemispheric Summer Insolation (NHSI). The frontal system, established through interactions between the Westerly jet and the Western North Pacific Subtropical High (WNPSH), separates moist tropical air masses and drier, colder extratropical air masses. The front is a primary rain-producing feature (Meiyu-Baiu rainband, or Changma for the Korean Summer Monsoon) in the EASM (Chen and Chang, 1980; Ding, 2004; Wang et al., 2008), and it controls the seasonality (Hu et al., 2019; Liu et al., 2022; Wang et al., 2008), rainfall distribution, intensity, and extent of the moisture advection from tropical oceans (Lin et al., 2024; Volonté et al., 2022).

Speleothem-based studies, complemented by climate simulations, have reported a connection between the EASM's characteristics (intensity and pattern) and frontal movements across different timescales and under extreme climatic conditions (Chiang et al., 2015; Liu et al., 2014; Volonté et al., 2022; Wang et al., 2008; Wu et al., 2018). For instance, Liu et al. (2014) demonstrate that frontal interactions with synoptic weather lead to heavier rainfall in northern China during high NHSI (boreal summer perihelion conditions). Later studies have reported wetter conditions in central-eastern China and spatially opposite rainfall signals across Asia during the last deglaciation (Termination I), in association with shifts in the position of the Westerlies (He et al., 2021; Zhang et al., 2018). A recent study by Xue et al. (2025) proposed that WNPSH intensity played a key role in shaping hydrological conditions across central China over various timescales.

On precessional timescales, a significant correlation exists between variability in Asian speleothem δ¹⁸O_c records and NHSI, indicating that the latter primarily controls the large-scale circulation of the Asian summer monsoon and its associated moisture transport. Chinese speleothem δ¹⁸O_c has commonly been interpreted as reflecting changes in large-scale monsoon circulation, integrated moisture transport, and upstream rainout processes, rather than local precipitation amount alone (Cai et al., 2015; Cheng et al., 2016; Kathayat et al., 2016; Wang et al., 2001). The interpretation of Chinese speleothem δ¹⁸O_c as a proxy for local precipitation amount remains debated because independent precipitation-sensitive proxies from East Asia, such as loess magnetic susceptibility and reconstructed seawater δ¹⁸O (Beck et al., 2018; Clemens et al., 2018; Sun et al., 2006; Zheng et al., 2022), show only a weak correspondence with the precessional cycle. Cheng et al. (2021) addressed the discrepancies and complexities of hydroclimate proxies in Eastern Asia, explaining that loess, marine, and cave records are complementary and preferentially capture specific aspects of Asian monsoon dynamics. This interpretation is further supported by Li et al. (2024), who found that the ratio of dithionite−citrate−bicarbonate extractable iron to total iron in Chinese loess deposits potentially records precipitation cycles on precessional timescales. However, uncertainties in interpreting the magnetic-susceptibility-based reconstructions remain, as they might reflect a smoothed precipitation signal or be dominated by other forcings (Cheng et al., 2021; Sun and Huang, 2006).

Climate model simulations generally show enhanced Asian summer monsoon precipitation under stronger NHSI (Bosmans et al., 2018), although the magnitude of the simulated response may differ from that inferred from paleoclimate evidence (Lee et al., 2019). Isotope-enabled modeling offers an opportunity to further evaluate the roles of atmospheric, land, and oceanic processes in a changing climate (i.e., transient simulations). They can help interpret precipitation-based proxies by tracing the stable isotopologues of water spatiotemporally (Battisti et al., 2014; Brady et al., 2019; Cheng et al., 2021). Such modeling studies have also shown that other processes may dominate on δ¹⁸O_p more than the Amount effect, including atmospheric circulation, moisture pathways, and upstream rainout (Liu et al., 2014; Maher and Thompson, 2012; Pausata et al., 2011). Subsequent studies have examined the complex relationships among orbital-scale EASM δ¹⁸O variability, NHSI, rainy seasons, and regional circulation patterns (Dai et al., 2021; Liu et al., 2022; Yang et al., 2024). Though the major features of simulated climate-isotopic signals across southern and eastern Asia agree with Chinese speleothem isotope records, some models suggest a pronounced out-of-phase precipitation response between the EASM rainband and the rest of the Asian monsoon on orbital timescales (Bosmans et al., 2018; Hu et al., 2019). Some isotope-enabled climate models also reveal large-scale differences in the response of δ¹⁸O_p to precessional forcing across Asia, with weaker anomalies over the eastern domain than over the western domain (Battisti et al., 2014; Tabor et al., 2018). A recent study elucidated the mechanisms by which the precession-driven insolation forcing reorganizes the large-scale Asian summer monsoon circulation, including the strength of the southwesterly monsoon flow, the position of the upper-level westerly jet, and the position and strength of WNPSH, thereby altering regional precipitation and the relative contributions of remote and nearby moisture sources to δ¹⁸O_p (Wen et al., 2024). These processes produce a grand dipole pattern of precipitation and δ¹⁸O_p variability across the Asian monsoon region, characterized by an antiphase precipitation pattern between South Asia and Japan. This pattern suggests that precipitation variability on precessional timescales may be weaker over nodal regions such as the Korean Peninsula (KP).

The KP lies at the boundary between northeast Asia and the Pacific, and moisture transport through the WNPSH strongly controls summer precipitation variability over the region (Ha et al., 2012; Wang and Lin, 2002). Under present-day conditions, precipitation isotope compositions over these regions are influenced by seasonal shifts in atmospheric circulation that alter moisture source regions, moisture transport pathways, and the associated upstream rainout history during transport (Araguás-Araguás et al., 1998; Jung et al., 2022; Kim et al., 2019; Zhu et al., 2023). Studies have suggested that variations in δ¹⁸O_p over KP are strongly influenced by seasonal changes in atmospheric circulation and moisture transport, as reflected in the pronounced seasonality of precipitation deuterium excess (d-excess) (Lee et al., 2001; Park et al., 2006). The seasonal cycle in d-excess primarily reflects changes in climatic conditions in the moisture-evaporation region, particularly relative humidity during evaporation. Further, d-excess can also be modified by changes in dominant air masses and moisture transport pathways, as well as secondary processes such as condensation and raindrop evaporation. Consequently, the seasonal contrast in precipitation d-excess, with generally lower values in summer and higher values in winter, has been associated with the distinct climatic conditions accompanying the EASM and East Asian winter monsoon, respectively (Yoshimura and Ichiyanagi, 2009).

While many studies (Clemens et al., 2018; Hayashi et al., 2017, 2021; Igarashi and Oba, 2006; Iwamoto and Inouchi, 2007; Wen et al., 2024) have extensively examined the eastern boundary of the EASM system, the KP has received less attention. We close this gap in our study by focusing on variations in δ¹⁸O_p over the KP and their underlying mechanisms, aiming to enhance our understanding of the pan-Asian monsoon system. Here, we explore the dynamics by tracking and quantifying the moisture sources and their precessional modulation of seasonal rainfall and isotopes over the KP. To this end, we consider the Korean stalagmites' 550 000-year δ¹⁸O_c records (discontinuous), which include a new δ¹⁸O_c record from the Eden Cave stalagmite (ED1) and previously published δ¹⁸O_c records from Gwaneum Cave stalagmite (GE1) (Jo et al., 2010, 2014). To link large-scale atmospheric circulation, moisture transport, and precipitation isotope variability over the KP on orbital timescales, we conducted a series of transient climate simulations, time-slice experiments, and water-tagging simulations using the isotope-enabled Community Earth System Model version 1.2 (iCESM1.2), focusing solely on the effect of orbital forcings. The transient simulation with an orbital acceleration (Lorenz and Lohmann, 2004) of 100 spans the past 130 000 years (130 ka). The iCESM1.2 time-slice simulations were conducted with fixed minimum and maximum values of the precession index, reflecting Northern Hemisphere summer perihelion and aphelion conditions, respectively. The results are further interpreted using previously conducted high-resolution CESM1.2 simulations (Chu et al., 2020; Huang et al., 2021; Raavi et al., 2023; Wengel et al., 2021) and a lower resolution transient CESM1.2 simulation, which includes orbital, greenhouse gas, and ice sheet forcings (Timmermann et al., 2022; Yun et al., 2023).

2 Materials and Methods

2.1 South Korean Speleothem Records

We present δ¹⁸O_c records from two Korean stalagmites obtained from Eden Cave (ED; Jo et al., 2014) and Gwaneum Cave (GE; Jo et al., 2010) (see locations in Figs. 1 and 2a), which grew intermittently. The δ¹⁸O_c records for ED1 and GE1 stalagmites were previously measured along with stable carbon isotope (δ¹³C_c) within the scope of Jo et al. (2010, 2014). However, the ED1 δ¹⁸O_c record has remained unpublished to date. Further details regarding the GE1 and ED1records and their covered time spans are provided in Table S1 in the Supplement.

Previous studies of the ED1 and GE1 stalagmites documented their growth fabrics, depositional histories, stable cave environment, and stable-isotope characteristics, thereby constraining the conditions under which their isotopic signals were preserved. The principal growth phase of GE1 was interpreted to have formed under conditions approaching isotopic equilibrium, based on its calcite textures and stable-isotope characteristics, whereas intervals associated with coralloid growth and limited dripwater supply were identified as distinct depositional phases (Jo et al., 2010). ED1, similarly, contains well-characterized growth intervals separated by hiatuses and changes in depositional fabric, supported by an extensive U-series chronology (Jo et al., 2006, 2014). Although certain intervals exhibit high δ¹⁸O_c variability, partially due to the low resolution of the records, the absence of a pronounced precessional cycle is notably observed in the δ¹⁸O_c profiles (Fig. S1 in the Supplement). From another Korean cave, Jo et al. (2017) identified clear millennial-scale resemblances among late Holocene δ¹⁸O_c records from both Korean and Chinese stalagmites, providing an in-depth discussion on the issue in question. These observations suggest that, despite potential site-specific modifications of the calcite isotopic signal, the well-characterized portions of these records can provide observational constraints for evaluating the modeled δ¹⁸O_p response over the KP.

2.2 Isotope-Enabled Model Simulations

For this study, we performed a series of fully coupled paleoclimate simulations with the iCESM1.2 (Brady et al., 2019), using a model resolution of approximately 1.9° × 2.5° in the atmosphere and land and 1° × 1° in the ocean and sea ice, respectively, providing simulated precipitation oxygen and hydrogen isotopes (δ¹⁸O and δ²H) (Brady et al., 2019; Hu et al., 2019; Tabor et al., 2018). We also combined simulated precipitation stable-isotope data to examine precessional variations in precipitation d-excess (δ²H − 8 × δ¹⁸O) (Dansgaard, 1964). To validate our iCESM1.2 present-day (iCESM-PD) simulations, we have used the available IAEA Global Network of Isotopes in Precipitation (GNIP) datasets for the Asian domain (accessible at https://nucleus.iaea.org/wiser, last access: 25 April 2023), published precipitation isotope records (Sinha et al., 2019), and data from the isotopes-incorporated Global Spectral Model – isoGSM2 (Yoshimura et al., 2008).

First, we conducted an idealized transient simulation covering the orbital forcing history of the past 130 000 years (iCESM-130ka), using fixed pre-industrial conditions of greenhouse gas (GHG) concentrations, ice sheets, and bathymetry. The transient simulation employs an orbital temporal acceleration (Lorenz and Lohmann, 2004; Timmermann et al., 2007) factor of 100, consistent with previous studies, which compresses the 130 000-year forcing history into 1300 model years. The acceleration choice is justified by the fact that we focus here on atmospheric processes primarily controlled by changes in land and ocean surface temperatures, with adjustment timescales that are orders of magnitude faster than the precessional cycle forcing. In our 130 ka-orbitally forced iCESM1.2 simulation, we analyzed the spatiotemporal patterns of annual-mean anomalies of precipitation amount and δ¹⁸O_p. We further compared our simulated precipitation patterns with the transient CESM1.2 simulation over the past three million years (CESM-3Ma) (Timmermann et al., 2022; Yun et al., 2023), which was forced by time-varying orbital, GHG, and ice-sheet forcings. CESM-3Ma simulation uses a lower horizontal resolution (3.75° in the atmosphere and 3° in the ocean), and an orbital acceleration of only 5, but does not include isotopes.

Next, we performed model sensitivity experiments with iCESM1.2 to evaluate the equilibrium climate responses to maximum (P_max) and minimum (P_min) precession conditions (iCESM-P_max and iCESM-P_min). These 100-year-long fully coupled time-slice sensitivity experiments branched off from year 900 of a 1000-year-long iCESM-PD control experiment. The two idealized simulations use a high eccentricity (Berger and Loutre, 1991) of 0.0493 to amplify the precessional effects. The simulations reach a quasi-equilibrium state for important atmospheric variables within 100 years. Much of our analysis from these two simulations is based on the June–July–August (JJA) seasonal averages of the last 40 model years. Furthermore, we compared iCESM P_min and P_max precipitation and circulation patterns with ultra-high-resolution simulations ($\mathrm{1}/\mathrm{4}$° in the atmosphere and $\mathrm{1}/\mathrm{10}$° in the ocean) conducted for Marine Isotope Stages (MIS)-5d (115 ka) and MIS-5e (125 ka) using the CESM1.2 (Chu et al., 2020; Huang et al., 2021; Raavi et al., 2023; Wengel et al., 2021) (hiresCESM-5d and hiresCESM-5e). This analysis provides further insights into the robustness of our results, enabling a more nuanced characterization of regional aspects.

To quantify the KP precipitation and δ¹⁸O_p responses to the P_min and P_max conditions, we performed water-tagging experiments using the iCESM to track changes in source regions and their relative contributions to the KP. The experiments were conducted using orbital parameters corresponding to 124 ka (P_min) and 119 ka (P_max), along with current atmospheric and oceanic conditions (including the initial isotopic composition) derived from the iCESM-PD simulation. We tagged a total of 14 boxes (7 each over land and ocean), based on the major contributors to the Asian region (Hu et al., 2019). The box latitude and longitude were chosen following Wen et al. (2024), although we have excluded regions from the calculations with minimal contribution (< 10 %) to the total rainfall and no effect on the anomalous δ¹⁸O_p. The regional pattern of the simulated δ¹⁸O_p is determined by the δ¹⁸O of the vapor from the tagged regions, followed by hydrological processes and rainout history (Tabor et al., 2018). Each tagging experiment was conducted over 40 model years, with the mean values from the last 20 years (both summer and annual) used for analysis. Further information on the tagging experiment is provided in the Supplement (Sect. S1).

3 Results and Discussion

3.1 Precipitation and δ¹⁸O: Model and Proxies

To extract the dominant modes of variability in Eastern Asia, we used precipitation amount and δ¹⁸O_p annual anomalies from the transient iCESM-130ka simulation. We calculated the first Empirical Orthogonal Function (EOF1) and the corresponding principal component (PC1) time series for the selected Eastern Asian region (22–45° N, 100–150° E), including eastern China, the KP, and Japan (Fig. 1). The leading EOF of annual mean precipitation variability exhibits a north-south dipole in the selected EASM domain (Fig. 1a), characterized by strong precessional variability and correlating with the NHSI at 30° N (Fig. 1b). The precipitation anomalies associated with the change in boreal summer insolation exhibit the strongest signal in an east-west elongated structure between 25 and 35° N. For Northern Hemisphere summer perihelion conditions (low precessional index), this area is characterized by anomalously dry conditions, whereas for northern-hemisphere summer aphelion conditions (high precessional index), it is characterized by anomalously wet conditions (Fig. S2). Most notably, Korea lies on the nodal line of the precipitation dipole pattern, suggesting a weak annual-mean precession signal over the KP (Fig. 1a). A similar precipitation pattern also emerges for the CESM-3Ma transient simulation (Yun et al., 2023) over the study area (Fig. S3), indicating that the effects of GHGs and ice sheet forcing are relatively minor in driving the EASM spatial precipitation patterns on precession scales. However, in this case, the nodal line lies slightly south of the KP. The ultra-high-resolution simulations confirm the overall large-scale response pattern in summer precipitation to changes in the precessional index, albeit with somewhat drier conditions over the KP, relative to the coarser-resolution CESM-3Ma and iCESM-130ka simulations (see comparisons in Fig. S2).

Figure 1 First empirical orthogonal function (EOF1) and corresponding principal component (PC1) of precipitation (a, b) and precipitation-δ¹⁸O (c, d) annual anomalies for the selected Eastern Asian region (22–45° N, 100–150° E). Annual anomalies were calculated by subtracting the long-term annual mean value (i.e., the climatological annual mean) from each annual mean. Annual means were obtained by averaging the monthly simulated values. Korean caves and cave locations for the China composite record are marked with yellow stars in panels (a) and (c). A red dashed box marks the Korean Peninsula (KP) region. The black dashed curve in panels (b) and (d) represents summer insolation at 30° N.

Contrasting the meridional dipole in precipitation, the leading mode of δ¹⁸O_p variability exhibits an east-west (zonal) dipole pattern over our selected study region (Fig. 1c). The δ¹⁸O_p dipole, which is characterized by relative ¹⁸O-depletion (enrichment) for western (eastern) Asia during Northern Hemisphere summer perihelion, has previously been presented as a notable feature over the broader Asian summer monsoon region on precession scales (Wen et al., 2024). In the western part of this dipole (China), we observe δ¹⁸O_p signals at pronounced precessional timescales; in the eastern part (Japan, Korea, and nearby seas), we generally find weaker variability that is antiphase with the west (Fig. 1c and d). For the KP region, the simulated precessional signal in δ¹⁸O_p is strongly suppressed (similar to the rainfall pattern), and no distinct precession cyclicity is evident in the simulated time series (Fig. S4).

The near absence of precessional variability in hydroclimate predicted by the models can now be tested using our speleothem isotope data from the KP (Fig. 2a). The Korean speleothem stable carbon isotope (δ¹³C_c) records aligns well with Chinese δ¹⁸O_c during overlapping periods (Fig. 2b), and the growth frequency of Korean speleothems exhibits clear signals that correspond to orbital-timescale hydrologic changes (Jo et al., 2014). However, Korean stalagmitesrecords show no apparent precessional signal (Fig. 2b). The weak covariation between δ¹³C_c and δ¹⁸O_c (ED1, R² = 0.16, p < 0.01; GE1, R² = 0.04, p = 0.02) suggests that kinetic effects are unlikely to be the primary control on their variability. Moreover, the precessional-scale variability in the Korean δ¹⁸O_c records is strongly suppressed relative to the relatively pronounced (∼ 4 ‰) variability observed in the Chinese composite δ¹⁸O_c record (Fig. 2b). This observational evidence is consistent with the muted response of δ¹⁸O_p over the KP to precessional forcing in the simulations (Fig. 1c). This raises the question of which physical processes govern the Korean δ¹³C_c on precessional timescales. Previous studies from the KP (Jo et al., 2010, 2011, 2014) suggested that glacial-interglacial changes in δ¹³C_c could be related to summer-temperature and rainfall-driven anomalies in terrestrial productivity. However, a more detailed investigation of the underlying processes driving precessional-scale variations in δ¹³C_c is beyond the scope of our comprehensive water-isotope analysis of the KP.

Figure 2 Korean Peninsula (red box, Fig. 1) with Korea cave locations (a). Speleothem stable isotope records (δ¹³C and δ¹⁸O) from Korean caves (Gwaneum (GE) cave and Eden (ED) cave) stalagmites (GE1 and ED1) and the Chinese caves composite-δ¹⁸O recordcomparison (b). The marine isotope stage (MIS) is indicated above the curves in panel (b). The black dashed curve in panel (b) represents summer insolation at 30° N.

Figure 3 Comparison of simulated Eastern Asian (EA) rainfall from the climate model simulations with prominent precipitation proxy records from the region over the last 550 ka. Summer Insolation (30° N) and eccentricity (high eccentricity and low/high insolation highlighted with pink/yellow bars) (a). Rainfall (mm d⁻¹) from iCESM (this study) and CESM 3Ma simulations, along with obliquity (b). Loess ¹⁰Be (central China) and EASM rainfall (Beck et al., 2018), along with CO₂ and palynological (Japan) rainfall reconstructions using the Lake Biwa pollen records (c). Lake Biwa-V (Hayashi et al., 2017) is based on vegetation response, and Lake Biwa-P (Tarasov et al., 2011) is based on the modern analogue technique. The marine isotope stage (MIS) is indicated above the top panel curves. For records' references, see Table S1. It should be noted that the Chinese ¹⁰Be record has been plotted after shifting the age by −6 ka (see Sect. S2), which then aligns well with other regional proxies (c) and model-simulated precipitation (Fig. S5).

Figure 4 Difference between precession minimum (Low) and maximum (High). Composites of summer mean precipitation (shaded) and sea level pressure (SLP, contour), from iCESM1.2 (a). Same as panel (a) with CESM Ultra High Resolution (UHR) simulation for MIS5e and 5d (b), and summer mean precipitation-δ¹⁸O, from iCESM1.2 (c). The dashed boxes shown are used to extract 130 ka time series from the model simulations: India, central China, South China, and Korea (left to right in increasing longitude order). The yellow curves on the map indicate major rivers in China. The red star represents the Korean cave region.

According to the transient iCESM-130ka and CESM-3Ma simulations, precipitation over Eastern Asia (Fig. 3b) exhibits precessional signals, which are out of phase with insolation. It is important to note that simulated precessional scale precipitation is suppressed during glacial maxima in eastern Asia, documenting the influence of the 80–120 kyr eccentricity forcing on the amplitude modulation of precession (Fig. 3a and b). There are several East Asian hydroclimate proxies, such as ¹⁰Be in loess from central China (Beck et al., 2018), pollen records from Japan (Hayashi et al., 2017, 2021; Tarasov et al., 2011) (Fig. 3c), and total nitrogen content in Lake Biwa (Iwamoto and Inouchi, 2007), which support and demonstrate the above-mentioned features in regional-simulated precipitation. Speleothem records from Japan (Kato et al., 2021; Mori et al., 2018) lack sufficient data points within the 0–130 ka simulated period, and the intermittent time span does not encompass high-insolation periods, which would be necessary to validate the simulated isotopic composition. Some mismatches between the pollen records and the model occur during MIS11 and low-CO₂ conditions, suggesting that vegetation dynamics in the region may, in part, be driven directly by CO₂, either through temperature or the CO₂ fertilization effect. However, uncertainties in the chronology of the palynological reconstruction preclude a definitive conclusion at this time. Moreover, the simulations indicate a stronger 41 kyr obliquity signal, particularly during periods of low eccentricity, which seems absent in the paleo proxy records (Fig. 3b).

3.2 Precession versus EASM Atmospheric Circulation

As discussed in the previous section, under low-precession conditions, a region of low precipitation in the eastern EASM domain zonally extends from KP across southern Japan into the central North Pacific. Our simulation shows that the region of low precipitation is co-located with an anomalously high sea level pressure (SLP; contour lines, Fig. 4a) and suppression of atmospheric convection (Fig. S6c, contour), in concordance with the intensification of the anomalous WNPSH (geopotential height, Fig. S6c). Supporting evidence comes from an integrated moist static energy budget analysis by Wen et al. (2024), which showed that the reduction in rainfall is linked to descending air motion. Our simulated precipitation pattern, meridional and zonal wind components, and strength of anomalous WNPSH are well supported by the difference between the hiresCESM-5e and hiresCESM-5d experiments (Figs. 4b and S6), even though the anomalies over the KP are slightly stronger compared to the ∼ 2° resolution iCESM.

In addition, the simulated western branch of the anomalously strong WNPSH (Fig. S6c, f) might help transport more moisture from the Pacific Ocean toward northern China and the northern part of KP, resulting in regional positive precipitation anomalies in the model. The rainfall pattern corresponds to an overall weakening and a slight northward shift of the Korean and Japanese summer monsoon systems, specifically the Changma and Meiyu-Baiu fronts, respectively. On the other hand, the KP hydroclimate reconstructed from speleothem proxies (viz., growth frequency) suggests signals of a wetter climate during interglacial periods (Jo et al., 2014), which are characterized by Northern Hemisphere summer perihelion. These observations highlight the importance of considering speleothem-specific characteristics, including site location, local orography, cave processes, and chronological precision, when comparing individual speleothem records with spatially resolved simulations, particularly for distinguishing hydroclimatic changes between closely spaced intervals such as MIS 5e and MIS 5d. It is worth noting that the Korean caves are primarily located in the mountainous Gangwon-do region of South Korea and range in altitude (Jo et al., 2014). It is possible that orbital forcing may also influence the precipitation intensity and the extent of the nodal line, relative to local topographic features (Fig. 4b). Thus, the apparent contrast between the speleothem growth-frequency records and the simulated muted response over the broader KP region may partly reflect spatial heterogeneity in the regional hydroclimate response.

We further explored the mechanisms underlying the precessional-scale muting of precipitation and δ¹⁸O for the KP and other parts along the nodal line, both for summer (JJA) and annual-mean conditions. Using the water tagging method, we estimated the contributions from different moisture sources to precipitation and δ¹⁸O_p at the KP in JJA (Fig. 5). The Pacific Ocean transports more moisture to the KP region during the Northern Hemisphere summer perihelion (Fig. 5a). This anomalous moisture predominantly originates from the subtropical North Pacific (SNP, 0.40 mm d⁻¹) due to the strengthening and northwest shift of the NWPSH. Land moisture (Fig. 5b) from the East Asia region also contributes to the anomalous P_min precipitation over the KP in JJA (EAS, 0.11 mm d⁻¹). However, during P_min less moisture from the northern North Pacific (NNP, −0.15 mm d⁻¹) and the equatorial Indian Ocean (EQI, −0.09 mm d⁻¹) arrives in KP, which compensates for the increase with a small cumulative change in summer precipitation of 0.28 mm d⁻¹ as a result. The annual-mean KP precipitation anomaly between P_min and P_max, shows only a small change, ranging from 0.023 mm d⁻¹ (SNP) to −0.15 mm d⁻¹ (EAS) in terms of moisture contributions (Fig. S7, Table S2).

Figure 5 Differences in average summer (June–August) precipitation (mm d⁻¹) on the Korean Peninsula between P_min and P_max from various moisture sources over (a) ocean and (b) land. Colored text refers to the major contributors (mean-state) to the Korean Peninsula (blue-positive and red-negative for P_min–P_max), as shown in Table S2, which was considered in the total calculations. Inset plot (lower right): the summer precipitation-weighted δ¹⁸O in Korea from major and nearby source contributors from Land (EAS, SAS) and the Ocean (SNP, NNP, EQI). Korea is marked by a yellow star. See the Supplement for the acronyms.

JJA δ¹⁸O_pover the KP, will be determined by both regional rainfall and water vapor transport contributions from the different source regions. Using the isotope tagging method, we calculated summer precipitation-weighted δ¹⁸O_p values (Sect. S1 and Table S2) at the KP (inset, Fig. 5). We find that the moisture transported from the SNP oceanic region to the KP is characterized by anomalously high δ¹⁸O values (i.e., enriched in ¹⁸O) during P_min relative to P_max by +0.22 ‰, with EAS land moisture also contributing higher δ¹⁸O values (P_min–P_max = +0.31 ‰). However, the moisture transported from the distant Indian Ocean (EQI) and the land (SAS) regions shows a compensating effect, with lower δ¹⁸O values (EQI + SAS = −0.38 ‰). Combining the positive δ¹⁸O contributions from the SNP and EAS regions with the negative contributions from the EQI and SAS regions results in a small net increase in δ¹⁸O, corresponding to higher JJA δ¹⁸O_p values (+0.22 ‰) over the KP during P_min. On annual timescales (inset, Fig. S7), this compensation is even more pronounced, leading to a slight net depletion in ¹⁸O (lower in δ¹⁸O_p values) over the KP by −0.09 ‰ (P_min–P_max). Furthermore, the tagging experiment reveals that, for annual means, the Indian region plays a negligible role in precessional-scale changes in precipitation and δ¹⁸O_p over the KP (Fig. S7, Table S2).

Overall, there is a tug-of-war between contributions from different oceanic and continental regions to KP δ¹⁸O_p, leading to the muting of the precessional signal. For instance, the anomalous WNPSH leads to a westward (eastward) shift of the ISM branch towards northern China under low (high) precession conditions, as documented by an analysis of the vertically integrated moisture transport (VIMT, Fig. S8). The schematic in Fig. 6 summarizes these findings by highlighting changes in circulation patterns and simulated δ¹⁸O_p values for P_min and P_max conditions.

Figure 6 Schematics of precession effects on East Asian summer monsoon circulation and spatial distribution of precipitation-δ¹⁸O (blue to dark orange: depleted (low δ¹⁸O) to enriched (high δ¹⁸O) in ¹⁸O in relative to low and high precession conditions). The thickness of the arrows represents the strength of the circulation systems. (a) In the low precession case, ISM strengthening (a more southerly component), Westerlies (northward), and WNPSH (northwest) result in low (high) δ¹⁸O_p over India to northern China (eastern Asia). The solid pink line indicates the EASM front location. (b) For the high precession case, the combined effects of the weakening of the ISM (more northeast component), the southward migration of the Westerlies, and the eastward movement of the WNPSH result in suppressed transport of Pacific moisture into the eastern domain of the EASM, result in high δ¹⁸O_p across India to East Asia, however leading to relatively low δ¹⁸O_p in eastern Asia. The interactions between two circulation systems (ISM and EASM) and the Westerlies create a nodal zone, which includes the Korean Peninsula (KP).

Given the oceanic-continental moisture contributions to the orbital-scale characteristics of precessional-scale changes in KP hydroclimate, we anticipate that the dynamics will also be well captured in d-excess, which reflects the source region (Araguás-Araguás et al., 2000; Jouzel et al., 2013; Pfahl and Sodemann, 2014) and evaporative conditions (Masson-Delmotte et al., 2005; Uemura et al., 2008) of the moisture source. Seasonal variations in d-excess, to some extent, can serve as an analog for precessional changes, as precessional forcing often modulates the amplitude of the seasonal signals in precipitation and isotopic composition, due to the strength and position of atmospheric circulation (Bosmans et al., 2018). As such, we examined d-excess variations on precessional timescales in our iCESM simulations.

3.3 Deuterium Excess in Simulations

d-excess, a secondary parameter derived from stable water isotopic ratios, provides additional insights into hydrological processes. Non-equilibrium conditions can lead to enrichment of ²H in the less stable phase (e.g., the gas phase) due to the higher diffusivity of hydrogen isotopes relative to oxygen. For instance, during evaporation from the ocean surface, the strong gradient in relative humidity accompanied by strong winds leads to higher d-excess in the vapor (Craig and Gordon, 1965; Uemura et al., 2008). It has been suggested that climatic influences on relative humidity, mediated by changes in atmospheric circulation, are captured in the d-excess proxy in precipitation-based paleoarchives (Pfahl and Sodemann, 2014) (e.g., ice core, speleothem). Studies have shown that the d-excess in paleowater can, in fact, be recovered from fluid-inclusion data in certain speleothems (Held et al., 2025; Matthews et al., 2021).

The precipitation d-excess is sensitive to various hydrological processes between the “source” and the “sink” (Xia et al., 2023), such as evaporation at the moisture source location (Craig and Gordon, 1965; Merlivat and Jouzel, 1979), moisture recycling from the land (Gat, 1996; Gat et al., 1994), raindrop evaporation (Risi et al., 2008; Stewart, 1975), and Rayleigh distillation (Xia et al., 2023). On precessional timescales, one crucial parameter is the near-surface relative humidity (Rh_s) at the source location. Over the tropical monsoon region, Rh_s changes may occur mainly in conjunction with wind anomalies, as the annual mean surface temperature changes are quite small on precession scales (Clement et al., 2004). We have analyzed the characteristics of simulated Rh_s at plausible source locations and d-excess variations at the sinks. The spatial pattern and gradient of surface relative humidity over the Asian region are reasonably well captured in the present-day simulation (Fig. S9). Brady et al. (2019) reported that the d-excess values simulated by iCESM1 exhibit a positive median bias of 3.3 ‰ and a poor correlation with GNIP observations. Our reinvestigation of simulated d-excess in the Asian region for present-day conditions shows that GNIP and model data align reasonably well (Fig. S10). A closer look at Fig. 4b in Brady et al. (2019) supports this conclusion. For the Asian region, we found good agreement (r = +0.52, p < 0.01) between the annual mean d-excess in GNIP observations and our iCESM-PD simulation (Fig. S11), with a smaller median bias of +1.75 ‰, compared to the 3.3 ‰ global mean bias. Moreover, the strong seasonality in the observed (GNIP) d-excess for East Asia, specifically for northern and eastern Asian sites, is well captured in iCESM-PD (Fig. S12), which is crucial in the current context because precessional shifts are likely to influence the amplitude of the seasonal cycle, rather than the annual mean.

Figure 7 Second empirical orthogonal function (EOF2, a) and corresponding principal component (PC2, b) for the summer mean surface relative humidity for the 130 ka transient simulation on iCESM1.2. Temporal evolution of surface relative humidity over oceanic source regions, Southeast Indian Ocean: blue curve (c) over blue dashed box (a); and Northwest Pacific Ocean: red curve (d) over red dashed box (a). The black dashed boxes in panel (a) represent the South Asian summer monsoon and the eastern extent of the EASM region. The dashed curve represents summer insolation at 30° N.

To refine our interpretation of d-excess as a proxy for precessional signals over Asia, we further investigate Rh_s variability in the Asian monsoon oceanic source domains, including the Indian and Pacific Oceans. In the calculated EOFs for the Rh_s from the iCESM-130ka simulation, EOF2 represents precessional-scale changes (Fig. 7). Over the southeast Indian Ocean and northwest Pacific Ocean (rectangular boxes in Fig. 7), we observe strong precessional-scale variability in Rh_s (Fig. 7b), with an amplitude change of approximately 3 %–4 % (Fig. 7c, d). The chosen blue and red boxes, which exhibit strong out-of-phase variability (Fig. 7a), respectively, closely correspond with the dominant moisture source regions to the South Asian (Tabor et al., 2018) and Eastern Asian (Wen et al., 2024) monsoon regions under P_min climate (black dashed boxes, Fig. 7a). Subsequently, our results indicate that the South Asian and East Asia (mainly Northeast China) continental domains also exhibit a substantial precessional-scale variability in Rh_s (Fig. S13a, b), with an amplitude of about 10 % change and in phase with the NHSI (Fig. S13c, d).

Moreover, our results show that the P_min and P_max conditions exhibit opposite signals in precipitation d-excess values between South Asia and East Asia (Fig. 8a). The d-excess is high over India and low over central China for P_min. Specifically, in the KP, d-excess is low (high) for high (low) NHSI. It should be noted that the KP d-excess differences for P_min–P_max are predominantly from the changes during the summer season; however, the Indian region's d-excess values are higher throughout the year for P_min relative to P_max (Fig. S14). Our tagging results also show that transported moisture from the Pacific Ocean and the East Asia continental regions contributes to lower d-excess at KP under P_min conditions (Fig. S15), likely due to relatively higher humidity over the source regions (Figs. 7d and S12d). To evaluate the Rh_s and d-excess relationship quantitatively, the translation of the source region's Rh_s change to the sink's precipitation d-excess has been estimated from the iCESM-130ka simulation. A regional Rh_s regression on d-excess shows about a 0.3 ‰–0.4 ‰ increase in summer precipitation d-excess values over South Asia with a 1 % decrease in source region (southeast Indian ocean) Rh_s (i.e., $\mathrm{\Delta }\left(\text{d-excess}\right)/\mathrm{\Delta }{\text{Rh}}_{\mathrm{s}}$ = ∼ −0.4) (Fig. S16a). And a 1 % increase in Rhs in the northwest Pacific Ocean can cause a 0.7 ‰–0.8 ‰ decrease in summer precipitation d-excess across Eastern Asia (i.e., $\mathrm{\Delta }\left(\text{d-excess}\right)/\mathrm{\Delta }{\text{Rh}}_{\mathrm{s}}$ = ∼ −0.8) (Fig. S16b).

Figure 8 Precipitation d-excess for Low (P_min)–High (P_max) precession conditions. Black dashed boxes represent regions: India, Central China, South China, and Korea (left to right), used for extracting the time series from the iCESM1.2 simulations.

3.4 Pan-Asian Monsoon Context

In this section, we refine our understanding of speleothem records from the KP in the broader context of hydroclimate changes across South to Eastern Asia in our transient simulation (Fig. S4). The iCESM-130ka simulation shows dominant anti-correlated precessional variability in precipitation between India (representative of Southern Asia) and central China (Fig. 9a and e). Furthermore, the amplitude of this variability is stronger (weaker) for high (low) eccentricity values, such as during MIS 5a–e and MIS3. On the other hand, precessional-scale variability in simulated KP precipitation, as well as in parts of southern China, is strongly muted (Fig. 9c and g), as discussed above. Moreover, over the KP, we found that precipitation exhibits an 11 kyr higher-order precessional harmonic, for instance, during MIS 5d and MIS 4 (Fig. 9g). This signal likely arises from the nonlinearity of the moisture transport equation, which multiplies precessional-scale winds with precessional changes in moisture, leading to semi-precessional variability.

Figure 9 Simulated annual mean variations in precipitation (dots) for the last 130 ka from an orbitally forced transient simulation conducted with iCESM1.2 for India (a, b), South China (c, d), central China (e, f), and Korea (g, h) (see selected region boxes in Fig. 8). Colors of dots show precipitation-δ¹⁸O (pink to green, right panel) and d-excess (blue to red, left panel) changes. Low insolation periods are marked by gray bars. The marine isotope stage (MIS) is indicated above the top panel curves.

In the Indian region, apart from the strong precessional signals in precipitation, we see precessional variability in d-excess and δ¹⁸O_p, with high (low) d-excess and depleted (enriched) δ¹⁸O_p values for high (low) insolation (Fig. 9a and b). Here, the change in d-excess values is predominantly controlled by the surface conditions at the ISM moisture source region, that is, high and low Rh_s values for low and high precession conditions, respectively, as discussed above. In central China, δ¹⁸O_p varies in phase with the anomalies over India (Fig. 9b and f). However, their d-excess variations are essentially anti-correlated (Fig. 9a and e), in particular for MIS5e and 5c, reflecting the complex interplay between the Indian and Pacific Ocean moisture sources and continental moisture effects. Our results support the previously suggested role of the ISM branch in controlling precipitation isotopic ratios across China (see Fig. S4a and b). Furthermore, we observe that KP precipitation and δ¹⁸O_p exhibit little precessional variability due to the compensating effects described above. In contrast, d-excess over the KP shows a more pronounced precessional signal (Fig. 9g). This d-excess variability responds to the balance between oceanic and continental moisture transport and changes in the surface conditions of source regions (i.e., Rh_s) on precessional scales. Moreover, the KP d-excess precessional-scale variability is a unique feature as compared to other Eastern Asian regions, with high values occurring during MIS 5d, 5b, 4, and 2. The positive feedback (kinetic effect) to the vapor d-excess could be even stronger for P_max relative to P_min when dry/cold winds pass over the ocean (Yoshimura and Ichiyanagi, 2009) (e.g., the Yellow Sea).

4 Summary and Concluding Remarks

Our study focuses on the muted precessional-scale hydroclimate variations in the KP region, as identified in our transient climate model simulations and δ¹⁸O_c records from South Korean caves. The muting of the δ¹⁸O_p signal can be attributed to compensating effects between oceanic and continental moisture sources. We further explored how the KP hydroclimate is linked to large-scale pan-Asian atmospheric variability. Our study illustrates that even along the nodal line of the oxygen isotope dipole pattern, where δ¹⁸O_cand the iCESM simulation indicate weak precessional variability in δ¹⁸O_p, the d-excess may still contain important precessional signals, which are associated with large-scale changes in Pacific, Indian, and continental moisture source region conditions, driven by the precessional changes in the Asian summer monsoon circulation and in the intensity of the WNPSH. Therefore, we advocate validating model simulations using speleothem fluid inclusions (paleowater) to shed more light on regional rainfall variability and large-scale atmospheric circulation changes across multiple timescales, even in regions like the KP where oxygen isotope signals are muted.

Our simulations exhibit another interesting feature that has not previously been recognized: simulated Korean precipitation shows evidence of semi-precessional hydroclimate variability. A few studies have suggested that such semi-precessional signals may result from a response to solar insolation in low latitudes, modulating the EASM. Other potential mechanisms include combination tones of the moisture transport, emerging from the product between precessional signals in specific humidity and atmospheric advection. However, the underlying mechanism for the semi-precessional variability over the KP remains unclear and may require further investigation.

The present study highlights the importance of a multi-proxy approach, utilizing additional speleothem proxies, such as trace elements, δ¹³C, and growth rate or frequency, in conjunction with fluid inclusions (d-excess). While reconstructions of paleo-water d-excess provide valuable information on past moisture source conditions and evaporation processes, suitable fluid inclusions are not always preserved or available for analysis. In such cases, oxygen-17 excess (¹⁷O-excess), which can be measured directly from speleothem calcite and is not sensitive to temperature variations, offers a promising complementary proxy (Sha et al., 2020, 2023). Future paired measurements of δ¹⁸O, d-excess, and ¹⁷O-excess in speleothems could therefore provide a powerful framework for testing the circulation and moisture-source mechanisms inferred in this study. In addition, multiproxy speleothem data, along with data from marine and lacustrine archives, may provide further insights into regional hydroclimatic variability in the EASM region, and help reconcile different hydroclimatic records and elucidate their underlying mechanisms.