Sediment phosphorus composition controls hot spots and hot moments of internal loading 1 in a temperate reservoir 2

14 Phosphorus (P) flux across the sediment-water interface in lakes and reservoirs responds to 15 external perturbations within the context of sediment characteristics. Lentic ecosystems 16 experience profound spatiotemporal heterogeneity in the mechanisms that control sediment P 17 fluxes, likely producing hot spots and hot moments of internal loading. However, spatiotemporal 18 variation in P fluxes remains poorly quantified, particularly in the context of sediment chemistry 19 as a controlling variable. We measured P flux rates and mobile sediment P forms along the 20 longitudinal gradient of a temperate reservoir every two months from February to October of 21 2020. Both aerobic and anaerobic processes mobilized sediment P throughout the year. High flux 22 rates at littoral sampling sites (8.4 and 9.7 mg P m day) occurred in late summer under oxic 23 conditions and mobilized labile organic P. High fluxes at the profundal site coincided with 24 hypolimnetic anoxia under ice cover and in mid-summer (11.2 and 17.2 mg P m day, 25 respectively) and released redox-sensitive P. Several high fluxes substantially skewed the flux 26

rate distribution, providing evidence of hot spots and hot moments of internal loading. We 27 further analyzed the ecosystem-level effects of elevated sediment P release by scaling measured 28 flux rates to representative areas of the lakebed and estimating P loads. We found that aerobic P 29 release from littoral sites had an outsized impact on total loads. Our findings demonstrate that 30 focusing solely on flux rates without ecosystem context can lead to misidentification of the 31 dominant biogeochemical mechanisms involved and ultimately impede eutrophication interface. Sediments hold a pool of phosphorus (P), which is the legacy of past external loading 50 and subsequent sedimentation (Søndergaard and others 2003, Walsh and others 2019). This P 51 may be retained in the sediments or mobilized and released into the overlying water (i.e., internal 52 P loading; Orihel and others 2017). Within the sediments there are many different P-containing 53 minerals, organic compounds, and surface complexes, which are vulnerable to different 54 mechanisms of internal loading (North and other 2015, Orihel and others 2017). As such, the 55 chemical composition of the sediment P pool is a pivotal slow variable that shapes how the rate 56 of P flux from the sediments (i.e., the fast variable) responds to fluctuations in external drivers 57 (Carpenter 2003). 58 The response of sediment P flux rates to changing dissolved oxygen availability, the 59 result of external drivers, is shaped by the forms of P present in the sediment. For example, if 60 redox-sensitive P forms (i.e., those associated with iron or manganese oxides) dominate the 61 sediment P pool, then hypolimnetic oxygen depletion will trigger P release due to reductive 62 dissolution of the host minerals (Mortimer 1941, Jensen and Andersen 1992). When a mixing 63 event or other external disturbance delivers dissolved oxygen to the sediment surface, oxidized 64 iron and manganese minerals will stabilize P, resulting in sediment P retention. The idea that 65 oxic conditions prevent P release and that internal loading primarily occurs under anoxia is a 66 persistent paradigm in limnology. However, there is ample evidence that internal P loading 67 occurs under a range of dissolved oxygen conditions based on sediment characteristics (Hupfer 68 and Lewandowski 2008). If sediments hold a large pool of labile organic P, then oxic conditions 69 are expected to mobilize and release P via microbial decomposition and mineralization (Joshi 70 and others 2015, Horppila and others 2017). In this case, an influx of dissolved oxygen to the 71 lakebed would stimulate rather than suppress internal loading through this mechanism. The 72 path, windspeeds at the reservoir are estimated to have exceeded 65 kph, fully mixing the water 165 column at all sites (Table 1) To quantify sediment P fluxes, we collected intact sediment cores and incubated them 184 under ambient conditions in the lab while measuring P exchange with the overlying water. For 185 each sampling site and event, three replicate sediment cores were collected using a gravity corer 186 (inner diameter 5 cm, length 50 cm), such that there were approximately 25 cm of sediment and 187 25 cm of overlying water. The sediment cores and the overlying water were sealed in clear, 188 acrylic core sleeves and transported at 4°C. In the laboratory, we exposed cores to temperature 189 and dissolved oxygen treatments corresponding to ambient conditions at each site. Temperature 190 treatments were achieved by securing the cores in either a water bath or an incubation chamber. 191 Dissolved oxygen levels were manipulated by bubbling either air mixtures or N2 through the 192 overlying water. A slow, consistent bubble rate was used to gently mix the water column without 193 disturbing the sediment surface. After the cores were placed in the incubation system, we 194 measured the height of the water column within each core tube to calculate the water volume. 195 Samples of the overlying water were collected 12, 36, 60, and 84 hours after the initial 196 incubation set-up. For each daily sampling, 50 mL of water was removed for analysis of TP. An 197 equivalent volume of hypolimnetic water, collected from 0.5 m above the sediment surface at the 198 corresponding site, was used to replace the volume removed. The replacement hypolimnetic 199 water was also analyzed for TP daily to account for changes in water column P due to sampling 200 and water replacement. Samples were preserved with concentrated sulfuric acid to pH 2 and 201 stored at 4°C before undergoing persulfate digestion and analysis for TP (Standard Methods 202 4500-P B.5, E). We monitored water temperature, dissolved oxygen, and pH (YSI ProDSS 203 Multiparameter Digital Water Quality Meter) daily to ensure that the overlying water remained 204 representative of ambient conditions in the reservoir at the time of sampling. 205 The change in TP in the overlying water was used to calculate daily, areal P flux rates 206 for each core. We first calculated the mass of P in the overlying water immediately following the 207 collection of the daily water sample as well as the mass of P in the replacement water. We then 208 determined how the addition of the replacement water changed the TP concentration of the 209 overlying water. The daily change in TP concentration was calculated as the difference between 210 this new TP concentration after the addition of replacement water and the TP concentration 211 measured in the water column the next day (see Supplementary Material for equations). The P 212 flux rate was then calculated as: 213 P flux rate (mg P m -2 day -1 ) = (Ct -C0) * V / A / d (Eq. 1) 214 Where Ct is the water column TP concentration on a given day, C0 is the TP concentration from 215 the previous day after the addition of replacement water, V is the total volume of water overlying 216 the sediment core, A is the area of the sediment surface, and d is the number of days between 217 measurements (Ogdahl and others 2014). Over a 4-day incubation, we calculated three daily P 218 flux rates for each sediment core. We took the mean of these temporal replicates to yield one P 219 flux rate per core, per incubation. Mean flux rate and standard error of the mean for the three 220 replicate sediment cores from each sampling site and event were then used to estimate P flux rate 221 through time at the various sites. 222 223

