A dual-biomarker approach for quantification of changes in relative humidity from sedimentary lipid D/H ratios

Abstract

Past climatic change can be reconstructed from sedimentary archives by a number of proxies. However, few methods exist to directly estimate hydrological changes and even fewer result in quantitative data, impeding our understanding of the timing, magnitude and mechanisms of hydrological changes.

Here we present a novel approach based on ${\mathit{\delta }}^{\mathrm{2}}$H values of sedimentary lipid biomarkers in combination with plant physiological modeling to extract quantitative information on past changes in relative humidity. Our initial application to an annually laminated lacustrine sediment sequence from western Europe deposited during the Younger Dryas cold period revealed relative humidity changes of up to 15 % over sub-centennial timescales, leading to major ecosystem changes, in agreement with palynological data from the region. We show that by combining organic geochemical methods and mechanistic plant physiological models on well characterized lacustrine archives it is possible to extract quantitative ecohydrological parameters from sedimentary lipid biomarker ${\mathit{\delta }}^{\mathrm{2}}$H data.

Introduction

Predicting future changes in the water cycle using state-of-the-art climate models is still associated with large uncertainties (IPCC, 2014). This is because we lack a mechanistic understanding of some of the key processes that influence the water cycle, in particular on regional spatial scales. A better mechanistic understanding of drivers and feedbacks within the hydrological cycle can be achieved from reconstructing past hydrological changes from sedimentary archives. Stable isotope ratios of meteoric water, expressed as ${\mathit{\delta }}^{\mathrm{18}}$O and ${\mathit{\delta }}^{\mathrm{2}}$H ($\mathit{\delta }$D) values are an excellent tool in this respect, because their variability is associated with changes in temperature and source water (Bowen, 2008; Gat, 1996). The isotope ratios of precipitation can be recorded in ice core (Alley, 2000), terrestrial and marine paleoclimate archives through a variety of proxies, such as carbonates (Kanner et al., 2013; von Grafenstein et al., 1999), silicates (Tyler et al., 2008) and lipid biomarkers (Sachse et al., 2012).

The interpretation of the stable isotope ratios from inorganic and organic proxies, despite their potential, often allows only a qualitative assessment of past hydrological changes, while quantitative reconstructions of hydrological changes from isotope proxy data, such as precipitation amount or relative humidity (RH), have been difficult to achieve. This is problematic as quantifiable data are necessary for identifying the mechanistic drivers of past hydroclimate changes as well as their continental-scale feedbacks and thresholds, for vegetation changes, for example. Moreover, quantitative data are needed to test the performance of state-of-the-art climate models in simulating past and future changes in the hydrological cycle.

The interpretation of isotope proxies is typically not quantitative because multiple drivers can influence meteoric ${\mathit{\delta }}^{\mathrm{18}}$O and ${\mathit{\delta }}^{\mathrm{2}}$H values, hampering the assignment of single quantitative relationships between a hydrologic variable and ${\mathit{\delta }}^{\mathrm{2}}$H values recorded in a geological archive (Alley and Cuffey, 2001). The increased understanding of the interplay between environmental and plant physiological factors affecting lipid biomarker stable isotope ratios over the last decade (Feakins, 2013; Kahmen et al., 2013a, b; Sachse et al., 2009; Smith and Freeman, 2006) has resulted in significant potential for quantitative paleohydrological approaches, exemplified by a reconstruction of seasonality in precipitation and bog surface wetness in a Norwegian peatland (Nichols et al., 2009). Here we take this a step further, combining lipid biomarker hydrogen isotope measurements and plant physiological modeling to constrain the influence of multiple drivers on ${\mathit{\delta }}^{\mathrm{2}}$H values recorded in organic material. Thus, we allow the extraction of quantitative information about changes in relative humidity from sedimentary archives.

Over the past decade, ${\mathit{\delta }}^{\mathrm{2}}$H values of lipid biomarkers from photosynthetic organisms have been increasingly used as proxies for reconstructing past changes in the continental hydrological cycle (Feakins, 2013; Rach et al., 2014a; Sachse et al., 2012; Schefuss et al., 2011; Seki et al., 2011). In particular $n$-alkanes are ubiquitous in marine and lacustrine sediments and can be preserved over geological timescales (Peters et al., 2007). $n$-Alkanes can be traced back to aquatic or terrestrial sources, where short-chain homologues ($n$C${}_{\mathrm{17}}$–$n$C${}_{\mathrm{21}})$ are primarily synthesized by algae and aquatic plants (Aichner et al., 2010; Ficken et al., 2000), mid-chain $n$-alkanes (e.g. $n$C${}_{\mathrm{23}}$–$n$C${}_{\mathrm{25}})$ by submerged aquatic macrophytes or mosses (Aichner et al., 2010; Ficken et al., 2000; Gao et al., 2011) and long-chain $n$-alkanes (> $n$C${}_{\mathrm{25}})$ predominantly by higher terrestrial plants as a protective leaf wax layer on the leaf surface (Bush and McInerney, 2013; Eglinton and Hamilton, 1967).

Algae and submerged aquatic plants directly use lake (or ocean) water as their hydrogen source for lipid synthesis. ${\mathit{\delta }}^{\mathrm{2}}$H values from $n$-alkanes from aquatic organisms (${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{aq}})$ are thus related to the ${\mathit{\delta }}^{\mathrm{2}}$H value of the water these organisms live in (Aichner et al., 2010; Sachse et al., 2004), offset by a biosynthetic fractionation (${\mathit{\epsilon }}_{\mathrm{bio}})$ between water and $n$-alkanes (Sachse et al., 2012) (Eq. 1). Laboratory culture studies (Zhang and Sachs, 2007) as well as field studies (Aichner et al., 2010; Sachse et al., 2004) have resulted in strong linear and nearly 1 : 1 relationships between source water and ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{aq}}$ (Sachse et al., 2012), but have shown that species-specific differences in ${\mathit{\epsilon }}_{\mathrm{bio}}$ do exist (Zhang and Sachs, 2007). $${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{aq}}={\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{precip}}+{\mathit{\epsilon }}_{\mathrm{bio}\left(\mathrm{aq}\right)}$$ Terrestrial plant leaf wax $n$-alkane ${\mathit{\delta }}^{\mathrm{2}}$H values (${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}})$ have also been found to be linearly correlated to the organism's source water ${\mathit{\delta }}^{\mathrm{2}}$H values, yet not in a 1 : 1 relationship (Sachse et al., 2012), indicating additional influences on ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ values. Recent greenhouse experiments and field studies have revealed that in particular the evaporative ${}^{\mathrm{2}}$H enrichment of leaf water shapes ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ values (Kahmen et al., 2013a, b). Soil water evaporation in the upper soil layers has been shown to be less significant for ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$, as plants usually access the deeper, isotopically unenriched, soil layers (Dawson, 1993). As such, ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ is affected mainly by the ${\mathit{\delta }}^{\mathrm{2}}$H value of plant source water (i.e. precipitation), the biosynthetic fractionation and leaf water deuterium enrichment (${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{e}})$ (Eq. 2). $${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{terr}}={\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{precip}}+{\mathrm{\Delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{e}}+{\mathit{\epsilon }}_{\mathrm{bio}\left(\mathrm{terr}\right)}$$ Systematic differences in ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ values have been observed for different plant types (especially between grasses and trees) (Diefendorf et al., 2011; Kahmen et al., 2013b), possibly indicating differences in either ${\mathit{\epsilon }}_{\mathrm{bio}}$ (Sachse et al., 2012) or the fraction of leaf water used for lipid biosynthesis (Kahmen et al., 2013b) or yet unidentified factors. As such, vegetation changes in sedimentary records have been suggested to affect ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ values and vegetation corrections have been proposed (Feakins, 2013).

Major model assumptions.

