Cover Sheet Title : Linking silicon isotopic signatures with diatom communities

13 The use of silicon isotope ratios (expressed as 30Si) as a paleolimnological 14 proxy in lacustrine systems requires a better understanding of the role of lake processes 15 in setting the 30Si values of dissolved Si (dSi) in water and in diatom biogenic silica 16 (bSi). 30Si of modern dSi ( SidSi) and bSi ( SibSi) in three lakes 17 in Lassen Volcanic National Park (LAVO), California (USA), and produced diatom 18 assemblage compositional data from the modern system and from sediment core 19 samples. In the modern systems, we observe the largest magnitude diatom Si isotope 20 fractionations yet reported, at -3.4 and Fragilaria dominated samples. Using 21 statistical approaches designed to condense multivariate ecological data, we can 22 deconvolve assemblage-specific Si isotope fractionations from the combined diatom 23 assemblage30Si data. For example, samples dominated by generally deeper water 24 Diatom silicon isotope ratios in Lassen Volcanic National Park 2 SibSi values (< . Conversely, samples dominated 25 by shallow water planktic or benthic periphyton SibSi values (> 26 30Si records from LAVO lakes reflect species specific Si 27 isotope fractionations and thus act as paleolimnological proxy for the aquatic-habitat of 28 bSi production. Silicon isotope analysis should be coupled with diatom community 29 composition data and other geochemical proxies for the most robust paleolimnological 30 interpretations. We also construct a Si mass-balance for Manzanita Lake based on 31 elemental fluxes. Despite a short residence time of ~4 months, it is an efficient Si sink: 32 about 30% of inflowing Si is retained in the lake sediments. An entirely independent Si 33 isotope-based estimate agrees remarkably well. Burial fluxes of bSi derived from 34 radiometrically dated sediment cores yield retention rates of about a factor of three 35 higher, which might suggest groundwater is an important term in the lake Si budget. 36 Diatom silicon isotope ratios in Lassen Volcanic National Park 3

deconvolve assemblage-specific Si isotope fractionations from the combined diatom 23 assemblage-30 Si data. For example, samples dominated by generally deeper water 24 Diatom silicon isotope ratios in Lassen Volcanic National Park 2 30 SibSi values (< -. Conversely, samples dominated 25 by shallow water planktic or benthic periphyton 30 SibSi values (> -26 30 Si records from LAVO lakes reflect species specific Si 27 isotope fractionations and thus act as paleolimnological proxy for the aquatic-habitat of 28 bSi production. Silicon isotope analysis should be coupled with diatom community 29 composition data and other geochemical proxies for the most robust paleolimnological 30 interpretations. We also construct a Si mass-balance for Manzanita Lake based on 31 elemental fluxes. Despite a short residence time of ~4 months, it is an efficient Si sink: 32 about 30% of inflowing Si is retained in the lake sediments. An entirely independent Si 33 isotope-based estimate agrees remarkably well. Burial fluxes of bSi derived from 34 radiometrically dated sediment cores yield retention rates of about a factor of three 35 higher, which might suggest groundwater is an important term in the lake Si budget. 36

Introduction 37
Diatoms, through their prodigious uptake of dissolved Si (dSi) and rapid 38 mineralization of biogenic Si (bSi), are central to freshwater Si cycling. Understanding 39 diatom uptake and recycling of this often limiting nutrient therefore shines light on the 40 continental Si cycle (Frings et al., 2014) and lake phytoplankton ecology (e.g. Kilham, 41 1971). Silicic acid is delivered to lakes from either surface water or groundwater, and 42 ultimately derives from the weathering of silicate minerals. Within-lake recycling of 43 biogenic silica is also important, as most freshwater systems are strongly undersaturated 44 with respect to bSi, particularly during periods of enhanced biological activity where 45 Si is removed from the water column. Given its role as a key nutrient in lake ecology, 46 and its coupling to the carbon cycle via the silicate-weathering feedback (Frings, 2019), 47 there is considerable incentive to be able to reconstruct aspects of the past Si cycle. As 48 integrators of catchment biogeochemistry, lake sediment archives are ideally placed to 49 achieve this, but estimates of bSi concentrations or accumulation rates can paint an 50 incomplete or misleading picture (Nantke et al., 2021). 51 Silicon isotope ratios (expressed 30 Si) are emerging as a powerful tool to 52 trace Si biogeochemistry. Differences in 30 Si between two phases result from isotope 53 fractionation, that can be expressed as a fractionation factor , defined as A-B = RA/RB, 54 where R is the ratio 30 Si/ 28 Si in phase A and B. Since A-B is typically very close to 55 unity it can also be presented in notation (Coplen, 2011) Si of the dSi initially supplied (e.g. of 85 the river flow into a lake), and is the fraction of initial dSi remaining (i.e. at , 86 no diatom growth has occurred, and at , all available dSi has been used). 87 indicates the instantaneously produced bSi, which is always one fractionation offset 88 from the coeval dSi. The cumulative product is given in Eqn. 1c. The alternative model 89 is a steady-dSi with 90 constant 30 Si into the system, of which a fraction is converted to bSi: 91 Eqn. 2a 92 Eqn. 2b 93 30 SibSi and 30 SidSi will evolve towards higher values, but along 94 different trajectories than that of the Rayleigh model. In this model, the bSi 95 instantaneously produced and the cumulative product have the same composition. Note 96 30 Si as a function of reaction completeness in a closed system but 97 where both forward and reverse reactions occur at equilibrium has the same 98 mathematical -99 -Two endmember possibilities for 100 the bSi produced are shared by both models. Firstly, at 0% utilization (i.e.
) an 101 infinitesimally small amount of diatom Si uptake will produce bSi offset by exactly one 102 fractionation from the source dSi. Secondly, at 100% Si utilization (i.e.
) the bSi 103 must have the same 30 Si as the source dSi ( ). The models differ in their 104 predict 30 SibSi 105 with increasing dSi utilization. This forms the basis of silicon isotope ratios as a 106 paleoproductivity proxy (De La Rocha et al., 1998). From this simple framework, the 107 source pool dSi isotope composition (i.e.
), and the Si isotope fractionation 108 associated with diatom production, emerge as key parameters. If these are constrained, 109 30 Si in terms of Si utilization (i.e. 110

). 111
There is a growing body of research using Si isotopes as proxies for diatom 112 productivity in freshwater systems (Chen et  Here, we explore the influence of diatom community composition and lake 132 meta-, and hypolimnion of each lake, and from the surface (0 m) in UMC and LMC. 226 (Epilimnion, metalimnion and hypolimion refer to the upper, intermediate and lower  227 thermal layers in a stratified lake). Samples for 30 Si analysis of ca. 100 ml were filtered 228 in the field through 0.4 µm polycarbonate (PC) filters (Whatman-Nuclepore) and stored 229 in acid-cleaned, opaque, polyphenylene ether (PPE) plastic bottles and kept cool in the 230 dark until return from LAVO. Chl-a and pheophytin were measured from 100 mL of 231 water via fluorometry (Turner Designs model 10AU Fluorometer) using methanol as a 232 solvent to determine algal biomass (Welschmeyer, 1994). Calibration was conducted 233 using a spectrophotometric method (Parsons et al., 1984) and a standard, Chl-a from 234 Anacystis nidulans (Sigma Corp.). Samples of lake diatoms for identification, 235 enumeration, and silicon isotope analysis were taken from: 1) surface tows (~300 m) 236 using a 20 µm plankton net; and 2) from substrate scrapes at the lake margin, and 3) 237 near surface inflow and outflow (ML only). The plankton tow and substrate scrape 238 samples were stored in acid-cleaned PPE bottles and kept on ice and then refrigerated 239 before sample preparation. 240

Sediment core samples 241
Lake cores were stabilized with Zorbitrol (sodium polyacrylate absorbent 242 powder) upon recovery in the field and then refrigerated before being sent to the 243 National Lacustrine Core facility (LacCore) in Minnesota, USA. At LacCore, the cores 244 were split and subsampled at a 0.5 cm or 1 cm resolution. Sediment samples from the 245 cores were freeze-dried and stored in opaque polyethylene containers before further 246 sample preparation and analysis. Age-depth models of the BL, ML, and WL cores were 247 established by 210  Laboratory (Northern Arizona University) using a wet-alkaline method of extraction 257 (Mortlock and Froelich, 1989). Briefly, ~20 mg of crushed freeze-dried sediment is 258 leached for 5 hr in 40 ml of 2M Na2CO3, then neutralized and analyzed for dissolved 259 Si concentration spectrophotometrically. 260

Diatom enumeration 261
Prior to cleaning, boil-and-burn mount slides were made from all samples to 262 examine the algal communities, including colony formation. Preparation of modern and 263 sediment samples for diatom enumeration followed protocols described in Battarbee et 264 al. (2001) and Stoermer et al. (1995). Diatom community enumeration was carried out 265 on samples also analyzed for 30 SibSi, including modern plankton, periphyton, and 266 sediment. Permanent slide mounts were made using ZRAX and were used to determine 267 relative abundance counts (n = 500). Samples were also examined during counting for 268 presence of chrysophytes and sponge spicules and presence was recorded in count data. 269 They were rare in all cleaned samples (<1% of counts) and are not discussed further. 270 Diatoms were identified and counted at 1000x using an oil immersion lens with DIC 271 (differential interference contrast) on an Olympus BX51 microscope. Grace, 2002). A 0.4 is considered (Peck, 2016), though this may 352 change as our understanding of this novel statistic improves (Peck, 2016). 353

Lake characteristics 355
A summary of catchment characteristics and lake physicochemical parameters 356 are reported in Table 1 for ML, WL, and BL. An in-depth discussion of catchment 357 characteristics, lake physicochemical parameters, modern diatom community 358 succession, and Chl-a values for Butte Lake (BL) is presented in Howard and Noble 359 (2018). Table 2 presents the modern water samples and phytoplankton silicon isotope 360 compositions and Chl-a concentrations (a proxy for phytoplankton biomass). All lake 361 characteristic data (e.g. oxygen, chlorophyll, temperature, conductivity) are given in 362 the online supplementary data file. Diatom relative abundance counts, sediment core 363 30 SibSi values are given in the online Supplementary Data. 364 Diatoms dominate the phytoplankton samples from all lakes, though BL and WL have higher mean Chl-a values than ML (Table 2). Both ML and WL appear to have a 366 subsurface Chl-a maximum in the summer. As with BL (Howard and Noble, 2018), 367 ML and WL have a consistent, planktic diatom succession following ice-out in the 368 spring, through to the fall, although the community composition differs between the 369 lakes. . 30 SibSi varies seasonally, and between 387 phytoplankton and periphyton (Table 3; Figure 5) 30 Si 388 values were higher than plankton dia 30 Si values (Table 3). In ML, Chl-a values suggest that the greatest productivity occurs in the meta-and hypo-limnia with the 390 exception of September 2014, when the highest productivity occurred in epilimnion 391 (Table 2). 392 In ML, Asterionella formosa is abundant in May and June, transitions to 393  Staurosira construens var. binodis (Howard and Noble, 2018). In addition to species 417 differences, total diatom productivity differed between seasons; total planktic diatom 418 biovolume decreased by ~50% between spring and summer 2014 (Howard and Noble, 419 2018). 420

Widow Lake 421
In contrast to ML and BL, the smaller WL does not exhibit sustained 422 stratification, and was nearly isothermal at all sampling dates with the exception of June 423 2012 ( Figure 3B). Vertical profiles of conductivity and DO also suggest that WL is a 424 polymictic system where frequent mixing is interspersed with short-term stratification. 425 dSi concentrations are relatively consistent across the sampling period, at 2.4 to 3 mg/L. 426 30 SidSi values of the three lakes 427 (Table 3). Widow Lake has a different diatom succession than its neighboring lakes; 428 Tabellaria flocculosa and Fragilaria tenera-nanana group dominate in spring, 429 followed by F. crotonensis and T. flocculosa in the summer, and F. tenera-nanana 430 group in the fall. 431

Manzanita Lake 433
Down-core 30 SibSi values, percent bSi, and diatom data are shown in Figure  434 5A. Around 1942 CE (Common Era), percent bSi is relatively low (<10%) and a 435 relatively high 30  Distinct trends are apparent both below and above the tephra layer ( Figure 5C identified for 4 of these 5 groups (Group 2 did not yield any significant indicator taxa 503 at the = 0.05 level) and are given in Table 3. 504 . 505

Patterns of Si utilization in LAVO lakes 507
Relatively high dSi concentrations (~16 mg/L Si at the UMC inflow, and ~12 508 mg/L at the LMC outflow; Table 3) compared to the other lakes likely reflect a large 509 drainage ratio, large capacity for surface water runoff (curve number = 63), and the 510 dominant rhyodacite bedrock (Table 1) (Table 3), lower than epilimnetic values and likely 522 reflect mixing of dSi pools from the metalimnion, hypolimnion, or groundwater into 523 outflow water. 30 Si suggests most biosiliceous production occurs in 524 the upper water column. Chl-a values are a general measure of productivity for all algae 525 groups, and in ML are higher in the meta-and hypolimnia compared to the surface 526 (Table 2). Higher Chl-a at depth in ML may result from non-diatom algal groups (e.g., 527 chlorophytes or cyanophytes). Alternatively, the values might reflect export or 528 migration of surface diatom production from the epilimnion to deeper depths. The  SidSi values and dSi concentrations vary slightly over the growing season in WL, 556 suggesting changes in dSi utilization and very weak or short-term stratification (Table  557 3, Figure 3B). 558

Distinguishing between fractionation models 559
30 SidSi and 30 SibSi as a function of Si utilized for ML, we can 560 attempt to distinguish between the two endmember fractionation models (see 561 introduction; Figure 7). The x-axis in this plot, i.e. the degree of dSi utilization, is hard 562 to constrain. We take as an approximation, where is the Si 563 concentration of inflows to ML via UMC at the time of sampling and is the Si 564 concentration of the epilimnion in ML. Accurately accounting for groundwater is not 565 possible, since no wells exist within the park (water supply for consumption is surface 566 water sourced), so we assume the groundwater has the same dSi concentration as UMC. 567 The groundwater dSi should be measured in future work. UMC is thus used as the 568 closest proxy available for background Si concentration of inflows in the absence of 569 biological Si isotope fractionation. Figure 7 shows that the overall range of Si utilization 570 (ca. 19 35%; see also section 5.4 below) is not large enough to be able to distinguish 571 Rayleigh-style system evolution (Eqn. 1) from a steady-state model (Eqn. 2). The 572 reality likely falls somewhere between these two endmember scenarios; for much of 573 the growing season, diatom growth rates exceed dSi supply, causing a transient 574 depletion in lake water Si inventory. This lack of balance between supply and removal 575 invalidates the steady-state assumption of Eqn. 2. On the other hand, the non-zero 576 supply of new or recycled dSi, invalidates the assumptions inherent in a Rayleigh model. Where a system behavior falls between these two models will depend on the 578 rate of bSi production relative to supply; when it is much greater, a Rayleigh model 579 (Eqn. 1) will be closer to the truth, but when they are more closely matched, a Steady-580 State model (Eqn. 2) may best capture the dynamics of the system. Overall, even simple 581 lake ecosystems cannot be condensed to the models typically used. Seasonal 582 imbalances between supply and demand, combined with variable Si isotope 583 fractionation factors (Section 5.2 below), suggest more nuanced models that capture 584 time-dependence, are required if we are to fully describe system behavior. 585

Diatom silicon isotope fractionations 586
Diatom Si isotope fractionation factors, degree of Si utilization, and degree of 587 system openness combine to define 30 SibSi and residual dSi 30 Si (Section 2). From our 588 LAVO data, we can get a snapshot of the magnitude of the isotope fractionation for 589 Because these functional groups can be related to specific habitats (see above), the 695 contribution of potential micro-habitat effects on diatom Si isotope fractionation for example, via micronutrient availability, growth rate, ecological interactions, etc. -697 rather than intrinsic species-specific effects cannot be strictly ruled out. 698

Si mass-balance for Manzanita Lake 699
Manzanita Lake is an efficient Si sink, removing approximately 30% of 700 inflowing dSi. This conclusion is based on three independent lines of evidence: 1) the 701 construction of elemental Si mass balance calculations (Figure 9, see below); 2) the 702 isotope ratio difference between inflows and outflows and 3) sediment accumulation 703 rates of bSi.  with the estimate based on element fluxes alone (ca. 35%, see above and Figure 9). One 714 implication is that any unaccounted for dSi inputs to the system (notably groundwater 715 or soil-interflow water, plus potential direct hydrothermal fluid recirculation, though 716 we have no evidence for this) has the same dSi concentration 30 Si as UMC waters. 717 Were this not the case, the good agreement between the two methods would be unlikely. 718 Down-core percent bSi and average mass accumulation rates (based on core 719 chronology data presented in Howard (2018)), produce an average bSi mass accumulation rate (MAR) for recent ML lake sediments of ~0.023 g/m 2 /yr (Figure 9). 721 This is approximately a factor of three higher than the estimate based on Si 722 concentrations. Two aspects may contribute to this discrepancy. Firstly, the fluvial dSi 723 flux does not represent 100% of the inputs of Si to ML. Other potential Si sources were 724 no dSi (see Figure 9). contribute Si to the system, though this is likely to be small (Frings et al., 2014). There 732 is no evidence for direct hydrothermal vent features in the lakes, though this cannot be 733 conclusively ruled out. Secondly, we do not attempt to correct for sediment focusing 734 the tendency of sediment to move to deeper areas of a lake. In other words, the bSi 735 MAR is likely overestimated here since the core derives from the lakes depocenter. 736 Work elsewhere has shown that single cores can distort estimates of mass-accumulation 737 rates (Dillon and Evans, 2001;Hilton, 1985;Likens and Davis, 1975). 738 Recent work has highlighted how not accounting for groundwater can give 739 extremely misleading or nonsensical lake Si mass-balances (Zahajská et al., 2021). One 740 benefit of silicon-isotope based estimates of lake Si retention is that they can help 741 reconcile the divergent Si burial rates derived from elemental mass-balance (0.008 742 g/m 2 /yr) with those from burial rates (0.023 g/m 2 /yr). Assuming other sources of Si are 743 limited (e.g. within-lake sediment dissolution), it suggests that groundwater is supplying about two-thirds of the lake dSi, corroborating the work of (Zahajská et al.,745 2021) who argue it should be considered in lake budgets more widely. One of the few 746 lake Si mass-balances that explicitly consider groundwater Si reported a groundwater 747 contribution of 70 % for an Argentine pampean lake (Miretzky and Cirelli, 2004). 748 Elsewhere, (Hofmann et al., 2002) have acknowledged their Si mass-balance for Lake 749 Lugano, on the Swiss/Italian border, may be biased by the lack of groundwater data. In 750 a French crater lake, (Michard et al., 1994) have shown that groundwater supplies ca. 751 90% of dSi to the system. We thus suggest that analogously to the attention submarine 752 groundwater discharge is receiving for its importance to ocean element and isotopic California. We generated data on lake characteristics, silicon isotope compositions, and 790 diatom abundances from modern and sediment core samples for Manzanita, Butte, and 791 Widow Lakes. These systems exhibit a range of diatom Si isotope fractionations, including the largest reported to date. There are strong, significant relationships 793 between diatom species composition and diatom silicon isotope composition that were 794 revealed by ecological dimension-reducing statistical approaches. Samples that are 795 dominated by specific diatom functional groups tend to cluster within a specific range 796 30 SibS values. These functional groups can be related to specific habitats. The 797 possibility of micro-habitat effects thus cannot be excluded as an explanation for 798 30 SibSi values between species, rather than species-specific Si isotope 799 fractionation factors sensu stricto. Data from LAVO lake cores suggest that diatom 800 species composition most closely covaries with 30 SibSi values over the past century. 801 While nutrient supply does play a role, the relationship between productivity and 802 30 SibSi is likely more indirect, mediated by the effects nutrient concentration and 803 stoichiometry have on the species composition of the sample. Finally, a Si mass-balance 804 for Manzanita Lake highlights the Si retention efficiency of lacustrine systems and 805 suggests an understudied role for groundwater in lake dSi supply. 806 Wes Rubio, Bud Schmidtbauer, and Sarah Luse for assisting in fieldwork. We thank 817

Acknowledgments
Nicole Fernandez and an anonymous reviewer for their constructive and valuable 818 comments on an earlier version of this manuscript, and Jérôme Gaillardet for editorial 819 handling. 820  The range of  866  30 SibSi values for each of the groups is found in table 3, as are the indicator taxa  867 associated with each group. Axis 1 explains 43% of the variance, Axis 3 explains 18% 868 of the variance. The statistical analyses indicate that approaches designed for the 869 analysis of noisy, ecological data can be used to deconvolve species-specific isotope 870 effects from bulk samples. 871   ML lake gravity core, data from Howard (2018)) is presented. Sediment-core derived 908 burial fluxes are a factor of three higher than lake mass-balance derived fluxes, which 909 may be due to groundwater dSi supply, to sediment focusing, or to a combination. 910    Table 1: Summary of Manzanita, Butte, and Widow lakes catchment characteristics and morphometry. For closed or semi-closed lakes (BL, WL), size metrics calculated based on GIS data representing full-lake conditions (see Howard and Noble 2018). Curve Number is an empirical parameter that predicts the likelihood of overland runoff for a rainfall event; Shoreline Development is shoreline length relative to a circle of the same area     Although not significant, top taxa for Group 2 are also given.