Accelerated phosphorous leaching during abrupt climate transitions in a temperate Atlantic ecosystem in Northwest Spain recorded by stalagmite P ∕ Ca variations

In natural ecosystems, phosphorus cycling regulates terrestrial productivity and may respond to climate variations. Seasonal to several year monitoring studies capture the short-term controls on P release but may miss longer term feedbacks. There is an important observational gap of the centennial to millennial scale response of the P cycle to climate oscillations. Cave carbonates such as stalagmites and flowstones, which precipitate from infiltrating groundwater, may record past changes in P loss on these timescales. Here, we examine trends in $\mathrm{P}/\mathrm{Ca}$ ratios in four coeval stalagmites from coastal caves in NW Iberia during two climate transitions, the Penultimate Glacial Maximum through the Last Interglacial (145 to 118 kyr BP) and an intermediate glacial climate state interrupted by an abrupt cooling event of the Greenland Stadial 22 (92 to 80 kyr BP). We conduct sensitivity tests with a model to assess the degree to which drip water pH and in-cave drip water chemical evolution could affect the stalagmite $\mathrm{P}/\mathrm{Ca}$ record. Both during the last deglaciation and during Greenland Stadial 22, we find large (3–10-fold) transient increases in stalagmite $\mathrm{P}/\mathrm{Ca}$ at the onset of abrupt cooling events and during the rapid recovery from some events. These increases are much larger than can be explained by variations in P incorporation due to drip water pH or in-cave chemical evolution and likely reflect significantly increased drip water $\mathrm{P}/\mathrm{Ca}$ ratios at the onset and end of abrupt stadial events. Two climatic factors may contribute to this increased leaching. First, soil temperatures may reach minimum values during these transition states, if the temperature minimum leads to increased thickness and duration of snow cover which raises soil temperatures. Minimum winter soil temperature suppresses microbial recycling of P. Second, the transitions into and out of stadial events may feature the highest frequency of freeze-thaw events which change the physical soil structure and lead to stronger spring flushing of P. Additionally, during cooling, reduced soil respiration rates may raise the pH of soil water and lead to increased mobility of P bound to soil minerals.

1 Introduction

Phosphorus availability is an important regulator of production in terrestrial ecosystems (Vitousek et al., 2010; Achat et al., 2016; Elser et al., 2007; Elser et al., 2000). Bedrock represents the ultimate source of P for terrestrial ecosystems. Dissolved inorganic orthophosphate is the primary form of P cycled through ecosystems because it is the main form assimilated by organisms (Blake et al., 2005). As P is released from bedrock minerals, it may be taken up in living biomass, stored in litter, or stored in soil adsorbed to Fe and Al (hydr)oxides in colloids or precipitated with cations (Hinsinger, 2001; Achat et al., 2016). The efficient capture of P in biomass means that organic litter is a major P reservoir and in a midlatitude forest ecosystem, water in the organic layer of the soil contains high P concentrations which are released from litter (Sohrt et al., 2019). These concentrations are reduced by an order of magnitude in the deeper mineral layers of the soil and further reduced in groundwater due to immobilization in soil oxides (Sohrt et al., 2019).

Both the weathering rate releasing P from bedrock (Walker and Syers, 1976; Filippelli and Souch, 1999), and the processes driving the loss of nutrients from the ecosystem (Buendía et al., 2010) can be climate sensitive, but also mediated by biological function (Oeser and Von Blanckenburg, 2020; Pastore et al., 2022). Current monitoring networks capture the short-term dynamics of P cycling but may not capture the decadal to centennial scale processes regulating P cycling. On million year timescales, potential changes in P cycling have been investigated using ocean sediments as the ultimate sink for terrestrial P (Filippelli, 1997; Dodd et al., 2021) and models of P cycling (Buendía et al., 2010). Yet, at intermediate timescales, in karst regions, cave carbonate formations such as stalagmites and flowstones, may also record changes in P cycling from the decadal to hundred-thousand-year timescales. Air filled caves intercept groundwater percolating through the non-saturated (vadose) zone, and the stalagmites and flowstones forming from this carbonate potentially preserve the history of P outflows in deep groundwater. Karst regions comprise 20 % of earths non ice-covered land surface (Goldscheider et al., 2020). While globally carbonate rocks have lower P than shale and most crystalline rocks (Porder and Ramachandran, 2013), the phosphate hosting phases in carbonates are similar to other bedrock classes in which phosphate is present in calcium phosphates as well as P bound to Fe and Mn-oxides (Dodd et al., 2021). Soil pH strongly influences the solubility and mobility of inorganic phosphorus by controlling both mineral dissolution and adsorption–desorption processes (Hinsinger, 2001). In acidic soils, phosphate is commonly bound to Fe- and Al-(oxyhydr)oxides, with adsorption maximized at low pH due to positively charged mineral surfaces; as pH increases, surface charge decreases, weakening adsorption and enhancing phosphate release. In alkaline or calcareous soils, phosphorus is often present as low-solubility calcium phosphate minerals such as apatite, whose dissolution is promoted as pH decreases from around 8 toward acidic values (Hinsinger, 2001). Microbial communities can accelerate these processes through localized acidification and the production of organic ligands that solubilize mineral-bound P, even under relatively stable bulk pH conditions (Pastore et al., 2022). Together, these mechanisms suggest that long-term or microbially mediated shifts in soil pH could mobilize mineral-derived P, providing an additional source of P to drip waters beyond organic matter contributions.

Although caves provide a ready access for sampling infiltrating groundwaters, few determinations of P concentration in cave drip waters have been reported. Drip water data from Trentino Italy reported in Fairchild et al. (2001) indicate typical total P contents of 3 to 8 µg L⁻¹, with increases to 30–50 µg L⁻¹ during infiltration events. Concentrations <10 µg L⁻¹ were found in drip water from a cave with mixed forest and pasture coverage (Kost et al., 2023). These mean concentrations are consistent with the 7 µg L⁻¹ mean groundwaters in temperate forest in Southern Germany under crystalline gneiss (Sohrt et al., 2019) a rock type with similar global mean P as carbonates (Porder and Ramachandran, 2013). In a tropical atoll with bedrock featuring high P concentrations, drip water P mean value was 280 µg L⁻¹ (Frisia et al., 2012). In a cave with agricultural land use above the cave, a seasonal drip had mean P concentrations between 5 and 40 µg L⁻¹, correlated with NO${}_{\mathrm{3}}^{-}$ and K⁺ concentrations applied in fertilizer (Jiménez-Sánchez et al., 2008).