Assumption	Explanation
${\mathit{\delta }}^{\mathrm{2}}$H${}_{\text{lake water}}={\mathit{\delta }}^{\mathrm{2}}$H${}_{\text{mean annual precipitation}}$	Stable hydrogen isotope composition of lake water equals mean annual stable hydrogen isotope compositions of precipitation (source water), as observed for small catchment lakes in temperate environments (Moschen et al., 2005)
${\mathit{\epsilon }}_{\text{terr-aq}}=$ leaf water evaporative ${}^{\mathrm{2}}$H enrichment	Difference between terrestrial- and aquatic-plant-derived $n$-alkane ${\mathit{\delta }}^{\mathrm{2}}$H values equals evaporative deuterium enrichment of leaf water (Kahmen et al., 2013b; Rach et al., 2014a)
${\mathit{\epsilon }}_{\mathrm{bio}}=$ constant	Biosynthetic fractionation is constant for aquatic as well as terrestrial source organisms on temporal and spatial scales of sedimentary integration (Sachse et al., 2012)
No significant delay (i.e. below sample resolution, i.e decades) of terrestrial $n$-alkanes transfer from source organisms into lake sediment	Due to the very small catchment of MFM with steep and wind-sheltered crater walls, we can assume an almost instantaneous transfer of $n$-alkanes and pollen from source organisms to lake sediment. It is likely that autumn leaf litter is the main $n$-alkane source to the sediment. This is supported by the similar sample to sample (i.e. decadal) variability in the lipid ${\mathit{\delta }}^{\mathrm{2}}$H values. If, for example, terrestrial leaf wax $n$-alkanes had a substantially longer residence time in the soils before being transported into the lake, then the decadal variability would be much smaller, as the soil would already deliver a more integrated signal into the lake
$e_{\mathrm{atm}}=$ constant	The atmospheric pressure is inferred from the altitude above sea level (0 m $=$ 1013 hPa), which remained unchanged. Short-term weather-related fluctuations (on the order of 100 hPa) do not affect the model outcome (see text).
$T_{\mathrm{leaf}}=T_{\mathrm{air}}$	Leaf temperature equals air temperature on the timescale of sediment integration (decades) (Kahmen et al., 2011b)
${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{wv}}=-{\mathit{\epsilon }}_{+}$	Atmospheric water vapor equals equilibrium isotope fractionation between vapor and liquid, as often observed for long-term (several years) time series in temperate climates (Jacob and Sonntag, 1991)
No significant influence from Péclet effect	Variations in the Péclet effect are minimal over time, in particular for angiosperm species (Kahmen et al., 2009; Song et al., 2013)
Concentration of produced $n$-alkanes from monocots and dicots are almost equal	Both of our vegetation correction approaches assume that palynological reconstructions are representative of leaf-wax-producing plants and that both monocots and dicots produce similar quantities of $n$-alkanes.

Since evaporative ${}^{\mathrm{2}}$H enrichment of leaf water only affects terrestrial plants and not aquatic organisms, changes in sedimentary ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ (Sachse et al., 2006) can be seen as a record of variations in terrestrial evaporative ${}^{\mathrm{2}}$H enrichment over time. Thus, by combining Eqs. (1) and (2) under the assumption that ${\mathit{\epsilon }}_{\mathrm{bio}}$ of both aquatic and terrestrial organisms was constant on the temporal and spatial scales of sedimentary integration, the difference between ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{aq}}$ and ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ values should mainly reflect the evaporative ${}^{\mathrm{2}}$H enrichment of leaf water (Eq. 3). Whenever referring to an isotopic difference between two pools (such as ${\mathrm{\Delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{e}}$), we employ the mathematically correct epsilon formula to calculate differences between two $\mathit{\delta }$-values (Sessions and Hayes, 2005). For simplicity we use the following expression: $${\mathrm{\Delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{e}}={\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{terr}}-{\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{aq}}.$$ Variants of this concept (Sachse et al., 2004) have been used to qualitatively interpret changes in evapotranspiration through the isotopic difference between ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ and ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{aq}}$ (i.e. expressed as ${\mathit{\alpha }}_{\text{TA/wat}}$, ${\mathit{\delta }}^{\mathrm{2}}$H C${}_{\mathrm{23}}$–C${}_{\mathrm{31}}$ and ${\mathit{\epsilon }}_{\text{terr-aq}}$; Jacob et al., 2007; Rach et al., 2014a; Seki et al., 2011). With recent progress in understanding of the determinants of ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ values and the existing mechanistic understanding of the processes governing leaf water evaporative ${}^{\mathrm{2}}$H enrichment (Craig, 1965; Kahmen et al., 2011b; Sachse et al., 2012), we propose a new framework – which we term the dual-biomarker (DUB) approach – to extract quantitative hydrological information, namely changes in relative humidity ($\mathrm{\Delta }$RH) from sedimentary records. To illustrate the power of this approach with paleohydrological data, we combine compound-specific hydrogen isotope measurements with plant physiological modeling on a previously published late glacial record of ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{aq}}$ and ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ from sediments of the lake Meerfelder Maar (MFM), Germany (Rach et al., 2014a).

Approach and model

The key assumptions of the DUB approach (Table 1) are that the difference between terrestrial- and aquatic-plant-derived $n$-alkane ${\mathit{\delta }}^{\mathrm{2}}$H values (${\mathit{\epsilon }}_{\text{terr-aq}})$ equals evaporative deuterium enrichment of leaf water (Kahmen et al., 2013b; Rach et al., 2014a) over the timescale of sediment integration (i.e. decades in our case) and that ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\text{lake water}}$ equals ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\text{mean annual precipitation}}$, a condition fulfilled for small catchment lakes in temperate environments without any major inflow. Also, the temporal delay in transfer of terrestrial $n$-alkanes from source organisms into lake sediment should be below the temporal resolution of the samples, which is fulfilled for sites with a very small catchment area and steep terrain, such as maar lakes. Furthermore, we assume that the biosynthetic fractionation (${\mathit{\epsilon }}_{\mathrm{bio}})$ is constant for terrestrial and aquatic source organisms on temporal and spatial scales of sedimentary integration (Sachse et al., 2012). We also assume that palynological data represent lake catchment vegetation so that those can be used to assess source organisms of aquatic and terrestrial $n$-alkanes (Rach et al., 2014a; Schwark et al., 2002). To assess the influence of vegetation changes on our reconstructions, we employ two different vegetation corrections based on palynological data, for which we assume that the concentration of $n$-alkanes produced by these different plants is equal to the pollen produced by them.

These assumptions and additional data are needed to parameterize the model; therefore, we emphasize that a robust application of the DUB model requires a good understanding of the paleolake system and its environment. As such, the DUB model should only be employed at a site that fulfills the conditions presented above and where a number of additional, well-constrained proxy data exist. As of now, this limits the application of the DUB model to precipitation fed, small catchment (ideally maar or crater) lakes in temperate regions.

${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{aq}}$ in such systems can be regarded as a direct recorder of growing season average precipitation ${\mathit{\delta }}^{\mathrm{2}}$H values and ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ values largely reflect leaf water ${\mathit{\delta }}^{\mathrm{2}}$H values, as has recently been demonstrated for greenhouse- and field-grown plants (Kahmen et al., 2013a, b). Leaf water in turn is a function of the plant's source water and leaf water evaporative ${}^{\mathrm{2}}$H enrichment. We argue that soil water evaporation is negligible as recently suggested by several observational studies and a global assessment (Jackson et al., 1996; Jasechko et al., 2013; Kahmen et al., 2013a) and that precipitation is the ultimate water source of aquatic organisms and terrestrial plants. In terrestrial plants, however, the source water becomes more enriched in deuterium due to plant transpiration before it is used for lipid biosynthesis. As such, the isotopic difference between ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ and ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{aq}}$ (${\mathit{\epsilon }}_{\text{terr-aq}})$ can be attributed to mean leaf water evaporative ${}^{\mathrm{2}}$H enrichment (${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{e}})$ (Sachse et al., 2004). Based on recent field and greenhouse studies we further assume that ${\mathit{\epsilon }}_{\text{terr-aq}}$ captures a growing season signal, probably biased towards the earlier summer months in temperate climate zones as the majority of leaf waxes are produced during leaf development, with suggested integration periods of somewhere between weeks (Kahmen et al., 2013b; Tipple et al., 2013) and several months (Sachse et al., 2015).

