Global analysis reveals climatic controls on the oxygen isotope composition of cave drip water.

The oxygen isotope composition of speleothems is a widely used proxy for past climate change. Robust use of this proxy depends on understanding the relationship between precipitation and cave drip water δ18O. Here, we present the first global analysis, based on data from 163 drip sites, from 39 caves on five continents, showing that drip water δ18O is most similar to the amount-weighted precipitation δ18O where mean annual temperature (MAT) is < 10 °C. By contrast, for seasonal climates with MAT > 10 °C and < 16 °C, drip water δ18O records the recharge-weighted δ18O. This implies that the δ18O of speleothems (formed in near isotopic equilibrium) are most likely to directly reflect meteoric precipitation in cool climates only. In warmer and drier environments, speleothems will have a seasonal bias toward the precipitation δ18O of recharge periods and, in some cases, the extent of evaporative fractionation of stored karst water.

Introduction

The oxygen isotope composition is by far the most widely reported climate proxy in cave deposits, or speleothems (e.g., stalagmites, stalactites and flowstones[1]). Multiple processes determine the oxygen isotope composition of speleothems (δ18Ospeleo), with the potential climate signal reflecting the source water (meteoric precipitation) δ18O (δ18Oprecip) and its relationship to local and regional climate. This signal is transferred to the cave through the vadose zone, where it may be mixed with existing waters and fractionated by evaporation. Finally, at the target (the speleothem), the δ18Ospeleo signal can be further altered by non-equilibrium fractionation processes and temperature-dependent fractionation during calcite precipitation. See refs. [1-3] for in-depth reviews of these processes and climate signal transformation.

Within the speleothem research community, it is widely acknowledged that a cave monitoring approach is necessary to fully understand, and constrain quantitatively, the extent that the climate signal is preserved in δ18Ospeleo (e.g., during transfer from the source to the target). The measurement of drip water hydrology[4], drip water geochemistry[5], cave environment[6] and calcite growth and geochemistry[7], as well as surface climate parameters, allows empirical relationships between the source and the target to be determined. With monitoring data, regression models between climate and speleothem proxy data can be developed[8], proxy interpretations can be evaluated[9], input data for forward or proxy system models can be generated[10-13] and the extent that speleothem calcite precipitates in isotopic equilibrium with its associated drip water can be assessed[7,14,15].

Recently, a new global database of speleothem carbon and oxygen isotope proxy records was compiled[16,17]. This archive includes 455 δ18Ospeleo records, with over 324 covering intervals within the last 21 ka[16,17]. Some regions have δ18Ospeleo records that span glacial–interglacial intervals (e.g. monsoon regions[18-20]), whereas other regions have records that are more complex (e.g. water-limited regions where δ18Ospeleo exhibits high magnitude and frequency variability[21,22]). In water-limited environments, potential mechanisms by which δ18Ospeleo can be modified during transit from the source, include evaporative fractionation of water δ18O in the soil; a shallow vadose zone or cave; selective recharge, whereby rainfall events with high amount or intensity have a distinct isotopic composition, typically low δ18O; non-equilibrium deposition during speleothem formation[23-27]. A fundamental research question is: what are the regional climate parameters where δ18Ospeleo values most faithfully preserve the source signal (δ18Oprecip)? Identification of such climatic regions, and speleothem samples, will have the greatest utility; for example, for research methodologies, such as data assimilation[28], which utilise proxy–climate model inter-comparison.

Interpretation of δ18Ospeleo proxy records would benefit from the best possible understanding of the climatic conditions under which oxygen isotope composition of drip water (δ18Odripwater) is most directly related to δ18Oprecip. Here, we compile cave monitoring data with the objective of understanding the modern-day relationship between δ18Oprecip and δ18Odripwater. We compile data sets where there are both cave δ18Odripwater data (1-year or longer data sets) and δ18Oprecip data (of equal duration, amount-weighted and collected close to the cave and similar altitude). The latter enables the amount-weighted precipitation oxygen isotope composition (δ18Oamountwprecip) to be compared with δ18Odripwater. By using a karst hydrology model developed for European climates, monthly modelled recharge amount is used to obtain an annual recharge-weighted δ18O (δ18Orechargewprecip) at European sites. This permits the first analysis of δ18Odripwater, δ18Orechargewprecip, δ18Oamountwprecip and climate parameters. The analyses show that drip water δ18O is most similar to the amount-weighted precipitation δ18O, when mean annual temperature is < 10 °C. The implications for speleothem palaeoclimatology are that speleothems (if formed near isotopic equilibrium) are most likely to directly reflect meteoric precipitation δ18O only in cooler climates.