Stalagmites may be expected to record variations in drip water P concentration because experimental precipitation of laboratory carbonates has shown that the P uptake into carbonates is proportional to its concentration in the water (Dodd et al., 2021). For calcite, at pH 8 the effective distribution coefficient relating $\mathrm{P}/\mathrm{Ca}$_solid to $\mathrm{P}/\mathrm{Ca}$_solution averaged 0.47, indicating P is slightly excluded from incorporation. However, there is evidence that this distribution coefficient is not constant. A 20 °C temperature increase leads to 50 % higher P incorporation (Dodd et al., 2021). Phosphate incorporation in carbonate is negatively correlated with solution pH, dropping by 3.5 fold as pH rises 6.5 to 8.5 during aragonite precipitation (Dodd et al., 2021); a similar decline in P incorporation with increasing pH is observed in calcite(Ishikawa and Ichikuni, 1981). For phosphate ion incorporation in speleothem calcite, a competitive isovalent (HPO${}_{\mathrm{4}}^{-\mathrm{2}}$) substitution for the carbonate ion has been proposed, which would be proportional to the solutional ratio of HPO${}_{\mathrm{4}}^{-\mathrm{2}}$ to CO${}_{\mathrm{3}}^{-\mathrm{2}}$. Because the second pK for carbonic acid is lower than the second pK for orthophosphoric acid, as pH increases from 7.3 to 8.3, there is a 2-fold decrease in the ratio HPO${}_{\mathrm{4}}^{-\mathrm{2}}$ to CO${}_{\mathrm{3}}^{-\mathrm{2}}$, which may lead to decreasing P incorporation at higher pH (Stoll et al., 2012). NMR studies of naturally formed stalagmites from two cave systems confirm that there is ubiquitous incorporation of phosphate groups in the crystal lattice during crystal growth (Frisia et al., 2012; Mason et al., 2007). In some circumstances, P may additionally enter stalagmites by occlusion of detrital particles containing P, or via formation of accessory phases such as monetite or ardealite which was found along two corroded “stromatolite-like” dissolution layers with microbially driven carbonate formation during a growth hiatus (Frisia et al., 2012).

The interpretation of variations in P concentrations in stalagmites has been diverse. On the decadal scale, comparison of P concentration with instrumental climate records from Cook Islands showed that for a stalagmite in which hydrological routing facilitates transmission of colloids, P concentrations are high and positively correlated with infiltration, but in another stalagmite from the same cave P concentrations were low and uncorrelated with infiltration, underscoring the influence of hydrological routing in this setting (Faraji et al., 2024). In Italian Bigonda Cave, over the early phases of deglaciation, lower flowstone P concentrations were attributed to periods of pioneer forest expansion increasing P sequestration, with higher P during the interglacial attributed to mature forest when P released by vegetation dieback was comparable to P sequestered by growth (Johnston et al., 2021). Alternatively, over deglacial transitions in a flowstone from Tana Che Urla cave in Italy, it was also proposed that higher P concentrations might reflect periods of increased soil development and colloid production facilitating P mobility (Regattieri et al., 2016). In contrast in tropical settings in Cuba, lower P concentrations were interpreted as evidence of increased vegetation and lower soil pH stabilizing P storage on Fe and Al colloidal phases in soils (Warken et al., 2019).

In the face of these diverse interpretations, the replication of trends in stalagmite P could clarify when regional processes or specific hydrological routing are responsible for variations, improving interpretation, yet few studies (with the exception of Faraji et al., 2024) have assessed the replication. Furthermore, it is unclear to what degree changes in P partitioning as a function of solution pH, as observed in experimental calcite growth (Dodd et al., 2021), may complicate using stalagmite P to infer changes in P concentrations in groundwater and drip water over the centennial to multimillennial timescales.

Here, we report a large dataset of P determinations in multiple stalagmites spanning the transition from the Penultimate Glacial Maximum (PGM) to the Last Interglacial (LIG) and spanning a stadial cooling event (Greenland stadial 22, GS22) around 87 kyr ago. By compiling records from multiple stalagmites in the same region, and conducting model calculations of P incorporation, we aim to identify which variations in speleothem P can be most confidently attributed to changes in P cycling, and those which may reflect in-cave processes and features specific to hydrological routing. The selected time periods are expected to span significant regional climate and ecosystem change. Pollen records show major transitions in the regional terrestrial ecosystems from glacial conditions dominated by shrubs and grasses and less than 20 % arboreal pollen, to interglacial conditions with up to 80 % arboreal pollen (Tzedakis et al., 2018). There is also a reduction in arboreal pollen from 40 % to <10 % at the onset of the GS22 event (Goni et al., 2008). Geochemical records also suggest major transitions in soil pCO₂ from minima of 800 ppm during stadial events to interglacial (Holocene) levels of 6000 ppm (Lechleitner et al., 2021). To assess the impact of changing P distribution coefficient in stalagmite records, we complete a set of sensitivity tests including both constant P distribution coefficient and a pH sensitive distribution coefficient, evaluating how they affect the water and stalagmite $\mathrm{P}/\mathrm{Ca}$ ratios during the degassing and calcite deposition process of speleothem formation. Our analyses reveal a reproducible trend of increased $\mathrm{P}/\mathrm{Ca}$ ratios of drip water at the onset of stadial cooling periods, which exceeds the magnitude of change in $\mathrm{P}/\mathrm{Ca}$ which could be attributed only to changes in P incorporation in calcite. This result is suggestive of transient periods of increased P loss from the terrestrial ecosystem of this region.

Figure 1 Map showing the study area (a) and the detailed location of the two coastal caves (b) from where the stalagmites originate as well as the climatology from coastal station of Llanes from 1981 to 2000 showing blue bars indicating the average monthly precipitation (mm) and a red line showing the average monthly temperature (°C) (c). Maps adapted from Stoll et al. (2015).

2 Methods

2.1 Cave Settings and Stalagmites