The major variables controlling leaf water isotope enrichment are well understood and mechanistic models that allow the accurate prediction or reconstruction of leaf water ${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{e}}$ values based on environmental and physiological input variables have been developed based on the Craig–Gordon evaporation model (Craig, 1965) (Barbour, 2007; Farquhar et al., 2007; Ferrio et al., 2009; Kahmen et al., 2011b) (Eq. 4). $${\mathrm{\Delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{e}}={\mathit{\epsilon }}_{+}+{\mathit{\epsilon }}_k+\left({\mathrm{\Delta }}^{\mathrm{2}},H_{\mathrm{wv}},-,{\mathit{\epsilon }}_k\right){\displaystyle \frac{e_{\mathrm{a}}}{e_{\mathrm{i}}}}$$ ${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{e}}$ is determined by the equilibrium isotope fractionation between liquid and vapor (${\mathit{\epsilon }}_{+})$, the kinetic isotope fractionation during water vapor diffusion from the leaf intercellular air space to the atmosphere (${\mathit{\epsilon }}_k)$, the ${}^{\mathrm{2}}$H depletion of water vapor relative to source water (${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{wv}})$, and the ratio of atmospheric vapor pressure and intracellular vapor pressure ($e_{\mathrm{a}}/e_{\mathrm{i}})$ and air temperature ($T_{\mathrm{air}})$. In addition, leaf temperature ($T_{\mathrm{leaf}})$, stomatal conductance ($g_{\mathrm{s}})$ and boundary layer resistance ($r_{\mathrm{b}})$ are essential secondary input variables for the prediction of $e_{\mathrm{i}}$ and ${\mathit{\epsilon }}_k$, respectively. Reformulating Eq. (4) allows the expression of $e_{\mathrm{a}}$ as a function of Craig–Cordon variables (Eq. 5). Since the atmospheric vapor pressure ($e_{\mathrm{a}})$ can also be calculated based on RH and saturation vapor pressure ($e_{\mathrm{sat}})$ (Eq. 6), we can merge Eqs. (5) and (6) to calculate RH and to estimate quantitative changes in RH ($\mathrm{\Delta }$RH) (Eq. 7).

$$\begin{array}{}& {\displaystyle }& {\displaystyle }{e}_{\mathrm{a}}={\displaystyle \frac{e_{\mathrm{i}}\left({\mathrm{\Delta }}^{\mathrm{2}},H_{\mathrm{e}},-,{\mathit{\epsilon }}_{+},-,{\mathit{\epsilon }}_k\right)}{{\mathrm{\Delta }}^{\mathrm{2}}{H}_{\mathrm{wv}}-{\mathit{\epsilon }}_k}}& {\displaystyle }& {\displaystyle }\mathrm{RH}={\displaystyle \frac{e_{\mathrm{a}}\times \mathrm{100}\mathit{\%}}{e_{\mathrm{sat}}}}& {\displaystyle }& {\displaystyle }\mathrm{\Delta }\mathrm{RH}={\displaystyle \frac{e_{\mathrm{i}}\left({\mathrm{\Delta }}^{\mathrm{2}},H_{\mathrm{e}},-,{\mathit{\epsilon }}_{+},-,{\mathit{\epsilon }}_k\right)\times \mathrm{100}\mathit{\%}}{e_{\mathrm{sat}}\left({\mathrm{\Delta }}^{\mathrm{2}},H_{\mathrm{wv}},-,{\mathit{\epsilon }}_k\right)}}\end{array}$$ Equation (7) illustrates that $\mathrm{\Delta }$RH can be inferred from a record of past changes in ${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{e}}$ (i.e. a record of ${\mathit{\epsilon }}_{\text{terr-aq}})$ if the additional variables $e_{\mathrm{sat}}$, $e_{\mathrm{i}}$, ${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{wv}}$, ${\mathit{\epsilon }}_{+}$ and ${\mathit{\epsilon }}_k$ can be constrained. In the following we discuss the model parameterizations necessary to apply the DUB approach to estimate quantitative changes in RH from sedimentary records.

Saturation vapor pressure $e_{\mathrm{sat}}$ (Eq. 8) as well as the equilibrium fractionation factor ${\mathit{\epsilon }}_{+}$ (Eq. 9) are a function of temperature (all given numbers and physically variable dependencies within the equations are transferred from the Péclet-modified Craig–Gordon model by Kahmen et al., 2011b and the original leaf water enrichment model; Craig, 1965; Dongmann et al., 1974; Farquhar and Cernusak, 2005; Farquhar and Lloyd, 1993). The atmospheric pressure term ($e_{\mathrm{atm}})$, which is also needed for calculation of $e_{\mathrm{sat}}$, describes (mean annual) atmospheric pressure as a function of the elevation above sea level (0 m $=$ 1013 hPa). $$\begin{array}{cc}{\displaystyle }& {\displaystyle }{e}_{\mathrm{sat}}={\displaystyle \frac{\mathrm{1.0007}+\mathrm{3.46}\times e_{\mathrm{atm}}\left[\mathrm{hPa}\right]}{\mathrm{1}\mathrm{000}\mathrm{000}}}\times \mathrm{6.1121}\\ & {\displaystyle }& {\displaystyle }\times \exp \left({\displaystyle \frac{\mathrm{17.502}\times T_{\mathrm{air}}\left[{}^{\circ },\mathrm{C}\right]}{\mathrm{240.97}+T_{\mathrm{air}}\left[{}^{\circ },\mathrm{C}\right]}}\right){\displaystyle }& {\displaystyle }{\mathit{\epsilon }}_{+}=\left[\exp \right({\displaystyle \frac{\mathrm{24.844}\times \mathrm{1000}}{{\left(\mathrm{273.16},+,T_{\mathrm{air}},\left[{}^{\circ },\mathrm{C}\right]\right)}^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{76.248}}{\mathrm{273.16}+T_{\mathrm{air}}\left[{}^{\circ },\mathrm{C}\right]}}\\ & {\displaystyle }& {\displaystyle }+\mathrm{0.052612})-\mathrm{1}]\times \mathrm{1000}\end{array}$$ For accurate estimates of $e_{\mathrm{sat}}$ as well as ${\mathit{\epsilon }}_{+}$, information on air temperature ($T_{\mathrm{air}})$ during the growing season is required. Estimates of past $T_{\mathrm{air}}$ variability can be derived from paleotemperature proxy data to estimate $e_{\mathrm{sat}}$ and ${\mathit{\epsilon }}_{+}$ (e.g. chironomids, Heiri et al., 2014, 2007; MBT/CBT, Blaga et al., 2013). In particular chironomid records, thought to represent spring and summer temperatures, provide an ideal proxy of past mean growing season temperatures in this respect (Heiri et al., 2007). Note that $e_{\mathrm{sat}}$ also depends on the atmospheric pressure (Eq. 8), which can be estimated from elevation above sea level, and is treated as a constant in the model. Leaf-internal vapor pressure $e_{\mathrm{i}}$, however, is a function of leaf temperature ($T_{\mathrm{leaf}})$. We assume for our calculations that $T_{\mathrm{air}}$ is a good estimate of a growing season average $T_{\mathrm{leaf}}$ and $e_{\mathrm{i}}$ can thus be calculated as $$e_{\mathrm{i}}=\mathrm{6.13753}\times \exp \left(T_{\mathrm{air}},\left[{}^{\circ },\mathrm{C}\right],\times ,{\displaystyle \frac{\mathrm{18.564}-\frac{T_{\mathrm{air}}\left[{}^{\circ },\mathrm{C}\right]}{\mathrm{254.4}}}{T_{\mathrm{air}}\left[{}^{\circ },\mathrm{C}\right]+\mathrm{255.57}}}\right).$$ We are aware that $T_{\mathrm{leaf}}$ can exceed air temperature in situations of extreme drought, when transpiration and evaporative cooling is reduced, or in bright and sunny conditions (Leuzinger and Korner, 2007; Scherrer et al., 2011). However, on cloudy days as well as on days with wind, $T_{\mathrm{leaf}}$ typically equals $T_{\mathrm{air}}$ (Jones, 2013). Given the spatial and temporal integration of leaves in sedimentary records (covering decadal to millennial timescales), it is thus unlikely that single drought events, in which $T_{\mathrm{leaf}}$ would exceed $T_{\mathrm{air}}$, dominate the overall relationship between $T_{\mathrm{leaf}}$ and $T_{\mathrm{air}}$. Recent studies also show that for temperatures between 15 and 20 ${}^{\circ }$C the $T_{\mathrm{leaf}}$ equals $T_{\mathrm{air}}$ on seasonal timescales (Kahmen et al., 2011b).

Another parameter affecting leaf water isotope enrichment is the ${}^{\mathrm{2}}$H depletion of water vapor relative to source water (${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{wv}})$. In temperate climates liquid water and atmospheric water vapor are often in isotopic equilibrium, especially when longer (annual to decadal) timescales are investigated (Jacob and Sonntag, 1991). We therefore assume that ${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{wv}}$ equals the equilibrium isotope fractionation between vapor and liquid ${\mathit{\epsilon }}_{+}$. $${\mathrm{\Delta }}^{\mathrm{2}}{H}_{\mathrm{wv}}=-{\mathit{\epsilon }}_{+}$$ In the model, ${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{wv}}$ can thus be replaced by $-{\mathit{\epsilon }}_{+}$ (Eq. 11).

The kinetic isotope fractionation (${\mathit{\epsilon }}_k)$ depends on the plant physiological variable stomatal conductance ($g_{\mathrm{s}})$ and boundary layer resistance ($r_{\mathrm{b}})$ (Eq. 12) (Kahmen et al., 2011b). $${\mathit{\epsilon }}_k={\displaystyle \frac{\mathrm{16.4}\times \frac{\mathrm{1}}{g_{\mathrm{s}}\left[{\mathrm{m}}^{\mathrm{2}},,\mathrm{s},,{\mathrm{mol}}^{-\mathrm{1}}\right]}+\mathrm{10.9}\times r_{\mathrm{b}}\left[\mathrm{mol},,{\mathrm{m}}^{-\mathrm{2}},,{\mathrm{s}}^{-\mathrm{1}}\right]}{\frac{\mathrm{1}}{g_{\mathrm{s}}\left[\mathrm{mol},,{\mathrm{m}}^{-\mathrm{2}},,{\mathrm{s}}^{-\mathrm{1}}\right]}+r_{\mathrm{b}}\left[{\mathrm{m}}^{\mathrm{2}},,\mathrm{s},,{\mathrm{mol}}^{-\mathrm{1}}\right]}}.$$ No direct proxies exist to reconstruct these plant physiological variables from sedimentary records, but paleovegetation data can be used to parameterize the model with biome-averaged values for $g_{\mathrm{s}}$ and $r_{\mathrm{b}}$ that are inferred from modern plants (Klein, 2014). We note that these plant physiological variables exert only minor control on the model outcome, expected to lie within the analytical error of ${\mathit{\delta }}^{\mathrm{2}}$H lipid measurements (Kahmen et al., 2011b); see also discussion below.

Schematic overview showing the functional relationships between model variables of the DUB approach. Grey boxes on top mark the input parameters while the box size corresponds to the sensitivity of each variable in the result (small box – low influence on $\mathrm{\Delta }$RH; larger box – higher influence on $\mathrm{\Delta }$RH).

[Figure omitted. See PDF]

The latest iterations of leaf water models also include a Péclet effect, which describes the ratio of convectional versus diffusional flow of water in the leaf (Eq. 4) (Kahmen et al., 2011b). However, we did not include the Péclet effect in our calculations because we assume that variations in the Péclet effect are minimal over time (Kahmen et al., 2009; Song et al., 2013), in particular for angiosperm species.

When combining Eqs. (9), (10), (11) and (12) with Eq. (7), we obtain a model for $\mathrm{\Delta }$RH (Fig. 1) that requires only four major input variables: ${\mathit{\epsilon }}_{\text{terr-aq}}$, air temperature ($T_{\mathrm{air}})$, as well as literature-derived values for stomatal ($g_{\mathrm{s}})$ and boundary layer conductance ($r_{\mathrm{b}})$, and one constant parameter (site altitude above sea level for atmospheric pressure, $e_{\mathrm{atm}})$ to calculate $\mathrm{\Delta }$RH:

$$\begin{array}{cc}{\displaystyle }\mathrm{\Delta }\mathrm{RH}=& {\displaystyle }{e}_{\mathrm{i}}^{\prime }\left(T_{\mathrm{air}}\right)\times \left({\displaystyle \frac{{\mathrm{\Delta }}^{\mathrm{2}}{H}_{\mathrm{e}}}{-e_{\mathrm{sat}}^{\prime }\left(e_{\mathrm{atm}},,,T_{\mathrm{air}}\right)\left({\mathit{\epsilon }}_{+}^{\prime },\left(T_{\mathrm{air}}\right),+,{\mathit{\epsilon }}_k^{\prime },\left(g_{\mathrm{s}},,,r_{\mathrm{b}}\right)\right)}}\right.\\ & {\displaystyle }& {\displaystyle }\left..+,{\displaystyle \frac{\mathrm{1}}{e_{\mathrm{sat}}^{\prime }\left(e_{\mathrm{atm}},,,T_{\mathrm{air}}\right)}}\right)\times \mathrm{100}\mathit{\%}.\end{array}$$ Since we use ${\mathit{\epsilon }}_{\text{terr-aq}}$ ($={\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{e}})$ as an input variable, which is representative of leaf water isotope enrichment above source water and not absolute ${\mathit{\delta }}^{\mathrm{2}}$H leaf water values, Eq. (13) predicts changes in RH ($\mathrm{\Delta }$RH) but not RH directly. In theory, Eq. (13) would also allow the calculation of RH values directly, if absolute ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{precip}}$ and ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{leafwater}}$ were available. The current lack of experimentally determined biosynthetic fractionation factors for the respective aquatic and terrestrial plants prevents this approach from being used, but future experimental research may result in robust estimates of ${\mathit{\epsilon }}_{\mathrm{bio}}$, potentially enabling the reconstruction of absolute RH values (Zhang et al., 2009).