Results

Global water oxygen isotope distributions

We find a strong positive correlation between δ18Odripwater and δ18Oamountwprecip. δ18Orechargewprecip provides a similarly strong correlation, but in this case with a slope and intercept indistinguishable from 1 and 0, respectively. Supplementary Data 1 presents the database of δ18Odripwater and δ18Oamountwprecip compiled from the literature and unpublished data comprising 163 drip sites from 39 caves on five continents. The location of the caves in comparison with modern mean annual temperature (MAT) and the global database of δ18Ospeleo records[17] are shown in Fig. 1. Climate regimes represented in the compilation include temperate maritime and semiarid monsoon, Mediterranean, montane and tropical, therefore including a wide range of MAT and aridity, as expressed by the ratio of precipitation to potential evapotranspiration (P/PET).

Fig. 1

Global map of sample sites, karst regions, surface temperature and speleothem records. Location of the cave δ18Odripwater samples (large circles). Global karst aquifer regions[41] are shown as coloured areas, with those with mean annual temperature < 10 °C (blue); 10 °C < mean annual temperature < 16 °C (green) and mean annual temperature > 16 °C (red). Dots show the locations with speleothem (δ18Ospeleo) records in the SISAL (Speleothem Isotopes Synthesis and AnaLysis Working Group) database[16,17]. a Europe, b Chinese monsoon region and c SE Australia. More information on the sites is presented in Supplementary Data 1

Figure 2a, b presents the global relationship between δ18Odripwater and δ18Oamountwprecip. The correlation is positive and strong (Spearman’s rank rs = 0.90, p < 0.00001), indicating that at a global scale, δ18Odripwater closely relates to δ18Oamountwprecip. The regression demonstrates that, at this scale, δ18Odripwater is greater than δ18Oamountwprecip where the latter is more positive, typically sites where MAT > 16 °C. Conversely, δ18Odripwater is less than δ18Oamountwprecip where the latter is more negative, typically at sites where MAT < 16 °C. Regional relationships between δ18Odripwater and δ18Oamountwprecip for Europe, China and Australia are quantified in Supplementary Fig. 1. At a regional scale, the correlation is positive, very strong and highly significant for the European region and moderately strong for China.

Fig. 2

Global oxygen isotope relationships. a Global relationship between δ18Odripwater and δ18Oamountwprecip. The global data set regression line is shown in red: δ18Odripwater = 0.64 (±0.25) + 1.10 ( ± 0.04) δ18Oamountwprecip (‰). Sites are coloured according to their mean annual temperature and symbols show their region: Europe (squares), China (circles), Australia (diamonds), United States (triangles) and other (stars). b Frequency histogram for the global data set for the difference between δ18Oamountwprecip and δ18Odripwater (Δawp-dw). c Relationship between the δ18Orechargewprecip, δ18Oamountwprecip and δ18Odripwater for the European data set. The amount-weighted data are shown in open black symbols, and the regression line is shown in gray: δ18Odripwater = 1.19 (±0.59) + 1.20 (±0.08) δ18Oamountwprecip (‰). The recharge-weighted data are shown by coloured symbols (as for (a)) and the regression line is shown in black: δ18Odripwater = 0.06 (±0.50) + 1.01 (±0.06) δ18Orechargewprecip (‰). The arrows show the direct effect of recharge weighting. d Frequency histogram for the European data set for the difference between δ18Orechargewprecip and δ18Odripwater (Δrwp-dw) and δ18Oamountwprecip and δ18Odripwater (Δawp-dw) for the European data. Typical analytical uncertainties for individual δ18O analyses are ± 0.2‰[42]

For cave drip water monitoring sites in Europe, we utilise a karst hydrology model[29] to determine the monthly recharge amount (see the ‘Methods’ section), and these monthly recharge values (see Supplementary Table 1) were then used to weight the δ18Oprecip in that month. At the European scale, the relationship between the δ18Odripwater and δ18Oamountwprecip is a strong positive correlation (Spearman’s rank rs = 0.90, p < 0.00001), similar to that observed globally (Fig. 2c, d), although over a more restricted range of δ18O. With recharge weighting, the correlation between the δ18Odripwater and δ18Orechargewprecip remains positive and strong (Spearman’s rank rs = 0.89, p < 0.00001). The intercept and gradient are indistinguishable from 0 to 1, respectively, indicating that after recharge weighting, at the European sites, δ18Odripwater can be explained by δ18Orechargewprecip.

