Millennial-length forward models and pseudoproxies of stalagmite δ18O: an example from NW Scotland

The stable oxygen isotope parameter δ18O remains the most widely utilised speleothem proxy for past climate reconstructions. Uncertainty can be introduced into stalagmiteδ18O from a number of factors, one of which is the heterogeneity of groundwater flow in karstified aquifers. Here, we present a lumped parameter hydrological model, KarstFor, which is capable of generating monthly simulations of surface water – ground water – stalagmite δ18O for more than thousand-year time periods. Using a variety of climate input series, we use this model for the first time to compare observational with modelled (pseudoproxy) stalagmite δ18O series for a site at Assynt, NW Scotland, where our knowledge ofδ18O systematics is relatively well understood. The use of forward modelling allows us to quantify the relative contributions of climate, peat and karst hydrology, and disequilibrium effects in stalagmite δ18O, from which we can identify potential stalagmiteδ18O responses to climate variability. Comparison of the modelled and actual stalagmite δ18O for two stalagmites from the site demonstrates that, for the period of overlapping growth, the two series do not correlate with one another, but forward modelling demonstrates that this falls within the range explicable by differences in flow routing to the stalagmites. Pseudoproxy δ18O stalagmite series highlight the potential significance of peat hydrology in controlling stalagmiteδ18O over the last 1000 yr at this site.


