Deciphering local and regional hydroclimate resolves contradicting evidence on the Asian monsoon evolution

Domain: earth_science
Published: 2023-09-14
Editorial outcome: major_revision

Abstract
The winter and summer monsoons in Southeast Asia are important but highly variable sources of rainfall. Current understanding of the winter monsoon is limited by conflicting proxy observations, resulting from the decoupling of regional atmospheric circulation patterns and local rainfall dynamics. These signals are difficult to decipher in paleoclimate reconstructions. Here, we present a winter monsoon speleothem record from Southeast Asia covering the Holocene and find that winter and summer rainfall changed synchronously, forced by changes in the Pacific and Indian Oceans. In contrast, regional atmospheric circulation shows an inverse relation between winter and summer controlled by seasonal insolation over the Northern Hemisphere. We show that disentangling the local and regional signal in paleoclimate reconstructions is crucial in understanding and projecting winter and summer monsoon variability in Southeast Asia.

Introduction
The Asian winter monsoon is commonly associated with cold and dry northwesterly winds across East Asia and northeasterlies in Southeast Asia1. While not associated with significant rainfall over most regions, the winter monsoon delivers substantial rainfall in some coastal regions, such as Vietnam, the Philippines, Southeast India, Sri Lanka and Japan, where it plays a critical role in agriculture, water resources and natural hazard risks related to flooding and landslides. These regions include some of the world’s largest food producers and exporters, causing not only the regional economy but also the already precarious global food trade to be vulnerable to changes in winter monsoon rainfall. In contrast to the very well-studied Southwest Summer Monsoon (SWM), there are no robust records documenting the long-term changes in Southeast Asian rainfall associated with the Northeast Winter Monsoon (NEM) under pre-industrial conditions. Therefore, the timing and mechanisms of hydroclimate variability in the region over longer timescales are not very well understood. Consequently, many climate models underestimate observed NEM rainfall by as much as 50%, leaving considerable uncertainty in future climate projections. These uncertainties are a harmful knowledge gap impacting planning for food production and geohazard management across this broad region1.

Modern meteorological observations show a strong coupling between the Asian winter and summer monsoons2. However, while several paleoclimate studies have found a positive correlation over glacial-interglacial timescales3–5, there is still no consensus on whether these two seasonal monsoon systems were similarly linked during the Holocene6–15. Several studies tried to address the co-evolution of the winter and summer monsoon in Southeast Asia over the Holocene, arguing for16 and against10,17,18 an inverse correlation, while others found that this relation changes depending on timescales13,19. Interestingly, solar insolation was argued to be the main driver of both a positive and negative correlation10,16,17. Thus, the nature and drivers of Holocene summer and winter monsoon variability and their co-evolution in Southeast Asia remain unclear and additional paleoclimate records are needed.

One of the issues often overlooked is the regional and local component of the Southeast Asian monsoon system. Both the winter and summer monsoon are characterized by large-scale circulation patterns common across the broad Asian monsoon region, which can be traced via oxygen isotopes in precipitation20. However, the timing of peak local precipitation is highly variable and oxygen isotopes fail to capture the local hydroclimate in many locations20,21. This is particularly true for coastal central Vietnam which does not receive substantial rainfall during the summer monsoon season like the rest of mainland Southeast Asia (MSEA; Fig. 1a, b). This is due to the Trường Sơn mountains, which act as an orographic barrier and create a rain shadow effect that limit the impact of southwesterly-sourced moisture during summer20,22,23. Instead the region receives over 75% of its annual rainfall between September to November from the western Pacific (Fig. 1c), and the effective rainfall (difference between rainfall and evapotranspiration) is positive from September to December only (Fig. 1d). Autumn and winter rainfall in Vietnam is related to the northeasterly winds of the NEM circulation sourced from the Pacific (Fig. 1b), unlike the SWM, which is sourced from the Bay of Bengal (Fig. 1a20). Therefore, paleoclimate records from central Vietnam allow us to reconstruct NEM variability and to test its relationship with the SWM. Central Vietnam is likely one of the most extreme examples for the decoupling of the local and regional signal: lows and highs in oxygen isotope ratios occur during summer and winter respectively, but peak rainfall is in autumn. However, this decoupling of local (rainfall amount) and regional (oxygen isotopes) components within the monsoon also occurs across other regions of Southeast Asia20,21,24,25. Deconvolving the correlation between the SWM and NEM can only be addressed if both components are considered and studied in detail across different timescales. We present a record of NEM variability during the Holocene using a speleothem from central Vietnam and distinguish between the regional circulation, reflected in oxygen isotopes, and local rainfall amount via carbon isotopes and trace elements. We apply a proxy system modeling approach to disentangle the regional component in our oxygen isotope record from the local signal. Using a Monte-Carlo principal component analysis, we investigate the correlation between NEM and SWM on a local and regional scale. Further, we compare the co-evolution of NEM and SWM during the mid-Holocene, a time with pronounced changes in local hydroclimate and use climate modeling to understand the mechanisms behind this shift in summer and winter monsoon.

Results
Multiproxy records of Southeast Asian monsoon precipitation
We use oxygen and carbon isotope ratios (δ18O and δ13C) and high-resolution trace elements (Mg, Sr, Ba, P, U, Al, Si) measured in stalagmite TD3 (Fig. S1), recovered from Thien Duong Cave, central Vietnam (17.5°N, 106.5°E) in 2015 (Fig. 1) to reconstruct the NEM rainfall and circulation patterns in Southeast Asia over the Holocene. In central Vietnam, δ18O in precipitation has a pronounced seasonal cycle, with low values during summer (June to August) and high values during the rest of the year, similar to other parts of MSEA20,26,27. This seasonal cycle in δ18O is controlled by the shift between moisture sources, which include the Indian Ocean during summer (low δ18O) and the Pacific Ocean (high δ18O) during the rest of the year20. Seasonal rainfall patterns in central Vietnam differ strongly from the observed cycles in δ18O, with a peak in rainfall between September and December (Fig. 1c). Therefore, central Vietnam offers the opportunity to reconstruct regional circulation patterns and local rainfall forced by the NEM. While the former is reflected in the speleothem δ18O signal, the latter is recorded by speleothem δ13C and Mg/Ca.

Here, we interpret δ13C and Mg/Ca variability as primarily driven by prior calcite precipitation (PCP), a process commonly associated with co-variance between δ13C and Mg/Ca. When PCP occurs, 12CO2 is preferentially degassed and magnesium, as well as other trace elements with partition coefficients smaller than 1, is preferentially excluded from calcite precipitated in the epikarst or on stalactite surfaces, leading to increased δ13C and Mg/Ca in the remaining solution which ultimately feeds the stalagmite28. PCP is enhanced during dry periods, as has been shown by previous studies29–32. Therefore, we interpret increased Mg/Ca and δ13C values in TD3 as indicative of relatively dry conditions at our cave site (Fig. S1a–e). This interpretation is further supported by Principal Component Analysis (PCA), which reveals that Mg/Ca, Sr/Ca and Ba/Ca co-vary in TD3, consistent with a common PCP control (Fig. S1f).

While several studies have interpreted speleothem trace elements and carbon isotopes as proxies for PCP and local water balance33, recently, it has been recognized that PCP can also alter the oxygen isotopic composition of cave waters28, and hence speleothem calcite. Similar to δ13C, δ18O is subject to Rayleigh-based isotope enrichment during calcite precipitation34 enhanced with increasing PCP, especially if PCP occurs at the cave ceiling28. Previously, it was believed that such an enrichment of the oxygen isotope ratio in drip water is sufficiently counteracted by HCO3- ↔ H2O buffering with the water reservoir35. However, more recent studies showed that this buffering is limited by lengthy exchange times between the bicarbonate and the water reservoir36 and that oxygen isotopes might be significantly increased by PCP. In order to distinguish between PCP-induced changes in δ18O and the initial rainfall δ18O signal, we estimate the extent of PCP which occurred throughout the reconstructed time period based on the Mg/Ca variability and a Rayleigh-based oxygen isotope evolution model37–40. We conduct the correction procedure with varying input parameters (temperature and pCO2) to estimate uncertainties of our approach (Fig. 2a). The simulation using 426 ppmv and 21°C (see Methods) is assumed to be the most accurate and is used as the main time series, whereas the remaining 12 time series are used as uncertainties. Details about cave monitoring and the δ18O correction can be found in the Methods section. We apply a Rayleigh-model to estimate the length of PCP for each time point with corresponding Mg/Ca values and use it to correct the speleothem δ18O record for the PCP-induced δ18O changes. We then test the impacts the correction has on our analysis using a synthesis of δ18O records across Southeast Asia. Because of the uncertainties inherent to our δ18O correction due to simplifying cave internal processes, we focus on large-scale trends in the time series.