Climate controls on selective recharge and partial evaporation

We provide empirical evidence from the global δ18Odripwater data set that increasing temperature and decreasing rainfall both increase the absolute difference between δ18Odripwater and δ18Oamountwprecip. Figure 3 explores the global relationship between climate parameters and the difference between amount-weighted precipitation and drip water isotopic composition (Δawp-dw = δ18Oamountwprecip − δ18Odripwater). It can be observed that there is a narrowing in the range of Δawp-dw when MAT is relatively low (<10 °C), the total annual P is high (>1750 mm), the annual PET is low (<800 mm) or the total annual P/PET values are high (>1.5). Linear single and stepwise multiple regression analyses on the global data set showed that the strongest correlation (Spearman’s rank) of the absolute value of Δawp-dw was with the ratio of mean annual temperature (MAT) to the total annual P:

Fig. 3

The global relationship between Δawp-dw and climate parameters. a Mean annual temperature (MAT), b total annual precipitation (P), c total annual potential evapotranspiration (PET), and d mean annual potential evapotranspiration relative to mean annual precipitation (P/PET). Colours represent different regions: Australia (black), China (green), Europe (blue), United States (cyan) and all other regions (magenta). Black vertical lines show the 0.3‰ criterion for determining the significant difference between δ18Oamountwprecip and δ18Odripwater

To further explore the relationship between Δawp-dw and these climate parameters, we define a threshold for |Δawp-dw| of more than 0.3‰ as a criterion for determining the significant difference between δ18Oamountwprecip and δ18Odripwater. This is chosen taking into consideration potential uncertainties in δ18O determinations of water and speleothem calcite (analytical uncertainties of 0.06–0.2‰, depending on measurement technique). Considering the climate parameter MAT, 91% of all drip waters with a MAT < 10 °C (n = 34) have a |Δawp-dw| of <0.3‰. Considering the P, then for a P threshold of 1750 mm, 61% of all drip waters (n = 31) have a |Δawp-dw| of < 0.3‰. These empirical observations agree with theoretical understanding that in warmer, water-limited climates, δ18Odripwater may be affected by evaporative fractionation of the water in the soil or shallow karst[22,30], or by selective recharge, with an isotopic composition dominated by those rainfall events or seasons that generate recharge[26,27]. However, we note that a combination of post-infiltration evaporative fractionation and isotopically depleted recharge could lead to observations of |Δawp-dw| < 0.3‰ for some sites with warm and dry climates.

Discussion

Our recharge modelling demonstrates the importance of selective recharge, and suggests that for a MAT < 16 °C, δ18Odripwater is best interpreted as δ18Orechargewprecip. The 1:1 linearity of the relationship between δ18Orechargewprecip and δ18Odripwater for European sites confirms the importance of selective recharge for this climate range (seasonal climates with a MAT ranging from 7.1 to 16.1 °C, Supplementary Data 1). Selective recharge is minimised at MAT < 10 °C. At these temperatures, the opportunity for soil and shallow karst evaporation is decreased, and karst water stores are more likely to be maintained, allowing mixing of recharge waters that buffer the isotopic impact of any individual recharge event. At a MAT < 10 °C, speleothems that have been deposited close to equilibrium would have the potential of recording past variations of δ18Oamountwprecip, plus a temperature signal from the fractionation during calcite precipitation.

Latitudes poleward of ~35° and high-altitude sites, where MAT < 10 °C (Fig. 1), would be most likely to contain a δ18Ospeleo record of amount-weighted precipitation (northern Europe, high-altitude and northern regions of the Asian monsoon, northern North America and New Zealand). In contrast, δ18Ospeleo records in regions of higher MAT are more likely to have |Δawp-dw| > 0.3 ‰ and would be sensitive to moisture balance changes, due to limited mixing with stored water, selective recharge and/or increased chance of evaporative fractionation of δ18O in the vadose zone. δ18Odripwater, and the associated δ18Ospeleo, can be more positive than amount-weighted precipitation (evaporative fractionation dominates), or either greater or less than amount-weighted precipitation (selective recharge dominates). Regions where this compound signal is most likely are predominantly in latitudes equatorward of ~35° (most of Africa, India, southern Asia, southern Europe, North America and Australia; Fig. 1). Modelling of δ18Orechargewprecip suggests that for seasonal climates with a MAT between 10 and 16 °C (the higher value being the upper bound of the European data set), selective recharge dominates these processes. At this range of MAT (and, we anticipate, at higher MAT), δ18Ospeleo may be a proxy for δ18Orechargewprecip and provides records of paleo-recharge. In addition, when considering δ18Ospeleo, any relationship between δ18Odripwater and climate could be additionally overprinted by non-equilibrium deposition.