Sediment P Content and Composition 224
At each sampling site and event, we collected an additional sediment core for analysis of 225 sediment P composition, total P content, and physical characteristics (see Supplementary 226 Material for sediment physical characteristics methods). We extruded the first 10 cm of the 227 sediment profile into an acrylic core sleeve, which was sealed immediately to maintain ambient 228 redox conditions. The top 10 cm of sediment is considered actively exchanging with the 229 overlying water as diffusive processes can occur in sediments this deep (Boström & Pettersson 230 1982, Forsberg 1989. Samples were transported and stored at 4°C until analysis, which began 231 within 18 to 36 hours of sample collection. Sediments were handled under N2-atmosphere in a 232 glove bag and thoroughly homogenized before removing three replicate subsamples from each 233 core. The replicates were analyzed for four mobile P species (loosely-bound, redox-sensitive, 234 aluminum-bound, and labile organic P) via sequential extraction following Lukkari and others 235 (2007). Dried sediments were used to quantify total P. 236 To begin the sequential P extractions, subsamples of fresh sediment equivalent to 0.5 g of 237 dry sediment were weighed into polyethylene centrifuge tubes. This same sediment pellet was 238 used throughout the sequential extraction procedure. All extractions were performed on an  Loosely-sorbed and pore water P were extracted in 0.46 M N2-purged sodium chloride 251 (NaCl) for one hour. One rinse in 0.46 M N2-purged NaCl was used, and the combined extract 252 solution was preserved for TP analysis (Standard Methods 4500-P B.5, E). Redox-sensitive P 253 species were extracted in a 0.11 M bicarbonate -0.1 M sodium dithionate (BD) solution for one 254 hour. This extraction included two rinses with BD solution and one NaCl rinse. The combined 255 extract supernatant was bubbled with compressed air for at least 90 minutes to remove dithionite 256 before being preserved for TP analysis. Labile organic P and P associated with aluminum oxides 257 were determined with an 18-hour extraction in 0.1 M sodium hydroxide (NaOH). One NaOH 258 rinse was used, followed by one NaCl rinse. The combined extract supernatant was analyzed in 259 two portions. First, a portion was filtered (0.45µm GF/C filters) and analyzed for SRP (Standard 260 Methods 4500-P E). This SRP concentration was used to calculate the sediment concentration of 261 aluminum-bound P. The remaining supernatant was digested and analyzed for TP. The labile 262 organic P fraction was determined as the difference between the total NaOH-extractable P and 263 the aluminum-bound P. Total sediment P concentrations were measured following a hot acid 264 digestion on dried, ground, and homogenized sediments. Three replicate subsamples (0.2 g) of 265 the dried sediment were combusted it at 550°C for 2 hours and then boiled on a digestion block 266 in 50 mL of 1 M HCl for 2 hours at 150°C. Following digestion, samples were diluted to 50 mL 267 with deionized water, and adjusted to pH 2 using 0.1 M NaOH before TP analysis. 268 269

Statistical Analyses 270
As there is no universal, quantitative method for delineating hot spots and hot moments of 271 biogeochemical fluxes in a distribution of flux measurements (Bernhardt and others 2017), we 272 used a variety of approaches to explore the distribution of P flux rates and determine whether the 273 observed spatiotemporal variation indicated hot spots or hot moments. We first identified 274 statistical outliers in the distribution, defined as flux measurements falling above or below 1.5 275 times the interquartile range. We further quantified the shape of the P flux distribution by 276 calculating skewness (m3; see Supplementary Material for equation), which identifies whether 277 the distribution is symmetric (m3 < 0.5) or if extreme flux rates, presumably due to hot spots or 278 moments, skew the distribution (m3 > 0.5). We further evaluated the influence of high flux rates 279 on the distribution by iteratively removing the highest flux rates and recalculating skewness. 280 Through this process, we determined how many of the highest flux rates would need to be 281 removed to produce a nearly symmetric distribution (Gakuruh 2017). 282 In order to evaluate how the composition of the sediment P pool varied across sites and 283 seasons, we performed a compositional data analysis. Compositional data analysis tests for a 284 difference of proportions among multivariate observations, allowing us to test differences in the 285 relative abundance of P fractions across different sites and events (Filzmoser and others 2018). 286 The compositional analysis was defined by the concentrations of loosely-bound (porewater and 287 surface sorbed), redox-sensitive (Fe-and Mn-bound), aluminum-bound, and labile organic P. 288 The sediment P concentrations were center logratio transformed prior to a principal components 289 analysis (PCA) on the covariance matrix, as not to bias the analysis to the most abundant 290 sediment P fractions. 291 To estimate the total sediment P load (kg day -1 ) in the reservoir for each sampling event, 292 we multiplied the mean P flux rate (mg m -2 day -1 ) at each sampling site by the lakebed area (m -2 ) 293 within representative depth contours corresponding to each site. The shallow site measurements 294 were assigned to the 1.2 -3.6 m depth contour, the intermediate depth site measurements were 295 assigned the 3.6 -5.6 m depth contour, and the deep site measurements were assigned the 5.6 -296 6.8 m depth contour. The threshold of 5.6 m was determined based on the historical mean depth 297 of hypoxia in the water column. We excluded areas shallower than 1.2 m, including the sediment 298 retention basins north of the reservoir, because these areas are mainly in secluded, wind-299 protected bays ( Figure 1). Our sampling sites, which were centrally located along a branch of the 300 reservoir, cannot reasonably be extrapolated to these shallow areas due to expected differences in 301 wave disturbance and the depositional environment (Kleeberg and others 2013). To calculate the 302 lakebed area within the representative depth contours, we used the length of each depth contour 303 from a bathymetric map produced by the Iowa Department of Natural Resources. The area 304 between depth contours is a trapezoid-shaped area of lakebed wrapping around the reservoir 305 basin. We calculated this area as: 306 Where length 1 and 2 are the lengths of the bounding depth contours, and h is the assumed 308 average distance between the depth contours across the lakebed, estimated using the Pythagorean 309 theorem: 310 Where depth 1 and 2 refer to the water depth of the top and bottom depth contours, and area 1 312 and 2 represent the planar areas at the top and bottom depth contours. This method assumes that 313 the lakebed follows a linear slope between the bounding depth contours, so it is likely an 314 underestimate of the true area of the sediment surface. After determining the lakebed area 315 corresponding to each sampling site, we multiplied this area by the flux rates and then summed 316 across all three load estimates for the total load for that sampling event. To propagate the 317 uncertainty of this estimate, we added the standard error of the mean values in quadrature.  rates differed between the deep site and more shallow sampling sites, and elevated P release rates 367 were observed under both oxic and anoxic conditions. 368 The shape of the distribution of sediment P flux rates over the course of 2020 provides 369 evidence of hot spots and hot moments of sediment P release in GVL. The distribution of P flux 370 rates was centered near 0 mg P m -2 day -1 with the majority of the rates falling between -10 and 371 10 mg P m -2 day -1 (Supplementary Figure S2). The first statistical moment, or mean, of the 372 distribution was 3.4 mg P m -2 day -1 . The distribution was moderately positively-skewed (third 373 standardized statistical moment m3 = 0.818) due to five high release rates (range 16.  internal loading is likely system-specific, but to assume that winter P fluxes are negligible risks 473 biasing estimates of annual internal loading. 474 Sediment core incubations are a common tool for measuring sediment P flux rates (Orihel 475 and others 2017); however, the approach also has limitations (Oghdal and others 2014). The 476 main assumption when using core incubations to quantify sediment flux dynamics is that the 477 conditions in the core are representative of the conditions in the lake. In an effort to meet this 478 assumption, we incubated cores at ambient temperature and oxygen conditions at the time of 479 collection, monitoring the temperature, dissolved oxygen, and pH daily in the cores. We also 480 used measurements from replicate cores and multiple days of incubation to estimate daily mean 481 flux rates at a site for a given sampling event in order to capture small-scale spatiotemporal 482 variability in the estimate while comparing across larger spatial and temporal scales. 483 Additionally, we limited our incubations to 3.5 days in order to minimize artifacts that can occur 484 in long-term incubations such as the depletion of organic matter or other key nutrients. Finally, 485 we compared our core incubation-based flux measurements to TP dynamics in the reservoir as 486 another way to verify that the qualitative patterns we observed in flux rates matched the changes 487 in TP concentration measured in the reservoir. These strategies combined provide confidence 488 that the broad-scale spatiotemporal patterns in sediment P flux we measured in GVL reflect the 489 dynamics occurring across sites and seasons in the ecosystem. 490 491

Sediment P Composition Controls Flux Response to Dissolved Oxygen 492
Over the course of the year, we measured internal P loading under a wide range of 493 dissolved oxygen concentrations at the sediment-water interface. High P flux rates occurred 494 under nearly anoxic conditions at the deep site and oxic conditions at the shallower sampling 495 sites. In order to understand why sediment P fluxes responded differently to dissolved oxygen 496 conditions across space and time, we measured the chemical composition of the sediment P pool. 497 We hypothesized that this slow variable shapes how fluxes respond to external drivers that alter 498 dissolved oxygen availability. Spatiotemporal variation in sediment P composition corresponded 499 to variation in sediment P flux rates due to two disparate mechanisms: oxic conditions liberated 500 P from labile organic materials and anoxic conditions mobilized redox-sensitive P species. The 501 dominance of aerobic versus anaerobic internal loading shifted over time and space in the study 502 reservoir, but both processes were important pathways for P recycling between sediments and the 503 overlying water. 504 At the deep site of the reservoir, change in redox-sensitive P concentrations mirrored P 505 flux rates over time. Redox-sensitive P declined steadily from winter to mid-summer, and then 506 sharply decreased (35% decrease) from mid-to late summer. The mid-summer sampling event 507 was a hot spot-hot moment of P release from the profundal sediments that coincided with 508 hypolimnetic anoxia at the deep site. The decline in redox-sensitive P from mid-to late summer 509 suggests that the high flux rates in mid-summer were the result of reductive dissolution of redox-510 sensitive P minerals under anoxic conditions. Concentrations of labile organic P in the sediments 511 generally increased over the course of the year at all sampling sites. However, declines were 512 measured at the deep and intermediate sites from winter to spring, suggesting that the aerobic 513 release observed in spring was due to the mineralization of labile organic P. From mid-summer 514 to late summer, labile organic P also declined at the shallow and intermediate sites, indicating 515 that the elevated rates of aerobic P release measured at these sites in late summer resulted from P 516 mineralization following decomposition of labile organic materials. 517 Our hypothesis that there are hot spots and hot moments of internal loading resulting 518 from interactions between sediment P composition and biogeochemical conditions at the 519 sediment-water interface was supported in our study reservoir. We found that anoxic conditions 520 can trigger the rapid mobilization and release of P from redox-sensitive P pools regardless of 521 water temperature. We also measured P release originating from labile organic materials under 522 aerobic conditions and found that these fluxes were enhanced following a storm disturbance. Our 523 findings underscore the importance of considering both aerobic and anaerobic pathways of 524 internal loading, especially in productive waterbodies. Our work supports the idea of a 525 "perpetual cycle of internal P loading" in hypereutrophic waterbodies, as proposed by Song and 526 Burgin (2017). The perpetual cycle describes a positive feedback loop that develops as lakes 527 become increasingly eutrophic. Increased algal production enhances inputs of detritus to the 528 sediments, producing a large pool of sediment labile organic P that is susceptible to aerobic 529 release (Baines and Pace 1994, Frost and others 2019). Sediment P release may then occur under 530 both anoxic and oxic conditions. High internal P loads sustain frequent algal blooms, the detritus 531 of which further fuels aerobic sediment P release via decomposition and mineralization. This 532 positive feedback loop likely introduces hysteresis to maintain waterbodies in a hypereutrophic 533 state. Clear evidence of sediment P release under both oxic and anoxic conditions in GVL 534 indicates that the reservoir has entered this proposed cycle in which both anaerobic and aerobic 535 internal loading will continue to fuel intense algal blooms. 536 537

Scaling Fluxes to the Whole Ecosystem 538
A hot spot-hot moment is defined as having a disproportionate influence on elemental 539 cycles at the ecosystem-scale (McClain and others 2003). As such, we scaled measured P flux 540 rates to representative areas of the lakebed to estimate P loads and determine how the fluxes we 541 classified as hot spots-hot moments actually influenced reservoir-wide internal loading. Hot 542 spots and hot moments of aerobic P release from the shallow and intermediate depth sites in late 543 summer produced the greatest total sediment P load to the reservoir over the study period. This 544 substantial P load resulted from both the high flux rates and the broad spatial extent of sediments 545 releasing P at this time. In contrast, high rates of P release from the deep site in winter and mid-546 summer did not result in large total P loads due to the small area of lakebed represented by the 547 deep site as well as sediment P retention at other sampling sites. The scaling results further 548 illustrate how even low rates of P release can result in elevated total P loads if sustained over the 549 entire lakebed. Specifically, we found that low rates of P release across all sampling sites in 550 spring resulted in a high total P load. 551 The biological relevance of these P loads varied with the timing and location of sediment 552 P release. The spring P load occurred at an ideal time and location to fuel early summer algal 553 blooms. Previous work has shown that GVL is strongly P-limited in early spring (Butts and 554 others, personal communication). Additionally, a large proportion of the lakebed releasing P at 555 this time was in direct contact with the mixed surface layer, so released P would be readily 556 available for algal uptake. The substantial P load from the shallow and intermediate sites in late 557 summer also occurred in contact with the mixed surface layer. In contrast, P loads from the deep 558 site in winter and mid-summer are less likely to have reached the euphotic zone and to be 559 available to algae (Tammeorg and others 2017). Overall, instances of aerobic sediment P release 560 resulted in the greatest total sediment P loads to the reservoir and released P at relevant times and 561 locations for algal uptake. 562 Focusing on rates alone may not be sufficient to understand how extreme values of 563 biogeochemical fluxes actually effect ecosystem structure and function. It is essential to evaluate 564 the biological relevance of observed rates and scale measurements to the whole system. At the 565 same time, scaling measurements from a single sampling station to an area of the lakebed is 566 fraught with uncertainty. However, we took a conservative approach in assigning representative 567 areas of the lakebed. We also took care in accurately describing the reservoir basin geometry and 568 calculating sediment surface area. As such, we have produced conservative estimates of sediment 569 P loads that can be used to evaluate the ecosystem effects of spatiotemporal variation in P fluxes. 570 Other studies have scaled discrete measurements of sediment P fluxes to larger areas of the 571 lakebed, generally based on waterbody surface area (Scicluna and others 2015, Noffke and 572 others 2016). Our approach to characterizing basin geometry allows for more accurate estimates 573 of sediment surface area and thus scaling at finer spatial resolutions based on water depth.  9.0 ± 5.0 1.9 ± 0.2 13.8 ± 0.9 1.9 ± 3.0 -1.   Table S1 details 807 surface and bottom water chemistry including total, volatile, and involatile suspended solids; 808 total and soluble reactive P; total nitrogen; and nitrate concentrations. Table S2 reports the 809 mean P flux rates by site and day of year, as presented in Figure 2 in the main manuscript. Table S3 demonstrates the effects of excluding high flux rates on the skewness of the P flux 811 distribution, as quantified by the standardized third statistical moment. We included a 812 visualization that compares 2020 epi-and hypolimnetic total P and soluble reactive P 813 concentrations to the P dynamics from previous years ( Figure S1). Figure S2 illustrates the 814 distribution of sediment P flux rates across sampling sites and seasons, with statistical outliers 815 noted. Figure S3 breaks down the sediment P flux rate data by site and sampling event to 816 identify the points in space and time associated with high flux rates. Figure S4

860
Where Wt is the weight of the aluminum weigh boat, Ww is the weight of the weigh boat and 861 fresh sediment sample, and Wd is the weight of the weigh boat and dry sediment.

864
Where Wa is the weight of the weigh boat and the ashed sediment after combustion.