Uncertainties and sensitivity tests

Uncertainties

The DUB approach contains different variables (Fig. 1) with specific error ranges that can be quantified. These quantifiable errors (i.e. analytical uncertainties during isotope measurement or paleotemperature determination as well as ranges of values) can be used to set up an error propagation function and finally to provide an error range for the results (e.g. Eq. A1, Appendix). However, in addition to these quantifiable uncertainties there are still some catchment-related non-quantifiable uncertainties (see Table 1 and Sect. 2), which can increase the error of the results and therefore need to be taken into consideration before applying them to a certain catchment and/or record. These unquantifiable uncertainties can however be minimized through the selection of a particular well-characterized lacustrine archive, fulfilling the conditions we outlined in Sect. 2.

Sensitivity tests

To evaluate the robustness of our DUB approach for predicting $\mathrm{\Delta }$RH in the context of uncertainties, we tested the sensitivity of the model to uncertainties in the four key input variables $T_{\mathrm{air}}$, ${\mathit{\epsilon }}_{\text{terr-aq}}$, $g_{\mathrm{s}}$ and $r_{\mathrm{b}}$. In these sensitivity analyses we used a leaf water model, in which all secondary variables ($e_{\mathrm{i}}$, $e_k$, $e_{+}$ and $e_{\mathrm{sat}})$ are coupled to the primary input variables $T_{\mathrm{air}}$, $T_{\mathrm{leaf}}$, $g_{\mathrm{s}}$ and $r_{\mathrm{b}}$ (Kahmen et al., 2011b). We performed this test under a range of dramatically different climatic and ecological settings reflected by the climate conditions of Lista (Norway), Koblenz (Germany), Genoa (Italy) and Perth (Australia), which differ in mean growing season temperatures and prevailing vegetation types. While the vegetation in Norway and Australia is dominated by conifers and Mediterranean shrubland, respectively, the prevailing vegetation in Germany and Italy is broad leaf tree species. As baseline values for the sensitivity tests we set $T_{\mathrm{air}}$ in the analyses to the growing season mean temperatures of each site, which were 9.4, 15, 17.2 and 20.4 ${}^{\circ }$C for Lista, Koblenz, Genoa and Perth, respectively (IAEA/WMO, 2006). Leaf water evaporative enrichment ${\mathit{\epsilon }}_{\text{terr-aq}}({\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{e}})$ was set to 25 ‰ (Lista), 35 ‰ (Koblenz), 45 ‰ (Genoa) and 55 ‰ (Perth), which reflects average growing season leaf water enrichment values for the tested environments (Kahmen et al., 2013a). Baseline data for plant physiological variables were biome-typical estimates that we obtained from the literature (Jones, 2013; Klein, 2014): stomatal conductance ($g_{\mathrm{s}})$ for Lista and Koblenz was set to 0.25 mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, while for Genoa and Perth the preset values were 0.45 and 0.35 mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, respectively (Klein, 2014). Boundary layer resistance ($r_{\mathrm{b}})$ for Lista and Perth was set to 0.5 m${}^{\mathrm{2}}$s mol${}^{-\mathrm{1}}$, while for Koblenz and Genoa this variable was set to 1.0 m${}^{\mathrm{2}}$s mol${}^{-\mathrm{1}}$ (Jones, 2013).

Sensitivity analyses for major model input variables (${\mathit{\epsilon }}_{\text{terr-aq}}$, $T_{\mathrm{air}}$, $g_{\mathrm{s}}$ and $r_{\mathrm{b}})$ on resulting $\mathrm{\Delta }$RH values tested for four different climatic and ecological environments (Norway, Germany, Italy and Australia). Bars represent the effect on model output ($\mathrm{\Delta }$RH) for each tested environment and its variation when the respective input variable will be varied by the marked value. Missing bars (i.e. for negative $g_{\mathrm{s}}$ and $r_{\mathrm{b}})$ result from a bigger (negative) variation than the preset value (below 0).

[Figure omitted. See PDF]

The temperature sensitivity tests were performed by increasing and decreasing the respective $T_{\mathrm{air}}$ values for a location by 0.5, 1, 2 and 5 ${}^{\circ }$C (encompassing reconstructed temperature variations during the last major abrupt climate shift in western Europe – the Younger Dryas (YD) period with about 4–6 ${}^{\circ }$C; Goslar et al., 1995; Heiri et al., 2007). ${\mathit{\epsilon }}_{\text{terr-aq}}$ (${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{e}}$) values were varied by $\pm$5, 10, 15 and 20 ‰ for each location, which corresponds to evaporative leaf water enrichment in the test areas (spring months) (Kahmen et al., 2013a). Plant physiological variables ($g_{\mathrm{s}}$ and $r_{\mathrm{b}})$ were varied by $\pm$0.1, $\pm$0.2, $\pm$0.4 and maximum $\pm$0.6 mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ and $\pm$0.6 m${}^{\mathrm{2}}$ s mol${}^{-\mathrm{1}}$. These tested variations in plant physiological variables cover the expected variation in $g_{\mathrm{s}}$ and $r_{\mathrm{b}}$ for the local vegetation at the sites described in the sensitivity analysis.

