Modern anthropogenic drought in Central Brazil unprecedented during last 700 years

Domain: earth_science
Published: 2024-02-26
Editorial outcome: minor_revision

Abstract
A better understanding of the relative roles of internal climate variability and external contributions, from both natural (solar, volcanic) and anthropogenic greenhouse gas forcing, is important to better project future hydrologic changes. Changes in the evaporative demand play a central role in this context, particularly in tropical areas characterized by high precipitation seasonality, such as the tropical savannah and semi-desertic biomes. Here we present a set of geochemical proxies in speleothems from a well-ventilated cave located in central-eastern Brazil which shows that the evaporative demand is no longer being met by precipitation, leading to a hydrological deficit. A marked change in the hydrologic balance in central-eastern Brazil, caused by a severe warming trend, can be identified, starting in the 1970s. Our findings show that the current aridity has no analog over the last 720 years. A detection and attribution study indicates that this trend is mostly driven by anthropogenic forcing and cannot be explained by natural factors alone. These results reinforce the premise of a severe long-term drought in the subtropics of eastern South America that will likely be further exacerbated in the future given its apparent connection to increased greenhouse gas emissions.

Introduction
Since the 1970s, observational records and models point to a widespread drying trend affecting many land areas across the globe1–5. Drought events impact societies and their countries’ economies in multiple ways, leading to agricultural losses, reduced energy production from hydropower, and affecting public water supply systems6,7. Intensive agriculture is highly dependent on water availability and, usually, is the first and most heavily affected sector impacted by drought, with important implications for food security. Accurate projections of future drought risk are crucial for policy formulation and implementation of adaptation and mitigation strategies and for effective disaster and risk reduction policy6. Regions within southern Amazonia and central–eastern Brazil are now characterized by a long-term hydrologic drought, which appears to be highly correlated with changes in surface temperature5,8–10. Central Brazil plays an important role in the country´s economy, due to its intensive agriculture, but it is also home to important ecosystems. The increased aridity observed in the region during the last decades has raised questions regarding the role of human-induced warming in driving changes in the hydrologic balance and how drought risk might be exacerbated by anthropogenic greenhouse gas emissions in the coming decades2,3,9,11.

In general, a drought can be defined as a period with a hydrologic deficit, when the water evapotranspiration exceeds the input from effective precipitation4,12. This type of drought is also commonly known as “hydrologic drought” and its effects are most apparent in reduced streamflow. Hence, in its simplest form, a hydrologic drought can be represented by a simple hydrologic balance model based on precipitation (P) and potential evapotranspiration (PET) according to the following formulation: Hydrologic Balance = P-PET1,13,14. Future projections of drought risk, however, are a major challenge for the climate community. At the local scale, temperature plays an important role in the hydrologic balance, with warmer temperature leading to an increase in evapotranspiration and a decrease in soil moisture. On a larger scale, changes in sea surface temperature also influence land precipitation, driving hydroclimate variability at interannual to multidecadal timescales15–17.

Reconstructing past climatic conditions in the region, prior to the period covered by instrumental records, is vital to document the full amplitude of natural climate variability and to quantify the contributions from external and internal forcings to regional hydrologic variability11. Speleothems can be very useful in this regard, providing high-resolution hydroclimate reconstructions with accurate age control spanning decades to hundreds of thousands of years. Drought events can affect the chemistry and the isotopic composition of speleothems in well-ventilated cave rooms. In such rooms, isotopic and geochemical exchange between drip water and cave atmosphere allow reconstructing past variations of the atmospheric conditions associated with local potential evapotranspiration (PET) and precipitation variability18–20. However, few paleoclimate studies have investigated the potential of speleothems from well-ventilated cave rooms to record past variations of the local environmental conditions21. Speleothem δ18O records from tropical regions have been used to reconstruct past monsoon activity and rainfall amount22 by relying on speleothems collected in isolated chambers deep inside caves, to avoid any atmospheric influence from outside the cave. However, in ventilated cave environments, the loss of water by means of evaporation during the dripping process significantly affects the oxygen isotopic composition (δ18Ocalcite)23. Hence, both a reduction in the amount of rainfall and an increase in evaporation will result in higher δ18Ocalcite values of the speleothem. Similarly, dry conditions related to reduced cave dripping and atmospheric relative humidity affect speleothem carbon isotopes in the same direction, resulting in higher δ13Ccalcite23,24. Reduced precipitation and decreasing humidity of the cave atmosphere also affect the trace element composition of the speleothems by inducing the so-called “Prior Calcite Precipitation” (PCP) effect. PCP simply refers to the fact that carbonates will precipitate before the drip water reaches the stalagmites25. The PCP increases under dry conditions, and the process leaves the remaining solution enriched in both 13C and in incompatible trace elements24–28. As a result of the PCP, the δ13C values and trace element ratios, such as Mg/Ca, Sr/Ca, and Ba/Ca, tend to increase in the forming stalagmite during dry periods. When combined, oxygen isotope, carbon and trace elements from near-cave-entrance speleothems provide a powerful tool that allows reconstructing variations of the hydrologic balance (precipitation-evapotranspiration) back in time, beyond the instrumental period, thereby yielding unique, annually resolved records of tropical hydrologic balance (e.g., Bernal et al.29).

To put currently observed changes in the hydrologic regime in a long-term context, here we present a multi-proxy reconstruction of annually resolved oxygen and carbon isotopes (δ18O and δ13C) and trace elements concentrations in speleothems from a cave located in central–eastern Brazil. Our results provide robust evidence that the current drying trend has no precedent in the last seven hundred years. By applying a detection and attribution approach, we demonstrate that anthropogenic forcing is the main driver of the current hydrologic regime, responsible for the observed hydrological deficit in the region through enhanced evapotranspiration combined with a negative precipitation trend.

Results
The study site is located in central–eastern Brazil, in a tropical Savannah-like biome in the northern part of Minas Gerais State (Fig. 1). The climate in this region is characterized by a strong precipitation seasonality, with an annual mean rainfall of near 950 mm. The wet season occurs during the South American Summer Monsoon (SASM) season, between October and April, and is followed by a pronounced 5-month dry season from May to September. The dry season accounts for only 3% of the total annual precipitation and can last as long as 6–7 months in some years based on historical records from local stations available in the Instituto Nacional de Meteorologia (INMET) in operation since 1915. Due to its prolonged dry season, central Brazil is highly dependent on a reliable wet season and accordingly very vulnerable to variations in precipitation.

