Wildfires and geochemical change in a subalpine forest over the past six millennia

The frequency of large wildfires in western North America has been increasing in recent decades, yet the geochemical impacts of these events are poorly understood. The multidecadal timescales of both disturbance-regime variability and ecosystem responses make it challenging to study the effects of fire on terrestrial nutrient cycling. Nonetheless, disturbance-mediated changes in nutrient concentrations could ultimately limit forest productivity over centennial to millennial time scales. Here, we use a novel approach that combines quantitative elemental analysis of lake sediments using x-ray fluorescence to assess the geochemical impacts of high-severity fires in a 6200 year long sedimentary record from a small subalpine lake in Rocky Mountain National Park, Colorado, USA. Immediately after 17 high-severity fires, the sedimentary concentrations of five elements increased (Ti, Ca, K, Al, and P), but returned to pre-fire levels within three decades. Multivariate analyses indicate that erosion of weathered mineral material from the catchment is a primary mechanism though which high-severity fires impact element cycling. A longer-term trend in sediment geochemistry was also identified over millennial time scales. This decrease in the concentrations of six elements (Al, Si, K, Ti, Mn, and Fe) over the past 6200 years may have been due to a decreased rate of high-severity fires, long-term ecosystem development, or changes in precipitation regime. Our results indicate that high-severity fire events can determine elemental concentrations in subalpine forests. The degree of variability in geochemical response across time scales suggests that shifting rates of high-severity burning can cause significant changes in key rock-derived nutrients. To our knowledge, these results are the first to reveal repeated loss of rock-derived nutrients from the terrestrial ecosystem due to high-severity fires. Understanding the future of fire-prone coniferous forests requires further documentation and quantification of this important mechanism linking fire regimes and biogeochemical cycles.