The sensitivity analyses showed similar results for all four tested environments (Fig. 2). This suggests a similar behavior of the model under very different climate and ecological conditions. The DUB model is most sensitive to changes in ${\mathit{\epsilon }}_{\text{terr-aq}}$ (i.e. ${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{e}})$ and $T_{\mathrm{air}}$, while the plant physiological variables ($g_{\mathrm{s}}$, $r_{\mathrm{b}})$ showed only minor effects on $\mathrm{\Delta }$RH (Fig. 2). Specifically, a change of $\pm$20 ‰ in ${\mathit{\epsilon }}_{\text{terr-aq}}$ (i.e. ${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{e}})$ resulted in a change of $\pm$20 % in $\mathrm{\Delta }$RH. A $\pm$5 ${}^{\circ }$C change in $T_{\mathrm{air}}$ resulted in a 3 % change in $\mathrm{\Delta }$RH. Varying $g_{\mathrm{s}}$ and $r_{\mathrm{b}}$ within the specified limits caused only changes in $\mathrm{\Delta }$RH of 0.01 to 0.5 % (Fig. 2), suggesting low model sensitivity to plant physiological variables. A sensitivity test with variations in atmospheric pressure ($e_{\mathrm{atm}})$ of $\pm$100 hPa led to changes in $\mathrm{\Delta }$RH of 0.05 %. The difference in calculated $\mathrm{\Delta }$RH for sites with low (e.g. Lista) and high (e.g. Perth) growing season mean temperature were smaller than the regional model sensitivity of the different input variables and is therefore negligible. Our sensitivity analyses shows that the most critical variables for estimating changes in RH with our model are ${\mathit{\epsilon }}_{\text{terr-aq}}$ and $T_{\mathrm{air}}$ (Fig. 2).

Application: reconstructing quantitative changes in $\mathrm{\Delta }$RH during the Younger Dryas (YD) in western Europe

In general, there are two approaches to validating a climate proxy. The most straightforward way is to test the proxy under modern hydroclimate conditions through variations in space or time and compare results with actual instrumental data, either along a modern climatological gradient or over the time period when instrumental data are available. The second possibility is the analysis of a longer time series during a period with otherwise known major changes in the parameter to be tested for.

For testing the DUB model, the first approach is not feasible. While highly resolved (ideally annual laminated) lacustrine sediments from temperate Europe covering the instrumental period (roughly the last 150 years) do exist, no major changes in RH have occurred during this time. Using only (nonlaminated) core top sediments (i.e. only one data point integrating the last decade) would not allow for testing the performance of the DUB approach, which aims to reconstruct relative changes in RH, not absolute data. Testing the DUB approach along a modern climatic gradient is also difficult because we cannot assume that the source of aquatic biomarkers (in our case $n$C${}_{\mathrm{23}})$ is always the same aquatic macrophyte in different lakes and ecosystems (Sachse et al., 2004), i.e. it is unlikely to encounter enough lake systems where the sources of aquatic biomarkers are comparable and cover a large enough aridity gradient.