To characterize drought variability in central–eastern Brazil we used the instrumental records of precipitation and temperature from local meteorological stations to calculate the P-PET index. The meteorological stations were compiled over an area of nearly 4°×4°, centered over the cave site (Supplementary Table S1). Local precipitation was calculated as the mean of monthly rainfall from all meteorological stations listed in Supplementary Table S1, obtained from the INMET and the Agência Nacional de Águas (ANA). The station-based precipitation data were compared with gridded precipitation data at 1° resolution from the Global Precipitation Climatology Centre (GPCC), version 2018 (Supplementary Fig. S1). The potential evapotranspiration (PET) was calculated based on the FAO-56 Penman–Monteith equation using regional gridded data from ground-based stations provided by Brazilian Daily Weather Gridded Data (BR-DWGD)30 (Fig. 1). In addition, we used measurements of potential evaporation from Piche evaporimeters from local meteorological stations to assess the local atmospheric evaporative demand. These records were used for comparison with speleothem proxies that are sensitive to local atmospheric evaporative conditions.

Since streamflow integrates the net hydrologic balance over a wider catchment area, we compared the calculated P-PET index with the z-scored streamflow data of local rivers obtained from the ANA (Fig. 1 and Supplementary Table S1). This approach permits us to assess how variations in precipitation and evapotranspiration demand relate to streamflow on a regional scale (Fig. 1). The average streamflow was calculated based on the mean monthly values of each fluviometric station. To determine the onset of the drying trend, we applied the change point detection algorithm31 to the streamflow time series (see “Methods”) (Supplementary Fig. S2). Next, we calculated the streamflow z-score of each station prior to averaging the monthly data. Finally, the reconstructed hydroclimate variability obtained from our speleothem record was compared with a gridded surface temperature reconstruction from the Paleo Hydrodynamics Data Assimilation (PHYDA) product32. PHYDA combines annually resolved proxy data with the physical constraints from the Community Earth System Model–Last Millennium Ensemble (CESM-LME) during the past 2000 years. Here we compared the PHYDA regional surface temperature (12–16° S; 41–47° W) with local observational data from the INMET meteorological stations. The two datasets show a strong coherence, and are highly correlated (r=0.73, P value<0.05) (Supplementary Fig. S3).

To reconstruct the local evapotranspiration demand we combine two speleothems from near the entrance of the Onça Cave (Fig. 1). This cave is located in the karst valley formed by the Peruaçu River (15°05'19.86“S–44°14'41.39“W), one of the few permanent tributaries of the São Francisco River, which is the fifth largest river of Brazil, extending along 2860 km north-south along central–eastern and Northeast Brazil.

The Onça Cave lies at the bottom of a karst canyon, below ~100 m of limestone, which serves as a water reservoir for the cave dripping water. The study site is located in a protected area under a karst relief characterized by cliffs, areas of exposed epikarst and dolines with steep-sided depressions, which effectively prevents the expansion of farms and agricultural activity in areas surrounding the cave. Therefore, land-use change is unlikely to have played a significant role in affecting the climate surrounding the cave or the geochemical composition of the speleothems. The cave displays a linear plan pattern and a wide-open entrance lying at the bottom of a steep-sided hill. The entrance measures about 50 m in length and nearly 10 m in height, giving access to the main cave room (Supplementary Fig. S4). The atmospheric conditions at the sampling site are typical of a well-ventilated environment, characterized by large seasonal amplitudes in relative humidity and temperature, reflecting the atmospheric and environmental conditions outside the cave (Supplementary Fig. S5). Variations in Onça Cave relative humidity reflect the relative humidity and evaporative potential from the exterior environment, as documented by the comparison between in-situ measurements inside the cave and data from a nearby local meteorological station (Supplementary Fig. S5). In the main cave room, the relative humidity varies from ~50 to ~100%, while the cave temperature oscillates between 17 and 25°C, with low and high values during winter and summer, respectively (Supplementary Fig. S5).

The record from Onça Cave is composed of two speleothems, named Onça2 and Onça4 (Supplementary Fig. S6 and Supplementary Table S2), providing an annually resolved chronology covering the last 718 years spanning from 1298 to 2016 CE. The chronology of Onça2 was constructed using a layer counting method, validated against 19 U-Th dates ranging from 2015±3 to 1760±4 years CE (“Methods” and Fig. S7). Conversely, the chronology of Onça4 was constructed based on U-Th ages, integrating a set of 19 ages ranging from 1852±2 to 1298±3 years CE (“Methods” and Supplementary Fig. S7). To reconstruct long-term changes in the regional hydroclimate and changes in the evaporative demand we combine trace elements ratios (Mg/Ca, Sr/Ca, and Ba/Ca) with oxygen and carbon isotopes. The role played by the evaporative demand in determining variations of the oxygen isotopes was assessed by comparing the δ18O of speleothems dripping from five different points with the isotope composition of the local rainfall. The comparisons allow us to detect non-equilibrium isotope fractionation related to evaporative effects in the cave environment (see “Methods” and Fig. 2). In addition, we performed an experiment using artificial substrates to analyze the seasonal geochemical effects of temperature and relative moisture on speleothem calcite from Onça cave. Five dripping sites (P1–P5) were monitored and sampled every month from July 2018 to October 2019. The monitored parameters include temperature, relative humidity, and calcite deposition rate (see “Methods” for details) (Supplementary Fig. S8). In addition, we tested the relationship between the geochemical proxies with instrumental timeseries from local meteorological stations33–35 by estimating a scaling factor using a least-squares regression between the geochemical proxies and the instrumental records that permits us to assess the relative contribution of the related environmental forcing to the proxy signal by testing the linear similarity among the time series (see “Methods”).