Introduction
Fire is an important Earth system process that transfers large amounts of carbon from the biosphere to the atmosphere (Bowman et al 2009). Increased wildfire activity in recent decades (Jolly et al 2015) and under global warming scenarios (Flannigan et al 2009, Moritz et al 2012 has caused concern about the consequences of fire on human health, economic vitality, forest function, and biogeochemical cycles (Bowman et al 2011, Smith et al 2016. However, little is known about the effects of fire on nutrients, which ultimately constrain ecosystem productivity and could serve as either a positive or negative feedback among changes in climate, forests, and fire regimes. Post-fire nitrogen (N) dynamics have been measured on short timescales (<1-8 years) (Smithwick et al 2005), but long-term changes in nutrient stocks remain almost completely unstudied in the context of disturbance regimes (Randerson et al 2009). Although several mechanisms of nutrient loss and accrual after high-severity forest fires have been proposed (Ryan 2000, Wright andHeinselman 2014), it has been difficult to study these mechanisms on appropriate temporal scales. Thus, the long-term biogeochemical consequence of wildfires is a key knowledge gap preventing full understanding of the impacts of this natural disturbance on ecosystem function.
In western North America, debate about the causes and consequences of changing fire activity has centered largely on the relative roles of climate change, land management (e.g., fire suppression, grazing, logging), and interactions with other large-scale disturbances, including drought and tree-killing insect epidemics (e.g., mountain pine beetle) (Littell et al 2009, Harvey et al 2014, Trumbore et al 2015. Lodgepole pine (Pinus contorta Douglas ex Loudon), a broadly distributed tree species, dominates subalpine forests that have historically experienced high-severity (i.e., stand-replacing) wildfires every one to several centuries (Romme andDespain 1989, Sibold et al 2006). Recent increases in large wildfires across the western US have been particularly pronounced in high-elevation forests (Westerling et al 2006), a pattern expected to continue in coming decades (Westerling et al 2011, Barbero et al 2015). Uncertainty over the consequences of increased fire activity on forest ecosystems is a pressing concern to ecologists and land managers alike (Turner 2010). Because biogeochemical processes evolve over timescales of decades to centuries (Smithwick et al 2005), it has been challenging to obtain appropriate data to characterize the biogeochemical impacts of wildfire and changing fire regimes. Lake-sediment records have the potential to overcome these challenges and provide high temporal resolution, long duration, and continuous measurements of both fire events and geochemical responses.
Previous research using lacustrine sedimentary sequences to study the biogeochemical impacts of wildfire demonstrated a significant impact of highseverity fires on N cycling for decades post-fire (Dunnette et al 2014). While these results broadly match and extend chronosequence-based estimates of how disturbance influences N stocks during forest recovery (Perakis et al 2015), the mechanisms of fire effects on N cycling have remained elusive. N cycling is linked through ecosystem processes with other nutrients, notably the base cations of P, Ca, Mg, and K. These base cations are essential for plant growth (Chapin 1980), such that P starts to limit terrestrial primary productivity in old growth tropical forests (Vitousek et al 2010), and Ca limits productivity in base-poor forests (Battles et al 2013). Base cations are contained in both rock and organic matter. The cycling of the base cations thus is fundamentally different from those of N because they are supplied solely from rock sources, and none has a gaseous form under ambient conditions. In addition, the loss pathways of base cations can differ from those of N because of differences in volatilization points and solubility. N volatilizes at a lower temperature than do base cations, and this temperature is easily reached during fire events (DeBano 1991). Alternatively, post-fire losses through erosion pathways may be more significant for base cations (than for N) due to inherently low supply rates. Thus, the response of N after fire events cannot necessarily be extended to other nutrients.
Here we evaluate how high-severity wildfires affect cycling of base cations, primarily with non-volatile loss pathways, using a high-resolution lake-sediment record spanning six millennia. Over this record, we quantified the concentrations of 10 major elements in bulk sediment samples using x-ray fluorescence (XRF) at a median temporal resolution of 4 year/samplesufficient to detect both short-term responses to fire events and long-term trends in response to climate and fire-regime variability over centuries to millennia. This novel technique to measure the concentration of nutrients and other elements in lacustrine sediments is combined with a recently refined approach that isolates high-severity fire events using sedimentary charcoal and magnetic susceptibility . By identifying high-severity fire events within the lake catchment (hereafter 'high-severity catchment fires'), we are working at the appropriate spatial scale to evaluate the biogeochemical impacts of individual fires. We apply this method to a subalpine watershed in the western US dominated by lodgepole pine and characterized by a high-severity fire regime over the late Holocene . We hypothesize detectable transfers of rock-derived nutrients from terrestrial to aquatic ecosystems over millennial timescales, as previously seen in the lacustrine sedimentary record.

Data and methods
Our study utilizes and expands on previous work by Dunnette et al (2014), who developed a high-resolution record of fire history and biogeochemical change for Chickaree Lake spanning the past 4200 years. We expand this record to extend to 6200 years BP, and conducted additional measurements spanning the entire record to quantify the elemental concentrations of 16 elements. Detailed methods for the development of the fire history record are found in Dunnette et al (2014). Below, we briefly describe the site, sample preparation, and fire-history methods, before describing the methods for elemental analyses in more detail.
2.1. Study site and fire history reconstruction Chickaree Lake (40.334°N, 105.847°W) is located in Rocky Mountain National Park, Colorado, USA (figure 1(A)), and is surrounded by an even-aged stand dominated by lodgepole pine (Pinus contorta Douglas ex Loudon) that established after a stand-replacing fire in 1782 C.E. (Sibold et al 2007). The sedimentary record from this small, deep lake (1.5 ha surface area; 7.9 m maximum depth) has a well-constrained chronology based on 210 Pb samples and 25 14 C dates (see figure 1 in Dunnette et al 2014). The ∼31 ha watershed has thin soils with sandy loam texture derived from granite, gneiss, and schist and a thin layer of decomposing litter (US Department of Agriculture NRCS Soil Survey Staff 2007). The historic fire regime is characterized by high-severity crown fires with average return intervals of 100-300 years, associated with severe seasonal drought (Sibold et al 2006). The climate in nearby Grand Lake (c. 8 km from the lake) is continental, with average January and July temperatures of −8.5°C and 14°C, respectively; average total annual precipitation is 483 mm, with an average annual snowfall of 3503 mm (Western Regional Climate Center 1940-2013 observations).
The fire history record for Chickaree Lake is based on a time series of 1486 continuous samples of macroscopic charcoal concentrations (>125 um), converted to charcoal accumulation rates using the estimated sediment accumulation rates from the chronology. The 0.5 cm thick samples on average include 4 years of sediment accumulation. Statistically significant peaks in the charcoal record were identified using the program CharAnalysis version 1.1 (available online at http://charanalysis.googlepages.com), and highseverity catchment fires were further defined as those charcoal peaks that coincided with significant peaks in magnetic susceptibility (MS) ( figure 1(B)). We used the same parameters and methods described in detail by Dunnette et al (2014).

Elemental concentrations with XRF
A total of 768 sediment samples was dried at 60°C, powdered, and analysed for a suite of biogeochemical proxies. The samples were 1-2 cm thick throughout the sequence, and half a centimetre thick for the 10 samples immediately preceding and following a fire event. Carbon (%) was analysed as described by Dunnette et al (2014), while all other analyses are new to this study and described below. We measured the concentration of the 10 major elements traditionally listed as oxides (Mg, Al, Si, P, S, K, Ca, Ti, Mn and Fe) and six trace elements (Cu, Zn, As, Pb, Rb and Mo) using a handheld wavelength-dispersive XRF spectrometry with a Bruker Tracer III SD. The major elements were measured for three minutes per sample at 15 kV, with a vacuum attached to the instrument to reduce background noise. Trace elements were measured for two minutes per sample at 40 kV, using a Ti-Al metal filter to increase the magnitude of the trace element spectra, and decrease the sensitivity to elements lighter than Ca. Data were measured as counts of photons emitted that are unique to each element, and were calibrated with the mudrock calibration (Rowe et al 2012), which is a combination of 200 standards of similar matrix (grain-size, homogeneity) and composition, to obtain a concentration value (mg/100 mg of dried powdered lacustrine sediment). Repeatability of this method was tested on a previous set of lacustrine sedimentary samples and analytical error was ±2% for both major and trace elements (Leys et al 2016).
To help interpret the source materials contributing to the bulk lake sediment samples, we measured the elemental concentration of targeted material from the sediment core and from within the watershed, including: lodgepole pine needles (n=6), Engelmann spruce needles (n=2), forest litter (n=4), decomposing wood (n=1), organic soil (n=9), terrestrial charcoal (n=1), surface lake sediment (n=4), and sedimentary charcoal (n=1). We also obtained existing elemental concentration measurements from the mineral soil (n=3) from a nearby site with similar bedrock (see appendix for more details).

Accuracy of XRF-derived concentrations
We tested the trends between count data, ratios to Ti, and calibrated data for all of the nine elements, excluding Ti (example shown with Ca, figure S1). The Pearson correlation coefficient between the count and the calibrated data is high (0.98 for Ca, appendix, table S1), indicating that the count and the calibrated data are similar in trends and magnitudes (figure S1). The calibration process transforms count to a concentration, facilitating comparison of multiple elements as well as calculating budgets, or rates of loss and supply. The protocol used for the count acquisition is based on more than 200 samples of regular matrix, which differs from the more commonly used scanning XRF methods, which are based on uncalibrated counts of wet sediment only. The comparison between counts obtained using our protocol and results from scanning XRF should be tested in the future.
We also calculated the Pearson correlation coefficient between calibrated elemental concentration and the ratio of count of elements to Ti count (e.g. Ca calibrated and Ca/Ti ratio). The comparison of the correlation between ratios of elements to Ti count and the calibrated elemental concentrations is overall lower than correlations between the count and calibrated element concentrations (e.g. equal to 0.41 for Ca/Ti compared to Ca calibrated) (table S1). Sometimes the ratio of elements is used to amplify signals that are not significant in the count data (e.g. Löwemark et al 2011, Chawchai et al 2016). We found that ratios of element to Ti counts and elemental calibrated concentrations shared the same long-term trends but differed in the magnitude of variability, and in some cases the directions of extreme values differed (figure S1). In our study, the magnitude and direction of the values, especially the extreme values, are fundamental to testing the response of elements to discrete fire events.
We quantified concentrations of eight trace elements and 19 rare Earth elements on 13 samples by inductively coupled plasma mass spectrometry (ICP-MS, conducted at the GeoAnalytical Lab at Washington State University) to estimate biases with the standard XRF calibration technique (figure S2, appendix). The two methods were strongly correlated (r 2 >0.80; table S2) for elements displaying concentrations close to or over 0.01 mg/100 mg, with a mean absolute percentage error of 4.8% (table S2). Although the ICP-MS method is the most accurate to quantify the absolute concentration of elements in a sample, it was not possible to conduct ICP-MS analysis at high resolution in this study due to high sample volume requirements. The high correlation between elemental concentrations obtained from calibrated XRF counts and the ICP-MS method demonstrates that using a handheld XRF with a mudrock calibration provides accurate measurements of 10 elements (Mg, Al, Si, P, S, K, Ca, Ti, Mn, and Fe).

Statistical analysis
We assessed the geochemical impacts of fire events using superposed epoch analysis (SEA, figure 2), following methods described by Dunnette et al (2014). The analysis was conducted on high-severity catchment fires (n=17) and a subset of lower severity/ extra local fires (n=13) using custom scripts written in MATLAB (Higuera and Dunnette 2014). All highseverity catchment fires were included in our analyses. Although four high-severity catchment fires occurred within 100 years of a previous high-severity catchment fire, our results were similar with those fires included or excluded from the SEA.
We used principal component analyses (PCA, Dray and Dufour 2007, R Core Team 2014) to summarize (i) the relationships among elements through time (figure 3(A)), and (ii) to compare lake sediment composition to plant remains, organic and mineral soils, charcoal, and surface lake sediments ( figure 3(B), appendix, table S2). Both SEA and PCA were performed using only the 10 major elements of highest interest for nutrient and ecosystem dynamics (Mg, Al, Si, P, S, K, Ca, Ti, Mn, and Fe).  facilitate a broader goal, beyond this work, of reconstructing elemental budgets in an ecosystem.

Patterns and mechanisms of change over years to decades
We demonstrate sensitivity of sedimentary geochemical proxies to wildfires on decadal timescales through analysis of proxy response to 17 high-severity catchment fires using SEA (figure 2). Concentrations of five elements-Ti, Ca, K, Al and P-increase significantly after fires, and remain elevated for an average of 15-30 years post-fire. In contrast, no significant geochemical changes were observed after low-severity/ extra-local fires ( figure S4). These data indicate that erosion is likely an important mechanism through which fire activity affects elemental composition of lacustrine sediments and catchment soils. Specifically, mineral material high in Ti, Ca, K, Al, was weathered from the source granitic rocks of the catchment and transferred to the lake after fire events, thereby representing a loss from the terrestrial system. These post-fire erosion events represent disturbances to both the terrestrial and aquatic ecosystems (Pickett and White 1985). The importance of erosional process is consistent with observations in the Greater Yellowstone Ecosystem where lodgepole pine forests experience stand-replacing fires that induce high tree mortality, combust the thin organic litter and upper soil layers, and expose mineral soil (Certini 2005). Although we see no evidence in the Chickaree Lake record of the fire-related debris flows, which can contribute significantly to sediment loads in fluvial systems (Meyer et al 1995), there seems to be sufficient precipitation to deliver mineral material to the lake sediments through overland transport in the catchment.
While concentrations of most elements consistently increase in sedimentary samples deposited immediately after high-severity catchment fires, other elements display two different patterns (figure 2). The concentration of S decreases after high severity catchment fires and remains low for three decades. Both S and N are considered sensitive to fire because they have low volatilization temperatures of 375°C and200°C, respectively (DeBano 1991, DeBano et al 1998). In contrast to other nutrients such as K, P, Mg, Ca, or Mn, S is readily volatilized from organic matter during combustion, particularly after highintensity fires which can exceed temperatures of 1000°C in forested systems (DeBano 1991). Thus, the decreases in S after high-severity catchment fires were most likely caused by volatilization, as opposed to erosion. This is consistent with the findings of Dunnette et al (2014), who concluded that combustion drove multi-decadal declines in terrestrial organic matter subsidies to the lake. Alternatively, the decreases in S concentrations could have also been caused by the dilution of S in post-fire material inputs. Si demonstrates no change in concentrations after high-severity catchment fires. The neutral response of Si indicates either an absence of fire impacts on Si, or that the two sources of Si-biogenic and minerogenic-vary in opposite directions and therefore result in no net change in lacustrine sedimentary concentration. Postfire nutrient pulses to the lake (e.g. Kelly et al 2006) could stimulate a short-lived increase in lake productivity, and sediment biogeochemistry suggests that diatom productivity has been relatively high in Chickaree Lake (Dunnette et al 2014). However, minerogenic Si from the parent material (i.e. granite and ilmenite) might be simultaneously decreasing relative to other sedimentary components.
When all elements measured over the entire sedimentary sequence are combined into a PCA, a coherent pattern emerges highlighting the similar behaviour of eight elements ( figure 3(A)). Axis 1 in the PCA explains 40% of total variation in the elements over time, with the strongest correlation to Ti (−0.47), followed by K, Al, and Ca (−0.42, −0.35, and −0.34, respectively). These four positive post-fire responders are followed by Fe, Mg, P, and Mn, which occur in lower concentrations, and are less correlated to Axis 1 (−0.32, −0.30, −0.27, and −0.25, respectively). Axis 2 is primarily driven by Si and Mn-non-fire-responsive elements. Samples following lower-severity/extra local fires are a subset of those from high-severity catchment fires, but centered on zero, indicating little geochemical change relative to average ( figure 3(B)). Sediment samples following high-severity catchment fires also occupy more ordination space and are shifted to negative values of Axis 1, indicating high values of Ti, K, Al, and Ca. Axis 1 therefore serves as an integrated geochemical proxy that reflects in part the impacts of high-severity catchment fires on elemental concentrations, generally reflecting mechanisms linked to soil erosion. Although MS is also an erosion proxy, used here to help identify high-severity wildfires, less than 1/3 of the variability in MS can be explained by Axis 1 values (i.e., Pearson's correlation of −0.53, table S3), indicating that Axis 1 contains unique information.
To better understand the composition of the material delivered to the basin after high-severity catchment fires and the origin of the sedimentary matrix, we conducted a second PCA that included all 768 sedimentary samples in addition to 31 watershed samples of organic and clastic material (figure 3(C)). Axis 1 explains 34% of the total variation, and is well correlated with Ti, Al, Fe, K, and P concentrations. Axis 2 explains 25% of the total variation, and positive values reflect high Mn, Ca, and Mg concentrations whereas negative values reflect high S concentration. Watershed samples display a wide range of geochemical characteristics, with plant tissue, organic soil, and mineral soil strongly separated by their elemental concentrations ( figure 3(D)). Organic matter and lake sediments are dispersed along Axis 1, whereas plant materials are dispersed along Axis 2. Samples from the lake sediments occupy a distinct area of ordination space. Although there is some component in the sediments that is not being measured with current samples (indicating a likely aquatic source), lake sediments in general are closest to processed organic matter, not plant tissues. Mineral soil samples have the most negative values on Axis 1, indicating that one consequence of high-severity catchment fires is the removal of upper soil layers, likely down to the mineral horizon.

Patterns and mechanisms of change over centuries and millennia
Although the geochemical patterns of the Chickaree Lake sedimentary record are dominated by highfrequency variability largely associated with highseverity catchment fires, the record also features millennial-scale trends that are important to understand, especially in the context of ongoing environmental change in subalpine forests of western North America. The main long-term trend is that the integrated geochemical metric, PCA Axis 1, increases toward present, indicating a long-term decrease in the concentration of elements associated with fireinduced erosion ( figure 4(A)). There are several potential causes of this change, which are not mutually exclusive and have different implications for the future of forest ecosystem function.
First, the millennial-scale decrease in elemental concentrations associated with fire-induced erosion coincides with an observed decrease in the frequency of high-severity catchment fires over the past several millennia ( figure 4(B)). This shift occurred in the context of a millennial-scale trend from a rain-dominated to snow-dominated precipitation regime, and gradual summer cooling over the past 6000 years These changes could be at least partially responsible for the decrease in high-severity catchment fires. Decreased fire activity c. 2000 yr BP is also consistent with lower biomass burning inferred from a network of sites in Rocky Mountain National Park (Higuera et al 2014), although we caution that the observed pattern at Chickaree Lake represents processes occurring within or very near the watershed, and the ability to focus on high-severity fires at Chickaree Lake is unique among regional fire-history records. Climate variability could have also directly altered weathering rates of the granitic parent material, with a warmer, rain-dominated regime leading to faster chemical weathering rates compared to a snow-dominated regime (Schaller et al 2010), potentially independent of fire activity. Second, changing geochemical elemental concentrations were occurring in the context of the longterm subalpine forest ecosystem development around Chickaree Lake, including vegetation dynamics, terrestrial organic matter accrual, and landscape stabilization. The pollen records from Chickaree Lake and other lakes in Rocky Mountain National Park show no evidence of large-scale vegetation changes over the past six millennia; rather, the Chickaree Lake watershed was consistently dominated by lodgepole pine and other subalpine forest species (figure 4(D); Dunnette et al 2014). Soil development and stabilization over millennia, on the other hand, could influence the quantity or makeup of geochemical material delivered to the lake after wildfires (Engstrom et al 2000), potentially altering the sensitivity of elemental concentrations to fire-induced erosion. A link between organic matter accrual or landscape stabilization and the sensitivity of subalpine forests to fire-induced ecosystem change would have important implications for ongoing and future environmental change. Specifically, such a link implies that ecosystems in earlier stages of development would be more vulnerable to fireinduced changes in elemental concentrations. Although the Chickaree Lake record suggests less fireinduced erosion towards present, it also highlights that any increase in resilience to high-severity fires could be reversed if the frequency or severity of fire increases, as seen after the most recent high-severity catchment fire around Chickaree Lake, in 1782 CE ( figure 1(B), Sibold et al 2006). Taken together, the geochemical record from the Chickaree Lake watershed suggests that elemental concentrations were shaped by interactions among climate variability, fire-regime variability, and long-term ecosystem development.

Conclusion
Our results indicate that high-severity fires in a lodgepole pine ecosystem result in a depletion of elements from watershed soils, including the nutrients Ca, K, Al, and S, for approximately three decades. The mechanisms driving these changes appear to be primarily erosion, and secondarily volatilization (of S). These multi-annual and decadal-scale changes are overlain on millennial-scale changes in elemental concentrations that likely reflect a suite of ecosystem processes, linked to climate and disturbance-regime variability as well as long-term ecosystem development. However, there are important remaining questions about the millennial-scale effect of these repeated disturbance events. To determine if terrestrial nutrient losses were significant enough to limit net primary productivity of lodgepole pine forests, future research should focus on quantifying the amount of base cation transfer to the lake during the past six millennia. The approach presented here takes an important first step toward matching fire reconstructions with quantitative geochemical reconstructions. Additional site-level studies across a broad spatial extent will be necessary to calculate the impacts of fire regimes on the mass-balance of elements. Ultimately, the characterization of feedbacks between biogeochemical cycling and fire regimes will improve predictions of future ecosystem function under scenarios of changing climate and disturbance regimes.
Appendix. Detailed material and methods information on the XRF-derived concentrations of elements, and the statistical analyses A.1. Vegetation reconstruction The pollen record from Chickaree Lake was used unmodified from Dunnette et al (2014).

A.2. Tests on XRF-derived concentrations
We tested the trends between count data, ratios on Ti, and calibrated data for all of the 10 elements (table S1). The correlation coefficient between ratios of elements to Ti count, and the calibrated elemental concentrations is overall lower than the count and calibrated element concentrations correlation coefficient (table S1). An example is shown with Ca in figure S1.
We quantified concentrations of eight trace elements and 19 rare Earth elements on 13 samples by ICP-MS (conducted at the GeoAnalytical Lab at Washington State University) to estimate biases with the standard XRF calibration technique (figure S2). Because ICP-MS requires a sample mass of at least 1 g, sediment samples were chosen throughout the core based on their weight, minimizing the number of contiguous samples needed to attain this mass. These eight trace elements represent a small portion of the sediment composition, and their concentrations are the lowest recorded (30 to 300 times lower than the base cation concentrations). Because of low sample mass in several elements (e.g., Th, U, and Nb), many samples are close to the detection threshold; thus, the calibration results for those elements are the least accurate.
The eight trace elements quantified both by ICP-MS and calibrated XRF were compared with correlation analysis. Pb, U, Th, Nb were present in the lowest concentrations, and also had the lowest correlations between ICP-MS and XRF-derived concentrations (table S2). Y, Zr, Sr, Rb displayed concentrations close to or over 0.01 mg/ 100 mg, and the two methods were strongly correlated (r 2 >0.80; table S2). The mean absolute percentage error calculated for these four elements is 4.8% (table S2).

A.3. Superposed epoch analysis (SEA)
To assess the impacts of wildfires on geochemical proxies, we followed methods described in Dunnette et al (2014) for SEA. A slightly modified version of these methods are provided here for reference, taken from Dunnette et al: 'We assessed the biogeochemical impacts of fire using SEA, a nonparametric method used to evaluate the average response to multiple events in a time series (Adams et al 2003). The analysis was conducted on high-severity catchment fires (n=11 [17 in the current study]) and a subset of lower severity/extra local fires (n=9 [13 in the current study]) using custom scripts written in MATLAB (available online at doi: 10. 6084/m9.figshare.988687). Response variables [of the 10 elements measured in this study (Mg, Al, Si, P, S, K, Ca, Ti, Mn, and Fe)] were sampled every 0.5 cm for c. 5 cm before and after each fire event. Outside of this window, subsamples from two consecutive 0.5 cm samples were analysed at c. 1 cm intervals. The highresolution sampling window represents c. 50-125 year (depending on local sediment accumulation rates). The SEA was performed specifically on samples from 50 year pre-fire to 75 year post-fire (a span selected based on trends during lodgepole pine stand development; Pearson et al 1987). Because many samples c. 50-75 year post-fire were not included in high-resolution sampling, composite patterns over this span may be more muted than if they were sampled at higher resolution.
Before SEA, samples were interpolated to the median sample resolution (for biogeochemical samples) of 5 year, and low-frequency trends were summarized with a 500 year locally weighted regression robust to outliers [ figure S3]. To minimize bias arising from long-term changes in response variables, we subtracted the low-frequency trends from the interpolated time series to obtain residual series. Residuals were averaged across events to produce composite time series (hereafter 'response series') showing each variable's mean response to fire events, in 5 year bins.
To assess the statistical significance of the response series, confidence intervals (CIs) were generated with a Monte Carlo randomization method. For each response variable, 10 000 time series were created by randomly shuffling samples in five-sample blocks to account for temporal autocorrelation. Each random time series was sampled via the same method used to generate the observed response series, producing random composites. The 0.5th, 2.5th, 97.5th, and 99.5th percentiles from the 10 000 random composites were used to construct 95% and 99% CIs.' In the current study, all high severity catchment fires were included in the SEA. Although there were four high severity catchment fires that occurred within 100 years of another high severity catchment fire, the direction and significance of SEA results did not differ if those fires were excluded from the SEA. Response of elements to high severity catchment fires, and low severity/extra local fires are presented in figure S4.

A.4. Principal component analyses (PCA)
PCA was conducted using the R software package (R Core Team 2014) to assess (i) the relationship among elements, and (ii) to compare lake sediment composition to plant remains, organic and mineral soils, charcoals and surface lake sediments. The first PCA was conducted on the 10 element concentrations (Mg, Al, Si, P, S, K, Ca, Ti, Mn, and Fe) of the 768 sediment samples from Chickaree Lake. The second PCA was run on 9 element concentrations of the 768 sediment samples from Chickaree Lake, and 31 samples from the watershed, as described in section 2.2 in the main text. The mineral soil samples correspond to the top 20 cm, the A horizon, and the C horizon of the US soil geochemical landscape site #12891, Grand County, Colorado. The elemental concentrations were obtained by ICP-MS technique. The PCAs were centred and scaled, with no axis rotation.
A.5. Correlations between biochemical proxies and the 10 rock-derived elements To test correlations between (i) previously published biogeochemical proxies (%N, %C, C:N, MS, charcoal accumulation rate, δ 15 N, and δ 13 C) by Dunnette et al (2014), and the 10 elements selected here (Mg, Al, Si, P, S, K, Ca, Ti, Mn, and Fe), and (ii) in among the 10 elements, we performed both Spearman and Pearson correlation analyses (table S3).
None of the Spearman and Pearson correlations are over 0.65 between the 10 elements and the biogeochemical proxies, indicating an absence of strong relationship between atmospheric source nutrients, and rock source nutrients. In addition, Ti and MS, both used in the literature as erosion proxies, are not strongly correlated (0.44 and 0.48 Spearman and Pearson correlation coefficients, respectively).
Among the 10 elements, there are strong positive correlations between the four elements responding positively to high-severity catchment fires, and explaining the PCA axis 1: Ti, K, Al, and Ca (table S3). All the negative correlations for both Spearman and Pearson tests involved Si or S in the pairwise comparison (table S3), indicating that S and Si are the two elements independent to the others in their source or loss pathways.