We interpret the corrected δ18O record of TD3 as reflective of changes in the moisture source related to the NEM20, and thus a measure of NEM strength. The corrected δ18O record of TD3 shows lower δ18O values during the mid-Holocene and higher values in the early and late Holocene (Figs. 2f and S6). This pattern agrees with previously reported speleothem records from Borneo, indicating decreasing δ18O values in the mid-Holocene41 (Fig. 2g). This trend has been interpreted as a result of enhanced convection over this region, broadly following autumn insolation in the Northern Hemisphere42. The δ18O records from northern MSEA and Sumatra instead show decreased values during the early Holocene and an increasing trend after (Fig. 2d, e and h), which is consistent with summer rather than autumn insolation.

The lower δ18O values observed in TD3 and the Borneo record suggest an intensification of the NEM during the mid-Holocene, in contrast to local hydroclimate proxies (Mg/Ca and δ13C) in TD3 and other paleoclimate records (Figs. S2; S6) which suggest pronounced drying during this time. The dry period started at 5 ka and ended abruptly at 3.2 ka, interrupted by two brief wet episodes at ~4.3 and ~3.7 ka. A dry period is also observed in other speleothem and lake records from Southeast Asia indicating widespread dry mid-Holocene conditions (Fig. S2). Specifically, Mg/Ca ratios and carbon isotopes in a speleothem record from Laos (Fig. S2a) show a pronounced local drought, beginning around 5 ka and lasting until around 3 ka43. Contemporaneously, a speleothem from southern Thailand exhibits a 2000-years-long hiatus, potentially indicating reduced water supply to sustain speleothem growth in the south of MSEA (Fig. S2c). While the δ18O record from Borneo suggests strong convection over this region from 5 ka to 3 ka, as with our δ18O record in central Vietnam (Figs. 2a and S2d), the δ13C record from Borneo reveals a period of increased values at 4 ka, potentially indicating locally dry conditions. Two speleothem records from the Maritime Continent indicate a shift towards dry conditions around 5 ka in Sumatra (Fig. S2e) and Flores (Fig. S2f). Dry conditions during this time have also been reported from lake sediments in northern Thailand44, Cambodia45,46, and southern China47. Our data (δ13C and Mg/Ca) shows that NEM precipitation declined during the mid-Holocene, while our δ18O record suggests a more distal moisture source, potentially reflecting enhanced NEM winds. Here, we aim at resolving this apparent contradiction. Firstly, we collate a dataset of proxy records across Southeast Asia and study their coherence with clear implications on common forcing mechanisms. As opposed to coherent local SWM and NEM rainfall variability, their regional trends are anti-correlated. Secondly, we underpin these findings on the interplay of regional vs. local drivers based on model simulations.

Holocene monsoon coherence across Southeast Asia
To investigate the spatiotemporal pattern of regional monsoon variability across Southeast Asia during the Holocene and to compare our proxy record in Vietnam with nearby records, we synthesized published speleothem δ18O records in Southeast Asia (Fig. S3) by using Monte-Carlo Principal Component Analysis (MC-PCA48, see Methods). To evaluate the impact of our PCP-corrected δ18O time series on the MC-PCA results and estimate the reliability of our approach we perform the MC-PCA on both the uncorrected (Fig. S4) and corrected (Fig. S5) δ18O record of TD3. Other δ18O records used in this analysis were not corrected for PCP, because the authors of the original studies found limited or no evidence for PCP, except for the Tham Doun Mai record in Laos43. This record shows no correlation between δ13C and δ18O, suggesting that either there was sufficient time to reach equilibrium after PCP or other processes might have dominated the isotopic signatures in this cave.

The leading Principal Component (PC1) of δ18O records explains 57% of the variance when TD3 is not corrected for PCP (Fig. S4a). PC1 exhibits a strong positive loading on the Laos (Tham Doun Mai Cave) (0.55) and China (Dongge Cave) (0.56) records, and a modest loading of 0.26 on the Sumatra (Tangga Cave) record. These results suggest a positive relation between these records which mean they co-varied (increased or decreased synchronously) over the Holocene. The Tangga Cave record and other speleothem records from India and East Asia showed a high coherence over glacial-interglacial timescales25, which appear less pronounced during the Holocene in our results. By contrast, PC1 displays a strong negative relation with Borneo (Bukit Assam Cave) (−0.48) and a modest negative relation with Vietnam (Thien Duong) (−0.29) (Fig. 3a). The orbital trend in the five δ18O records is aligned with summer and autumn insolation (Fig. 3g): Here, we compare the long-term trends in the δ18O records with temporal changes in insolation and find that records loading strongly positive on PC1 follow summer insolation most closely, while those with a negative correlation to PC1 follow autumn insolation. The long-term trends of each δ18O time series used for the PCA shows significant strong correlations with seasonal insolation for each site (Fig. 3g), further evidencing this link. We conclude that seasonal insolation is a main driver of δ18O across Southeast Asia for both the SWM and NEM. δ18O in all records becomes more negative when insolation is increased, however, autumn and summer insolation peak at different times during the Holocene, potentially causing the apparent anti-correlation (Figs. 3 and 4). Modern observations and model simulations indicate that decreased δ18O in rainfall over Southeast Asia is linked to enhanced convection, more distal moisture source, and increased rainout along the trajectory20,21,49. Thus, we propose that higher insolation leads to increased rainout and the moisture travel distance from the source, which results in decreasing δ18O values in both the SWM and NEM.

Interestingly, the majority of our PCP-corrected time series (presented as the cubic fits, Fig. 3g) closely follows insolation, however, those corrections with high cave air pCO2 (3000 and 2000 ppmv) are deviating the most. Values much higher than 1600 ppmv (highest modern cave air pCO2 is ~1600 ppmv) are unlikely to have been reached over the Holocene, due to the cave being well-ventilated and changes in atmospheric CO2 over the Holocene only varied between ~240 to 270 ppmv (by 30 ppmv)50. While unlikely, it should be noted that an increase in cave air pCO2 to 2000 and 3000 ppmv alters the long-term trend of δ18O and reduces the correlation between the cubic fit of δ18O and SON insolation to r=−0.69 (for 3000 ppmv, 21°C), compared to r=−0.92 (for 426 ppmv, 21°C). Considering the uncertainties in our correction approach and recognizing that cave environment processes are simplified, we find that our results are a promising step in accounting for PCP influence on speleothem δ18O, which can enable the separate investigation of local and regional dynamics in our proxy record.

The overall positive loading on the PC2 of δ18O indicates that all records share, to some extent, a common forcing. Indeed, PC2 exhibits a strongly-to-moderately positive loading on Sumatra (0.80), Vietnam (0.36), and Borneo (0.42), while weak relationship exist with Laos (0.06) and China (Dongge) (0.10) (Fig. 3b). Thus, after removing the trends contained in PC1, PC2 is more consistent with the local hydroclimatic pattern (Fig. 4c), however, also contains some of the autumn insolation trends reflecting the Vietnam and Borneo time series (Fig. 4b). When comparing the time series of PC2 of δ18O and PC1 of δ13C, we also find that PC2 of the uncorrected δ18O time series shares more similarities with PC1 of the carbon records (Fig. 4c). This highlights that our correction indeed has the potential to remove the local hydroclimatic (PCP) signal. The mismatch between local and regional signals of the Southeast Asian monsoon is likely one of the limiting factor in deconvolving summer and winter monsoon variability across this region in past studies.