The significant negative trend seen in the local streamflow since the 1980s (Fig. 1b) points to a prolonged hydrologic drought prevailing over central–eastern Brazil. In 2010 the local streamflow dropped to its lowest levels since the start of the monitoring around 1940. As shown in Fig. 1, the good correspondence between the negative trend observed in the streamflow of local rivers and the regional precipitation indicates that rainfall plays an important role in contributing to the hydrologic deficit. The persistently low precipitation since the 1990s, has contributed to increased hydrologic stress in recent years. However, precipitation is not the only determining factor for hydrologic drought. Indeed, the negative trend recorded in the z-score values of the streamflow (slope: −0.054±0.007) is comparatively higher than the one observed in the precipitation time series (slope: −0.039±0.013). Conversely, the calculated P-PET index shares nearly the same trend as the streamflow (−0.052±0.011) (Fig. 1b). The significance of the regression coefficients calculated for each variable are confirmed by the F-test, presenting P values < 0.001. From 1979 to 2016, streamflow of local rivers decreased by 20% per decade in terms of m3/s. During the same time period local precipitation decreased only 7% per decade (70 mm.decade−1) while the P-PET presents a decreasing trend of 18% (125 mm.decade−1), indicating that evapotranspiration changes represent an important component of the current drought. These comparisons suggest that the negative trend observed in the river discharges starting in the 1980s is partially a result of the increased evapotranspiration demand with increasing local temperatures. A recent detection and attribution study from Chagas5 shows that the reduced streamflow over central–eastern South America is aligned with decreasing P-PET. Furthermore, a strong trend towards aridification over central–eastern South America after 2010 is also consistent with model projections that point to a severe decrease in soil moisture, even under the moderate RCP 4.5 emission scenario, as a consequence of increased evapotranspiration3,5,36. Yet, we cannot completely rule out a contribution from land-use changes to the rapid decrease observed in regional streamflow. As pointed out by Chagas5, the accelerated reduction of streamflow over central–eastern Brazil from the year 2000 onward point to intensive agriculture as a relevant driver of the water deficit.

The long-term trends in δ18O and δ13C recorded in the Onça Cave speleothems bear a close resemblance to the observed increase in aridity documented in the instrumental record (Fig. 2), suggesting increasing temperature and local evaporative conditions as important drivers for the decadal variability recorded in the speleothem geochemical proxies from Onça Cave (Fig. 2c). For oxygen isotopes, temperature and relative humidity represent competing forces that drive values in opposite directions23,37–39. In remote and isolated cave chambers, the constant temperature and high and stable relative humidity prevent substantial isotopic changes, allowing the speleothems to provide reliable records of rainfall isotope variability38,39. However, an increase in temperature during calcite deposition leads to a reduction of the δ18O values of the forming calcite due to a reduction of the calcite-water isotope fractionation factor (αcalcite-water). Empirical measurements point to a rate of approximately −0.2‰/°C40–42. However, our δ18O record presents a trend in the opposite direction, with a ~0.7‰ increase, despite a warming of about 2° since 1970. Conversely, an increase in evaporative conditions also increases the δ18O values of the calcite formed in speleothems. Thus, the trend seen in the Onça2 δ18O values after the 1970s is primarily a result of kinetic fractionation related to the highly evaporative conditions and not a direct consequence of temperature variations. This assessment is supported by the isotope monitoring of Onça Cave dripping water, which clearly shows that the δ18O and δD values of the drip water plot on an evaporation line (Fig. 2). As can be observed in Fig. 2a, the δ18Owater and δDwater values of the dripping solution from Onça Cave define a meteoric water line with a slope of about 5.4, which is significantly less than the meteoric water line of the local rainfall (slope=7.4), suggesting that evaporation is affecting the isotopologues of the cave dripping water. Under such conditions, the speleothems potentially reproduce the relative variations of the evaporative demand, recording remarkably high positive excursions during periods characterized by a water deficit. The dripping rate also plays an important role in modulating the isotope effects produced by changes in the cave atmosphere environment. For instance, the evaporative effect is strongly controlled by the dripping rate, since the time of re-equilibration between oxygen isotopes and cave atmosphere controls the extent of the isotope fractionation effect23,43,44. Thus, as the dripping rate decreases, the longer exchange time between dripping solution and cave atmosphere enhances the evaporative effect on δ18O and δ13C values of the speleothem23.

Finally, the increase in the δ18O values of the Onça Cave speleothem could also result from a decrease in the local rainfall amount45–47. However, the positive trend observed in the δ18O from Onça2 stalagmite starting around 1970 (Fig. 2) is too high to be explained by rainfall oxygen isotope variability alone. Indeed, estimates of changes in rainfall oxygen isotope values based on monitoring studies47 document a positive δ18O trend of 0.007±0.003 ‰.yr−1. Meanwhile, the Onça2 stalagmite shows a trend that is twice as large (0.015±0.002 ‰.yr−1) (Supplementary Fig. S9).

In general, the trace elements also share a common trend with the isotope record, characterized by an increase in the elemental ratio values since the 1970s (Fig. 3). These trends are opposite to the P-PET (Fig. 1a) and in phase with evaporation trends (Fig. 2b), which suggests a close link between decreased water recharge and processes like the PCP controlling the geochemical composition of the Onça Cave speleothems. How the cave atmospheric conditions determine the evaporative forcing and thus drive proxy variability can also be observed in the results from the calcite deposition experiment carried out at Onça Cave (Supplementary Fig. S8). In general, Mg/Ca ratios in farmed calcite vary out of phase with cave relative humidity, i.e., high values occur during periods of reduced relative humidity in the cave (Supplementary Fig. S8). These results confirm the PCP as the driver of trace element variability.