We evaluate stalagmites from two cave systems (Fig. 1) in northern Spain. La Vallina Cave is hosted in Carboniferous limestone of the Barcaliente formation and Cueva Rosa is hosted in the Escalada Formation. Modern vegetation above the cave sites is variably affected by anthropogenic disturbance, such as planting Eucalyptus crops, and maintaining pasture (including occasional burning for pasture regeneration). The climax vegetation of this area is temperature deciduous forest (Roces Díaz et al., 2015). La Vallina cave (43°24^′36^′′ N, 4°48^′24^′′ W; 70 m a.s.l. (above sea level) is on a northward-facing hillslope at a distance of 2.5 km from the coast, with rock thickness above the cave ranging from 5–30 m. Modern vegetation includes native Quercus ilex and Quercus rober as well as recently grown eucalyptus and pasture (see Kost et al., 2023, for overview) Cueva Rosa (43°26^′37^′′ N, 5°08^′25^′′ W; 121 m a.s.l.) features 50 m of bedrock above the cave, and a lower gallery is transited by a perennial cave stream which periodically floods some of the upper galleries (González-Lemos et al., 2015).

The annual precipitation at the cave sites is 1250 mm, based on the averages from 1970 to 2009 (data on annual rainfall from http://idebos.bio.uniovi.es/GeoPortal/Atlas/Pan19702009.html, last access: 28 April 2013). A water deficit emerges in summer months due to lower precipitation and increased evapotranspiration (for a review, see Kost et al., 2023). Mean cave temperatures are 12 °C (Kost et al. 2023), close to the annual average air temperature, with monthly air temperatures at the coastal site ranging from 9 to 20 °C (Fig. 1c).

Figure 2 Comparison of the range of $\mathrm{P}/\mathrm{Ca}$ in stalagmites within the 145 to 118 ka time frame and the 92 to 80 ka time frame and in the bedrock from La Vallina (LV) and Cueva Rosa (CR). Blue indicates stalagmites and bedrock from La Vallina and red-orange indicates stalagmites from Cueva Rosa. Stalagmite data are illustrated if they fall below the $\mathrm{Al}/\mathrm{Ca}$ threshold of 0.04 mmol mol⁻¹. Line indicates the median, box indicate the 80th and 20th percentiles and whiskers indicate the 95th and 5th percentiles.

Drip water monitoring in La Vallina cave found modern P concentrations to be below the analytical detection limit of 0.32 µmol L⁻¹ (Kost et al.2023). The median bedrock $\mathrm{P}/\mathrm{Ca}$ (determined as described in Kost et al. (2023) at La Vallina is 0.19 mmol mol⁻¹ and at Calabrez it is 0.08 mmol mol⁻¹ (Fig. 2). Because these estimates derive from opportunistically collected bedrock samples, we do not have a way to estimate the weighted mean bedrock composition at either cave setting.

Table 1 Stalagmites examined in this study, sampling details, and the source of age model and $\mathrm{Mg}/\mathrm{Ca}$ results from previous studies: (1) Stoll et al. (2022), (2) Stoll et al. (2023), (3) Stoll et al. (2013), (4) Stoll et al. (2015). Cave LV indicates La Vallina and CR indicates Cueva Rosa. Bedrock BAR is Bar and Esc is Escalada Formation.

* indicates contiguous drill samples

Stalagmites GAE, GLO, and ROW span the 92 to 90 kyr interval and their carbon isotope record, corrected for in-cave degassing and prior calcite precipitation processes, record the cooling of the Greenland Stadial 22 (GS22) event as a positive shift in the δ¹³C_init (Stoll et al., 2023) (Table 1). Stalagmites GAR, NEI, GAE, GLD, and BEL span the transition from the Penultimate Glacial Maximum (PGM) to the Last interglacial (LIG); isotope and $\mathrm{Mg}/\mathrm{Ca}$ records from GAR, NEI and GAE and GLO have been described (Stoll et al., 2022, 2015) and the δ¹³C_init record from GAR discussed (Kaushal et al., 2025) and the constraints for calculation of δ¹³C_init for GAR, GAE, GLD and BEL presented elsewhere (Stoll et al., 2023). In several of these stalagmites, previously published isotope and $\mathrm{Mg}/\mathrm{Ca}$ records span longer time intervals than the $\mathrm{P}/\mathrm{Ca}$ determinations presented here, because some previously published $\mathrm{Mg}/\mathrm{Ca}$ were measured by ICP-OES for which we did not obtain reliable P measurements. As detailed elsewhere, trends in δ¹³C_init (and the index fCa, the fraction of calcium remaining in solution after prior calcite precipitation) may be calculated more reliably than the absolute values (Stoll et al., 2023), and here we focus on the interpretation of trends in these parameters. Sample chronology from U-Th dates derives from published studies (Stoll et al., 2013, 2022, 2023) (Table 1). Age models are plotted in Fig. S1. Median growth rates of the studied samples range from 4 to 36 µm yr⁻¹. Age models for GAR include sectors with layer counted growth rates, and in counted sectors layer thickness suggests growth rates varying by 5 to 10-fold on the decadal to centennial scale. In this sample, layer thickness decreases during cold stadial events, and annual layers become too thin for counting between 134.1 and 132.3 ka, as well as during the LIG (Stoll et al., 2022). Age models for other stalagmites are based on the interpolation schemes of BCHRON (Haslett and Parnell, 2008) and STALAGE (Scholz and Hoffmann, 2011) or linear interpolation when the density of dates is low. With the exception of GAR, the available chronological resolution does not permit us to confidently identify variations in growth rates between $\mathrm{U}/\mathrm{Th}$ dates, and in most samples growth rate variations cannot be identified at timescales shorter than several thousand years (typical interval between dates Fig. S1).

2.2 Analytical Methods

Stalagmite samples were drilled using a Sherline drill which digitally monitors sampling position. Five of the presented stalagmite datasets are based on drilling resolution that continually sampled all of the deposited stalagmite along the growth axis, using drilling resolution of 0.25 to 1 mm per sample (Table 1). These include most of the presented records covering the PGM-LIG and one of the four records covering the GS22. For these samples, each drilled increment typically integrates 30 to 230 years (Table 1). The other presented records sampled 1 mm of powder but at sample increments ranging from 2 to 5 mm, meaning that only 20 % to 50 % of the growth conditions were sampled. The typical interval between samples reflects 150 to 400 years for these datasets (Table 1). During GS22, each drilled sample from GAE, NEI and ROW integrates 30 to 70 years based on the median growth rates. For a 1.5 cm increment in GAE with higher sample resolution of 0.5 mm, this portion of the record has also been downsampled to reflect a constant sample spatial resolution. The impact of variable stalagmite growth rate on the variance in $\mathrm{P}/\mathrm{Ca}$ is assessed in stalagmite GAR, which features a layer-counted age model from 136 to 128 ka (Fig. S2).

An aliquot of 330 to 400 µg of powders was dissolved in 350 µL of 2 % HNO₃. $\mathrm{P}/\mathrm{Ca}$ and $\mathrm{Al}/\mathrm{Ca}$ were analyzed with an Agilent 8800 QQQ-ICP-MS at ETH Zürich using the intensity ratio calibration (de Villiers et al., 2002) and in-house multielement standards of varying trace element to Ca ratios but with Ca concentrations matched to samples (400 ppm). We analyze ³¹P mass shifted in O₂ reaction mode (+16 mass units) to reduce ¹⁵N¹⁶O polyatomic interferences. We report data as the stalagmite $\mathrm{P}/\mathrm{Ca}$ ratio. Because Ca is the major ion which is accounting for stalagmite growth, increases in $\mathrm{P}/\mathrm{Ca}$ cannot be only due to higher growth rate because a higher growth rate also means a greater accumulation of Ca per unit time.

2.3 Model of $\mathrm{P}/\mathrm{Ca}$ evolution in stalagmites

The incorporation of P in stalagmite calcite may vary due to changes in the pH of dripwater and the evolution of dripwater chemistry during calcite precipitation. Water infiltrating through soil and root zones is acidified by respired CO₂ and dissolves limestone bedrock; higher soil CO₂ leads to lower water pH and greater dissolution. When infiltrating water reaches the lower pCO₂ cave environment, CO₂ begins to degas from the drip water, leading to an increase in pH and the precipitation of CaCO₃. CaCO₃ precipitation can occur on the cave ceiling before the drip water falls to the top of the stalagmite, and this “prior calcite precipitation” affects the pH and concentration of ions in the water forming the stalagmite. This process is simulated by a number of models.

To evaluate the potential effects of prior calcite precipitation and changing drip water pH on P incorporation, we carry out a set of calculations accounting for the change of Ca, pH, dissolved inorganic carbonate (DIC), and dissolved P in drip water. We employ CAVECALC (Owen et al., 2018) to simulate the evolution of the DIC from bedrock dissolution through progressive degassing for a set of 5 initial soil pCO₂ concentrations which lead to a range in initial pH of drip water (Table S1). For each set of conditions, from Cavecalc, we export the sequence of evolution of Ca, DIC, HCO${}_{\mathrm{3}}^{-}$, CO${}_{\mathrm{3}}^{-\mathrm{2}}$, and pH from dissolution to complete degassing when the water is equilibrated with the cave atmosphere. We report the progressive Ca loss using the index fCa, where fCa is the fraction of initial Ca remaining in solution (see Stoll et al., 2023). From a stalagmite sample, one indicator of the soil pCO₂ can be attained from the stalagmite carbon isotope ratio corrected for in-cave fractionation processes (δ¹³C_init), because δ¹³C_init becomes more negative at higher soil pCO₂. We use a previous set of Cavecalc simulations (Lechleitner et al., 2021), to provide the δ¹³C_init value of calcite corresponding to open system dissolution at the given soil pCO₂.

Table 2 Summary of model scenarios for P incorporation in calcite.

We seek to evaluate the incorporation of orthophosphate in calcite and do not simulate the incorporation of P associated with detrital phases or other accessory minerals. For comparison with the $\mathrm{P}/\mathrm{Ca}$, we evaluate also the impact on $\mathrm{P}/\mathrm{Mg}$ since the drip water concentration of the latter, highly incompatible element is negligibly affected by calcite precipitation. We evaluate each of these 5 CAVECALC runs for four different possible models for the behavior of P incorporation in calcite (Table 2). For the calculations, we define the $\mathrm{P}/\mathrm{Ca}$ ratio of drip water at initial dissolution (hereafter initial dripwater $\mathrm{P}/\mathrm{Ca}$) to reflect the median $\mathrm{P}/\mathrm{Ca}$ of the host bedrock of La Vallina Cave which is 0.2 mmol mol⁻¹; this ratio is consistent with the modern drip water concentrations in La Vallina being below 0.32 µmol L⁻¹ and could reflect a steady state value where P inputs from rock dissolution are comparable to P loss through groundwaters. We set the initial drip water $\mathrm{Mg}/\mathrm{Ca}$ ratio at 2 mmol mol⁻¹ for all simulations and employ a DMg of 0.03. Following Dodd et al. (2021) we define P incorporation coefficient (D) relative to Ca, which matches the analytical measurement of $\mathrm{P}/\mathrm{Ca}$ ratios, although the actual P incorporation mechanism may reflect substitution for the carbonate ion rather than Ca⁺². We calculate the evolution of the drip water P concentration iteratively through the progressive calcite precipitation, evaluating various assumptions for the incorporation of P in calcite. In two models, we specify that the P incorporation is proportional to the solution ratio of HPO${}_{\mathrm{4}}^{-\mathrm{2}}$ to CO${}_{\mathrm{3}}^{-\mathrm{2}}$ and therefore evolves during degassing, starting with a maximum D of 1.46 in one scenario and a maximum D of 0.75 in another (Table 2) (D, the partition coefficient, is the ratio of P concentration in calcite to that in drip water). In two other models, we specify a constant P incorporation (D=1, and D=0.5). The partitioning proportional to the solution ratio of HPO${}_{\mathrm{4}}^{-\mathrm{2}}$ to CO${}_{\mathrm{3}}^{-\mathrm{2}}$ is consistent with a mechanism of P incorporation dominantly substituting for CO${}_{\mathrm{3}}^{-\mathrm{2}}$, and is consistent with the observed pH dependence of P incorporation in carbonates (Dodd et al., 2021). A constant D could be indicative of P incorporation as both HPO${}_{\mathrm{4}}^{-\mathrm{2}}$ and H₂PO${}_{\mathrm{4}}^{-\mathrm{1}}$, the latter substituting HCO${}_{\mathrm{3}}^{-}$ (potentially analogous to borate incorporation in calcite). The incorporation of HCO${}_{\mathrm{3}}^{-}$ as well as CO${}_{\mathrm{3}}^{-\mathrm{2}}$ in growing calcite remains discussed (Andersson et al., 2016; Huang et al., 2021) hence we include this possibility in our sensitivity analysis. There are limited constraints on the effective P incorporation coefficient in calcite. The experiments of (Dodd et al., 2021) for calcite suggest D=0.47 at pH 8. An additional constraint is provided by the maximum the drip water P concentration (the analytical detection limit of 0.32 µmol L⁻¹ and 61 ppm average Ca concentration (Kost et al., 2023) and the average $\mathrm{P}/\mathrm{Ca}$ of 0.13 mmol mol⁻¹ for the actively growing stalagmite SNO in La Vallina Cave for which the negligible correlation of $\mathrm{P}/\mathrm{Ca}$ and $\mathrm{Al}/\mathrm{Ca}$ suggests limited detrital P contribution (Sliwinski et al., 2023). For this single, rapidly growing stalagmite, we estimate a minimum D of 0.6 but we do not have an estimate for the pH during precipitation of this sample.

2.4 Comparison with coupled climate model

To evaluate the changes in snow cover and soil humidity under evolving boundary condition, we employ a previously published simulation (Romé et al., 2022). The simulation features a quasi-idealised glacial climate state with an oscillating AMOC strength, which is triggered by a constant meltwater flux of 0.084 Sv corresponding to the GLAC-1D (Ivanovic et al., 2016) ice sheet history at 17.8 ka BP. When comparing periods of a strong (∼16 Sv) to a weak (∼8 Sv) AMOC state in this simulation, the patterns significant for an AMOC slowdown emerge (i.e., temperature reduction in the North Atlantic realm and Southward shift of ITCZ), allowing us to study the range of plausible conditions above the study site.

3 Results

3.1 Range of stalagmite $\mathrm{P}/\mathrm{Ca}$

The median $\mathrm{P}/\mathrm{Ca}$ ratio ranges from 0.1 to 0.2 mmol mol⁻¹ across 9 datasets from 7 different stalagmites (GAE and NEI provided data during two distinct periods) (Fig. 2). The PGM-LIG section of GAE features numerous analyses with higher $\mathrm{P}/\mathrm{Ca}$ than other samples, despite very slow growth and higher signal smoothing than other coeval stalagmites such as BEL and GAR, and features a positive correlation with $\mathrm{Al}/\mathrm{Ca}$ ratios in the acid soluble fraction (Fig. S3). GAE in this time span also features intervals of anomalously high ²³²Th in the series of dated samples, indicative of intervals with appreciable detrital content (Stoll et al., 2022, 2013). These high P and high Al data points originate from the base of the stalagmite and from an interval of condensed and potentially interrupted growth in the last interglacial (121 to 124 ka BP). For examination of subsequent trends, we evaluate the $\mathrm{P}/\mathrm{Ca}$ ratios which have $\mathrm{Al}/\mathrm{Ca}$ ratios <0.04 mmol mol⁻¹, the interval below which $\mathrm{Al}/\mathrm{Ca}$ and $\mathrm{P}/\mathrm{Ca}$ exhibit no correlation in GAE. This filter eliminates few points in stalagmites other than GAE (Fig. S3). In recent actively growing stalagmites in Zoolithen cave, detrital contribution to stalagmite P was inferred from high correlations with Al (Riechelmann et al., 2020). This removes singular high P layers such as growth hiatus which may be enhanced by accumulation of (aluminum-bearing) silicates and the impact of microbially dominated carbonate phases following intervals of stalagmite dissolution (Frisia et al., 2012).

Figure 3 Range of $\mathrm{P}/\mathrm{Ca}$ in stalagmites from penultimate glacial (GL, 145 to 135 ka), deglacial (DE, 135 to 127 ka) and last interglacial (IG, 127 to 118 ka). Line shows the median, box shows 80th and 20th percentile and whiskers indicate the 5th and 95th percentile.

We divide the PGM-LIG datasets into three intervals, the last interglacial (127 to 118 kyr BP), the penultimate deglaciation (135 to 127 kyr BP) and the penultimate glacial (145–135 kyr BP). Among samples from the same interval, there are substantial differences in the mean $\mathrm{P}/\mathrm{Ca}$ ratios, with BEL and GAE exhibiting lower median $\mathrm{P}/\mathrm{Ca}$ than GAR or NEI (Fig. 3). These differences are not due to contrasting bedrock $\mathrm{P}/\mathrm{Ca}$, since the contrasting BEL and GAE are from samples from the same cave. Although the median $\mathrm{P}/\mathrm{Ca}$ in bedrock from Cueva Rosa is half that of La Vallina, the stalagmite from Cueva Rosa (NEI) does not feature lower mean $\mathrm{P}/\mathrm{Ca}$.

Figure 4 Time series of δ¹³C_init and $\mathrm{P}/\mathrm{Ca}$ for stalagmites GAE, NEI, GLO, and ROW spanning the 92 to 80 ka time interval of GS22. The $\mathrm{P}/\mathrm{Ca}$ symbols are color coded according to fCa as illustrated in the legend. For GAE, in the $\mathrm{P}/\mathrm{Ca}$ record three black circles illustrate downsampling to the same 5 mm resolution as the remaining record. Green vertical lines highlight the onset of increases in $\mathrm{P}/\mathrm{Ca}$.

3.2 Temporal evolution of stalagmite $\mathrm{P}/\mathrm{Ca}$

3.2.1 Greenland Stadial Event 22

In stalagmites GAE, GLO, and ROW, and NEI, a positive shift in δ¹³C_init marks the transient cooling event of GS22, as described in previous studies (Stoll et al., 2023). Prior to the cooling event, $\mathrm{P}/\mathrm{Ca}$ values are low in all the stalagmites; all four stalagmites feature an increase in the $\mathrm{P}/\mathrm{Ca}$ ratio synchronous with the shift in δ¹³C_init marking the onset of the cooling event (Fig. 4). $\mathrm{P}/\mathrm{Ca}$ in NEI increases nearly 5-fold, in GLO the peak is >5 fold higher than the baseline, $\mathrm{P}/\mathrm{Ca}$ nearly triples across the onset in GAE, and $\mathrm{P}/\mathrm{Ca}$ nearly doubles in ROW. The significantly lower amplitude of the δ¹³C_init anomaly in ROW (2 ‰ vs. 4 ‰ in the other sampled stalagmites) suggests that peak cooling may not have been sampled, potentially due to the interplay of slowed growth rate and lower geochemical sampling resolution used for this stalagmite. The low temporal resolution of sampling may also contribute to the lower sampled amplitude of the $\mathrm{P}/\mathrm{Ca}$ peak compared to the faster growing NEI or the continuously sampled GLO (Table 1). Following this increase at the onset of the cooling event, $\mathrm{P}/\mathrm{Ca}$ in GLO drops to values comparable to or below pre-event levels by midway through the cold event. $\mathrm{P}/\mathrm{Ca}$ in Rowena decreases slightly at the end of the cooling event, but without reaching the values prior to the event. In both NEI and GAE, a subsequent peak in $\mathrm{P}/\mathrm{Ca}$ occurs at about 82 to 84 ka BP.

Figure 5 Time series of δ¹³C_init and $\mathrm{P}/\mathrm{Ca}$ for stalagmites GAE, BEL, and GAR spanning the 145 to 118 ka time interval. The $\mathrm{P}/\mathrm{Ca}$ symbols are color coded according to fCa, as illustrated in the legend. Green vertical lines highlight the onset of increases in $\mathrm{P}/\mathrm{Ca}$. The sum of pollen from Mediterranean and Eurosiberian forest species on the Iberian margin is illustrated, from Tzedakis et al. (2018).

3.2.2 Penultimate Glacial Period to Last Interglacial

GAR provides the highest resolution and longest record from the penultimate glacial to last interglacial. In GAR, at the onset of a significant stadial cold event manifest in δ¹³C_init at 140 ka, $\mathrm{P}/\mathrm{Ca}$ abruptly increases and decreases at the end of the stadial event (Fig. 5). A series of transient increases occur between 138 and 130 kyr ago, and the onsets of several of them are coincident with the onset of centennial to millennial scale cooling events manifest as positive excursions in δ¹³C_init. Comparison of drilled resolution and average $\mathrm{P}/\mathrm{Ca}$ in 200 years age bins confirms that these trends are not artifacts of changing sample resolution (Fig. S2). The $\mathrm{P}/\mathrm{Ca}$ decreases at the culmination of the warming at 130 kyr, followed by another transient increase in $\mathrm{P}/\mathrm{Ca}$ which initiates during a subsequent set of stadial events around 129 ka. $\mathrm{P}/\mathrm{Ca}$ during the latter part of the last interglacial are highly variable.

BEL records a first abrupt increase in $\mathrm{P}/\mathrm{Ca}$ at a stadial event at ≈141 ka which is likely equivalent to the 140 ka event in GAR. Two extreme positive anomalies coincide with stadial cold events manifest as positive excursions in δ¹³C_init between 136 and 130 ka. No stadial cold events are recorded at 129 ka, potentially because stalagmite growth ends before the stadial events, an interpretation consistent with age model uncertainty. GAE features a pronounced increase in $\mathrm{P}/\mathrm{Ca}$ with the onset of stadial cold events between 136 and 130 ka, although the amplitude of change is lower in this slower growing stalagmite which averages twice as many years per sample as in BEL or GAE (Table 1). Albeit with lower resolution, GLD shows increased $\mathrm{P}/\mathrm{Ca}$ at the onset of a positive δ¹³C_init excursion, and decreasing $\mathrm{P}/\mathrm{Ca}$ with the deglacial warming (Fig. 6). In contrast, in NEI, with the onset of stadial cooling (manifest as positive excursions in δ¹³C_init), $\mathrm{P}/\mathrm{Ca}$ concentrations decrease significantly and recover during the warming (Fig. 6).

Figure 6 Time series of δ¹³C_init and $\mathrm{P}/\mathrm{Ca}$ for stalagmites GLD and NEI spanning the 145 to 118 ka time interval. The $\mathrm{P}/\mathrm{Ca}$ symbols are color coded according to fCa, as illustrated in the legend. Green vertical lines highlight the onset of increases in $\mathrm{P}/\mathrm{Ca}$.

3.3 Relationship of $\mathrm{P}/\mathrm{Ca}$ with fCa and δ¹³C_init

Soil pCO₂ and δ¹³C_init are inversely correlated due to contrasting isotopic compositions of respired CO₂ and atmospheric CO₂, an effect simulated by CAVECALC (Lechleitner et al., 2021).

Figure 7 For the stalagmites covering the PGM-LIG period, the $\mathrm{P}/\mathrm{Ca}$ vs. δ¹³C_init, with symbols indicating the fCa. Lines illustrate additionally the modeled variation of $\mathrm{P}/\mathrm{Ca}$ with δ¹³C_init for a drip water of constant initial $\mathrm{P}/\mathrm{Ca}$ (0.2 mmol mol⁻¹) for the case of fCa 0.95 and 0.52. Panel (a) illustrates the results of model 1 where D_max=1.46 and panel (b) illustrates the results of model 3 where D=0.5. The gray rectangle highlights maximum $\mathrm{P}/\mathrm{Ca}$ for BEL, GAR, and GLD.

Most stalagmites spanning PGM to LIG exhibit a wide range in $\mathrm{P}/\mathrm{Ca}$ for a given δ¹³C_init (Fig. 7). Models of P incorporation into calcite show that while soil pCO₂ affects drip water pH and may potentially affect P incorporation into stalagmite calcite, this partitioning effect and concomitant evolution of the dripwater through degassing, would give a much smaller range of variation in $\mathrm{P}/\mathrm{Ca}$ than observed (Figs. 7, S4–S5, Supplement). Stalagmites BEL, GAR, and to a lesser extent GLD, attain peak $\mathrm{P}/\mathrm{Ca}$ concentrations at intermediate δ¹³C_init. For this intermediate δ¹³C_init, the range in $\mathrm{P}/\mathrm{Ca}$ is >2 fold (up to 4-fold in GAR) larger than could be predicted by variations in fCa according to a range of reasonable models either model 1 (D_max=1.46) or model 3 (D=0.5) (Fig. 7). At the same time, in GAR, the minimum $\mathrm{P}/\mathrm{Ca}$ values trend higher with decreasing δ¹³C_init, consistent with the trend predicted by some models (Figs. 7, S4; Model 2 (D_max=1.46) for fCa = 0.95) This trend may be attributable to variations in P incorporation. In NEI, samples of moderately high fCa (0.8 to 0.87) closely follow the modeled trend in $\mathrm{P}/\mathrm{Ca}$ for model 2 (D_max=1.46), thus some of the variation in $\mathrm{P}/\mathrm{Ca}$ may be attributable to variations in P incorporation.

A similar comparison of model and data for the GS22 interval also indicates that the range of $\mathrm{P}/\mathrm{Ca}$ variation exceeds that which is simulated to result from changing P incorporation at a constant dripwater $\mathrm{P}/\mathrm{Ca}$ (Fig. S6).

4 Discussion

4.1 Inference of changing drip water $\mathrm{P}/\mathrm{Ca}$ vs. changing P incorporation

4.1.1 GS22 event

Several lines of evidence suggest that during the GS22 interval, the examined stalagmites have grown from drip water which experienced temporal variations in the initial $\mathrm{P}/\mathrm{Ca}$ ratio, and that variations in the in-cave evolution of the drip water (manifested by fCa changes) are not the only factor varying $\mathrm{P}/\mathrm{Ca}$ in the stalagmites. For the GS22 event, in GLO, GAE, NEI, and ROW samples, the increased $\mathrm{P}/\mathrm{Ca}$ at the onset of the GS22 cold event is expressed in samples of similar fCa and is not likely to reflect enhanced P incorporation in calcite due to a higher fCa and concomitant lower pH. Furthermore, the stadial event is interpreted as a decrease in soil pCO₂, which would lead to higher pH drip waters and lower incorporation according to models 1 and 2 or have no impact in P incorporation according to models 3 and 4. Additionally, the higher $\mathrm{P}/\mathrm{Ca}$ samples do not correspond to samples of higher $\mathrm{Al}/\mathrm{Ca}$ ratios (Fig. S7). Thus, we infer that the increase in $\mathrm{P}/\mathrm{Ca}$ reflects an increase in drip water $\mathrm{P}/\mathrm{Ca}$ across the onset of the stadial event.

It is more complex to distinguish whether trends during the end of GS22 may reflect variations in fCa rather than in initial drip water $\mathrm{P}/\mathrm{Ca}$. At the end of the stadial event, the fCa in GAE increases significantly at the same time as the $\mathrm{P}/\mathrm{Ca}$. Across this transition, the $\mathrm{P}/\mathrm{Ca}$ is negatively correlated with δ¹³C_init. According to models 1 and 2, both the increase in fCa and increase in soil pCO₂ could increase the $\mathrm{P}/\mathrm{Ca}$ in stalagmite even if initial drip water $\mathrm{P}/\mathrm{Ca}$ were constant. In GLO at the end of the GS event, a similar increase in fCa and lower δ¹³C_init coincides with a $\mathrm{P}/\mathrm{Ca}$ increase which may also be affected by changing P incorporation. However, in NEI this further increase in $\mathrm{P}/\mathrm{Ca}$ after the δ¹³C_init recovery is expressed in samples of similar fCa, suggesting that there may have been a second increase in drip water $\mathrm{P}/\mathrm{Ca}$.

4.1.2 PGM-LIG transition

Similar to the GS22 event, the transient $\mathrm{P}/\mathrm{Ca}$ peaks characterizing the onset of many stadial events during the penultimate deglaciation are expressed in samples of similar fCa. These peaks significantly exceed the magnitude of variation in $\mathrm{P}/\mathrm{Ca}$ which can be expected due to altered P partitioning according to any model, and therefore strongly suggest periods of significantly increased $\mathrm{P}/\mathrm{Ca}$ ratios in drip water.

On the other hand, comparison of glacial to interglacial $\mathrm{P}/\mathrm{Ca}$ concentrations is complicated by changes in P partitioning which could accompany the increase in soil CO₂ during interglacials causing decreased initial drip water pH; changes in partitioning could also be driven by the systematic decrease in fCa from the glacial to the interglacial in the majority of samples. GAE features similar fCa during intervals of the glacial and interglacial; this sample features similar $\mathrm{P}/\mathrm{Ca}$ range and variability among both the glacial and interglacial samples. GAR and NEI exhibit slightly higher $\mathrm{P}/\mathrm{Ca}$ during warmer intervals, but with a slope that would be consistent with a pH effect on partitioning according to model 1. Furthermore, a 25 % increase in P incorporation might be expected to accompany a 10 °C warming (Dodd et al., 2021). For these modest differences in observed $\mathrm{P}/\mathrm{Ca}$, a better constraint on the controls on P partitioning would be required to accurately interpret initial drip water $\mathrm{P}/\mathrm{Ca}$ ratios.

4.2 Periods of high P flushing during stadial climate oscillations

The stalagmite records provide evidence for periods of elevated initial drip water $\mathrm{P}/\mathrm{Ca}$ ratios. Because greater bedrock dissolution would proportionally increase both the P and the Ca (and DIC) concentrations in drip water leaving the $\mathrm{P}/\mathrm{Ca}$ ratio unchanged, the periods of elevated $\mathrm{P}/\mathrm{Ca}$ cannot be explained only by greater bedrock dissolution. Rather, they likely reflect increased release of P from non-carbonate reservoirs of P. Independent data on the evolution of vegetation from pollen abundances in coastal sediment cores indicate that there is not a consistent relationship between tree cover and P export. While the abrupt cooling at the onset of GS22 led to decreased vegetation productivity, including a reduction in Atlantic forest pollen from 40 % to <10 % recorded in marine sediment cores off the NW Iberian margin (Goni et al., 2008), during the penultimate deglaciation, temperate tree pollen progressively increased from < 10 % during the PGM to nearly 60 % by 129 ka (Tzedakis et al., 2018) (Fig. 5). Thus, because there is no consistent sign of change in vegetation during the pulses in $\mathrm{P}/\mathrm{Ca}$, we suggest that changes in P storage in aboveground living biomass are not the dominant cause of the peaks in P export.

During the PGM to LIG transition, our age models suggest these pulses may have lasted about 1 kyr. Most events fall during the onset of the cold stadial events indicated from rapid positive shift in δ¹³C_init, but in GS22 the end of the stadial may also feature elevated $\mathrm{P}/\mathrm{Ca}$ and one event recorded by GAR occurs during a weak stadial in between two cold extremes, while in the different cave setting of NEI the pulses were more linked to stadial recoveries than onsets.

Figure 8 Comparison of monthly simulated soil temperatures and snow cover under simulations of weak AMOC typical of stadial periods and simulations of strong AMOC typical of interstadial periods. Results from simulations of Romé (2024), Romé et al. (2022).

Cold stadial interludes are widely attributed to weakened AMOC and therefore we consider the simulations during weakened AMOC to provide indicative information for how stadial climate in the study region differed from glacial climate conditions. GCM climate model simulations suggest hydrological balance (precipitation-evapotranspiration) in this region was similar during cold events of significantly weakened AMOC as during glacial climate with strong AMOC (Romé et al., 2022). However, the infiltration regime may have changed in response to winter temperature and snowcover. During simulated weak AMOC conditions, monthly average top soil temperatures remained below freezing several months of the year, which does not occur in simulated interstadial simulations, and base soil temperatures are colder (Fig. 8). Furthermore, the winter snowcover is simulated to be significantly higher during the weakened AMOC conditions than during the preceding glacial. In combination with the below freezing temperatures, the increased snowcover may suppress infiltration during winter months of stadial periods and lead to large infiltration peaks during spring snowmelt. Changes in the winter temperature distribution also likely affect the intensity and severity of freeze-thaw cycles, a parameter which could be evaluated in future paleoclimate model studies exporting hourly resolution data from paleoclimate simulations.

Modern monitoring studies suggest several aspects of changing cold season climate which contributed to enhanced leaching of P from soil and biomass. Soil microbial cycling is strongly temperature sensitive and soil microbes store large amounts of inorganic and labile organic phosphorus compounds (Gao et al., 2021). Strong decreases in winter soil temperatures lead to lysis of microbial cells, and a decreased soil microbial population may suppress microbial recycling of P, leading to greater loss to groundwaters. In climate chamber experiments, modern soils from cold, snow rich sites featured an up to 400 % increase in soluble phosphate when subjected to high magnitude freeze-thaw cycles (Kreyling et al., 2020). Across a global dataset, freeze-thaw cycles increased dissolved total P and dissolved inorganic P by 312 % and 115 %, respectively. This effect may partly reflect the P release from the death and lysis of soil microorganisms, but it may also arise because ice crystals disrupt the bonding of soil aggregates and accelerate organic matter release (Gao et al., 2021). Finally, frozen topsoil and accumulation of snow change the hydrological cycle by storing water and preventing infiltration during the period of frozen topsoil. If the rate of melting exceeds the infiltration rate when topsoil remains frozen, spring snowmelt may disproportionately enhance surface runoff, and historically the high rates of nutrient loss during snowmelt have been attributed to overland erosion, especially in agricultural catchments (Costa et al., 2020). When areas formerly continuously covered by winter snowfall are subjected to intermittent melt and rain-on snow events, increased frequency of cold season flushing may accelerate nutrient loss, at least via surface waters (Seybold et al., 2022). But, depending on the infiltrability of the soil, meltwater can also infiltrate into the soil, accelerating nutrient flushing (Granger et al., 1984).

Figure 9 Schematic showing components of the soil-karst cave system relevant to the changes in P cycling detected in stalagmites. Panel (a) illustrates continuous infiltration of rainfall through soil and epikarst in warm climates, while panel (b) illustrates the potential for more intense spring snowmelt infiltration and mobilization of soil P through frost cracking during freeze-thaw cycles during onset of abrupt stadial cooling events. Red arrows illustrate P influx, with arrow thickness being proportional to relative magnitude of P influx.

In addition to the direct climatic effects on freeze-thaw frequency and hydrology, the climatic transitions into stadial cold events likely also affected soil pH. More positive δ¹³C_init values during the stadials indicate reduced soil respiration and lower soil CO₂ concentrations, which would raise the pH of soil water. Thus, the elevated $\mathrm{P}/\mathrm{Ca}$ at the onset and termination of GS events could also reflect pH-driven mobilization of mineral-bound P, if intermediate pH conditions destabilized a fraction of the soil P pool. At the onset of the stadial, the positive pH shift may have enhanced the dissolution of Ca-phosphate minerals in calcareous soils and reduced P adsorption to Fe- and Al-(oxyhydr)oxides in more acidic soils, thereby increasing the flux of dissolved P to drip waters. This pathway, acting alongside changes in organic P cycling, may therefore also contribute to the magnitude of the observed $\mathrm{P}/\mathrm{Ca}$ excursions during these climate transitions.Unlike the recent observations of nutrient export or experimental treatment which deal with controls of nutrient release over weeks to months, we have detected events indicating enhanced groundwater P export during hundreds of years to a thousand years. We suggest that during the transitions into stadial cold periods, a combination of more frequent or intense freeze-thaw cycles, suppression of microbial P recycling, and more intense spring infiltration events during snowmelt all increased the P flux in groundwater infiltrating into the cave (Fig. 9). This period of increased flux likely represents a transient imbalance between the P supply from bedrock weathering and the P loss in groundwaters. During several events, during the coldest period of the stadial, P fluxes returned to the lower levels characterizing the period prior to the onset of cold conditions. Potentially, this decline could reflect the depletion of soil P pools following a sustained imbalance in the cycle, given estimates for a mean residence time of P in the soil pool of ≈500 years (Wang et al., 2010). Alternatively or in addition, it may reflect climatic conditions which were less favorable to P export, such as persistent cold with fewer freeze-thaw cycles or soil warming due to the insulation of snowcover which allowed a more effective soil microbial recycling of P. Finally, during some events, enhanced P export also occurs during the climatic recovery from the cold stadial period. Potentially, these transitions may have featured greater frequency of freeze-thaw cycles, or the climatic recovery may have facilitated an increase in bedrock dissolution rates which replenished P stocks sufficiently that P could be exported during these leaching events.

5 Conclusions

Our results suggest that stalagmites can be used to investigate past changes in P flushing into groundwater, but that variations in P partitioning need to be taken into account. Our simulation of P incorporation into stalagmites, evaluating a number of possible partitioning behaviors, confirms that significant variations in stalagmite $\mathrm{P}/\mathrm{Ca}$ can arise despite a constant initial drip water $\mathrm{P}/\mathrm{Ca}$. Future interpretation of $\mathrm{P}/\mathrm{Ca}$ would be significantly enhanced by correction of these effects if the appropriate model for P incorporation can be independently ascertained. Ideally, the appropriate model for P incorporation in speleothems could be assessed with in-cave sampling of dripwaters along a pathway of degassing and calcite precipitation (e.g. as done for carbon isotopes (Mickler et al., 2004; Mickler et al., 2019)), or in laboratory analogues (Hansen et al., 2017; Hansen et al., 2019; Day and Henderson, 2013). Additionally, the influence of other factors such as growth rate or fabrics on P incorporation remains to be explored.

Despite the potential for variation in P partitioning, we find peaks in $\mathrm{P}/\mathrm{Ca}$ for many investigated stadial cooling events which are of a magnitude which exceeds range which could occur from constant dripwater $\mathrm{P}/\mathrm{Ca}$ and variation in P incorporation alone. The reproducibility of several peaks in $\mathrm{P}/\mathrm{Ca}$ concentrations in multiple stalagmites from the same cave system suggests that there are processes affecting drip water P concentration which are not unique to a particular hydrological routing but rather reflective of widespread critical zone processes. Our finding of sustained peaks in P loss in groundwaters during the abrupt climate transitions into and out of stadial cold events, suggests that some processes observed on experimental timeframe of weeks to months may also alter natural P cycling in ways that can be sustained for centuries to millennial. Thus, stalagmites may offer a new opportunity to evaluate the longer term feedback on nutrient cycling in response to climate change.