In order to investigate the local components of the Southeast Asian monsoon, we now apply the MC-PCA on δ13C speleothem records and one lake record to compensate for the δ13C from Dongge cave which is not publicly available. Results of the PCA for speleothem δ13C records and the titanium content of Lake Huguang Maar14 from southern China reveal that all records are strongly to moderately positive correlated with PC1 (Fig. 3e). The strongest loading between PC1 and the proxy records is found for Laos (0.62), followed by Sumatra (0.46), and Borneo (0.41) (Fig. 3e and S7). The winter monsoon records—e.g. the stalagmite from Vietnam (TD3) and the Lake Huguang Maar record from southern China—show a weaker, but also positive relation with PC1 (Vietnam: 0.27 and southern China: 0.37). This suggests that the first PC derived from stalagmite δ13C records and a lake record from Southeast Asia likely share a common forcing throughout the Holocene, despite the different modern rainfall regimes unique to each site (Fig. 1c, d20). However, the records more dominated by winter monsoon circulation also show a stronger loading on PC2, with 0.57 for Vietnam and 0.43 for southern China (Huguang Maar). Additionally, Sumatra shows a strong negative relation with PC2 of −0.59, which can be explained by local hydroclimate, characterized by minimum rainfall in boreal summer and strong influence of the Indo-Australian summer monsoon25. The Borneo record loads most heavily on PC3 (−0.81). This may be due to the local hydroclimate with year-round rainfall or, alternatively, to the dominance of non-hydrologic controls on speleothem δ13C in Borneo (Bukit Assam Cave)51. Nevertheless, the first three PCs combined explain 77% of the variance, thus covering most of the variability present in all five records. Crucially, the dry conditions observed in the proxy records between ~5 to ~3 ka are evident in both PC1 and PC2, suggesting strong coherence between the sites during this time (Fig. S7). Furthermore, these results highlight that the dry period observed in TD3 and other paleoclimate records, affected most of Southeast Asia, with reduced rainfall confined not just to the summer monsoon season, but also to the autumn/winter season. In the next step we investigate the forcing mechanisms of this dry period and their effects on the regional and local aspects of the Southeast Asian monsoon.

Monsoon coherence during the mid-Holocene
The mid-Holocene is characterized by large-scale changes in tropical climate conditions and mean ocean states—e.g. reduced ENSO variance, pronounced droughts or monsoon intensification (Figs. S8 and S9a). Based on our results, Southeast Asian speleothem records suggest substantial local drying during this period, whereas the regional NEM intensity (reflected in speleothem δ18O) seem to have increased (Figs. 2 and S2). We investigate this apparent mismatch by using three state-of-the-art climate models and synthesize results of mid-Holocene idealized sensitivity experiments from EC-Earth52, iCESM53, and the University of Toronto version of CCSM4 (UofT CCSM4;54) climate models. These model experiments simulate mid-Holocene climate following PMIP protocol, but additionally account for changes in vegetation cover across the Saharan region, a forcing which has recently been suggested to affect local SWM rainfall43. For all modeling experiments, we compare results from two sensitivity experiments for the mid-Holocene (6 ka). (1) MH simulations with mid-Holocene insolation, greenhouse gases, vegetation cover and dust load based on the PMIP3 (EC-Earth) or PMIP4 (iCESM and UofT CCSM4) protocols. (2) MH simulations have a prescribed vegetated ‘Green Sahara’ (GS), additionally, EC-Earth and iCESM have a reduced dust load. Further details on the varying methods and parametrizations incorporated in these different modeling experiments can be found in Pausata et al. (2017; EC-Earth)52, Tabor et al., (2020; iCESM)53, and Chandon and Peltier et al., (2020; UofT CCSM4)54. In summary, we compare runs simulating the vegetation over the Sahara as pre-industrial with simulations using a Green Sahara, to isolate the effect of dust loading (MH ORB -MH GS) during the mid-Holocene.

Starting with the regional Southeast Asian monsoon, we use iCESM simulations for the mid-Holocene, to examine the impact of changing mid-Holocene climate on Asian rainfall δ18O variability (Figs. 5 and 6). During summer, rainfall throughout southern China and MSEA exhibits significantly higher δ18O values (Fig. 5a). Simulated autumn (SON) δ18O values in rainfall across southern China/northern MSEA and Sumatra also increase, which is in agreement with proxy observations over the Holocene (Figs. 3a, c and 5b). While absolute values cannot be compared between the model and the proxy records, the spatial patterns concur: during the mid-Holocene the δ18O proxy records can be grouped by trends: low values peak during the mid-Holocene in Borneo and Vietnam (Figs. 3g and 5b), whereas records from southern China/northern MSEA and Sumatra show continuously decreasing trends from the early to the late Holocene (Figs. 3g and 5a). The same spatial pattern can be seen in iCESM (SON) δ18O values, indicating an increase in δ18O values in northern MSEA and Sumatra only. These results suggest that the inverse relation between SWM and NEM evident over the entire Holocene (Fig. 3) was not only present but likely enhanced during the mid-Holocene due to vegetation changes in the Sahara.

In terms of rainfall, all three models show widespread drying in Southeast Asia during the Mid-Holocene, with EC-Earth showing the largest changes (Fig. 6). The EC-Earth model indicates decreased summer and autumn rainfall across much of Southeast Asia in a scenario of reduced Saharan vegetation and higher dust load relative to the Green Sahara state, contributing to the overall decline in annual rainfall (Fig. 6). While EC-Earth agrees with the dry conditions inferred from paleoclimate proxies across Southeast Asia, UofT CCSM4 suggests wetter summer conditions during the mid-Holocene in northern Southeast Asia, contrasting proxy observations from Laos and southern China (Fig. S2). Here, we compared the changes in the absolute local proxy record before (8–5 ka BP) and after (5–3.5 ka BP) the end of the Green Sahara, where all proxy records used in the MC-PCA suggest drier conditions during 5 to 3.5 ka BP (Fig. 4a, b). During autumn the UofT CCSM4 model simulates reduced precipitation in northern MSEA and Sumatra and wetter conditions in parts of Borneo and no significant change in annual mean precipitation over MSEA. Combined, the mean of all three models suggests drying during autumn and a lower annual mean in the north of MSEA and southern Southeast Asia (Fig. 6 bottom panel).

While there is some disagreement among the models, the EC-Earth simulations generally agree with our interpretation of reduced rainfall in MSEA. Specifically, reduced Saharan vegetation and increased dust forcing result in an ‘El Niño-like’ mean state across the Pacific (Fig. S9b). Such an El Niño-like mean state is associated, with weakened easterly trade winds leading to warmer central/eastern Pacific SSTs, and increased easterlies in the eastern Indian Ocean resulting in lower SSTs off Sumatra (Fig. S9b). Upwelling in the South China Sea and the East China Sea could potentially inhibit convective activity over the moisture source region during the peak of the northeast winter monsoon season in Vietnam (Fig. S9b), contributing to the much lower average rainfall observed in our proxy records. Our simulations and proxy synthesis show that rainfall during the mid-Holocene was reduced in summer, autumn and year-around (annual mean), not simply during the summer monsoon as previously suggested43.

Interestingly, the mid-Holocene spatial δ18O pattern simulated in iCESM agrees well with spatial δ18O patterns of proxy records in the PCA. However, iCESM results indicate no significant changes in precipitation amount, despite pronounced hydrological changes observed in the proxy records during the mid-Holocene. Thus, δ18O and precipitation share no direct link in iCESM, a result which is also found in the deviation between δ18O and δ13C records from MSEA. These are significant results because climate models still have limits to accurately simulate local rainfall in the tropics, due to the inability to properly account for convection. The simulation of δ18O appears to be more reliable in our results than rainfall. This highlights that local rainfall and the regional monsoon signal preserved in δ18O need to be considered in both climate model simulations and proxy reconstructions when investigating summer and winter monsoon relations in Asia.

Here, we addressed the long-term co-evolution between summer and winter monsoon in Asia over the Holocene by reconstructing NEM rainfall variability. We employed speleothem-based proxies from central Vietnam to infer local hydroclimatic and regional monsoon changes over the last ca. 8 ka. We showed that SWM and NEM rainfall in Southeast Asia strongly covaried over most of the Holocene, which is in good agreement with climate models simulations (Fig. 6). As the main driver for monsoon coherence between sites in Southeast Asia, we identify eastward shifts of the Walker Circulation as regulators of rainfall over the region (Fig. S9b). The tight positive coupling of the local rainfall from SWM and NEM stands in contrast to the regional components, which shows an inverse relation between SWM and NEM driven by changes in insolation in the Northern Hemisphere over the Holocene. Thus, their main forcing mechanisms are internal, via ocean-atmosphere feedbacks and external via insolation, respectively. Our results highlight the need for careful distinction between the local and regional aspects of the Southeast Asian monsoon, and inform on the choice of proxies when aiming to reconstruct monsoon variability.