The combined effect of reduced rainfall and exacerbated evaporation lead to a decrease in the meteoric water recharge that may favor the PCP in response to decreasing water levels along the fractures and conduits of the cave vadose zone, which is reflected as higher Mg/Ca in the speleothems29,48–50. The role of the pluvial recharge in driving trace element variability in the Onça speleothems becomes apparent when comparing the detrended time series of Mg/Ca from Onça2 with local precipitation (Supplementary Fig. S10). The two variables are negatively correlated with an antiphased decadal variability, which reinforces the concept of hydrologic recharge leading to PCP as the main driving force of trace element variability (Supplementary Fig. S10). Nevertheless, a further increase in the trace element values can be attributed to the increase in the potential evapotranspiration that significantly reduced the meteoric water recharge in the last 40–50 years (Fig. 3). As we demonstrate in Fig. 3, variations in the trace elements from the Onça stalagmite are consistent with changes in evaporative demand, which ultimately, are driven by rising temperature (Fig. 1). To assess the correspondence between the observed trend in isotope and trace element ratios with the observed contemporaneous change in temperature, precipitation, and the evaporative demand, in Fig. 3 we regressed the z-scored time series of δ18O, δ13C, and trace element ratios against z-scored time series of these suspected climatological drivers using an ordinary least square regression. The significance of the regression coefficients calculated for each variable are confirmed by the F-test, presenting P values < 0.001. For this test we relied on the same change point detection algorithm by Killick31 to first determine the onset of the common trend (Supplementary Fig. S2), which occurs around 1942 (see “Methods”). We then compared this trend with the long-term temperature record from CRUTS 4.06 (Fig. 1c). To extend the PET to 1942, we used the Thornthwaite (1948) equation but applied a monthly coefficient correction factor that allows for an approximation to the PET based on the Penman–Monteith equation51.

The scaling factor calculated using ordinary least square regression highlights that the δ18O, δ13C, and trace element ratios share similar trends with the evaporative demand, (P-PET or evaporation), at multidecadal timescales, which ultimately is temperature-dependent. Based on the comparison of the linear trend, we can state that roughly 70% of the increase in these geochemical ratios can be attributed to the increase in the evaporative demand alone. This analysis is based on z-score data and the scaling factors thus represent units of standard deviation. On the other hand, temperature may drive significant variations in Mg/Ca in speleothem through changes in the partitioning of Mg uptake into calcite, leading to high Mg/Ca values in warm climates52–54. However, variations in Mg/Ca in the dripping solution, mainly associated with PCP and incongruent dissolution of dolomite, override the effects of temperature by more than an order of magnitude54.

Simulations from the Coupled Model Intercomparison Project Phase 6 (CMIP6)55 show that global temperature trends since the mid-1970s cannot be reproduced in climate models without considering anthropogenic forcing. Similarly, the timing of the observed increase in aridity recorded in the speleothem from central–eastern Brazil, starting in the 1970s, suggests an anthropogenic contribution (Fig. 4). To investigate the main drivers of changes in the regional hydrologic balance and to assess the anthropogenic contribution, we perform a Detection and Attribution (D&A) analysis33–35. The method is based on the estimation of scaling factors (regression coefficients) obtained from the regression between the geochemical proxies Mg/Ca and δ18O values of the speleothem and PET and P-PET simulated by an ensemble of CMIP6 models from the Detection and Attribution Model Intercomparison Project (DAMIP) (see “Methods”). The test was performed using data from natural and anthropogenically forced historical simulations (1850–2014 CE) covering a region of 4°×4° (13°–17° S; 42.5°–46.5° W) centered over the cave site; the same region used when calculating the hydrologic balance using instrumental records. For this purpose, we selected three different experiments: (i) a historical experiment (GHG+NAT) that includes natural forcing represented by time-evolving externally imposed volcanic aerosol and solar irradiance variability, as well as anthropogenic greenhouse gas and aerosol forcings; (ii) a natural historical experiment (NAT)—this simulation is forced only with solar irradiance and stratospheric aerosol forcing from volcanic eruptions, and (iii) a historical anthropogenic forcing experiment (GHG)—forced with the well-mixed greenhouse-gas changes only34. Our analysis is based on a multi-model ensemble, composed of 8 different CMIP6 earth system models (Supplementary Table S3). The times series were smoothed using a 5-year running mean in order to focus on the multiannual variability.

The relative contribution of the individual forcing components to the increased aridity observed in the Mg/Ca and δ18O speleothem records can be assessed in the whisker plots and in the time series for PET and P-PET calculated based on temperature and precipitation data extracted from the ensemble of the different model experiments (Fig. 4). The scaling factors, represented by the beta coefficient, calculated based on the z-score data indicate the extent to which the trends in Mg/Ca and δ18O values are explained by the individual forcing. They show the magnitude by which the dependent variable (Mg/Ca and δ18O) changes in response to a one-unit increase in the predictor variable (forcing). The distribution of the beta coefficients across the eight models for both the GHG+NAT and GHG simulations are significantly different from zero, indicating a detectable influence of anthropogenic and combined natural-anthropogenic forcing on Mg/Ca and δ18O in the simulations. On the other hand, the scaling factors calculated using natural forcings simulations only (NAT) are not significantly different from zero, regardless of the variable considered. These results show that Mg/Ca and δ18O at our cave site are responding to anthropogenic climate change, and that natural forcing alone cannot explain the observed trends in the speleothem records. The spread of the scaling factors obtained in the comparison between δ18O and P-PET derives from the differences in the simulated precipitation across the models. In general, simulating precipitation at such a regional scale is much more challenging for climate models than simulating temperature. Problems related to boundary conditions and misrepresentation of physical processes associated with precipitation-forming mechanisms result in significant inconsistencies between simulated and real precipitation, particularly in convective environments and on local scales56,57. However, the trend expressed by Mg/Ca and δ18O and the increased aridity starting around the 1970s is consistent only with the forcing in the anthropogenically forced experiments (GHG+NAT and GHG), which confirms the notion of a major influence of the anthropogenically induced temperature increase on the potential evapotranspiration as recorded in our speleothems (Fig. 4).

The role of temperature in driving hydroclimate variability in this region is evident in the 20th century, but the lack of longer and precisely dated and high-resolution temperature records covering eastern South America complicates a comparison prior to the 20th century. To circumvent this issue, we compared the regional hydroclimate variability recorded in the speleothems from central–eastern South America with the gridded surface temperature reconstruction from the PHYDA product32.