Introduction
Paleoclimate reconstructions from stalagmite and peat proxy archives in NW Scotland have been the focus of research for the last decade (Proctor et al., 2000(Proctor et al., , 2002;;Baker et al., 2002;Charman et al., 2001;Trouet et al., 2009), benefiting from the frequent occurrence of annually laminated stalagmites in the Uamh an Tartair cave system (Baker et al., 1993).Of additional significance is the correlation between rainfall amount in the region and the winter North Atlantic Oscillation index (wNAO; Hurrell, 1995).Previous speleothem studies have focused on analyses of the annual growth rate record; however, recently a high-resolution (annual to decadal) 1000-yr stable oxygen isotope record has been obtained from an annually laminated stalagmite (Baker et al., 2011).This sample was deposited for ∼ 1000 yr until 1996 AD and is from the same stalagmite that has previously provided annual resolution climate reconstructions of local rainfall and wNAO from variations in annual growth rate (Proctor et al., 2000;Trouet et al., 2009) and a proxy archive of severe winters (Baker et al., 2002).Cave climate and hydrological monitoring demonstrates that deposition close to isotopic equilibrium for δ 18 O is highly likely, and hence stalagmite δ 18 O directly reflects the composition of infiltrating water during periods of hydrologically effective precipitation (Fuller et al., 2008).However, Baker et al. (2011) failed to find any correlation over either the instrumental period (the 20th century) or that of reconstructed climate (proxy and instrumental; 17th century to present) between stalagmite δ 18 O and precipitation, temperature or atmospheric circulation, preventing a simple paleoclimate interpretation of the stalagmite δ 18 O proxy.Hence, a major goal of this paper is to advance our understanding of the δ 18 O proxy record by comparing stalagmite δ 18 O with pseudoproxy δ 18 O series generated using a hydrological forward modelling approach.
Recent developments in cave drip water hydrological modelling have quantified the uncertainty in stalagmite δ 18 O time series that arises as a result of variations in the hydrological routing of water from the surface to the cave (Baker and Bradley, 2010;Bradley et al., 2010;Wackerbarth et al., 2010).This permits an alternative approach to the assessment of stalagmite δ 18 O series in the form of forward models and the generation of pseudoproxies (Mann and Rutherford, 2002;Mann et al., 2005;Moberg et al., 2008).Forward modelling is a process in which the climate series is taken as an input to a process-based model, as opposed to the target of an empirical regression equation (Sturm et al., 2010).Driven by input climate data, forward models take an integrated Earth systems approach to predict proxy series.As such, they are dependent on good quality (i.e.long and continuous) measured or derived series, often including δ 18 O of precipitation.Where measured series are unavailable, a pseudoproxy series can be generated by a forward model, in which the model input parameters have a mean, a variance and an autocorrelation of that to be expected over the model period.At the time of writing, this approach has been postulated as a suitable test of stalagmite and other proxy series, but this has been rarely implemented to-date (for a review and example see Christiansen et al., 2009).More generally, attempts have been made to quantify the uncertainty in proxy climate reconstructions using pseudoproxies that replicate the statistical properties of the proxy series (Mann et al., 2005;Christiansen et al., 2009).For example, Moberg et al. (2008) used the statistical time series properties of different proxies (used in the Moberg et al. (2005) reconstruction of the last 1000 yr climate) to generate a pseudoproxy reconstruction of climate over the time period, the pseudoproxy series having the same variance as the original proxy records.Such approaches may be valid for proxy climate series where their statistical properties can be readily replicated, but given that non-linear hydrological behaviour of drip water flow within karstified aquifers can lead to non-linear responses (Genty and Deflandre, 1998;Baker and Brunsdon, 2003), we prefer a modelling approach that can replicate the actual processes of karst groundwater flow and which can be inferred from field evidence.
Stalagmite forward models published to date have focused on modelling stalagmite or major ion chemistry over 10 1 -10 2 -yr timescales.The first karst forward model, described by Fairchild et al. (2006a), comprised a two-layer hydrochemical model which was used to capture flow and geochemical variability in the soil and karst at Brown's Folly Mine (UK Jurassic oolitic limestone) in SW England.A linear reservoir approach was used in which matrix and macropore soil components routed water to a single output draining a single reservoir.In this instance, the focus of the model was solute delivery rather than water isotopes over the 10 1yr timescale, although the expected oxygen isotope output from this model (derived from Wallingford rainfall input) was presented in Fairchild et al. (2006b, Fig. 14c).Baker and Bradley (2010) also used a linear reservoir approach to model speleothem δ 18 O over century-length simulations.In this case, a linear reservoir of total storage volume S was considered to have variable capacity, initial volume of water (St i ) and outflow rate (v).Drip waters were presumed to be fed from the integrated reservoir volume.Further developments of this approach include Baker et al. (2010) who applied a single reservoir model to an Ethiopian stalagmite that was assumed to be fed by waters following two distinct flow paths.Here, the first pathway was assumed to represent diffuse seepage through the karst matrix, while the second pathway comprised fracture-fed flow, on the basis of the observed characteristics of δ 18 O and the growth rate variability of the speleothem being modelled.Water flows in each pathway were parameterised via specified hydrological thresholds, and both seepage and fracture-fed components of flow were mixed in the reservoir, and drip-water was modelled as an "overflow" type (sensu Tooth and Fairchild, 2003) to reflect observed seasonal drying of the drip water.
This approach was continued by Bradley et al. (2010) who developed a monthly lumped parameter karst hydrology model (KarstHydroModel) specifically to attribute the uncertainty arising as a result of the dynamics of karst water movement and storage in the soil, epikarst and karst.Waters percolating through the karst were envisaged to move between five discrete water stores: (i) a surface store (that retained precipitation when air temperatures were < 0 • C); (ii) the epikarst; (iii) a principal water store (Store 1); (iv) an underflow store; and (v) an overflow store.Water movement between stores occurred at varying rates depending upon hydrologically effective precipitation in any month and the store capacity.In parallel to the hydrological model, the isotopic composition of each water store was modelled as a function of precipitation δ 18 O, the δ 18 O of the water store, and the quantity of water in the store.The theoretical ranges in drip water and stalagmite δ 18 O were determined for three regions of contrasting climate (NW Scotland, Ethiopia and Gibraltar) using actual climate data (monthly mean temperature; total precipitation and evapotranspiration and mean precipitation δ 18 O) and a common set of model parameters, over the last ∼ 45 yr.
In this paper, we present a combination of the forward modelling, cave monitoring and pseudoproxy approaches to improve our understanding of the 1000-yr δ 18 O record preserved in stalagmite SU-96-7 from Assynt in NW Scotland (Baker et al., 2011) as well as a previously unpublished record from stalagmite SU-03-2 from the same cave (Fuller, 2007).We refine the recently developed KarstHy-droModel described above (Bradley et al., 2010), to forward model at monthly resolution stalagmite drip water δ 18 O over 10 3 -yr timescales for the first time.Model parameterisation is chosen to fit our process-based understanding of water movement, observations of cave climate and cave hydrological and geochemical monitoring data (undertaken occasionally between 1995(undertaken occasionally between -2002(undertaken occasionally between and intensively between 2003(undertaken occasionally between and 2005 AD; AD;Fuller et al., 2008).For the last 1000 yr, multiple pseudoproxy δ 18 O series are then generated using a variety of climatologically realistic datasets.The aim of the paper is to compare these pseudoproxy series against the δ 18 O records from stalagmite SU-96-7 and SU-03-2, and hence improve our understanding of the different processes that determine the stalagmite δ 18 O series.

Site description and methods
Stalagmites SU-96-7 and SU-03-2 were sampled from Uamh an Tartair (UAT; translated as "Cave of the Roaring") at Assynt, NW Scotland, in 1996.Site descriptions of cave environment, surface climate, and overlying soil and vegetation have been published previously (Proctor et al., 2000;Charman et al., 2001;Fuller et al., 2008), are detailed in Fuller (2007), and are summarised here and in Fig. 1.UAT is situated in the upper catchment of a peatland-dominated river basin.The geology at the site is Cambro-Ordovician dolomite, which has been completely cemented, tectonized and intruded by a sill, allowing water movement only through fractures or micro-fractures, i.e. with no matrix flow.Directly above UAT, the soil cover comprises mineral soils and localised shallow peats that have accumulated in small pockets generally < 100 m 2 in area and with a maximum depth not exceeding ∼ 1 m.Several sink holes occur above the cave which focus surface drainage, the largest of which approximately overlies the "The Grotto", where SU-96-7 and SU-03-2 were sampled.The peatland is ombrotrophic and radiocarbon dates demonstrate that it is relatively young, and that 80 cm of peat has accumulated above the site over the last ∼ 2300 yr (approximately 1 cm 30 yr −1 ; Charman et al., 2001).A sampling well, installed in the peat immediately above the cave and monitored over 2004-2005 AD, indicated that the depth to the water table is highly variable over time and spans the entire depth of the peat profile over the hydrological year.Typically, water levels rise to the surface of the peatland in winter but fall to the base of the peat deposit in summer.Thus, the available peat water level monitoring data suggest that periods of soil moisture deficit are likely to occur annually but only in the warmest summer months (JJA) and that subsequent rainfall produces an autumn flush of dissolved organic material to the groundwater.Although not determined, the peat at UAT is likely to have a high specific yield; Boelter (1965) found the specific yield (i.e. the drainable porosity) of surficial moss peats in a comparable peatland in Minnesota to be between 0.79 and 0.52 falling to 0.22 with increasing depth and decomposition.Hence, the organic deposits represent a substantial and almost continuous store of water at this site.
Based on records for the period 1961-1990, the climate at the site is oceanic with > 1900 mm annual rainfall, with 250-270 rain days (> 0.2 mm in a 24-h period), 4-6 snow days (fresh snow lying at 09:00 LT), and an average of 77 % cloud cover over the year.Surface and cave temperature monitoring over 2003-2005 AD demonstrate a mean temperature of 7.2 • C, and within the Grotto of Uamh an Tartair cave, air temperature has an annual range of 4.7 • C, due both its proximity to the base of the overlying peat, as well as possible heat transport from the small stream within the lower levels of the chamber (Fuller et al., 2008).Water levels respond rapidly to rainfall due to a combination of a nearby upstream sink and a downstream constriction.Despite the presence of the stream, the Grotto is poorly ventilated, with no detectable air movement.Drip water hydrology and isotope data from a selection of drip water sites including that supplying SU-96-7 and SU-03-2 were collected over the same time period; mean drip water δ 18 O of all sites was identical to the weighted mean of monthly precipitation samples.The temporal variability of drip water δ 18 O was notable for the convergence of all monitored drip site δ 18 O in autumn 2004 and autumn 2005 (Fuller et al., 2008), (Fig. 1).This would indicate a flush of homogenous water from the overlying peat each autumn, although a variability of drip water δ 18 O of 0.3-0.4‰ at these times suggests a more complex hydrological routing may be occurring.The timing of drip waters feeding stalagmite SU-96-7 was observed to range between 25 and 35 min per drip over the period 2003-2005 AD, while drip waters for stalagmite SU-03-2 ranged between 2 and 25 min per drip.
The annual growth rate record from SU-96-7 provides a proxy archive of rainfall and the North Atlantic Oscillation, with higher growth rates in drier conditions that are also associated with a lowering of the water table and increased production of soil CO 2 (Proctor et al., 2000) (Fig. 1).During particularly severe winters, annual growth laminae exhibit a "doublet" structure, which is interpreted as a second flush of organic matter following thawing of the peat (Baker et al., 2002).Most recently, Mariethoz et al. (2012) analysed the growth increment (acceleration) of 11 annually laminated stalagmite records, including SU-96-7, and demonstrated that all exhibit a "flickering" in growth around a quasi-equilibrium value.This inability of the stalagmite to maintain continuous growth acceleration is interpreted hydrologically as demonstrating that stalagmite growth rate is controlled by a storage or reservoir, which is perturbed by variable annual inputs.Interestingly, SU-96-7 exhibits one period of non-flickering in growth, which occurs at 1615 AD and lasts ∼ 30 yr, which could be interpreted as a change in hydrological routing to the stalagmite.
Stalagmite SU-03-2 is also annually laminated, with continuous fluorescent laminae preserved from 1679 AD to the date of sampling (2003 AD); the annual lamina series has been previously published in Baker et al. (2008).Drip water supplying the stalagmite has also been analysed (Fuller et al., 2008;Baker et al., 2011).The drip water hydrology, isotope and inorganic element geochemistry demonstrates that this stalagmite has, in comparison to SU-96-7, a higher mean discharge (0.5 µl s −1 vs. 0.05 µl s −1 ); more variable drip rate (56 % relative standard deviation); and, higher dissolved calcium and magnesium concentrations (inferred from electrical conductivities of 547 vs. 469 µS cm −1 ).This suggests that stalagmite SU-03-2 has a different flow pathway (higher dissolved calcium, more variable discharge) and may potentially have a greater sensitivity to within-cave fractionation effects (correlation between discharge and δ 13 C of dissolved inorganic carbon; Baker et al., 2011).To obtain a high-resolution climate record, calcite samples were obtained from stalagmite SU-03-2 using a micromill at the University of Innsbruck.Samples were milled at 100 µm intervals and approximately 150 µm depth, creating trenches 2 mm wide along the growth axis.To enable the milled samples to be aligned to the lamina chronology, every 10th sample was milled to a wider trench.Sample SU-03-2 was milled down the central growth axis from the top to within ∼ 2 mm of a visible growth hiatus.The calcite powders were analysed by conventional acid-digestion methods as outlined in Spötl and Vennemann (2003).C and O isotope values are reported on the VPDB scale.The 1 σ error is ±0.06 ‰ and ±0.08 ‰ for δ 13 C and δ 18 O, respectively.The annual fluorescent laminae in the milled section of the stalagmite were reimaged using identical techniques as previously published (Proctor et al., 2000(Proctor et al., , 2002) ) and each 100 µm isotope analysis aligned to a lamina year for the last 331 yr.

Comparison of stalagmite δ 18 O series and climate correlations
The previously unpublished δ 18 O series for stalagmite SU-03-2 from Fuller ( 2007) is presented in Fig. 2, along with the published record for SU-96-7 (Baker et al., 2011).The δ 13 C series are also given for both stalagmites, as is the δ 18 O vs. δ 13 C relationship for each sample.Comparison of δ 18 O series (Fig. 2a) demonstrates that mean δ 18 O in SU-03-2 is lower than SU-96-7 over the period of overlap of the last ∼ 330 yr (−4.8 ± 0.4 vs. −4.5 ± 0.4, statistically different at 99 % confidence interval, Student's t-test), and the two stalagmites exhibit a similar δ 18 O range (−4.0 to −5.8 ‰ vs. −3.8 to −5.4 ‰), despite having different temporal sampling resolutions.Figure 2a demonstrates that a divergence in δ 18 O occurs between the two stalagmites from ∼ 1630 AD to ∼ 1720 AD, and ∼ 1800 to ∼ 1860 AD, with δ 18 O > 1 ‰ higher in SU-96-7 in both instances.δ 13 C also diverges to lower values at the same periods (Fig. 2b).Comparison of the two δ 13 C series demonstrates that mean δ 13 C for SU-03-2 is lower than SU-96-7 (−11.5 ± 0.6 ‰ vs. −10.6 ± 0.5 ‰), statistically different at the 99 % confidence interval (Student's t-test).δ 13 C has a greater range in SU-03-2 (−9.6 to −13.0 ‰ vs. −9.8 to −12.0 ‰), indicative of the higher sampling resolution and likely revealing the influence of kinetic effects due to drip rate variations over inter-annual timescales.Plotting δ 13 C vs. δ 18 O for both samples (Fig. 2c) demonstrates a correlation that is similar in strength and gradient for both SU-03-2 (r = 0.56, gradient = 0.37±0.03)and SU-96-7 (r = 0.69, gradient = 0.44±0.03).A correlation between δ 13 C vs. δ 18 O could indicate the dominance of a fractionation process controlling these isotopes.Isotopic fractionation during speleothem formation can be caused by factors such as cave ventilation, the degassing process, possible equilibriation between solution and cave atmosphere and evaporation (for a review, see Fairchild and Baker, 2012, 273-281), and these processes are increasingly well understood and modelled (Romanov et al., 2008;Mühlinghaus et al., 2007Mühlinghaus et al., , 2009;;Dreybrodt and Sholz, 2011).Baker et al. (2011) demonstrated from modern cave monitoring of drip water chemistry that kinetic fractionation is relevant in SU-96-7 for δ 13 C, but not δ 18 O; given the lower mean δ 18 O in SU-03-2 compared to SU-96-7, we consider the same will apply for SU-03-2.We investigated the strength of the relationship between δ 18 O and various climate parameters.For the period 1880 to 2003 AD, we investigated the correlation with monthly, seasonal and annual temperature and precipitation, using the local adjusted instrumental series of Proctor et al. (2000), updated to 2003 AD by Fuller et al. (2008).We also calculated monthly, seasonal and annual water excess (precipitation minus evapotranspiration) following Thornthwaite (1948).All correlations with temperature, precipitation and water excess were weak; the strongest correlation was with total annual water excess (r = −0.27)(Fig. 3a).A negative correlation between δ 18 O and effective recharge is typical of what would be expected from a precipitation δ 18 O "amount effect".However, the correlation is not stationary with time: split analysis of the instrumental data demonstrates that this relationship was strongest in the earlier half of the instrumental series.
We undertook further analyses using the 20th Century Reanalysis Project data provided by NOAA/OAR/ESRL PSD, Boulder, Colorado, USA (from their Web site: http://www.esrl.noaa.gov/psd/;Compo et al., 2006Compo et al., , 2011)), which allowed us to investigate additional climate relationships with atmospheric pressure and zonal and meridional winds.Again, correlations were weak, with the strongest correlations being between both local sea-level pressure and 500 mb geopotential height and stalagmite δ 18 O. Figure 3b presents the spatial correlation between δ 18 O and sea level pressure, demonstrating a positive relationship between sea-level pressure and δ 18 O.This is in agreement with the west-east δ 18 O gradient observed across European in Holocene stalagmites (McDermott et al., 2011) and reflects the importance of Atlantic vs. continental sourced moisture for the isotopic composition of precipitation in the region.Figure 3c shows the correlation between δ 18 O and sea-level pressure in the region 60-65 • N, 0-10 • W. Interestingly, correlations with sea-level pressure were less strong when other reanalysis products of shorter time duration (post 1948 AD) were used, indicating that the correlation is not stationary over time.
Finally, we also considered the relationship between δ 18 O and averaged climate data to investigate possible water mixing and storage in both soil and groundwater.Figure 3d demonstrates that the best correlations between δ 18 O and annual water excess and mean annual sea-level pressure occurred with an average of the previous 3 to 6 yr of the climate series, with correlations of up to ∼ 0.5 between sea level pressure and δ 18 O.More complex mixing models (Baker and Bradley, 2010;Jex et al., 2010) did not improve any of the correlations with sea-level pressure or water excess, nor did a multiple regression of water excess and sea-level pressure against δ 18 O.
Overall, comparison with instrumental and reanalysis product data suggests that δ 18 O in stalagmite SU-03-2 does contain a climate signal that can be related to rainfall source (sea-level pressure) and amount (water excess), but that this explains < 15 % of the data at an annual timescale and < 25 % at a 3-6 yr timescale, and is also non-stationary over time.A similar lack of climate correlation was observed in the δ 18 O record of stalagmite SU-96-7 (Baker et al., 2011), suggesting that, for this region and this time period at least, stalagmite δ 18 O is not sensitive to climate variability.Additionally, stalagmites SU-03-2 and SU-96-7 have different mean δ 18 O and exhibit different time series over their period of overlap during the last ∼ 330 yr, which requires further investigation.Therefore in this paper, we develop a forward modelling approach to generate stalagmite δ 18 O pseudoproxy series, which we can use to investigate sources of uncertainty in stalagmite δ 18 O.

KarstFor forward model
The setting of UAT differs in a number of respects from the generalised context envisaged in the karst hydrology forward models described by Baker and Bradley (2010), Baker et al. (2010) and Bradley et al. (2010), and the model configuration was revised to correspond with our understanding of the hydrology of the UAT cave system.In particular, stalagmite SU-96-7 is characterised by a slow but continuous drip rate.This indicates the importance of water storage within the epikarst, which is very possibly within fine fractures, the abundance of which is apparent from visual observation of a dense network of soda straw stalactites in the cave.Both stalagmites also contained continuous annual fluorescent laminae, inferring continuous hydrological connectivity to the overlying blanket bog (Proctor et al., 2000).Moreover, for discrete annual fluorescent laminae to be preserved, this suggests the seasonally distinct movement of infiltrating waters from the overlying peat, and hence water flow to the epikarst is likely to be maintained by the substantial storage capacity of the surficial peat.The combinations of these lines of evidence suggest that the stalagmite is predominantly fed by a small epikarst water store, probably associated with dissolutionally enhanced micro-fractures.This store receives inflows from the overlying peat, which transport fluorescent organic matter produced by microbial degradation in the peatland during periods of summer soil moisture deficit.
Given the detailed information available for UAT, the karst hydrology model described by Bradley et al. (2010) was modified to derive the model structure summarised in Fig. 4 and Table 1.The hydrological system is envisaged to comprise two principal water stores: the epikarst, and a surface peatland, with water storage in both stores determined by a water balance equation.The sole water inflow to the peatland is precipitation, with outflows of evapotranspiration (ET), and two flow pathways (F1 and F2) representing infiltration to the karst.F1 describes a flush from the peat store which occurs every autumn, two months after the observation of seasonal peak evapotranspiration (i.e.PET t < PET t−1 < PET t−2 > PET t−3 ).F1 occurs annually, but only during one autumn month at a time when a significant proportion of monthly rainfall is not lost via evapotranspiration, and autumn rainfall is flushed through to the cave system, prior to the peatland water table rising to the surface.F2 occurs whenever there have been two consecutive months with surface temperatures below freezing, permitting a snowmelt flush, as observed in stalagmite SU-96-7 by the presence of occasional fluorescent double laminae.The threshold for F2 was set to 0.4 • C, which generated the same number of flushes over the 1000-yr simulations as observed as lamina doublets in stalagmite SU-96-7 (n = 43).For both flows (F1 and F2), a volume of water equivalent to 10 % of the peat water store flows from the peat store to the epikarst store.
F1 and F2 are the sole water inflows to the epikarst while water outflows (F3) from the store every month, at a rate of ∼ 1 % of total water storage within the store.In our model, the principal hydrological control on stalagmite drip water δ 18 O is the relative size of the epikarst store (envisaged as a dissolutionally enlarged fracture or proto-cave; EPXS-TOR) in relation to the proportion of water draining this store (through micro-fractures to other stalagmites or stores; F3).Where the storage volume is large in relation to the drainage rate, δ 18 O variability will be low and drip water supply will be continuous, whereas a high drainage rate relative to storage volume would result in high δ 18 O variability and discontinuous water supply.Knowing that stalagmites SU-96-7 and SU-03-2 both contain a continuous annual lamina sequence and variable δ 18 O, we can constrain the storage volume to drainage ratio.Using data from a general circulation model (GCM) as the input series (see Sect. 4.2), we can demonstrate that this ratio has to be within the range 0.005 to 0.2.At a ratio of > 0.2, the epikarst store occasionally drains, which would result in discontinuous growth and hiatuses, which are not observed in either SU-96-7 or SU-03-2.At a ratio of < 0.005, all short-term variability in δ 18 O is removed and only a long-term trend is visible, which is also not apparent in stalagmite SU-96-7 or SU-03-2.Using these constraints, we take a median value for F3 of 0.08 × EPXSTOR.The δ 18 O composition of each store is determined by a mixing equation as a function of precipitation δ 18 O and the isotopic composition of the store at the previous time-step.The isotopic composition of water in the peat store is not adjusted to account for possible fractionation of peat water due to evapotranspiration since evaporation at this site is low in proportion to the total wetland water storage.Seepage waters from the epikarst are assumed to provide drip waters to stalagmites at UAT, and our model generates continuous drainage from the epikarst store (F3) as well as an overflow from the same store (F4).The latter is not considered further here, as an overflow routing would likely lead to discontinuous water supply and therefore stalagmite deposition with hiatuses, which is not observed in stalagmites SU-96-7 or SU-03-2.Stalagmite drip water δ 18 O is output as (1) the equivalent to the EPXSTOR composition, which can be conceptualised as hydrological routing via basal flow from the peat and the epikarst store; (2) mixtures of epikarst and precipitation δ 18 O, the equivalent of including a preferential flow component through the shallow peat such as from peat cracks or peat erosion (both of which are observed in the overlying soil) and fractures within the limestone.Stalagmite δ 18 O is estimated as a function of epikarst δ 18 O and calcite fractionation (mean temperature over the previous 7months).The most appropriate formulation of this fractionation equation for speleothem deposition is the subject of research debate (see, for example, Tremaine et al., 2011 andFairchild andBaker, 2012).Here we use the widely used equation of Kim and O'Neil (1997).
In its initial form, the model was developed by Bradley et al. (2010) using ModelMaker (Walker and Crout, 1997) but, as part of the changes in the model structure described here, the model (here termed KarstFor) has been coded in Fortran to reduce computation time and increase availability to the academic community.The model code is available in the Supplement.

Model validation
The model was validated by comparing measured stalagmite δ 18 O data from SU-03-2 and SU-96-7 with simulated series generated by KarstFor for the period 1960 to 2005 AD.Model input data comprised climate data (temperature, precipitation, potential evapotranspiration) taken from Proctor et al. (2000), and rainfall δ 18 O from the IAEA station at Valentia, Ireland, 800 km to the SW of the site.The IAEA rainfall data series was incomplete, and for months where δ 18 O is missing (predominantly between 1965 and 1977 AD), data were infilled by calculating a δ 18 O value from the relationship between δ 18 O and precipitation amount for the whole data series.The infilled Valentia series was then adjusted to provide a local estimate of precipitation δ 18 O using the correlation described by Fuller et al. (2008) from an analysis of 18 months of data.Given the potential uncertainty in the precipitation δ 18 O input series using distant and incomplete IAEA data, global precipitation isotope fields were obtained from IsoGSM, which is a historical isotope simulation computed using a GCM nudged with atmospheric reanalysis data (Yoshimura et al., 2008).The GCM is an up-to-date version of the Scripps ECPC global spectral model, with a horizontal resolution of T62 and 17 vertical levels.By applying spectral nudging towards the meteorology captured by the NCEP/DOE reanalyses, the spatial and temporal patterns of precipitation isotopes simulated by IsoGSM closely agree with observational data (the Global Network of Isotopes in Precipitation (GNIP) database).Thus, IsoGSM can provide realistic isotope information from regions where insitu observations are not available.Precipitation isotope data were extracted from the IsoGSM grids for Assynt (58 • N, 355 • E) and Valentia (51.93 • N, 349.75 • E), using bilinear interpolation.The data spanned the period 1979-2007 AD.A comparison between the IsoGSM time series at the nearest IAEA monitoring station (Valentia, Ireland) and NW Scotland, with our monthly precipitation series, is shown in Fig. 5. Figure 5a shows the good agreement between IsoGSM and IAEA GNIP data for Valentia, and Fig. 5b shows a similarly strong agreement between IsoGSM δ 18 O and our limited data for UAT (Fuller et al., 2008).This demonstrates that the use of IsoGSM data for NW Scotland is likely to be most representative of actual rainfall δ 18 O, rather than using the relatively distant Valentia series, as demonstrated by the contrasting gradients in the δ 18 O Valentia vs UAT using IsoGSM and adjusted IAEA GNIP data (Fig. 5c).
Four stalagmite δ 18 O series were generated (as shown in Fig. 4): (i) stal 1 was derived directly from the epikarst store δ 18 O; (ii) stal 2 was from an "overflow store" that received inflow from the epikarst when the latter > 2500 mm; (iii) stal 3 comprised 50 % epikarst δ 18 O, 25 % precipitation δ 18 O, 25 % precipitation δ 18 O in previous month; and (iv) stal 4 was 75 % epikarst 25 % precipitation δ 18 O.Calcite fractionation was described in each series as a function of mean monthly temperature over the previous seven months. Figure 6 presents a 12-month mean output for stal 1, stal 2, stal 3, and stal 4 together with stalagmite SU-96-7 and SU-03-2, although, as discussed earlier, overflow routing is not expected to be relevant to SU-96-7 and SU-03-2.Figure 6a uses IsoGSM as the precipitation isotope input series, and Fig. 6b uses the adjusted IAEA precipitation series.Comparison of model output between the two input series shows a good agreement, except over the period 1977 to 1981 AD when the GNIP-based output modelled series was ∼ 1 ‰ higher than the IsoGSM-based series.This occurs just after a heavily interpolated part of the GNIP-based series: the IsoGSM series also has a closer fit to the two stalagmite series over this time period.Combined with the good agreement with our limited field data (Fig. 5b), we therefore consider the IsoGSM-based model to be more reliable and only consider this input series for model validation.

Millennial-length input data series
Given the ability to model stalagmite δ 18 O for 10 3 -yr periods, it was then necessary to develop input climate data over longer time periods to enable 1000-yr forward model simulations.a.We generated multiple random 1000-yr monthly series for precipitation, temperature and δ 18 O of precipitation (Fig. 5).The precipitation and temperature series were both constrained to have a mean and normally distributed variance identical to that of the modern instrumental data series of Proctor et al. (2000).Evapotranspiration was estimated from temperature by the Thornthwaite method (Thornthwaite, 1948), while precipitation and temperature series were generated independently to avoid introduction of autocorrelation into the series.Global precipitation isotope fields were obtained from IsoGSM, described previously, and random series were generated with the same mean and normally-distributed variance.This series is hereafter called RANDOM.
b.We generated multiple adjusted random series by weighting both temperature and precipitation by the best available paleoclimate data.Temperature was adjusted using the Northern Hemisphere mean annual temperature of Mann et al. (2008), which was applied equally to all months of the year of the temperature series.Figure 7a demonstrates that, as expected, this results in input series that typically have a mean temperature that is up to 0.5 • C cooler than RAN-DOM for much of the last 1000 yr, with the increase in temperature over the last 100 yr particularly apparent.The adjusted series has the same variance as the RANDOM series, but now contains the autocorrelation structure of the Mann et al. (2008) series.This adjusted temperature series was then used to recalculate potential evapotranspiration.An adjusted precipitation series was derived using the wNAO reconstruction of Trouet et al. (2009), with winter (December-March) precipitation increased by 2.27 % for every unit change in wNAO.The Trouet et al. (2009) reconstruction only covers the period from ∼ 1050 AD, and consequently no adjustment was possible for the first 48 yr of the reconstruction and so the original randomised series was used for this period.Figure 7b demonstrates the small effect on mean monthly precipitation that this adjustment has.For δ 18 O of precipitation, the δ 18 O series was again generated using the mean and standard deviation of δ 18 O time series obtained from the nearest grid square using IsoGSM model output.Evapotranspiration was again estimated from temperature using the Thornthwaite method.These series are called MANN-TROUET.Given our relatively poor understanding of how precipitation δ 18 O relates to climate dynamics in the region over the last 1000 yr (Trouet et al., 2012), we do not adjust the δ 18 O input data for changes in temperature and precipitation, instead generating δ 18 O series with the same mean and normallydistributed variance as used for the RANDOM series.
c.A three-member ensemble of transient simulations of the last 2000 yr was conducted using the CSIRO Mk3L climate system model version 1.2 (Phipps et al., 2011).This is a coupled atmosphere-land-sea ice-ocean GCM, with horizontal resolutions of 5.6 × 3.2 • in the atmosphere and 2.8 × 1.6 • in the ocean.The model was forced with changes in the Earth's orbital parameters, atmospheric greenhouse gases (MacFarling Meure et al., 2006), solar irradiance (Steinhilber et al., 2009) and volcanic emissions (Gao et al., 2008).Monthly temperature and precipitation data were extracted for the model grid point closest to the study site.Temperature was adjusted by −0.7 • C to reflect the difference between the grid square mean and the local instrumental temperature series, and precipitation adjusted by ×2.12 to reflect local instrumental precipitation (both using the instrumental series of Proctor et al., 2000) in a region where there is a large precipitation gradient due to orographic effects.Evapotranspiration was again estimated from temperature using the Thornthwaite method.Figure 7a demonstrates that the GCM series has greater temperature variability than the RAN-DOM and MANN-TROUET series and a similar autocorrelation to the MANN-TROUET series.However, while precipitation variability is similar to the RAN-DOM and MANN-TROUET series, there is now autocorrelation in the series, as indicated in Fig. 7b.The precipitation δ 18 O series was again generated using the mean and normally distributed variance of δ 18 O time series obtained from the IsoGSM model output, without adjustment for temperature or precipitation.These three series are called GCM.
Model output comprises monthly water volume and isotopic composition in the peat and epikarst stores and stalagmite δ 18 O.The latter was calculated from modelled drip water δ 18 O following Kim and O'Neil (1997), with temperature-dependent fractionation applied to match the observed seasonal variability of cave air temperatures within the cave.

Forward model output
Pseudoproxy time series are presented in Figs. 8 and 9. Figure 8 presents rainfall input and model output from a 10-yr period, where model output comprises that using 3 GCM series and 10 MANN TROUET and 10 RANDOM series.Months 9000 to 9120 represent a 10-yr period from 1750-1760 AD, and model output is shown for the water isotopic composition of the epikarst store, and the stalagmite isotope composition for stal 1, stal 3 and stal 4 assuming a 2-yr sampling interval (representative of stalagmite SU-96-7).Figure 8 demonstrates that the epikarst store removes the monthly variation in δ 18 O as observed in the input rainfall series (Fig. 5b), with annual recharge of the epikarst store.Isotopic variability of the epikarst store is < 0.5 ‰ over the 10-yr period, with no difference between the MANN TROUET, RANDOM or GCM input series.This variability mostly derives from the variability between model simulations rather than climate variability within an individual simulation.Modelled stalagmite isotopic composition demonstrates increasing variability with increasing proportions of preferential flow, as would be expected.Stal 1 exhibits similar variability as the epikarst store, with variability increasing to ∼ 1 ‰ for stal 4 and 1.3 ‰ for stal 3.
Figure 9 presents model output for the months 0 to 12 000, covering the last 1000 yr, as well as data from stalagmites SU-96-7 and SU-03-2.Stalagmite SU-96-7 and SU-03-2 δ 18 O have both been sampled at 100-µm resolution equating to time intervals of between ∼ 1 and ∼ 10 yr depending on growth rate.Consequently, Fig. 9 compares the pseudoproxy series smoothed to 12 and 24 month moving averages to reflect the different stalagmite sampling resolution.Figure 9a compares stalagmite series with model output for stal 1 (derived from drip water from the epikarst store) and smoothed with a 24-month moving average, which model validation suggested was most appropriate for stalagmite   9b provides a comparison with model output stal 4 (25 % preferential flow) and 12-month smoothing, which validation suggested was more appropriate for stalagmite SU-03-2. Figure 9 shows that, with the exception of the first ∼400 months of the simulation run, when the KarstFor model was stabilising, pseudoproxy δ 18 O was relatively constant and with a variability of < 0.5 ‰ for stal 1 and < 1.0 ‰ for stal 4. No Figure 9b, where drip water series stal 4 comprises 75 % routed through the peat and epikarst stores and 25 % preferential flow, shows that the stalagmite SU-03-2 δ 18 O falls within the range of pseudoproxy δ 18 O over the ∼ 350 yr of stalagmite growth.However, if the drip water is sourced only from the peat and epikarst stores with no preferential flow, as inferred for stalagmite SU-96-7, pseudoproxy δ 18 O variability is significantly less than that observed.In particular, isotopically heavy δ 18 O occurred persistently between 1620 and 1880 AD in the SU-96-7 stalagmite series, falling well outside the range of all the pseudoproxy series.For this period, stalagmite δ 18 O is between 0.5 and 1.0 ‰ higher than the highest stalagmite and pseudoproxy δ 18 O in the last 100 yr and during the validation period.The early part of the stalagmite δ 18 O shift to higher values is coincident with a period of severe winters recorded as lamina doublets (Baker et al., 2002) as well as the change in stalagmite growth rate properties (Mariethoz et al., 2012).
Both stalagmites SU-96-7 and SU-03-2 contain lowfrequency variability within their δ 18 O series: this is not reproduced by the pseudoproxy series due to the independent nature of δ 18 O, temperature and precipitation in the input climate series.Although the GCM series contains some low frequency variability in temperature and precipitation, this is not large enough to imprint on the low-frequency component of pseudoproxy δ 18 O through its temperature dependent fractionation during calcite formation.The presence of low-frequency variability within the stalagmite δ 18 O series therefore remains enigmatic.It is probable that this has a climatic origin, either directly (e.g.due to variations in air mass source and trajectory) or indirectly (e.g.changes in fractionation processes that are not captured by the KarstFor model).With better dynamical understanding of the relationship between δ 18 O and temperature and precipitation in the region, or the development of 1000-yr isotope-enabled δ 18 O GCM simulations, comparison of stalagmite and pseudoproxy series in the spectral domain will be possible.Our forward modelling approach helps understand the potential sources of uncertainty in stalagmite δ 18 O at UAT, but a definitive interpretation of the stalagmite series remains elusive.One possible hypothesis is that there is no change in flow routing over the last 1000 yr, such that flow routing never involves preferential flow, and therefore precipitation δ 18 O over the period 1620-1880 AD would have to be isotopically heavier than the modern period.Such a hypothesis would be supported by the replication of higher δ 18 O values in both stalagmites, which is not clearly the case.An isotopically heavier rainfall composition would be consistent with a pervasive negative NAO state within the Little Ice Age (Trouet et al., 2012;Graham et al., 2011).However, it would also have to be outside the range of that modelled by the IsoGSM modelling over the modern period.A > 1 ‰ shift in mean precipitation δ 18 O would be required, which would be of a similar order of magnitude to that observed and modelled in the 8.2 ka event (LaGrande and Schmidt, 2009;Dominguez-Villar et al., 2009).Further independent proxy evidence, such as replicate speleothem series from a different cave in the region, or other independent proxies of δ 18 O in the region, as well as comparison with isotope enabled GCM output for the last 1000 yr, would be required to test this hypothesis.
An alternative hypothesis is that the preferential flow model is more realistic (Fig. 9b), in which case the SU-96-7 δ 18 O series can be explained without a change in precipitation δ 18 O over time.However, this would require systematic change(s) in hydrological routing to explain the change in stalagmite δ 18 O to heavier isotopic composition.This could be due to changes in the nature of peat accumulation over time.Pollen evidence indeed demonstrates changes between Calluna-and Sphagnum-dominated peat over time at the site (Charman et al., 2001); and changes in hydrological routing could also be achieved by increased preferential flow through cracking of the peat.This might lead to rapid infiltration of recent precipitation to water stores in the epikarst store and/or in karst fractures, and could permit recharge of isotopically heavier rainfall in summer, or in low intensity precipitation events.Such waters might otherwise be stored in the peat or lost via evapotranspiration.Alternatively, cracking or erosion of the peat might have increased the connectivity between surface water pools on the peat surface.Repeat analysis of waters sampled from four nearby pools on the peat surface between December 2003 and October 2004 revealed δ 18 O average compositions that ranged between −6.28 and −6.70 ‰ (n = 20), with the highest value measured of −3.0 ‰ (Fuller, 2007).Surface peat pool waters are all isotopically heavier than the mean cave drip waters recorded over the same period (−7.14 ± 0.14 ‰, n = 86), demonstrating the influence of evaporation on the surface water body.The heavy isotopic signal that occurs over the period 1620-1880 AD could therefore be a non-linear response to severe winters, with the cracking of the peat due to freezing, or freeze-thaw action on the shallow epikarst, enabling a periodic flush of isotopically enriched water.One would also expect a change in hydrogeochemical response at this time, and the period of more steady inter-annual growth rates, identified by non-flickering in the growth of SU-96-7 at ∼ 1615 AD (Mariethoz et al., 2012), can be interpreted as a change in hydrological routing to the stalagmite providing further evidence that this mechanism is a likely explanation for heavier stalagmite δ 18 O at this time.It could also explain the lack of correlation between δ 18 O and other climate series over the last 300 yr as attempted by Baker et al. (2011).An independent proxy for precipitation δ 18 O, such as that extractable from peat sphagnum (Daley et al., 2009), could help test our alternative hypothesis.
A final alternative hypothesis is that the periods of heavier stalagmite δ 18 O, which fall outside the range of the pseudoproxy series, are the result of non-equilibrium deposition at these times.Baker et al. (2011) demonstrate a kinetic fractionation effect in δ 13 C from cave monitoring, but modern cave drip water δ 18 O values are identical to those expected from equilibrium deposition.Changes in climatic conditions over the last 1000 yr could lead to increased isotope fractionation during calcite deposition.For example, SU-96-7 has a slow drip rate and SU-03-02 a high drip rate variability: therefore both might be sensitive to degassing and equilibriation fractionation effects, but to different extents.However, although disequilibrium can never be completely ruled out, one would expect this to have resulted from significant changes in the saturation state of the waters, or humidity of the cave, which should be reflected in the crystallization style or other aspects of speleothem chemistry.In the absence of such evidence, we believe that a climate forcing is the most probably explanation.et al. (2011) demonstrated that δ 18 O for European Holocene speleothems, using multiple samples binned into 50-yr resolution slices, systematically decreases as a function of distance from Atlantic Ocean, which can be explained by a Rayleigh distillation model of progressive rainout of easterly advected moisture-bearing air masses.They concluded that speleothem calcite δ 18 O values were not unduly affected by site-specific or variable disequilibrium effects, at least on multi-decadal timescales.The UAT site forms part of that compilation, and here we investigate the δ 18 O of two stalagmites at this site, sampled at 1-10 yr resolution.At this temporal sampling resolution, both stalagmites demonstrate only weak correlations with instrumental and recent climate parameters, and between-sample replication is poor.At such temporal resolution, evidently siteand sample-specific factors become increasingly important.To better understand the interplay between climatic and local hydrological effects of stalagmite δ 18 O at UAT, we have presented for the first time a lumped parameter model for the generation of 1000-yr hydrological simulations of groundwater storage and isotopic composition relevant to the generation of stalagmite δ 18 O pseudoproxy series.Using the Karst-For model, we have generated the first 1000 yr pseudoproxy simulations of stalagmite δ 18 O for the UAT site.A good agreement between pseudoproxy mean δ 18 O and stalagmite mean δ 18 O over the modern validation period demonstrates that the calcite was precipitated close to isotopic equilibrium.We demonstrate the importance of preferential flow in determining δ 18 O variability; however, with model optimisation not possible, we are not able to fit an absolute proportion of epikarst and preferential flow for each stalagmite δ 18 O series.However, our forward modelling approach allows the quantification of stalagmite δ 18 O variability that is related to flow routing.Comparison of pseudoproxy and stalagmite SU-96-7 δ 18 O series demonstrates that stalagmite δ 18 O over the period 1620-1880 AD falls outside that of the model validation period and of the pseudoproxy series.Either recharge δ 18 O was isotopically heavier at this time than the modern period, or this difference reflects changes in karst and peat hydrology over time, possibly due to peat cracking or erosion.Useful further research could include the use of isotope-enable GCM output for the last 1000 yr to better constrain pseudoproxy δ 18 O variability; additional stalagmite records for the last 1000 yr to better constrain climate versus karst and peat hydrological variability in δ 18 O; analysis of δ 18 O proxies from overlying peat sphagnum, to obtain an independent record of precipitation δ 18 O, and the development of integrated hydrology and isotope fractionation models.

Fig. 1 .
Fig. 1.UAT site description and 2003-2005 AD monitoring data summary.From top to base: site location; monthly precipitation, monthly soil moisture excess/deficit; UAT drip water δ 18 O; comparison of mean drip water δ 18 O and monthly precipitation δ 18 O samples.From Fuller et al. (2008).

Fig. 4 .
Fig. 4. Conceptualisation of the KarstFor model (left) with only baseflow from a homogenous peat store (right) with preferential flow through the peat and epikarst stores via peat cracks and bedrock fractures.

A.
Baker et al.: Millennial-length forward models and pseudoproxies of stalagmite δ 18 O