Methods
Uranium-thorium dating
Hand drilled powder samples from TD3 weighing 250–300 mg were prepared for U-Th dating at the Massachusetts Institute of Technology following Edwards et al. (1987)55. Samples were dissolved in nitric acid and spiked with a 229Th–233U–236U tracer, followed by isolation of U and Th by iron co-precipitation and elution in columns with AG1-X8 resin. The isolated U and Th fractions were analyzed using a Nu Plasma II-ES multi-collector inductively coupled plasma mass spectrometer (MC-ICP-MS) equipped with an Aridus 2 desolvating nebulizer. Reported errors for 238U and 232Th concentrations are estimated to be ±1% due to uncertainties in spike concentration; analytical uncertainties are smaller than 1%. Ages were corrected for detrital 230Th assuming an initial 230Th/232Th atomic ratio of (4.4 ± 2.2) × 10−6 and the age was reported in years before present, where the present is defined as 1950 C.E. Decay constants for 230Th and 234U were taken from Cheng et al. (2013)56 and, the decay constant for 238U of 1.55125 × 10−10 yr−157 is used. The depth-age model was calculated in COPRA58 based on 21 U-Th dates, from which 6 have been excluded by stratigraphic order method.

Trace elements
Trace element measurements were made using a RESOlution-SE Compact 193 nm excimer laser ablation (LA) system in tandem with an Agilent 8900 Inductively Coupled Plasma Mass Spectrometer (ICP-MS) at Waikato University, New Zealand. Analyses were conducted by pulsing the laser at 15 Hz with a 50 μm spot size and energy density of 5 J/cm2 at a scan speed of 26.118 μm/s. Measurements were made inside the milled track to ensure that stable isotope and trace elemental data have the same depth-age relation. Ultra-high purity helium is used as the carrier gas to deliver the ablated sample from the LA system to the ICP-MS. The ICP-MS was optimised to maximum sensitivity daily using the NIST 612 glass standard. Background counts (gas background, measured with the laser off) were collected for 40 seconds between samples, allowing for complete washout. NIST glass standards (610, 612) were analysed after each sample track to account for any instrumental drift. Data processing was performed using Iolite (v3.32;59). Background counts are subtracted from the raw data and all data were standardised to NIST 612 glass certified reference values. NIST 610 glass was utilized as a secondary standard. GeoReM database60 is utilized for NIST glass certified values.

Stable isotopes
In total 371 low-resolution (1mm) and 853 high-resolution samples (200μm) were drilled and analyzed. Carbonate stable isotope analyses were performed using a Delta V Advantage isotope ratio mass spectrometer (IRMS) (Thermo Fisher Scientific) interfaced with a GasBench II and a CTC GC PAL autosampler at ETH in Zürich (high-resolution) and Northumbria University (low-resolution). Results were normalized against the international NBS18 and NBS19 standards and in-house standard Wiley (δ18O = −7.2, δ13C = −0.41), and reported in permil (‰) relative to VPDB. Carbonate was dissolved using phosphoric acid and the sample weight was between 70 to 100 µg.

Mineral determination
Fourier-transform infrared spectroscopy (FTIR) was used to determine carbonate mineralogy at Northumbria University, Newcastle. 3 mg of stalagmite powder was analyzed on a Perkin Elmer Spectrum RX I FTIR over a spectrum of 640-4200 cm−1 with 32 replications. The detection limit was about 5 wt% and the measurement had a precision of approximately ±10% resolution is 0.8 cm−1. Before each run the background spectrum was measured and automatically subtracted from the sample spectrum. TD3 shows a typical calcite FTIR spectrum61, with pronounced peaks at 712 cm−1 and 872 cm−1.

PCP removal
Similarity of Mg/Ca and δ13C together with MC-PCA analysis reveals that Mg/Ca of TD3 is a proxy for PCP. PCP-influenced bicarbonate of drip water is not representing the initial rainfall δ18O signal because it is modified by Rayleigh-induced isotopic enrichment. Mg/Ca can provide information on the extent of PCP and its variation over time and, thus, can be used to remove the PCP signal from the speleothem δ18O. With this approach it is feasible to obtain a record of initial bicarbonate δ18O of cave drip water, which in turns allows determination of the δ18O of infiltration/rainwater.

To estimate PCP length from the Mg/Ca time series it is necessary to calculate calcite precipitation from a solution, which we approximate by using an exponential approach62. In addition, we use a solely Rayleigh-enabled isotope evolution model37–40,63 to simulate δ18O evolution while CO2 degasses and CaCO3 precipitates. In the present study, we extend this approach by implementing the behaviour of the Mg/Ca ratio with evolving CaCO3 precipitation. Skiba & Fohlmeister, 202340 found that changes in drip interval result in only small deviations in the isotope composition of CaCO3 in comparison to the potential effect of variations in the duration of PCP. As the evolution of the Mg/Ca ratio of a solution is treated similar to the isotope evolution in the model (both are described by a Rayleigh approach), the influence of drip interval changes on the Mg/Ca ratio of speleothem CaCO3 can also be considered small compared to those induced by changes in the duration of PCP. Thus, we neglect changes in isotopic composition which are caused by variations in drip interval here.

To simulate calcite precipitation, the evolution of Mg/Ca and the isotope composition over time, we use a modern cave temperature of 21°C64 which is assumed to have remained relatively constant throughout the last 8000 yrs of speleothem growth. This is a reasonable assumption as cave temperatures mainly reflect mean annual air temperatures and the Holocene is known to have had a relatively constant temperature history in the (sub-)tropics. Furthermore, the model needs additional parameters such as an initial Ca2+ concentration of drip water, Ca2+initial. This parameter depends mostly on soil or karst pCO2 and the equilibrium Ca2+ concentration, Ca2+eq, with respect to the cave air pCO2. During monitoring of Thien Duong cave in spring (2015, 2016) and summer (2016), cave air CO2 values of ~420 ppmv to ~1580 ppmv have been measured64. This corresponds to a range of Ca2+eq values of ~0.7 to ~1.1 mol/m3. We use the minimum value of 0.7 mol/m3, which corresponds to a pCO2 of 426 ppmv. This value represents the lower end member of observed cave air pCO2, where carbonate precipitation, i.e. speleothem growth, is enhanced compared to periods of elevated cave air CO2 levels. When using a Ca2+eq of 0.7 mol/m3, we need to set the Ca2+initial to a minimum value of ~2.73 mol/m3 to explain the range of Mg/Ca observed in the speleothem data. This value corresponds to an equilibrium pCO2 value of ~25,300 ppmv, which is one order of magnitude larger than observed in present-day soil air pCO2 (i.e. 3400 and 1600 ppmv in spring 2015 and 2016, respectively;64). However, there is growing evidence that the air in karst systems can have much higher pCO2 values than soils32,65–68. Therefore, we consider the Ca2+initial value of 2.73 mol/m3 or more as realistic and simulate two cases: 1) Ca2+initial of 2.73 mol/m3, i.e. in equilibrium with a pCO2 of ~25,300 ppmv and 2) Ca2+initial of 3.18 mol/m3, i.e. in equilibrium with a pCO2 of ~40,000 ppmv. Similar to temperature, we assume that modern CO2 dynamics did not change substantially over the last 8000 yrs and keep the Ca2+initial and Ca2+eq parameters constant throughout the simulation.

During calcite precipitation, the Mg/Ca ratio of the solution changes according to the partition coefficient, D(Mg), between the solution and precipitated calcite (e.g., Johnson et al., 200669). We use a D(Mg) of 0.015, calculated for T=21°C33 which is based on laboratory experiments performed under cave analogous conditions.