The most prominent feature recorded in the Onça2 stalagmite is the anomalous increase in the Mg/Ca values pointing to a prominent drying trend after 1970, which is consistent with the other geochemical proxies such as the δ18O and δ13C records (Fig. 5). Considering just the Onça2 stalagmite, it is readily apparent that the Mg/Ca, δ18O, and δ13C values over the last decades have no counterpart since 1729. Moreover, the combined δ18O and δ13C records from the Onça2 and Onça4 stalagmites show that the drying trend observed from the 1970´s onward is unprecedented in the past 720 years (Fig. 5). Notably, the observed drying trend is clearly detached from the decadal variability prior to the 20th century. From 1300 to 1850 CE the isotopic excursions are characterized by abrupt short-lived departures from the mean state, lasting ~10 years. Those events have no counterpart with the climate variability of last two-centuries and may be associated with abrupt climate changes related to variations in external forcing, such as large volcanic eruptions or solar irradiance changes. Conversely, in the 20th century, the isotopic variability is superimposed on a long-term isotopic enrichment. This notion is supported by results from a wavelet analysis, which indicates that multidecadal variability with a periodicity of ~60 years dominated the record from 1300 to 1800 CE (Supplementary Fig. S11). When combined, the two speleothem records from Onça Cave reinforce the idea that temperature is playing a major role in the hydrological balance of central–eastern Brazil during the last decades by modulating the potential evapotranspiration, thereby further exacerbating the drought conditions in the region.

In conclusion, the speleothem record from central–eastern Brazil points to an increased risk of future drought in the savannah biomes, primarily forced by the evaporation trend driven by anthropogenic global warming. The unprecedented upward trend seen in the geochemical proxies results from a combination of factors, including surface warming, increasing evaporative potential, and reduced precipitation, reflecting environmental conditions not seen in at least the past 720 years. The risk of an exacerbated future hydrologic drought is likely, given the role played by GHG forcing as shown in our study. This risk may be further exacerbated by the projected enhanced rainfall seasonality (e.g., Liang58), which will enhance the water deficit during the long dry season. Biomes such as the tropical savannah, where the dry season usually lasts five months or more, and which are home to vast areas of crops, native forests and cities with large populations, will therefore be at severe risk of long-term future drought. Our results also point toward temperature as an important factor to consider in the hydrological balance over the tropics, due to its ability to exacerbate droughts13,59,60.

Methods
Speleothem chronology
The Onça2 stalagmite has a robust chronologic control, based on 19 U-Th ages, covering the time period from 1760±4 to 2015±3 years CE, with an age uncertainty varying between 1 and 5 years (Supplementary Figs. S6 and S7). The Onça2 stalagmite also shows an annual layer structure that consists of sub-annual pairs of light and dark laminae, that enabled us to apply layer counting as an additional chronologic control. The accuracy of the chronology based on layer counting was assessed in comparison with the U-Th ages by applying the StalAge algorithm, which uses a Monte-Carlo simulation to construct an age model based on the best adjustment of linear trends through sub-sets of three adjacent data points61. In general, the chronology based on layer counting is in good agreement with the U-Th dates along the speleothem growth axis. Overall, the speleothem includes 255 layers, consistent with the time interval obtained from the U-Th ages (U-Th sample On2-0.1=2015±3 years and On2-17.8=1760±4 years CE), defining a temporal range of 256±7 years, Supplementary Fig. S7 and Supplementary Table S2). The Onça4 stalagmite record is composed of 22 U-Th ages, spanning from 1298±2 to 1852±2 years CE (Supplementary Fig. S7). The chronology of Onça4 was constructed using the StalAge algorithm optimal adjustment of linear trends between ages using a Monte-Carlo simulation61. The results of the U-Th ages are presented in Supplementary Table S2. The dating work was performed at the Institute of Global Environmental Change, Xi’an Jiaotong University, Xi’an, China.

Streamflow data
We selected a total of 13 fluviometric stations from ANA (Supplementary Table S1, see also data availability) covering the main tributaries of the São Francisco River in a 4°×4° area (13°–17° S; 43°–47° W) centered at the cave site. The limited number of stations before 1970 prevents us from assessing the spatial coherence of the overall fluviometric record. Therefore, we limit our analysis to the time period between 1970 and 2018 (Supplementary Table S1). The fluviometric stations were selected according to the following criteria: (i) at least 25 years of data, with less than 5% of missing data; (ii) with data starting before 1990; and (iii) with data ending after 2010. To calculate the mean streamflow, we first standardized each annual time series (z-score) between 1976 and 2016 before creating a regional average based on the median during the hydrologic year (from September to August). Since we are interested in analyzing the effect of natural forcings, i.e., precipitation and evaporation, we removed the record from the São Francisco River, which is strongly affected by urban occupation and water reservoirs.

Iscam
To construct a single composite for the isotope and trace element records, we applied the iscam algorithm from Fohlmeister62. This age-depth modeling software uses a Monte-Carlo approach on absolute age determinations to find the best correlation between different climate proxies from adjacent archives. The iscam method was applied to the composites of δ18O and trace elements from the Onça2 and Onça4 stalagmites (Supplementary Fig. S12). Age uncertainties at 68%, 95%, and 99% significance levels are obtained from an evaluation of a set of 2000 first-order autoregressive processes (AR1) for each record, which have the same statistical characteristics as the individual records. This method allows significantly reducing the age uncertainty within the overlapping periods and it can be tested if the signal of interest is indeed similar in all the records62.

Carbonate oxygen and carbon isotopes
The oxygen and carbon stable isotope analysis of speleothem calcite was performed at the Stable Isotope Laboratory at the University of São Paulo. Oxygen and carbon isotope ratios are expressed in δ notation, the per mil deviation from the VPDB standard according to the following equations:

δ18O = ½((18O/16O)sample / (18O/16O)VPDB) − 1) × 1000 (1)

δ13C = ½((13C/12C)sample / (13C/12C)VPDB) − 1) × 1000 (2)

The collected calcite powder was analyzed with an online, automated carbonate preparation system linked to a Finnigan Delta Plus Advantage mass spectrometer. Samples were collected along the growth axis of the speleothem with a sampling resolution varying from 0.125 to 0.4 mm. This sampling approach provided a near-annual resolution.

Water oxygen and hydrogen isotopes
The analyses of δ18O and δD from cave dripping water and rainwater were performed at the Centro de Pesquisas de Águas Subterrâneas at the University of São Paulo (IGc‐USP) using a Laser absorption spectrophotometer of the brand PICARRO L2130i. Data were processed by Laboratory Information Management System (LIMS) for Lasers software63,64. Values are reported with an analytical precision of 0.09‰ for δ18O and 0.9‰ for δD relative to Vienna Standard Mean Ocean Water.