Our meta-analysis reveals that the oxygen isotope composition of drip water is primarily determined by the oxygen isotope composition of the recharge water δ18O. At a global scale, we show that the extent to which δ18Odripwater is representative of δ18Oamountwprecip is primarily related to the mean annual temperature and annual precipitation, which determines the extent to which δ18O is further altered by soil and karst processes. To confidently interpret the δ18Odripwater as a specific climate parameter, the relationship between recharge δ18O and climate needs to be understood for specific sites. For sites and regions, characterised by lower temperatures (MAT < 10 °C), where Δawp-dw is likely to be closest to zero, we show that the oxygen isotope composition of drip water is most directly related to the isotopic composition of local rainfall. These regions could produce δ18Ospeleo proxies (if the speleothems are deposited close to equilibrium), where δ18Ospeleo could be used to provide a signal of past δ18Oamountwprecip and cave air temperature (due to the temperature-dependent fractionation during calcite formation), useful for proxy–model assimilations. In these cooler climates, where water in karst stores and fractures is more likely to be well mixed, one would also expect greater agreement in δ18Odripwater between drip sites within a cave. In regions with higher temperatures (MAT > 16 °C), δ18Ospeleo is less likely to represent δ18Orechargewprecip, and instead can contain a compound signal that reflects selective recharge and evaporative fractionation. Such records are of palaeoclimatic value, and may yield a proxy for δ18Orechargewprecip, but are more likely to show greater heterogeneity between coeval records and therefore require a drip-specific interpretation.

Important Quaternary δ18Ospeleo records have been produced from around the world, and in the context of this analysis of modern conditions, we can make several conclusions. Firstly, many palaeoclimate studies interpret the relative changes in δ18Ospeleo over time, and in many cases, monitoring data are not available to guide the interpretation. The climatic controls made here can be used to help guide the interpretation of those records. This is particularly relevant over periods of significant climate change (e.g. glacial–interglacial transitions) and where the climate control on the difference between δ18Oamountwprecip and δ18Odripwater may change over time. A map of the cave sites at modelled last glacial maximum (LGM) surface temperatures is provided in Supplementary Fig. 2, and suggests that a change in the temperature control on the δ18Oamountwprecip–δ18Odripwater relationship is mostly observed in mid-latitudes, and most ubiquitously in the LGM in southern Europe. Secondly, in the Chinese monsoon region, the cooler northern sites are most likely to have δ18Odripwater similar to δ18Oamountwprecip, as reported previously[30]. However, given that monsoon rainfall requires a land–ocean temperature gradient, there is a trade-off between caves at cooler locations that have δ18Odripwater closest to δ18Oamountwprecip, and those in regions with the strongest monsoon signal. The latter are more likely to experience evaporative fractionation and selective recharge, and therefore less likely to be similar to δ18Oamountwprecip (but may reflect δ18Orechargewprecip). This trade-off would apply to all monsoon regions. At the modern monitoring sites in the Chinese region, where MAT > 10 °C and annual P < 2000 mm (Fig. 3), δ18Oamountwprecip does not correlate with MAT or the total annual P, but δ18Odripwater does positively correlate with both (Supplementary Fig. 3). This appears to be due to the combined overprinting of increasing selective recharge and evaporative fractionation over this range of MAT and offers new insights into the interpretation of δ18Ospeleo in the region. Thirdly, even in regions of exceptionally high rainfall, such as Mulu (Malaysian Borneo) and parts of India, δ18Odripwater can be higher than the δ18Orecharegewprecip[31], probably due to the continuous high temperatures, leading to the partial evaporation of vadose water. Analysis of speleothems at caves at higher elevations should help mitigate this effect. Finally, δ18Ospeleo records from regions with high aridity and temperatures should not be expected to preserve a record of δ18Oprecip. Our meta-analysis confirms the modern monitoring observations[25], which indicate that δ18Ospeleo in these regions would be an archive of alternating palaeo-aridity and palaeo-recharge and supports the interpretation of δ18Ospeleo as a palaeo-recharge and palaeo-aridity proxy for the last glacial maximum in arid southern Australia[22].