To account for the possible variations in cave air temperature and pCO2, we simulate 12 additional PCP corrections using modern 21° ± 2, which is the Holocene range of cave air temperatures reconstructed from a cave in Borneo70 and pCO2 ranging from 260 to 3000 ppmv. As the spatial resolution of the Mg/Ca record is higher than the δ18O record, we need to down sample the Mg/Ca record to the coarser timescale of the δ18O in order to use both proxies in our calculations. For down sampling, we use Gaussian kernels (Gaussian kernel width = 0.5 mm, which is the maximum integrated depth during sampling of carbonate powder for isotope measurements) and calculate smoothed Mg/Ca values for each time step for which δ18O is available. In absence of values for the initial Mg/Ca ratio of the solution, we identify the smallest Mg/Ca value from the smoothed record and assume that this value represents the initial Mg/Ca value. In other words, we assume that the Mg/Ca value at t = 0s, when no PCP has occurred (Mg/Ca0 = 8.3 mmol/mol) was constant through time and forms the starting point for the coupled calcite precipitation—element ratio—isotope composition simulation. After the drip water comes in contact with cave air, calcite precipitates with a decreasing trend while Mg/Ca and δ18O of the solution are increasing (Fig. S10).

From the speleothem Mg/Ca values, we can estimate the duration of PCP necessary to obtain each Mg/Ca value based on estimated Mg/Ca0 to calculate the PCP-induced δ18O change. For example, if the duration for PCP is ~650 s (t1 in Fig. S10) and Ca2+initial is 2.73 mol/m3, an increase in Mg/Ca of the solution by ~15 mmol/mol (from Mg/Ca0 = 8.3 mmol/mol to the Mg/Ca value of the speleothem sample of 23.3 mmol/mol) is required. This corresponds to an increase in δ18O of ~4.1‰, which is applied as a correction to the δ18O speleothem record. The two model set-ups 1) Ca2+initial of 2.73 mol/mol3 and 2) Ca2+initial of 3.18 mol/mol3) lead to similar results.

Monte-Carlo Principle Component Analysis (MC-PCA)
To calculate the spatiotemporal coherence between ‘n’ records across Southeast Asia, an age uncertainty-aware Principal Component Analysis (PCA) is applied. The approach largely follows Deininger et al. (2017)48, using Monte Carlo sampling of age model realizations to adjust for the age uncertainties in all records71. Firstly, all speleothem age models have been recalculated in COPRA using the initial published U-Th ages. With COPRA, 2000 distinct time axes compatible with age uncertainty are sampled employing a Monte-Carlo sampling scheme and subsequently transferred into uncertainty of the proxy values58. This results in a single ‘error-free’ time axis and 2000 distinct proxy realizations which are compatible with the U-Th dating uncertainties. PCA requires constant sampling intervals in order to compute and decompose the covariance matrix of normalized proxy time series. Thus, all records are resampled with a constant sampling interval of 75 years and an overlap of 80% between consecutive windows to ensure an equal sampling resolution, adapting to the lowest resolved record. As proposed by Deininger et al. (2017)48, Gaussian Kernel-based smoothing is applied in this step to align the temporal variability between the differently resolved proxy records. Smoothing windows are centered at the resulting regularly sampled dates. Smoothed proxy realizations are z-normalized, yielding zero mean and unit standard deviation. Finally, PCA is applied to each of the 2000 sets of smoothed and z-normalized proxy realizations. Since every proxy realization can be regarded as a shifted and stretched/compressed version of the underlying ‘real’ proxy time series, each set of proxy realizations yields an individual characterization of spatiotemporal coherence represented by m ⩽ n modes. Each mode consists of a principal component (PC) and n EOFs (empirical orthogonal functions, also called ‘loadings’). EOFs are normalized such that they are equivalent to the linear correlation between the respective PC and a proxy realization. As discussed by Deininger et al. (2017)48, a smoothing and a flipping effect can be observed in the MC-PCA. In order to obtain meaningful estimates of the median and confidence limits for the PCs and EOFs, the flipping effect requires additional processing of the ensemble of results. Firstly, an intrinsic sign ambiguity exists for the singular value decomposition used to numerically solve the eigendecomposition of the data covariance matrix72. Small deviations in one of the ‘n’ proxy records can suffice to result in a ‘random’ sign flip of the resulting PCs and EOFs, irrespective of the MC-based procedure. This will in some cases result in a bimodal distribution of the EOFs for a given PC. Secondly, a sign flip can result from age uncertainty whenever a set of proxy realizations from the age model ensembles is analyzed that is comprised of anti-phased instead of phased signals as discussed in Deininger et al. (2017)48. In this work, k-means clustering is applied for each ensemble of PCs to distinguish flipped and non-flipped solutions73. In this context, another effect was observed when age uncertainty of any of the included proxy records is large: apart from the flipped and non-flipped solution for the 1st PC, a third distinct solution can be identified which exhibits strong similarity to the 2nd PC. This results from an overlap between the variance explained by each of the modes and may occur between any given pair of ith and i+1th PC if dating uncertainties are large. While in most cases (e.g. 1800 out of 2000) the 1st mode is represented by strongly correlated (flipped and non-flipped) PCs, in fewer cases variance may be maximized by what is largely identified as the 2nd mode. Such a scenario is identified successfully by the proposed k-means clustering of PCs with k=3 and should be considered in future application of the MC-PCA method. In this work, dating uncertainties were found to be sufficiently low such that this effect does not result in a significant number of ‘swapped’ modes. Consequently, EOFs and PCs of the flipped modes identified by k-means clustering are sign-flipped and a single median PC/EOF is computed from the respective ensemble of solutions for each mode.

Age uncertainty is propagated by computation of the 2σ-error bars of the respective PC and explained variance. The statsmodel package was used to compute PCA solutions for each set of proxy realizations and the sklearn package was used for k-means clustering. The MC-PCA is applied to δ13C and δ18O records respectively form Tham Mai Cave (Laos)43, Thien Duong Cave (Vietnam) this study, Tangga Cave (Sumatra)74, Bukit Assam Cave (Borneo)41. Additionally, a lake record from Lake Huguang Maar in China14 was included. However, COPRA is specifically built for speleothem records, therefore the original age model without age uncertainties in all 2000 realization was used. For the δ18O MC-PCA, Dongge Cave record was included instead of the Lake Huguang Maar lake record, but no δ13C record was accessible. For both applications, three modes were computed motivated by their interpretability.

Climate model simulations
For this work, three climate simulations were adapted to understand the forcing mechanisms of the mid-Holocene. For more details on the methodology used, refer to Pausata et al. (2017)52 and Gaetani et al. (2017)75 for EC-Earth, to Chandan and Peltier (2020)54 for the UofT CCSM4 experiments and to Tabor et al. (2020)53 and Thompson et al. (2021)76 for the iCESM simulations.

EC-Earth
Numerical simulations of the mid-Holocene were run using the EC-Earth version 3, with the atmospheric model based on the Integrated Forecast System (IFS cycle 36r4) (European Centre for Medium-range Weather Forecasts), including the H-TESSEL land model. The simulations were run at T159 horizontal spectral resolution, corresponding to ca. 1.125° by 1.125° and a vertical resolution of 62 vertical levels. As an ocean model version 3.3.1 of the Ocean General Circulation Model-NEMO was used77. The model solves the primitive equations discretized on a curvilinear horizontal mesh with a horizontal resolution of ca. 1° × 1° and 46 vertical levels. The model is coupled every model hour with the Louvain-la-Neuve Ice Model-LIM370 at the surface, having the same horizontal resolution as NEMO.

Two experiments were conducted with both set to an orbital forcing corresponding to values at 6 ka, which were computed internally following Berger (1978)78. Greenhouse gases accord with the PMIP3/CMIP5 protocol and other boundary conditions were set pre-industrial values following to the PMIP/CMIP5 protocol. The first experiment, the MH ORB experiment, uses the dust climatology of the MERRAero data as the Community Atmosphere Model, based on the long-term monthly mean (1980–2015)79, simulating a desert over the Sahara region. In the second experiment (MH GS), the vegetation type over the Sahara domain (11° to 33°N and 15°W to 35°E) is set to shrub and the dust amount is reduced by up to 80%52 simulating a vegetated Sahara, as proposed by paleoclimate data80. Changing vegetation cover from shrub (MH GS) to desert (MH ORB) corresponds to an increase in the surface albedo from 0.15 to 0.30 and a decrease in the leaf area index from 2.6 to 0.2. The dust changes from the MH GS to MH ORB experiment correspond to an increase in the global total aerosol optical depth of 0.0281. To avoid abrupt transitions in dust concentrations, a smoothing filter was used and the dust reduction was applied evenly above 150 hPa, because, at these heights, aerosol particles are uniformly distributed. No accurate vegetation reconstruction for the mid-Holocene is available at present, thus the model cannot provide a realistic reconstruction, but give insights on potential feedback mechanisms during this time.