The monitoring of δ18O and δD at the Onça Cave dripping sites started in February 2018 and continued until November 2019. A total of five dripping sites for monthly sampling were installed. The water samples were manually collected and stored in 8 ml high-density polypropylene bottles. To avoid isotope equilibration between water and atmosphere the flasks were completely filled and stored in a refrigerator at 5°C.

A local rainfall isotope monitoring study was performed in the region from November of 2011 to October of 2017. The isotope monitoring station, here referred as “Januária-INMET Station” (15°26'53“S; 44°21'59“W) is located about 40 km south of the cave site, at the INMET meteorological station (OMM: 83386) in Januária City. Local rainfall was collected weekly using a modified version of the rainfall water collector developed by Gröning65. The device was improved here by including one more layer of heat dissipation. These samples were collected in 8–30 mL high-density polypropylene bottles. Monthly rainfall δ18O values for local precipitation (Supplementary Fig. S9) were calculated based on the weighted mean of total precipitation between samplings.

Detection and attribution analysis
Detection and Attribution (D&A) analysis was performed using z-scored time series according to the following equation:

z-score = (xi − μ) / σ (3)

where xi is data at the time i, μ is the mean value, and σ is the standard deviation of the time series. For this analysis we used a 5-year running mean time series which mitigates the influence of shorter-term variability and highlights the underlying multidecadal trends. The ensemble of each experiment was calculated using the median of the eight models. To overcome the problems related to inter-dependency among observations (non-independent and non-identically distributed data) and retain the multidecadal signal and the long-term trend we performed linear regressions using 20 Monte-Carlo simulations based on the minimax algorithm from Xu and Xiaihua66 (Supplementary Figs. S13 and S14). In this approach, no white noise filtering was applied. For the regressions between speleothem geochemical proxies and CMIP6 simulations we used the period between 1900 and 2015. The resulting regression coefficients (β) provide a scaling factor that is used to compare the trends between time series. A standard t-test was applied to test whether the regression coefficients were significantly different from one another at the 95% level (Supplementary Table S4)67. This analysis shows that the scaling factor for the NAT scenario is significantly different from the GHG and the GHG+NAT scenario, and that we cannot reject the hypothesis that the scaling factor for NAT equals zero, while the scaling factors for GHG and GHG+NAT are distinct from zero. The proxies, i.e. δ18O or Mg/Ca were used as the response variable (x) and the different forcing simulations (e.g., temperature in the GHG simulation) were used as the explanatory variable (y). To obtain values of PET from CMIP6 runs we used the annual mean temperature output. Absolute values of PET were then calculated using the Thornthwaite (1948) equation and the monthly coefficient from Aschonities51 that allows an approximation to the PET based on the Penman–Monteith equation.

Change point detection algorithm
To identify the exact point where the trend along a time series changes, we performed a regime-shift test using a change point detection algorithm31,68. The tests was applied to the streamflow and speleothem proxy time series (Supplementary Fig. S2). In the streamflow time series, we assumed two change points, based on a visual inspection, given the two clear peaks around 1980 and 2004 (Supplementary Fig. S2). The algorithm identified two change points around the times 1979 and 2004, separating three periods with different trends; the first trend representing a flat and stable streamflow while the second and the third trend show a consistent decrease toward lower values. After determining the desired number of change points, the algorithm proceeds iteratively by randomly selecting points in the time series. The final estimate is obtained by choosing points where the sum of the square errors is minimized. In other words, for each desired number of change points (n), the algorithm iteratively defines n+1 trends in each iteration and compares the sum of square differences between the signal values and the least-squares fit (error sum of squares). The combination of n points that minimizes the error sum of squares represents the best estimate for change point detection31,68. To detect changes in a trend, the F-test is employed on the linear coefficients obtained pairwise along the time series.

Calcite monitoring
The calcite monitoring was performed using concave glass dishes of 11 cm diameter, also known as watch glass. Each watch glass was placed under the monitored dripping site. The experiment was replaced during every cave trip from July 2018 to October 2019, resulting in a nearly monthly sampling resolution. To monitor the calcite deposition, the watch glasses were weighed prior to deployment and after collection, yielding a mass of calcite of about 200 mg month−1. A common sampling surface with a radius equal to 5 cm around the center of deposition was defined between the glasses. The calcite was collected through manual scraping using a scalpel.

The calcite sample preparation was performed by weighing 10 mg of the homogeneous sample and adding 5 mL of 2% (v/v) HNO3 solution, which was stirred and filtered with 0.22-μm PVDF filter. Each sample was diluted 100 times for Ca quantification, so 100 μL of each solution was taken for a total volume of 10 mL (in HNO3 2% v/v/ solution). Quench water samples were collected at the sampling site, filtered, and acidified to 2% v/v HNO3. The elemental quantification of Ca (315.8 nm) and Mg (279.5 and 280.2 nm) were analyzed via ICP-OES (iCAP 6000 Series, Thermo Scientific, Bremen, Germany), in radial view, in partnership with the Spectrometry, Sample Preparation and Mechanization Group (GEPAM) at the Institute of Chemistry (University of Campinas, Brazil, SP) coordinated by Prof. Dr. Marco A. Z. Arruda.

Trace elements
Mg/Ca, Ba/Ca, and Sr/Ca ratios in the Onça stalagmite were obtained using Cavity Ring-Down Spectroscopy (CRDS) at the Laboratorio de Estudios Isotópicos at the Centro de Geociencias of the Universidad Nacional Autónoma de México (UNAM) following the methodology29 using a Resonetics M-50 laser ablation workstation with a large-format two-volume ablation chamber previously described in ref. 69 and ref. 70. We projected a 40 μm ArF excimer laser (193 nm, 20 ns FWHM) onto the sample surface with a fluence of 3 J/cm2 and a 5 Hz frequency. The laser was rastered along the growth axis at a speed of 0.5 mm/min within a He atmosphere, allowing to obtain a nearly continuous record of elemental variability in the stalagmite with sub-annual resolution. The resulting aerosol was admixed with Ar and transported to an iCap-Q ICP-MS where signal intensity ratios were measured. The analysis of the samples was bracketed by the analysis of SRM NIST 612, which allowed us to calculate elemental ratios (mmol/mol) using the reported composition by Jochum et al.71 or the NIST glasses.