(a) ${\mathit{\delta }}^{\mathrm{2}}$H values of aquatic plants (${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{aq}}$, blue line) and higher terrestrial plants (${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$, green line; Rach et al., 2014a). (b) Terrestrial evapotranspiration (${\mathit{\epsilon }}_{\text{terr-aq}}$, orange line) during the Younger Dryas at MFM (Rach et al., 2014a). (c) Original chironomid-based temperature reconstruction from Hijkermeer (NL) (Heiri et al., 2007) (black line with X as data points) and interpolated temperature data for DUB approach (purple dots). (d) Variability in $\mathrm{\Delta }$RH during the YD cold period at MFM. The data are normalized to mean Allerød level (12 847–12 680 years BP). The bold line marks the moving average.

[Figure omitted. See PDF]

Therefore, we decided to employ the second approach, i.e. test the proxy during a period of known and significant changes in RH, such as the YD cold period (Rach et al., 2014a). The YD, the last major abrupt climatic shift in younger Earth's history (between 12 680 and 11 600 years BP), was characterized by a significant atmospheric temperature decrease of 4–6 ${}^{\circ }$C (Goslar et al., 1995; Heiri et al., 2007), a relocation of atmospheric circulation patterns (Brauer et al., 2008) and major hydrological changes (i.e. significantly drier conditions) and ecological variations (propagation of grass and reduction of tree vegetation) in western Europe (Brauer et al., 1999a; Litt and Stebich, 1999; Rach et al., 2014a). The relocation of atmospheric circulation patterns during northern hemispheric cooling led to drier conditions in western Europe. This forced changes in the regional vegetation composition (Brauer et al., 1999a, 2008; Rach et al., 2014a). For this period a high-resolution record of changes in ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{aq}}$ and ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ from a lacustrine archive, which fulfills the requirements outlined above (i.e. precipitation fed, a very small catchment, available palynological and other climate proxy data; Brauer et al., 1999a; Litt and Stebich, 1999), Lake Meerfelder Maar (MFM) in western Germany, exists. The presence of annual varves and a high temporal sampling resolution (decades) allows the evaluation of the timing of climatic and ecosystem changes – an ideal setting to illustrate the power of the DUB approach. A detailed description of the record and the available proxy data is given in Rach et al. (2014a). Briefly, the annually laminated sediments of MFM covering the YD period contain abundant aquatic ($n$C${}_{\mathrm{23}})$ and higher terrestrial ($n$C${}_{\mathrm{29}})$ lipid biomarkers ($n$-alkanes) (Fig. 3a). Based on the pollen record, the $n$C${}_{\mathrm{23}}$ alkane can be related to the aquatic submerged plant Potamogeton sp. and the $n$C${}_{\mathrm{29}}$ alkane to leaves originating from the terrestrial angiosperm trees Betula sp. and Salix sp. with input from grasses (Brauer et al., 1999a; Diefendorf et al., 2011). For the DUB approach we use the isotopic difference between ${\mathit{\delta }}^{\mathrm{2}}$H values of the $n$C${}_{\mathrm{29}}$ and of $n$C${}_{\mathrm{23}}$ alkanes (${\mathit{\epsilon }}_{\text{terr-aq}})$ (Fig. 3b) as a measure for leaf water ${}^{\mathrm{2}}$H enrichment (${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{e}})$.

Model parameterization for the MFM application

Temperature

Since no paleotemperature proxy data are directly available for MFM, we use a high-resolution chironomid-based temperature reconstruction from a nearby location, lake Hijkermeer in the Netherlands (Fig. 3c), ca. 300 km north of MFM (see the Appendix). The Hijkermeer record is interpreted as a record of mean July temperatures for western Europe with a mean error of about 1.59 ${}^{\circ }$C (Heiri et al., 2007). Since leaf wax synthesis most likely occurs during the early part of the growing season (spring and summer) (Kahmen et al., 2011a; Sachse et al., 2015; Tipple et al., 2013), the Hijkermeer record might slightly overestimate spring temperatures. However, when reconstructing $\mathrm{\Delta }$RH during the Younger Dryas, it is important that paleotemperature data capture the changes in temperature before and during that period, rather than absolute temperatures.

Plant physiological parameters

We estimated plant physiological variables ($g_{\mathrm{s}}$ and $r_{\mathrm{b}})$ based on literature data from the prevalent catchment vegetation inferred from available MFM pollen records (Brauer et al., 1999a; Litt and Stebich, 1999). These suggest that Betula sp. and Salix sp. were the dominant $n$C${}_{\mathrm{29}}$-producing taxa but that grasses became more abundant during the YD (Brauer et al., 1999a; Litt and Stebich, 1999). Reported $g_{\mathrm{s}}$ values for these species growing under humid to arid conditions today range from 0.1 to 0.5 mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ and boundary layer resistance ($r_{\mathrm{b}})$ values range from 0.95 to 1.05 m${}^{\mathrm{2}}$ s mol${}^{-\mathrm{1}}$ (Klein, 2014; Schulze and Hall, 1982; Schulze, 1986; Turner et al., 1984). We therefore used mean values, i.e. 0.3 mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ for $g_{\mathrm{s}}$ and 1.0 m${}^{\mathrm{2}}$ s mol${}^{-\mathrm{1}}$ for $r_{\mathrm{b}}$ as input variables for our modified model. We used the variance in $\pm$0.2 mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ for $g_{\mathrm{s}}$ and $\pm$0.1 m${}^{\mathrm{2}}$ s mol${}^{-\mathrm{1}}$ for $r_{\mathrm{b}}$ to calculate the error range of $\mathrm{\Delta }$RH. We note the low sensitivity of the DUB model outcome to variability in these variables (see Fig. 2, Appendix), as such that $\mathrm{\Delta }$RH changes of less than 0.1 % result from varying $g_{\mathrm{s}}$ values by 0.4 mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ or $r_{\mathrm{b}}$ values by 0.1 m${}^{\mathrm{2}}$ s mol${}^{-\mathrm{1}}$ (Fig. 2).

Estimation of uncertainty

The estimation of uncertainty for $\mathrm{\Delta }$RH is based on a linear error propagation (Eq. A1 in the Appendix) using specific error ranges for the individual input variables. For each input variable we used their individual reported or estimated error (i.e. for chironomid interfered temperature reconstruction: $\pm$1.5 ${}^{\circ }$C), for ${\mathit{\epsilon }}_{\text{terr-aq}}$ the analytical uncertainty (standard deviation) of the respective biomarker ${\mathit{\delta }}^{\mathrm{2}}$H measurements and for $g_{\mathrm{s}}$ and $r_{\mathrm{b}}$ the observed range of plant physiological parameters between different species ($g_{\mathrm{s}}$: 0.1–0.5 mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, $r_{\mathrm{b}}$: 0.95–1.05 m${}^{\mathrm{2}}$s mol${}^{-\mathrm{1}}$). The resulting average error for $\mathrm{\Delta }$RH estimation during the investigated interval is 3.4 % (see above and in the Appendix).

Model results for the YD period at MFM

Applying the DUB approach to the late glacial MFM record we can, for the first time, estimate the magnitude by which RH changed during a distinct period of abrupt climatic change in the past. Our quantification revealed substantial changes in RH on the order of 30 % (Fig. 3d) during the late glacial period, some of which occurred on multi-decadal timescales. To better illustrate these changes we normalized our results to the mean of the period between 12 847 and 12 680 BP (mean Allerød) (Fig. 3d), which is thought to have been warmer and moister than the YD (Hoek, 2009).

In particular, at the onset of the YD at 12 680 years BP, $\mathrm{\Delta }$RH decreased by 13 % $\pm$ 3.4 % over 112 years compared to mean Allerød level (Fig. 3d). During the YD (from 12 680 to 11 600 years BP) $\mathrm{\Delta }$RH values were on average 5 % $\pm$ 3.4 % lower compared to the mean Allerød level. Furthermore, in our high-resolution dataset we observe a division of the YD into two distinct phases: the first part of the YD (12 610–12 360 years BP) was characterized by low but relatively constant $\mathrm{\Delta }$RH (variability between $-$8 and $-$13 % and a mean of $-$10 %, compared to Allerød), whereas the variability in $\mathrm{\Delta }$RH increased after 12 360 years BP and ranged between $-$19 and $+$2 % with a mean of $-$8 % compared to Allerød mean values (Fig. 3d). Towards the termination of the YD we reconstructed a strong increase in $\mathrm{\Delta }$RH (up to $+$20 % above the Allerød level) over only 80 years. This increase started about 100 years before the YD – Holocene transition at 11 600 BP (Fig. 3d), indicating that hydrological changes lead major ecosystem changes, which formed the basis for the definition of the YD–Holocene boundary (Brauer et al., 1999a, b). The onset of the Holocene was characterized by substantial variability in $\mathrm{\Delta }$RH, with a strong increase followed by a decrease to mean Allerød levels 150 years after the transition. The reconstructed magnitude of changes, i.e. a ca. 9 % reduction in RH during the YD constitutes a shift from an oceanic to a dry summer climate, comparable to the difference in mean annual RH between central and southern Europe today (Center for Sustainability and the Global Environment (SAGE), 2002; New et al., 1999). The overall temporal pattern of reconstructed $\mathrm{\Delta }$RH changes is in good agreement with proxy data from western Europe (Bakke et al., 2009; Brauer et al., 1999a, 2008; Goslar et al., 1993), which indicate a shift to drier conditions due to a southward displacement of the westerly wind system channeling dry polar air into western Europe (Brauer et al., 2008; Rach et al., 2014a).

Our approach reveals for the first time that substantial changes in RH of up to 20 % can take place over very short timescales, i.e. several decades, leading to substantial changes in terrestrial ecosystems. While other proxy data reveal qualitative trends in aridification, our approach can be used to identify hydrological thresholds. Applied to high-resolution records, such as annually laminated lake sediments, the DUB approach can even be used to derive rates of hydrological changes and compare those with associated ecological changes (i.e. pollen records).

(a) Reconstructed $\mathrm{\Delta }$RH variability during the YD period (light grey shading), without vegetation correction (black line, $\mathrm{\Delta }$RH), with vegetation correction assuming a constant offset between C3 dicots and C3 monocots (blue line, $\mathrm{\Delta }$RH${}^{\ast }$), and with vegetation correction assuming different leaf water sensitivities among grasses and trees (red line, $\mathrm{\Delta }$RH${}^{\ast \ast }$). The shaded area marks the error range for $\mathrm{\Delta }$RH${}^{\ast \ast }$. (b) Relative distribution of trees and grasses in the catchment of MFM during the YD from pollen studies (Brauer et al., 1999a; Litt and Stebich, 1999). (c) Occurrence of Artemisia pollen in the catchment of MFM during YD (Brauer et al., 1999a; Litt and Stebich, 1999). Arrows highlight the contemporaneous major changes in $\mathrm{\Delta }$RH and Artemisia.

[Figure omitted. See PDF]

The effect of vegetation change on ${\mathit{\epsilon }}_{\text{terr-aq}}$ and the estimation of $\mathrm{\Delta }$RH

Numerous studies have established that vegetation changes can also affect the sedimentary leaf wax ${\mathit{\delta }}^{\mathrm{2}}$H record since significant differences in the net or apparent fractionation (${\mathit{\epsilon }}_{\mathrm{app}})$ between source water and lipid ${\mathit{\delta }}^{\mathrm{2}}$H values exist among different plant types, in particular between monocot and dicot (all grasses) plants (Kahmen et al., 2013b; Tipple et al., 2013). Since the YD period at MFM was characterized by an increased amount of grass, we tested how vegetation changes may affect $\mathrm{\Delta }$RH reconstructions through the DUB approach. For this we have developed two approaches to correct ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ values, based on either a constant offset between monocot and dicot ${\mathit{\epsilon }}_{\mathrm{app}}$ (Sachse et al., 2012) or a lower sensitivity of grass-derived leaf wax ${\mathit{\delta }}^{\mathrm{2}}$H values to leaf water isotope enrichment (Kahmen et al., 2013b). Both approaches assume that palynological reconstructions are representative of leaf wax producing plants and that both monocots and dicots produce similar quantities of $n$-alkanes.

We used available palynological data to quantify the relative distribution of major tree vegetation (Betula and Salix) and grasses over the investigated period (Fig. 4b), expressed as the fraction of trees and grasses, $f_{\mathrm{trees}}$ and $f_{\mathrm{grass}}$, assuming that leaf waxes and pollen share a similar transport pathway in this small, constrained crater catchment.

Correction – case 1 – constant difference in ${\mathit{\epsilon }}_{\mathrm{app}}$ between monocots and dicots

The first vegetation correction for reconstructed leaf water enrichment (${\mathit{\epsilon }}_{\text{terr-aq}}^{\ast }$) is based on the assumption of a constant offset in biosynthetic isotope fractionation (${\mathit{\epsilon }}_{\mathrm{bio}})$ between trees and grasses. Observational evidence shows that leaf wax lipid ${\mathit{\delta }}^{\mathrm{2}}$H values (${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}})$ from C3 monocots are on average 34 ‰ more negative than from C3 dicots (non-grasses) when growing at the same site (Sachse et al., 2012). This value is based on an observed mean difference between apparent isotope fractionation (i.e. the isotopic difference between source water and leaf wax $n$-alkanes, ${\mathit{\epsilon }}_{\mathrm{app}})$ values of C3 dicots ($-$111 ‰) and C3 monocots ($-$141 ‰) within a global dataset (Sachse et al., 2012).

The difference between monocot and dicot $n$-alkane ${\mathit{\delta }}^{\mathrm{2}}$H could potentially affect our modeled $\mathrm{\Delta }$RH values, especially since a 23 % increase in grass abundance in the MFM catchment during the YD has been suggested by pollen studies (Brauer et al., 1999a; Litt and Stebich, 1999).The causes for these differences in ${\mathit{\epsilon }}_{\mathrm{app}}$ have been hypothesized to be due to species-specific differences in biosynthetic fractionation (Sachse et al., 2012) or temporal differences in leaf wax synthesis during the growing season (Tipple et al., 2013). Both scenarios would result in a more or less constant isotopic offset between monocots and dicots growing under the same climatic conditions.