Methods

Data compilation

δ18Odripwater data were compiled from a literature search and unpublished data. To minimise uncertainties that could be introduced into our analysis, we have chosen to only include sites where both of the following two criteria were met. Firstly, δ18Oprecip was collected at or close to the sites (<20 km) and at a similar altitude, monthly integrated samples for at least 1 year, with an amount-weighted annual mean (δ18Oamountwprecip) value reported. Secondly, δ18Odripwater was collected over the hydrological year, for at least 1 year, with at least bimonthly sampling frequency. Monitoring results had to have at least 1 year of both δ18Odripwater and δ18Oprecip data, with overlapping time periods. We therefore have not included sites where δ18Oprecip is a derived parameter, e.g. from isotope-enabled GCM output or based on empirical relationship with distant Global Network of Isotopes in Precipitation (GNIP) stations. Average drip water age is unknown for all sites, and it is possible that for some locations, the δ18Odripwater integrates δ18Oprecip prior to the monitoring period.

For each site, the local MAT and the total annual P were taken from the publications, and PET was taken from the WorldClim Global Climate Database[32,33]. For one study[29], the total annual precipitation was not provided, and output from the gridded data set was used instead. The P/PET was calculated from the local P and gridded PET.

Climate comparison

δ18Odripwater and δ18Oamountwprecip data were compared with the following climate characteristics: mean annual temperature (MAT), total annual precipitation (P), potential evapotranspiration (PET) and the precipitation over PET ratio or  aridity index (P/PET). PET and the P/PET were taken from the global aridity and PET database[32,33], where PET is modelled at ~1-km resolution, using data from the WorldClim Global Climate Database using mean monthly extraterrestrial radiation, and mean monthly temperature and range (using the equation of ref. [34]). Sites are classified as humid where P/PET > 0.65; semi-arid and dry sub-humid at 0.2 ≤ P/PET ≤ 0.65; arid and hyper-arid at P/PET < 0.2. The difference between the δ18Oamountwprecip and δ18Odripwater was determined for each drip site (Δawp-dw).

As cross-checks on the gridded database, we compared for all caves local P and gridded P (Eq. 2) and local T and gridded T (Eq. 3), and for the Australian caves, we compared gridded PET with the mean PET (1960–1990 AD) calculated from the Australian Water Availability Project (AWAP) database[35,36]:

For the Australian sites, the gridded PET calculated by the two products agreed within 7% for all sites, except Golgotha Cave, where the AWAP PET was 30% higher than that calculated by WorldClim. The difference in PET at this site did not change the P/PET classification (using WorldClim: 1.06; using AWAP 0.82), and the WorldClim data are used for consistency.

Statistical analyses were undertaken using Microcal Origin. Correlations are Spearman’s rank-correlation coefficients (rs). Probability values (p), are conservatively determined using the lowest degrees of freedom (df), based on the number of cave sites (global: n = 39; Europe: n = 16; China: n = 10; Australia: n = 5), rather than the number of unique drip waters. Regression equation slope and intercept uncertainties are the standard error.

Karst hydrological model

To estimate recharge, we use a large-scale karst groundwater recharge model that was previously developed for European and Mediterranean climates[29,37,38]. The model simulates karstic groundwater recharge at a 0.25° × 0.25° resolution at a daily resolution for a 10-year period from 2002 to 2012, which we consider long enough to provide representative average values of monthly recharge. All relevant karstic and non-karstic processes, such as infiltration of rainfall and snowmelt, evapotranspiration, downward percolation from the upper soil layer to a lower soil/epikarst layer and vertical percolation from the epikarst layer towards the groundwater, are considered within its structure. The epikarst, which is a typical vadose-zone feature of karst systems, allows the dynamic separation of focused and diffuse groundwater recharge[38,39]. For the weighting of recharge, output from the epikarst is used: the epikarst in the model is a series of N = 15 linear storages with variable capacities (V [mm] and V [mm]) and with variable storage constants (K[d]), which are distributed by a Pareto function, with a shape parameter a [−]. Parameter estimation provided ranges of values for V V, K and a for the humid, mountain, Mediterranean and desert karst landscapes[29]. Here, we use the average recharge volumes (over all simulations obtained with the parameter sets within these confined ranges), and weight the δ18Oprecip in each month by the fraction of the total annual epikarst recharge that occurred in that month.