Analyses of Mg/Ca ratios from calcite deposited on glass substrates were performed in the Department of Chemistry of the Universidade de Campinas using ICP-OES. Near 10 mg of calcite powder was sampled from glass substrates and diluted in 15 mL of HNO3 2%. The relative standard deviation between the measurement levels yields an analytical precision of ~1.1% for Ca and 0.74% for Mg.

Data availability
The speleothem U-Th ages generated in this study are provided in the Supplementary information. The geochemical records, including δ18O, δ13C and trace elements records from Onça2 and Onça4 stalagmites and composite from Onça2 and Onça4 stalagmites are provided in Supplementary Data 1.xls. The instrumental records for regional precipitation, evaporation, potential evapotranspiration and streamflow are presented in the Supplementary Data 2.xls file. The data of potential evaporation calculated to 13–17°S; 42.5–46.5°W were obtained from Brazilian Daily Weather Gridded from Xavier30 and are available at https://sites.google.com/site/alexandrecandidoxavierufes/brazilian-daily-weather-gridded-data?authuser=0. The streamflow data are also available at http://www.snirh.gov.br/hidroweb/. The cave monitoring data generated in this study, including isotope data from cave drippings, modern calcite deposition rate and Mg/Ca ratio from farmed calcite and cave temperature and relative humidity are presented in the Supplementary Data 3 (tabs: Cave dripping isotope data; Monitoring Mg_Ca and dep rate; Cave temp. and RH). Local rainfall isotope data from Januária Station generated in this study, including the records from GNIP Station are presented in Supplementary Data 4. The rainfall isotope data from GNIP station: Belo Horizonte CDTN and Brasília Airport are also available at https://www.iaea.org/services/networks/gnip. The surface temperature reconstruction from PHYDA is presented in the in the Supplementary Data 5. The Potential evapotranspiration (PET) and hydrologic balance (P-PET) calculated from CMIP6 experiments GHG+NAT, GHG and NAT are presented in the Supplementary Data 6 (tabs: “CMIP6_GHG+NAT”, CMIP6_GHG and “CMIP6_NAT”, respectively). Source data are provided with this paper.