Assuming a mean isotopic difference of $-$34 ‰ between trees and grasses (Sachse et al., 2012), we calculated a vegetation-weighted correction value ($-\mathrm{34}\times f_{\mathrm{grass}})$ for each data point. This value is then subtracted from ${\mathit{\epsilon }}_{\text{terr-aq}}$, and results in the vegetation-corrected ${\mathit{\epsilon }}_{\text{terr-aq}}^{\ast }$ value (Eq. 14). Similar approaches for a pollen-based vegetation reconstruction have been recently proposed and applied (Feakins, 2013; Wang et al., 2013). $${\mathit{\epsilon }}_{\text{terr-aq}}^{\ast }={\mathit{\epsilon }}_{\text{terr-aq}}-(-\mathrm{34}\times f_{\mathrm{grass}})$$

Correction – case 2: different sensitivity to leaf water isotope enrichment in dicot vs. monocot leaf wax ${\mathit{\delta }}^{\mathrm{2}}$H values

The second vegetation correction (${\mathit{\epsilon }}_{\text{terr-aq}}^{\ast \ast }$) is based on the assumption that the isotopic difference between monocot and dicot leaf wax $n$-alkanes is not constant, but dependent on environmental conditions (Kahmen et al., 2013b). Previous greenhouse studies imply that the difference in ${\mathit{\epsilon }}_{\mathrm{app}}$ between dicots and monocots is variable depending on a change in humidity conditions (Kahmen et al., 2013b). In a high-humidity climate chamber treatment (80 % RH), monocots and dicots showed similar values for ${\mathit{\epsilon }}_{\mathrm{app}}$ ($-$220 and $-$214 ‰, respectively) whereas in a low-humidity treatment, ${\mathit{\epsilon }}_{\mathrm{app}}$ for monocots was substantially lower compared to dicots ($-$205 and $-$125 ‰, respectively) (Kahmen et al., 2013b), a finding that is in disagreement with the two hypotheses proposed above. Rather, the latter study hypothesized that grasses use a mixture of enriched leaf water and unenriched xylem water for lipid synthesis (Kahmen et al., 2013b). This hypothesis would imply that leaf wax $n$-alkane ${\mathit{\delta }}^{\mathrm{2}}$H values of monocots do not record the full magnitude of the evaporative leaf water enrichment signal, but only a fraction (Sachse et al., 2009). A recent greenhouse study on grass-derived $n$-alkane ${\mathit{\delta }}^{\mathrm{2}}$H values of a broad spectrum of C3 and C4 grasses support this idea (Gamarra et al., 2016). Gamarra et al. suggest that the differences between $n$-alkane ${\mathit{\delta }}^{\mathrm{2}}$H values from grasses and $n$-alkane ${\mathit{\delta }}^{\mathrm{2}}$H values from dicotyledonous plants are caused by an incomplete transfer of leaf water ${\mathrm{\Delta }}^{\mathrm{2}}$H to the $n$-alkanes. As such, a sedimentary record of $n$-alkanes derived partly from grasses would also underestimate mean ecosystem leaf water enrichment. Under dry conditions this fraction was estimated to be ca. 18 % for C3 grasses, based on one grass species (wheat) that was studied (Kahmen et al., 2013b). The data from Gamarra et al. show that for C3 grasses only 38–61 % of the leaf water evaporative ${}^{\mathrm{2}}$H-enrichment signal (depending on the species) was transferred to leaf wax $n$-alkane ${\mathit{\delta }}^{\mathrm{2}}$H values. To work with a conservative value and not to overestimate a potential leaf water enrichment signal in grass-derived $n$-alkane ${\mathit{\delta }}^{\mathrm{2}}$H values, we decided to use the data from Kahmen et al. (2013b) for the wheat C3 grass. As such, our correction approach rather underestimates changes in relative humidity and represents as such the lower limit of reconstructed changes.

Under the assumption of different sensitivities to leaf water isotope enrichment of $n$-alkane ${\mathit{\delta }}^{\mathrm{2}}$H values in monocot and dicot plants (Kahmen et al., 2013b), we developed a correction for ${\mathit{\epsilon }}_{\text{terr-aq}}$ based on the experimentally determined mixing ratio between leaf water and unenriched xylem water in wheat, a C3 grass (Kahmen et al., 2013b), by essentially weighing the fraction of grass cover with a factor of 0.18 (Fig. 4b) (Eq. 15): $${\mathit{\epsilon }}_{\text{terr-aq}}^{\ast \ast }=\left(f_{\mathrm{trees}},\times ,\mathrm{1},+,f_{\mathrm{grass}},\times ,\mathrm{0.18}\right)\times {\mathit{\epsilon }}_{\text{terr-aq}}.$$

Comparison of results from uncorrected (${\mathit{\epsilon }}_{\text{terr-aq}})$ and corrected (${\mathit{\epsilon }}_{\text{terr-aq}}^{\ast }$, ${\mathit{\epsilon }}_{\text{terr-aq}}^{\ast \ast }$) values

Results from the raw ($\mathrm{\Delta }$RH) and both vegetation-corrected scenarios ($\mathrm{\Delta }$RH${}^{\ast }$ and $\mathrm{\Delta }$RH${}^{\ast \ast }$) are within the calculated error range of 3.4 % of $\mathrm{\Delta }$RH (Fig. 4a) during the Allerød and the early Holocene, but diverge by up to 10 % during the YD, when C3 grass vegetation was estimated to have increased from 28 to 52 % in the catchment of MFM (Fig. 4b). Vegetation-corrected results (case 1, Fig. 4a) showed on average a 7 % stronger decrease for $\mathrm{\Delta }$RH${}^{\ast }$ and only a 2 % stronger decrease for $\mathrm{\Delta }$RH${}^{\ast \ast }$ compared to uncorrected results. As such, $\mathrm{\Delta }$RH${}^{\ast \ast }$ values (case 2) are within the error range of uncorrected $\mathrm{\Delta }$RH during the entire record.

Interestingly, both correction approaches, but in particular case 2, place the relatively large variability in uncorrected $\mathrm{\Delta }$RH at the onset and the termination of the YD, where abrupt vegetation changes occurred. For example, uncorrected $\mathrm{\Delta }$RH changes were predicted to be up to 35 % during the termination of the YD, corresponding to the modern gradient between western Europe and the semidesert areas in northern Africa (Center for Sustainability and the Global Environment, SAGE, 2002). Vegetation-corrected $\mathrm{\Delta }$RH${}^{\ast \ast }$ values were on the order of 20 %, seemingly more reasonably representing local late glacial changes (Fig. 4a).

Our analysis shows that vegetation changes have the potential to affect the DUB approach estimates, but a lack of mechanistic understanding of the causes of the differences in ${\mathit{\delta }}^{\mathrm{2}}$H${}_{\mathrm{terr}}$ between tree and grass vegetation (Sachse et al., 2012) makes an assessment of the validity of either (or any) correction approach difficult. Tentatively, the lower variability in $\mathrm{\Delta }$RH${}^{\ast \ast }$ within the YD as well as the less-pronounced shift in particular at the onset and termination of the YD (Fig. 4a) provides a more realistic scenario. However, as of now, we regard the differences in predictions as the error of quantitative predictions from the DUB approach. This uncertainty is larger during periods characterized by vegetation changes and in our case maximum differences in prediction of $\mathrm{\Delta }$RH between the Allerød and the YD are on the order of 11 % (mean Allerød vs. mean YD difference between $\mathrm{\Delta }$RH and $\mathrm{\Delta }$RH${}^{\ast }$).

Comparison of reconstructed $\mathrm{\Delta }$RH with other proxy data

We can further demonstrate the validity of our approach by direct comparison to other hydroclimate proxies from the MFM record. For example, a classical palynological marker for more arid conditions is Artemisia pollen (D'Andrea et al., 2003). In the MFM catchment a prominent increase in the occurrence of Artemisia has been used to infer drier conditions during the YD (Fig. 4c) (Brauer et al., 1999a; Bremer and Humphries, 1993; D'Andrea et al., 2003; Litt and Stebich, 1999). When comparing the abundance of Artemisia pollen as a percentage (note that the Artemisia abundance data are not part of the vegetation corrections discussed above) to the DUB $\mathrm{\Delta }$RH record, we observed striking similarities over the whole of the study period (Fig. 4a, c). Inferred wetter conditions during the second phase of the YD, or centennial-scale excursions to higher $\mathrm{\Delta }$RH (such as between 12 280 and 12 170 years BP) go in line with lower Artemisia pollen abundance after 12 100 BP. In fact, both independent datasets show an inverse, statistically significant relationship ($p$ < 0.001) (Fig. 5a–c), with high Artemisia pollen abundance during periods of low $\mathrm{\Delta }$RH values (Fig. 4a, c). The correlation between $\mathrm{\Delta }$RH and Artemisia is higher for vegetation-corrected $\mathrm{\Delta }$RH${}^{\ast }$ and $\mathrm{\Delta }$RH${}^{\ast \ast }$ (Fig. 5b, c) than uncorrected $\mathrm{\Delta }$RH; in particular, $\mathrm{\Delta }$RH${}^{\ast \ast }$ the variance of the dataset is greatly reduced (Fig. 5c), providing support for the hypothesis that vegetation changes could have affected the record.