University of Toronto CCSM4 (UofT CCSM4)
In addition to the EC-Earth, the University of Toronto version of the NCAR CCSM4 model (UofT CCSM4) is used to perform similar experiments82. Peltier and Vettoretti (2014)83 modified the ocean component of the model, namely the Parallel Ocean Program version 2, to be more appropriate for paleoclimate simulations. It was first applied by these authors to the study of glacial climates and subsequently used for studying the mid-Pliocene warm period in Chandan and Peltier (2017, 2018)84,85 and the mid-Holocene in Chandan and Peltier (2020)54. A detailed description of the modifications to the ocean component can be found in Peltier and Vettoretti (2014)83 and in Chandan and Peltier (2020)54. Similar to the EC-Earth experiments, a mid-Holocene simulation was performed using orbital parameters at 6 ka and concentrations of trace gases following PMIP4/CMIP6 protocol86. The UofT CCSM4 MH GS simulations follow the vegetation cover recommended by Otto-Bliesner et al. (2017)86 and Pausata et al. (2016)81, similar to the EC-Earth MH GS run. However, the Guineo-Congolian tropical forest belt was additionally extended northward. The change in surface albedo resulting from MH GS to MH ORB is 0.18 to 0.30, similar to the EC-Earth run.

iCESM
We use a water isotopologue tracking-enabled version of the Community Earth System Model (iCESM1;87), with a fully coupled configuration of iCESM1 with 1.9 × 2.5° atmosphere (CAM5) and land (CLM4), and nominal 1° ocean (POP2) and sea ice (CICE4). Similar to the other model runs, two experiments are performed:
A mid-Holocene configuration with orbital forcing corresponding to values at 6 ka and greenhouse gases (264.4 ppm CO2, 597 ppb CH4, 262 ppb N2O) as in PMIP4 (MH ORB; Otto-Bliesner et al., 2017), and a mid-Holocene configuration with prescribed Saharan vegetation (MH GS) (Fig. S11). In MH GS, vegetation coverage types come from the preindustrial default in CLM4 with agricultural land reverted to natural plant functional types. In MH GS, PI zonal average land conditions in Africa at ~11°N is extended, including vegetation/soil types and leaf/stem area indices, to all land in the Sahara and Arabian Peninsula. The addition of vegetation in the mid-Holocene reduces mean annual aerosol optical depth from 0.50 to 0.05 and surface albedo from 0.30 to 0.15 over the region of North Africa with modified vegetation. Further information about the simulation can be found in Tabor (2020)53.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_2": "Dr. Wolf et al. present a new continuous Holocene speleothem record (oxygen isotopes, carbon isotopes, and trace elements) from central Vietnam. They argue that this record can disentangle the large-scale circulation reflected by oxygen isotopes and local rainfall changes reflected by carbon isotopes and trace elements. They further use a Monte-Carlo principal component analysis to investigate the spatiotemporal changes in δ18Op and δ13C records across Southeast Asia. Besides, they use mid-Holocene climate simulation results to explore the potential mechanism of the mid-Holocene anomalies in δ18Op and δ13C observed from proxy records. They concluded that the synchronous rainfall change in Southeast Asia is forced by changes in the Pacific and the Indian Ocean, whereas the inverse winter and summer regional atmospheric circulation is controlled by seasonal insolation over the Northern Hemisphere.\nWhile I find the new Holocene speleothem record to be a useful indicator of the Northeast Winter Monsoon in Southeast Asia and potentially provides an opportunity to disentangle large-scale circulation and local rainfall, the discussion about the spatiotemporal change of monsoon and precipitation in Southeast Asia and the underlying mechanism of NEM change during the Holocene are not well supported and sometimes unconvincing. Below, I highlight my major concerns and suggestions.\n1. The features of δ18O, δ13C, and trace elements time series described in the paper are not very clear to me. For example, I think the “lower δ18O values during the mid-Holocene and higher values in the early and late Holocene” (Line 158-159) are not evident, given the combination of changes in high-frequency variability and millennial-scale changes in mean value. I would suggest the authors use statistical methods to identify the timing of step shifts in δ18O, δ13C, and Mg/Ca. Some potential options include the environmental time series change point detection (EnvCpt) method (Beaulieu and Killick, 2018, Journal of Climate), and regime shift detection (Rodionov, 2004, GRL). These methods can detect the change point in a time series based on changes in its mean value and variance.\nBesides, I think the cubit fitted curves of the δ18O record from TM (Fig. 2E) and TD (Fig. 2F) could be problematic as they extend outside the range of the original data (e.g., 0-1ka),\nThe mismatch between δ18O and δ13C and trace elements in the central Vietnam speleothem record is the main motivation of this work, and the different long-term patterns in δ18O records across the Southeast Asian revealed by the cubic fit is key evidence for the inverse relation between SWM and NEM. So, I think it is important to use appropriate statistical methods to detect shifts and the timing of shifts in the δ18O, δ13C, and trace elements time series.\n2. The mechanisms of the orbital-scale reverse response in SWM and NEM are unclear to me. The authors found positive loading on the PC1 of δ18O in southern China and Sumatra and negative loading on the PC1 of δ18O in Borneo and Vietnam. Based on the significant correlation between boreal summer and autumn insolation, they conclude that “Holocene insolation changes are the primary drivers of the long-term trend of δ18O speleothem records across Southeast Asia” (Line 239-242). However, the significant correlation between PC1 and autumn insolation is insufficient to prove NEM recorded in Vietnam and Borneo mainly follows autumn insolation because the temporal evolution of the PC1 of δ18O (long-term increasing trend from 8ka to 1ka) doesn’t capture the main features in the δ18O records from Vietnam and Borneo (more depleted δ18O during the mid-Holocene). In addition, the authors didn't clearly explain how summer and autumn insolation would influence the SWM and NEM, respectively.\n3. It’s very interesting to see SWM and NEM show reverse response (reflected by the PC1 of δ18O) on orbital timescale, whereas coherent rainfall pattern across Southeast Asia was found on millennial timescales (reflected by the PC1 of δ13C). It’s worth to further investigating the underlying mechanisms of the different regional circulation (δ18O) and local rainfall (δ13C) on different timescales. But I don’t think the mid-Holocene time-slice model simulation results alone are sufficient to explain the difference in the response of δ18O and rainfall on different timescales.\n4. It’s not clear to me why the authors use “MHORB-MHGS” to discuss the mechanisms of δ18O and precipitation changes during the mid-Holocene observed from proxy records.\nThe MHORB simulations incorporate mid-Holocene insolation, greenhouse gases, vegetation cover, and dust load, and MHGS simulations have a prescribed vegetated ‘Green Sahara’ (GS), additionally” (Line 287-290). It’s not clear to me what climate forcings would “MHORB-MHGS” indicate and why the authors use the difference between the ORB and GS runs to identify the mechanisms driving mid-Holocene anomalies. Another issue related to this is the argument about reduced Saharan vegetation could contribute to the mid-Holocene lower average rainfall observed in proxy records (Line 331-338, Fig. S9B). It’s not clear to me why “MHORB-MHGS” can specifically reveal the influence of Saharan vegetation. I would suggest the authors use the MHORB minus preindustrial (or early Holocene) to investigate the contribution of orbital forcing on mid-Holocene δ18O and precipitation anomalies, and (MHGS - MHORB) minus preindustrial (or early Holocene) to evaluate the influence of Green Sahara.\nOther comments:\nFig. 1: (1) Please change “mean daily precipitation per month [mm/day]” in Fig. 1A and 1B to “mean monthly precipitation rate [mm/day]”. (2) Please change TH to TD (Line 108).\nFig. 2: Please use appropriate statistical methods to detect shifts in δ18O, δ13C, and Mg/Ca. Please make sure the cubic fit in Fig. 2E and 2F can reveal the long-term trend in δ18O from TM and TD, respectively.\nLine 162-166: Please show changes in autumn and summer insolation in Fig. 2 to back up this statement.\nFig. 3: (1) Please add the time series of PC2 of δ18O orig, δ18O corr, and d13C. Please reverse the color bar for EOF, so that the warm color indicates positive loading, and the cold color indicates negative loading. (2) In Fig. 3A-3F, please change “PC1” and “PC2” to “EOF1” and “EOF2”.", "reviewer_3": "The authors present a new isotope and trace element record from speleothems from Vietnam. They claim that it is the first reconstruction of winter monsoon for the Holocene period, which is important for understanding the continent's paleoclimate. The data has high quality, and the comparisons with climate models were also addressed to reinforce their interpretation. The paper is well written, and this novel has the potential to be published in Nature Communications. However, I see a major problem with their correction on the d18O record. This correction compromises the quality of the data and the interpretation on the manuscript. See my concerns about this correction below. Apart from this, I have only minor comments to improve the manuscript.\n1) The correction is based on the approach presented by Deininger et al. (2021), which shows that PCP can also affect the d18O in stalagmites. This is something that so far has been ignored in the literature. However, the implication of the PCP on d18O presented in Deininger et al. is mainly theoretical and still needs further validation before being applied to speleothem data records. Furthermore, no calibration studies with speleothem were done before.\n2) For this correction, the authors assume temperature and CO2 in the soil to be constant during the Holocene period, which is unrealistic. Temperature changes were small but not constant. Slight elevations of temperature associated with drier periods can significantly increase evaporation with all associated effects on the isotopes. But the big problem is with regard the CO2. The monitoring of soil CO2 over two years already shows variability between 430 ppmv and 1580 ppmv. This variability was probably more significant during the Holocene. Furthermore, long periods of drier or wetter conditions can affect the vegetation and soil above the cave, changing the CO2 concentration in the soil, which is the primary source of the cave's CO2. The proxies presented in this paper present high variability, indicating significant rainfall variability above the cave, which will also affect the CO2 concentration in the soil.\n3) The authors use d13C and Mg as a “pure” proxy for PCP, which is unrealistic. I agree that both are proxies for PCP, but PCP is not the only force on these proxies. In using Mg to correct d18O, the authors are not only removing the PCP signal but also adding noise to the d18O signal.\n4) In this correction, the author assumes that all correlation between d18O and Mg and d13C is due to PCP. However, other factors provide a similar correlation between these proxies other than PCP. Therefore, in removing this dependency from the d18O, the author also removes additional information. For example, during a strong monsoon, the “amount effect” provides more negative values to the d18O, but the high rainfall during this event also decreases the PCP leading to lower values of Mg and d13C. Therefore amount effect and rainfall variability also drive d18O, d13C, and Mg in the same direction. Therefore, using the Mg/Ca and correcting the d18O, the authors are not just removing the PCP effect on the d18O but also the amount effect.\n5) The authors argue that they are removing the effect of PCP on the d18O. But the corrected record has much higher variability than the original record, making me believe that this \"correction\" adds more variability and noise to the records than filter PCP.\nMinor comments:\nFig.1: How the amount-weighted d18O was calculated? Were monthly averages used?\nFig2: I suggest including the name of the caves in the figure or at least on the captions. It isn't enjoyable to go to another figure to see the name of the records presented in this figure.\nFig. 2: it is difficult to distinguish the scales at the y-axes between D and E.\nThe authors also measured the Sr/Ca and Ba/Ca with laser ablation. Why were these data not shown in the paper?\nLine 258: The authors say “This suggests that the first PC derived from stalagmite δ13C records and a lake record from Southeast Asia likely share a common forcing throughout the Holocene...”.\nAccording to the authors, d13C in the stalagmites is connected to the local hydrology due to the PCP process. However, PCP is absent in lakes. d13C in lakes is related to the isotopic composition of the organic matter/vegetation. If d13C from stalagmite and lakes correlates, this is suggestive that the d13C from speleothems is also affected by the d13C from the vegetation (something ignored in this manuscript) and not that the d13C from the lakes is affected by PCP, which is impossible. How do the authors explain it?\nBest regards,\nValdir F. Novello"}, "author_response": "We thank the reviewer for the suggestion; we have updated the cubic fits using the polyfit and polyval functions in Matlab to ensure a more robust fitting to the time series boundaries. Additionally, we performed several change point and trend analyses which have been added as a figure to the supplementary material (Fig. 2 and 3 shown below; Fig. S6) using EnvCpt and BEAST. Both approaches support our interpretation and show that indeed TD3 δ18Ocorr values were lower during the mid-Holocene and higher during the early and late Holocene, as shown in panel 3 of Fig. 1 and 3 in this document. They also show that there are significant shifts towards increased values in the δ13C and Mg/Ca during the onset and end of the drought.\n[...]\nTo address Reviewer#2 concerns and improve clarity, we have added a new panel to Fig. 3 (G), which shows the newly preformed cubic fits of δ18O records used in the PCA, and compare them with insolation closest to the record’s location, considering the main season for peak rainfall. The summer monsoon records (Dongge and Tham Duon Mai) show similar trends to JJA insolation, while Thien Duong and Bukit Assam records share most similarities with SON insolation. The record from Tangga cave matches best with mean JJA and SON insolation, potentially due to the influence of the Australian monsoon. Here, we propose that an increase in insolation is strengthening the rainfall and associated winds, leading to longer rainout and more depleted δ18O values reaching the cave sites. The anti-correlation of the δ18O in rainfall associated with summer and winter monsoon thus results from the anti-correlation between insolation at these locations. We have added the following text to the main manuscript (Lines 218-230):\n[...]\nWe generally agree with the reviewer, however, isotope-enabled transient model simulations are, to our knowledge, not yet available for the time period covered by our proxy record. We would also like to highlight that the mid-Holocene simulations were not aimed at understanding drivers for the entire Holocene, as stated in lines 285-291: “Crucially, the dry conditions observed in the proxy records between ~5 to ~3 ka are evident in both PC1 and PC2, suggesting strong coherence between the sites during this time (Fig. S7). […] In the next step we will investigate the forcing mechanisms of this dry period and their effects on the regional and local aspects of the Southeast Asian monsoon.”\n[...]\nThe MHORB has the vegetation prescribed as pre-industrial as is in the standard PMIP protocol; however, such vegetation is unrealistic to represent the Green Sahara. On the other hand the MHGS include a more realistic representation of the Green Sahara. Hence, we use the MHORB-MHGS to isolate the effects of vegetation and dust forcing during the mid-Holocene and investigate the transition period leading to a dry Sahara but still with higher summer insolation than the pre-industrial. By plotting the difference between the MHORB run (standard protocol; pre-industrial Sahara) and MHGS (changes in vegetation and dust loading), we can isolate the effect of dust loading only on the spatial pattern in δ18O.\n[...]\nTo our knowledge, there is indeed no calibration study available, which investigated a direct link between δ18O and Mg/Ca during carbonate precipitation. However, in individual studies it was shown, that δ18O values are increasing with ongoing precipitation and there is a common agreement, that the Mg/Ca ratio is doing the same. See, e.g., the recent study by Stoll et al., in Climate of the past, which is still under review, but freely available), in which they also apply a first step to correct for the effect of PCP. We also oppose the comment of the 'mainly theoretical' nature of this effect. This process has been shown to be important in many laboratory experiments (e.g., Day and Henderson, 2012, Polag et al., 2010, Hansen et al., 2019, Mickler et al., 2004), where the stable isotope composition is always increasing along the path of carbonate precipitation.\n[...]\nWe are grateful for this suggestion and have conducted further simulations for our time series using varying parameter settings, including a range in cave air temperatures based on reconstructions from Borneo (Løland et al., 2022), which identified a range of ± 2°C in the cave over the last 8 kyrs. We also add a range of CO2 estimates to simulate potential changes in pCO2 ranging from 260 to 3000 ppmv. In the table below we list the effect of these parameters on δ18O during 6-4 ka. This period is likely to have seen the highest amounts in PCP over the last 8 kyrs. We find that changes in cave air temperature have only small effect on mean δ18O (± 0.05‰), whereas cave air CO2 variability affects the CO2 up to ~-1‰ (see table below). We consider our initial simulation with 426 ppmv and 21°C the ‘main’ simulation and the 12 other the ‘uncertainties’.\n[...]\nWe agree with the reviewer, but would like to highlight that we only focus on the Holocene trend and do not use / interpret the short-term variability of the record, due to the limitations of our PCP correction. Therefore we focus on the local proxies (δ13C and Mg/Ca) and the PCA simulations for our mid-Holocene interpretation. Our statistics (PCA and trend analysis) aim to filter any noise our uncertainties related to the local proxies. Further, we base our conclusions on a suite of 5 speleothem records from SE Asia to provide a more robust interpretation for our findings.\n[...]\nThere is strong evidence that the amount effect is not controlling δ18O at our cave sites (Wolf et al., 2020, Yang et al., 2016, Pausata et al., 2011). Data from the local GNIP station shows that there is no correlation between monthly δ18O and rainfall amount (for the last 5 years; Fig. 1 main manuscript; Wolf et al., 2020). Instead, it has been shown that the reversal of winds, changing the location of moisture source, drives the seasonal and interannual variability of δ18O in central Vietnam. Thus, a correlation between δ18O and δ13C /Mg based on the amount effect is not to be expected in our record.\n[...]\nAgain we would like to highlight that we are not using variability or the amplitude of individual peaks in our interpretation, but rather focus on the long-term trends, which is extracted by the cubic fits and PCA. Therefore, despite the limitations of our corrections, such as estimated cave parameters, we have confidence in our conclusions on the co-evolution of δ18O in SWM and NEM rainfall in Southeast Asia.\n[...]\nWe apologise for this confusion. The lake record is using the Ti content as a proxy for winter monsoon strength not δ13C, noted in lines (235-236 and 256-257). However, we recognise that the figure title was misleading and have renamed it to δ13C + Ti."}, {"round": 2, "reviewer_comments": {"reviewer_2": "Dr. Woft et al. have addressed most of my main concerns that were detailed in my earlier review. In particular, I think the incorporation of trend/change point detection and the comparison between the cubic fits for d18O records and seasonal solar insolation strengthen the manuscript. However, there are still some aspects that are unclear from my perspective.\nThe simulated influence of the Green Sahara on rainfall and d18O, and its comparison with proxy records, is still not clear to me. The authors have clarified that MHORB has prescribed pre-industrial vegetation, while MHGS uses a more realistic representation of the Green Sahara. Therefore, I assume that \"MHORB - MHGS\" represents the influence of \"reduced Saharan vegetation and increased dust.\" If that is the case, Fig. 6 (PRECT: MHORB - MHGS) demonstrates that the reduced Saharan vegetation and increased dust during MH would lead to drier conditions in Southeast Asia, as described in Lines 355-362. In other words, the mid-Holocene Green Sahara (a vegetated and less dusty Sahara) would lead to wetter conditions in Southeast Asia, which seems to contradict the widespread mid-Holocene drought in paleoclimate proxy records (EOF1 and PC1 of d13C+ Ti). Addressing this point is important for establishing the linkage between the Green Sahara and monsoon rainfall in Southeast Asia.\nRelated to the above point, please further explain how the spatial pattern of MHORB-MHGS d18O shown in Fig. 5B supports the different Holocene trends in SWM and NEM d18O (Lines 324-331). The MHORB-MHGS δ18O values increase only in northern MSEA and Sumatra, corresponding to reduced Sahara vegetation and increased dust (MHORB-MHGS). But as shown in Fig. 3G, all proxy records across Southeast Asia show lower d18O during the mid-Holocene relative to preindustrial.\nThe discussion on the linkage between insolation and changes in SWM and NEM has been improved in the current version. I would suggest the authors further test the mechanism of the insolation influence on SWM and NEM by comparing MHORB d18O and preindustrial d18O in iCESM (isolating the influence of solar forcing).\nOther comments:\nLine 227-230: More evidence is needed here to support the claim that changes in solar insolation during the Holocene can influence d18Op in the SWM and NEM by altering moisture source and rainout along the trajectory.\nLine 242: Please delete the quotation mark at the end of this sentence.\nLine 249: Please add \"(LH)\" after \"Lake Huguang Maar.\"\nLine 254: I would suggest replacing \"Information contained in the remaining trends (PC2)\" with \"The overall positive loading on the PC2 of d18O\" to make this sentence clearer.\nLine 257-258: Please change “Fig. 3E” to “Fig. C4” or change “PC2” to “EOF2”.\nLine 316: Please add \"during the Mid-Holocene\" to the end of this sentence to avoid potential confusion.\nLine 372: I would suggest changing \"accurate\" to \"reliable.\"\nFig. 4:\n(1) Please correct the title of Fig. 4C to \"PC2 of d18O and PC1 of d13C.\"\n(2) Please explain the meaning of the dashed lines in Fig. 4A and 4B in the figure caption.\n(3) Please add \"of d18Op\" after \"The uncorrected PC2.\"", "reviewer_3": "I'm still skeptical about the correction on the d18O values. But the authors presented good arguments to justify these corrections in the manuscript. Furthermore, all my other minor suggestions were incorporated. Therefore, I recommend the publication in Nature Communication without further changes.\nBest regards,\nValdir F. Novello"}, "author_response": "Thank you for this comment! Yes, it is correct that MH -MH represents the reduced Saharan vegetation and increased dust at the end of the MH, which lead to drier conditions in Southeast Asia. Thus, our simulations (Fig. 6) show that at the end of the Green Sahara increased dust load led to drier conditions across MSEA. This is well represented via the widespread drought (ca. 5-3.5 ka BP) in our proxy records in Figs. 4A and 4B and S2. Both the mean age of the onset of the GS termination (Fig. 1 below) and the start of the pronounced mid-Holocene droughts observed in our proxy data start around 5 ka. Here, we compare the conditions after the GS around 5 ka in the simulations to the drought also starting at 5 ka in our proxy records. We apologize for the confusion and have updated the main text as: ”Here, we compared the changes in the local proxy record before (8-5 ka BP) and after (5-3.5 ka BP) the end of the Green Sahara, where all proxy records used in the MC-PCA suggest drier conditions during 5 to 3.5 ka BP.” Lines 347-349.\n[...]\nThere was an incorrect label for the figure reference in this section. We apologize and have changed Fig. 3G to Fig. 3A and C in Line 324. Here, we aimed at comparing the spatial pattern: negative loading on Vietnam and Borneo vs positive loading on China, Laos and Sumatra, compared to the model simulation 18Op SON, which has a similar spatial pattern. Plotting MH -MH removes the effects of insolation, thus we cannot compare the temporal trends in the proxies to the simulation and instead compare the spatial pattern of the simulations to the spatial pattern in the timeseries (Fig. 3A and C in the main text). Since the end of the GS causes an increase in 18O over northern MSEA and Sumatra only we conclude that the pattern of an inverse relation between SWM and NEM is enhanced by the dust loading via an additional increase in 18O for SWM records, while NEM is at its lowest values of the Holocene.\n[...]\nThank you for this comment. Please see the figure below, which was added to the supplement as Fig. S11. We can see that by considering the changes in vegetation and dust the differences in isotopic composition between the mid-Holocene and pre-industrial are enhanced and in better agreement with proxy records than when considering orbital forcing only. Unfortunately, not all proxy records extent towards modern times and thus we cannot directly compare the simulated changes with the paleoclimate archives.\n[...]\nOur simulations do show a change in 18Op in the SWM and NEM (Fig. 5, main text); however, our simulations do not include moisture tracking, so we cannot verify it through our experiments. Nevertheless, the influence of the insolation on speleothem 18O has been shown in many speleothem records from East Asia and other monsoon regions, as well as that the effects of moisture source and rainout dynamics control modern 18Op in MSEA, please see references 22, 23, 58 in main text and other studies (e.g. Zhang et al., 2018; Hu et al., 2019 (using iCESM and speleothems); He et al., 2021). The role of moisture source dynamics affecting 18Op in iCESM has also been identified for the mid-Holocene in Tabor et al. 2020 (62) and Thompson et al. 2021 (83). Thus, we propose that moisture source dynamics and rainout changed based on insolation, which is reflected in our proxy record. Such changes are further enhanced by vegetation and dust changes which further alter the intensity of the Southeast Asian monsoon system. We agree that additional studies are required to back up this statement, however, this is beyond the scope of our paper."}]

Decision letter