Code availability
The “minmax algorithm” was performed using MATLAB (R2021a). The MATLAB function used in the minmax algorithm is provided in the Source Data file.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Strikis et al. present a new mulfiproxy record of rainfall and evaporafion from stable isotope (d18O,\nd13C) and trace element (Mg/Ca, Sr/Ca, Ba/Ca) measurements from two speleothems located in central-\neastern Brazil. Overall, the paper is quite well-wriften and structured, and the conclusions are well\nthought out. The record itself is also of high quality. The problem, however, is that many assumpfions are\nmade re: the controls on the speleothem proxies. Specifically, aftempts by the authors to quanfify P-PET\nusing speleothem proxies that may have different temporal resolufions to observafions and model\nsimulafions. In addifion, efforts to calibrate (or assess ‘linear congruence’) the speleothem proxies using\nobservafions (via linear regression) can pose problems when factoring in the inherent seasonal biases\nassociated with each proxy. For example, the d18O may be more sensifive to summer monsoon rainfall\nwhile the interannual trace element variability may be explained by dry season hydroclimate. These\nissues should be addressed in the manuscript.\nComments:\n(1) The calibrafion of the geochemical proxies (which have not been shown to be recording\nannual/seasonal variability) against climate model output using linear regression. Aside from the fact\nthat the different geochemical proxies have their own seasonal biases, the authors make liftle menfion\nof the fact that the speleothem proxies may be (likely) recording mulfi-year variability given karst mixing.\nHence, conducfing a linear regression of this type from data of (likely) varying resolufion is problemafic\nand warrants further discussion in the manuscript.\n(2) The outdated equafion used to calculate PET (using both observafions and model output) without\nconsiderafion of its uncertainfies given that it relies heavily on temperature alone. It is well known that\nevaporafive flux is also influenced by factors other than temperature, such as wind speed. Why did the\nauthors only use this rather outdated equafion to calculate PET and not test other more updated\nequafions? By doing so may have an influence on their findings. There is a rich body of work from the\nhydrology community looking at global and regional changes in PET from, for example, satellite data and\nusing more updated equafions. A deeper dive into this literature /research may benefit the study.\nFor a clarificafion of these points and others, please refer to the annotated pdf-version of the\nmanuscript.", "reviewer_2": "This paper presents an analysis of speleothem records in central Brazil that the authors argue relate to\nchanges in the watercycle in the region. The authors conduct a careful analysis of the role of\nprecipitafion and evaporafion in the data and reconstruct a changing ability to meet evaporafive demand\nfrom the data. A detecfion and aftribufion analysis then concludes that the driver of this change is\nanthropogenic. Overall it looks like there are some interesfing findings here, but the approaches could be\nexplained more clearly and i have some concerns about what exactly was done (which isnt everywhere\nclear either).\nthus i have the following comments\n1) i like the careful analysis of the role of the various contributors to hydrological changes and the use of\nmulfiple proxies. i dont quite understand why the individual proxies (13O and 18O) have such strong\nvariability while the combined anomaly which is interpreted as temperature anomaly (figure 5) has very\nrobust and consistent variability. can this be more clearly explained? and is the detecfion aftribufion\nanalysis conducted on the components - it seems like this from figure 4 - which is interpreted as P-PET.\nwhat is the uncertainty in the calibrafion of one to the other - does this need to be considered? I doubt\nthat the variability in the individual components is consistent with climate model variability.\n2) the authors find a change in hydrological variables that is only reproduced if all forcing simulafions are\nused to compare to the data. This is quite convincing as done, but the details are slim in the ms on the\\ Aftribufion approach - has a residual consistency test been conducted? (figure 5 suggests very strange\nvariability that is not simply red in the individual DeltaO 13 and 18 records this should not really agree\nwith model simulated variability?) it would be good to see a comparison of variability prior to confirming\nconclusions. Also, what analysis method was used, was exisfing code used following one of the recent\npapers (Ribes e.g.) or was a simple regression conducted? more should be said here and there definitely\nneeds to be a comparison of simulated control and residual variability of the approach to ensure that the\nmethod is self consistent. also the scaling between proxy data (black in fig 4) and model data and the\nrafionale of that comparison should be explained.\n3) the spafial scale at which the detecfion and aftribufion analysis is conducted is not clear to me. figure\n1 shows the locafion of the cave and a regional stafion measurement of temperature yet the variability is\nmuch larger than that of the wider CRUTEM region suggesfing that there are strong regional effects. the\nregional water cycle and precipitafion regionally can change for a lot of reasons, not all of them climafic.\ni am not fully convinced that the change seen in the stafion data and proxies is enfirely climafic based on\nwhat i can see in the paper. Also, the history prior to the 1980s shows quite a different evolufion of\nstreamflow vs rainfall and P-PET which suggests that other influences may be at work (figure 1) - could\nthat be discussed more clearly? generally, it is tricky to compare local precipitafion related records to\nlarger scale climate - can you rule out nonclimafic influences like land use change?\ndetail comments:\nwhat is the ANA in l 129? ( i think its explained in the supplement)\nl 196 reading later in the ms i think the comparison of trends involves all variables in z scores. however\nthe figure shows some in other units please clarify in text that the slopes are in z scores.\nl 301 i think you have said prefty much the same about use of a D+A method a bit further up in the ms,\nthis is a befter place for it. maybe remove above. It needs to be clear what spafial representafion of P-\nPET you use in your comparison.\nThe supplement gives the models but not the ensemble sizes. it would be good to menfion these.\nfigure 4 capfion: how is the detecfion and aftribufion analysis conducted - Hist vs Hist-nAT but then what\nis the role of GHG? and what is nat-GHG in the capfion? are you using 3 fingerprints but then there is\noverlap between GHG and Hist? overall the detecfion and aftribufion approach could use a lot more\nclarity including size of region of model data used and scaling between P-PET and Onca2. (unless that is\ndrawing on figure 3 which might be the case? could this be menfioned?)\nfig. 5: why does the combinafion of the speleothems in the composite remoe the strong variability in the\nindividual ones - wouldnt the temperature anomaly be then befter suited for aftribufion (although less\nexcifing) due to much more homogeneous variability? unless the climate model P-PET variability is\nsimilar to that in figure 5...which i doubt\nl 450: i am not convinced by analyzing P-PET from 1979 on and all else from 1900 this looks like\ncherrypicking to me.\nsome liftle typoes e.g. detected in l 179, a few more too so maybe careful proofreading by a nafive\nspeaker may be useful"}, "author_response": "We are grateful to the reviewers for their thoughtful and encouraging comments on our\nmanuscript NCOMMS-23-00133. We addressed the reviewer comments below in bold blue text.\n[...]"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I am safisfied with the revisions and commend the authors for their efforts in addressing my comments\nand suggesfions. Thus, I can now recommend the paper for publicafion.", "reviewer_2": "Thank you for considering my comments carefully and i find the revised ms much clearer (and thanks for\nthe z-scores in some literature (non hydro) this is call standardized but as long as its explained it is a good\nchoice.\nThere are some small issues remaining that i would like to see addressed:\n1) the significance tesfing is sfill not quite clear to me. You relate to Allen + Teft's paper that is based on\nan F-test, but it will only work on pre-whitened data (i.e. method with opfimizing - it is not menfioned\nthat this was done?). So the method text sfill doesnt quite give enough informafion, and sounds more\nlike a generic F test is conducted which would not be ok as the method will need to detect the role of\nforcing against decadal and mulfidecadal variability. Any method used needs to take that into account.\nmost approaches just use samples from either model control simulafions, or the within-ensemble\nvariability in the forced runs providing the fingerprints (e.g. all GHG+nat runs individually minus the\nensemble mean, which gives a sample of variability, corrected by a factor for the variance loss from\nsubtracfing the ensemble mean which is sfill noisy (I think its root(n/(n-1)). please clarify the uncertainty\nesfimate in the methods secfion - there is not much there and the rebuftal reference gives the wrong\nline numbers.Alternafively, maybe you used available software which would be helpful (and fair) to say,\nplease clarify this. l 475 menfions just a standard F test. the analyses you made to support robust\nresiduals is quite nice and could go into a supplement to support reliable results.\n2) thanks for clarifying the forcing simulafions. i am not quite clear if the anthropogenically forced\nsimulafions include aerosols and other anthropogenic forcings - that is usually the convenfion. can you\nclarify in text and fix? normally, HIST includes all forcings including other anthropogenic than ghg, such\nas aerosols, ozone, land use... while nat include solar and volcanic. this could use a quick check. (e.g, l\n325)\n3) this is a new quesfion sorry - i think its addressed but not clear from text: you illustrate nicely how all\nyour fimeseries show similar trends in z scores. However, this could happen for a lot of reasons, while in\nyour case they also covary strongly as obvious from the figure (but not quanfified at least not in capfion)\n- i recommend to highlight that too. that way you have many degrees of freedom supporfing your proxy.\nMinor comments:\n- Abstract, l 95 and on: would prefer not to call detecfion and aftribufion 'formal detecfion and\naftribufion' as i dont think thats a useful disfincfion and it may also get confused with approaches using\nsome of the available software and fairly sophisficated methods while here i think a mulfi regression only\nis used (which is fine for fimeseries parficularly and more transparent!). If you have used available\nsoftware and a prewhitening/opfimal fingerprint approach please say so.\nl 74 has a typo or strange format in the reference 20\ni suggest you refer to this paper in the introducfion which draws related conclusions but for different and\nlarger regions: Marvel, K., B.I. Cook, C.J.W. Bonfils, P.J. Durack, J.E. Smerdon, and A.P. Williams, 2019:\nTwenfieth-century hydroclimate changes consistent with human influence. Nature, 569, no. 7754, 59-65,\ndoi:10.1038/s41586-019-1149-8."}, "author_response": "We are grateful to the reviewers for their thoughtful and encouraging comments on our manuscript\nNCOMMS-23-00133. We addressed the reviewer comments below in bold blue text.\n[...]"}]

Decision letter