Correlation plots of normalized reconstructed $\mathrm{\Delta }$RH vs. Artemisia population. (a) Uncorrected $\mathrm{\Delta }$RH values vs. Artemisia. (b) Vegetation-corrected $\mathrm{\Delta }$RH values ($\mathrm{\Delta }$RH${}^{\ast }$) vs. Artemisia. (c) Vegetation-corrected $\mathrm{\Delta }$RH values ($\mathrm{\Delta }$RH${}^{\ast \ast }$) vs. Artemisia.

[Figure omitted. See PDF]

Conclusions

We present a novel approach for quantifying paleohydrological changes (i.e. changes in relative humidity) combining sedimentary lipid biomarker ${\mathit{\delta }}^{\mathrm{2}}$H values from aquatic and terrestrial lipids with mechanistic leaf water isotope modeling. This dual-biomarker approach relies on the observation that aquatic and terrestrial organisms within the catchment of small lakes from temperate climate zones use distinct water sources, namely lake (i.e. precipitation) and ${}^{\mathrm{2}}$H-enriched leaf water as a source for their organic hydrogen. By taking advantage of the mechanistic understanding and available models of leaf water isotope enrichment in terrestrial plants, we show that it is possible to extract quantitative information about changes in relative humidity from sedimentary records.

Parameterizing and applying the DUB model to a lacustrine lipid biomarker ${\mathit{\delta }}^{\mathrm{2}}$H record from western Europe, we find strong and abrupt changes in RH at the onset and the termination of the YD occurring within the lifetime of a human generation. Specifically, our approach showed that shifts in RH of up to 13 % $\pm$ 3.4 % occurred within only 112 years. This dramatic change corresponds to shifts in average biome RH from oceanic to dry summer climates. Our quantification showed that dry conditions prevailed during the Younger Dryas period, with RH being between 8 and 15 % lower on average compared to the Allerød, depending on how the possible effect of vegetation changes is accounted for. The pattern and magnitude of our RH reconstruction agree well with other proxy data, such as the increase in the abundance of specific taxa adapted to dry conditions (e.g. Artemisia) during that time period.

Our analyses shows that the DUB approach is capable of quantifying past hydrological changes in temperate environments, when additional proxy data, especially on vegetation distribution and paleotemperature, exist. We suggest that this approach can be particularly valuable in the future for the validation of climate models and to better understand uncertainties in predictions of future hydrological change under global warming. However, we stress that the DUB approach relies on a number of assumptions and is currently limited by our incomplete understanding of processes affecting the transport and deposition of, in particular, terrestrial biomarkers from their source to the sedimentary sink. To minimize the arising uncertainties, this approach should only be applied to small catchment lake systems which are fed by precipitation in temperate climate zones, when biomarker sources can be constrained by paleovegetation data (such as palynological records). It is particularly crucial to constrain the aquatic biomarker source, but in principle any aquatic lipid biomarker (macrophyte and algal) could be employed. Our reconstruction provides reasonable values of RH changes, which are in agreement with ecosystem changes in the region, during the YD cold period. As such, the present approach provides a first step towards quantitative paleohydrological reconstructions.

The lipid-derived stable hydrogen isotope data are available on the PANGEA database (Rach et al., 2014b). Chironomid-based temperature data for the YD period are from Heiri et al. (2007). Vegetation data are from Brauer et al. (1999a) and Litt and Stebich (1999).

Error propagation

The uncertainty estimation ($\mathrm{\Delta }f$, Eq. A1) for the reconstructed $\mathrm{\Delta }$RH variability is based on a linear error propagation, which is the most conservative method for error estimations. This method does not require the same kind of considered errors and therefore provides the possibility of combining different kinds of errors with their specific ranges (i.e. measuring error, counting error). The individual error ranges of the independent variables in our approach arise from different sources such as analytical errors (chironomid interfered temperature reconstruction: $\pm$1.5 ${}^{\circ }$C), observed variations of plant physiological parameters between different species (stomatal conductance: 0.1–0.5 mol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, boundary layer resistance: 0.95–1.05 m${}^{\mathrm{2}}$ s mol${}^{-\mathrm{1}}$), and standard deviation of ${\mathit{\delta }}^{\mathrm{2}}$H measurements of terrestrial and aquatic $n$-alkanes.

The specific uncertainty for ${\mathit{\epsilon }}_{\text{terr-aq}}^{\ast \ast }$ was preliminarily determined by a separate error propagation using the (analytical) standard deviation of the triplicate measurements of the sedimentary $n$-alkane ${\mathit{\delta }}^{\mathrm{2}}$H values as well as the plant-derived $n$-alkane ${\mathit{\delta }}^{\mathrm{2}}$H measurements by Kahmen et al. (2013b). The results of these separate error estimations were integrated into the general error estimation of $\mathrm{\Delta }$RH${}^{\ast \ast }$.

In contrast to the linear error propagation, a less conservative method (Gaussian error propagation) requires a similarity of the errors, i.e. all errors are measurement or counting errors, which is not the case in this study. The mean error when using the Gaussian method is however only 3.2 % and therefore only 0.2 % smaller than the calculated error using the linear propagation. $$\begin{array}{cc}{\displaystyle }\mathrm{\Delta }f=& {\displaystyle }\left|{\displaystyle \frac{\partial \mathrm{RH}}{\partial {\mathit{\epsilon }}_{\text{terr-aq}}^{\ast \ast }}}\right|\times \mathrm{\Delta }{\mathit{\epsilon }}_{\text{terr-aq}}^{\ast \ast }+\left|{\displaystyle \frac{\partial \mathrm{RH}}{\partial r_{\mathrm{b}}}}\right|\times \mathrm{\Delta }{r}_{\mathrm{b}}+\left|{\displaystyle \frac{\partial \mathrm{RH}}{\partial g_{\mathrm{s}}}}\right|\\ & {\displaystyle }& {\displaystyle }\times \mathrm{\Delta }{g}_{\mathrm{s}}+\left|{\displaystyle \frac{\partial \mathrm{RH}}{\partial T_{\mathrm{air}}}}\right|\times \mathrm{\Delta }{T}_{\mathrm{air}}\end{array}$$

Temperature data

The temperature data used for the DUB model parameterization of the MFM case were taken from Heiri et al. (2007) and constitute reconstructed summer temperatures based on chironomid analyses from Hijkermeer (NL) (Heiri et al., 2007), which, to our knowledge, constitutes the late glacial paleotemperature record closest to the MFM site (distance 311 km). However, the dataset of Hijkermeer consists only of 37 data points between 13 000 and 11 000 BP with a temporal resolution varying between 26 and 167 years per sample. Therefore, we determined a new equidistant time series for the temperature data, fitting data volume and temporal resolution of our ${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{e}}$ record from MFM (106 data points with an 8 to 33-year resolution). For calculating the equidistant time series we used the “interpl” method with the specification “linear” in MATLAB (version R2010b).

Vegetation data

Information about late glacial vegetation cover in the catchment area of MFM is based on palynological analyses (Brauer et al., 1999a; Litt and Stebich, 1999). We also used pollen percent data for determining the vegetation distribution between trees and grasses for each data point. For using these vegetation data in our model it was necessary to determine an equidistant time series according to the age model of our ${\mathrm{\Delta }}^{\mathrm{2}}$H${}_{\mathrm{e}}$ values. For calculating these time series we also used the “interpl” method with the specification “linear” in MATLAB (version R2010b).

OR conducted model modifications, calculations and contributed to writing the paper. AK provided the basic leaf water enrichment model and was responsible for the plant physiological part and contributed to writing the paper. AB was responsible for lake coring, provided the chronology and stratigraphy for the Younger Dryas hydrological reconstruction, and contributed to writing the paper. DS conceived the research, acquired financial support and contributed to writing the paper.

The authors declare that they have no conflict of interest.

Acknowledgements

This work was supported by a DFG Emmy-Noether grant (SA1889/1-1) and an ERC Consolidator grant (no. 647035 STEEPclim) to Dirk Sachse. It is a contribution to the INTIMATE project, which was financially supported as EU COST action ES0907 and to the Helmholtz Association (HGF) Climate Initiative REKLIM Topic 8, rapid climate change derived from proxy data, and has used infrastructure of the HGF TERENO program. The article processing charges for this open-access publication were covered by a Research Centre of the Helmholtz Association. Edited by: Erin McClymont Reviewed by: Sarah Feakins and two anonymous referees