Integrating Geochemical and Geophysical Datasets at the Three Sisters Volcanic Complex, OR, USA by Annika Dechert A dissertation accepted and approved in partial fulfillment of the requirements for the degree of Doctor of Philosophy In Earth Sciences Dissertation Committee: Josef Dufek, Chairperson Nathan Andersen, Core Member Katharine Cashman, Core Member James Watkins, Core Member John Christian, Institutional Representative University of Oregon Spring 2025 2 © 2025 Annika Dechert All rights reserved. This work, including text and images of this document, is licensed under a Creative Commons Attribution-Non Commercial 4.0 International License 3 DISSERTATION ABSTRACT Annika Dechert Doctor of Philosophy in Earth Sciences Title: Integrating Geochemical and Geophysical Datasets at the Three Sisters Volcanic Complex, OR, USA Modern volcanology aims to comprehensively understand volcanic eruption processes and mitigate risks for stakeholders to minimize loss of life. Volcanology has advanced significantly in recent decades, improving our understanding of past and present magmatic systems due to technological advances, precision in laboratories, increased computing power, and more experimentation. However, the community advocates for more interdisciplinary volcanology projects, urging geochemists, geophysicists, physical volcanologists, and experimental petrologists to collaborate in investigating the inner workings of volcanoes. My Ph.D. research enters this arena, exploring the magmatic evolution of the Three Sisters volcanic complex in the central Oregon Cascades. First, I examine magmatic storage and timescales at South Sister through zircon geochronology. This work establishes that the two oldest rhyolites are over 10 kyr younger than previously recorded, hypothesizes that the late Pleistocene and Holocene rhyolite sequences originate from distinct reservoirs, and indicates that the increased rate of rhyolite eruption during the late Pleistocene should inform future volcanic hazard monitoring and preparation at South Sister. Second, I investigate questions of magmatic evolution related to the modern-day system of the Three Sisters using a field Bouguer gravity survey and forward gravity model through magma dynamics simulations. This work suggests a magma flux range of 10-3 – 10-1 m3/m2/yr to recreate the past eruptive volumes at the Three Sisters. It highlights that the magmatic system does not have a significant density contrast with the crust, limiting the Bouguer gravity anomaly. Finally, I employ a zircon growth model to investigate the zircon record in magmatic simulations corresponding to the eruptive history of the Three Sisters. These simulations yield a zircon age record that resembles the natural record, suggesting a single magma body with cooling edges and eruptible melt pockets. Overall, this interdisciplinary study of the Three Sisters magma system updates the eruption timeline, indicates a lack of significant mobile melt in the shallow crust, and supports a focused magma 4 system for zircon crystal growth. Enhancing our understanding of this system aids hazard analysis and monitoring of this high-threat volcanic complex. This dissertation includes previously published and unpublished coauthored material. 5 CURRICULUM VITAE NAME OF AUTHOR: Annika Dechert GRADUATE AND UNDERGRADUATE SCHOOLS ATTENDED: University of Oregon, Eugene, OR, USA Occidental College, Los Angeles, CA, USA DEGREES AWARDED: Doctor of Philosophy, Earth Science, 2025, University of Oregon Bachelor of Arts, Geology, 2019, Occidental College AREAS OF SPECIAL INTEREST: volcanology geochemistry geophysics PROFESSIONAL EXPERIENCE: Graduate Employee, University of Oregon, 2019-2025 National Science Foundation INTERN, U.S. Geological Survey Cascades Volcano Observatory, 2024 GRANTS, AWARDS, AND HONORS: National Science Foundation Non-Academic Research Internships for Graduate Students (INTERN), Forward Geophysical Modeling of the Three Sisters Volcanic Complex, OR, USA, 2024. Earth Science Scholarship Award, Springfield Thunderegg Rock Club, 2024. Department of Earth Science Good Citizen Award, University of Oregon, 2023. Graduate Student Research Grant, Gravity modeling at the Three Sisters Volcanic Complex, OR, USA, Mazamas Organization, 2022. 6 National Science Foundation Graduate Research Fellowship Honorable Mention, Numerical modeling of zircon petrochronology in the context of magma dynamics, 2021. Student Research Grant, Central Oregon Geoscience Society. 2021 Student Research Grant, Community Foundation for Southwest Washington USGS Volcano Science Center Jack Kleinman Grant for Volcano Research, Numerical Modeling of Zircons at the Three Sisters, OR, 2020. University of Oregon Graduate School First Year Fellowship, University of Oregon, 2019-2020. PUBLICATIONS: Dechert, A. E., Andersen, N. L., Dufek, J., & Jilly, C. E. (2024). Zircon constraints on the eruptive sequence and magma evolution of rhyolites at South Sister volcano, Oregon. Geochemistry, Geophysics, Geosystems, 25, e2024GC011680. https://doi.org/10.1029/ 2024GC011680 Andersen, N. L., Dechert, A., Ruth, D.C.S., Sas, M., Chouinard, J., & Dufek J. (in press). The transition from melt accumulation to eruption triggering recorded by orthopyroxene Fe- Mg diffusion timescales in late Holocene rhyolites, South Sister volcano, Oregon Cascade Range. Geochemistry, Geophysics, Geosystems. Le Mével, H., Andersen, N.L., Dechert, A.E., Dufek, J. (submitted). The magmatic-hydrothermal system of the Three Sisters volcanic cluster, Oregon, images from field gravity measurements. Journal of Geophysical Research – Solid Earth. 7 ACKNOWLEDGMENTS To say it takes a village is an understatement. So many people have helped support me through this Ph.D. process, and I am forever grateful. Thank you to my parents, Jerry and Becky Dechert, for encouraging me to chase my dreams and putting up with my requests to visit volcanoes during family vacations. A big shout-out to friends who have cheered me on for the last few years, including Kat Donovan, Ellie Amann, Kara Guttridge, Olava Goerges, Nicole Abib, Lissie Conners, Yara Rossi, Kate Scholz, and Christina Cauley. I so appreciate all the coffee breaks, dog walks, movie nights, and camping trips that have made my grad school years even more special. Thank you to the Dufek lab groupmates, including Allison Kubo, Paul Regensburger, Johan (Yoshi) Gilchrist, Chris Harper, Gabe Eggers, Erin Fitch, and PJ Zrelak, for humoring me when I bring rock samples to lab group meetings and helping me debug code. Finally, thank you to all the friends, family, graduate students, research staff, and professors who have encouraged me through this endeavor. Thank you to Marla Trox and the other department staff who helped me navigate the university system. I appreciate the Department of Earth Science at the University of Oregon for supporting my research, conference attendance, DEI efforts, and field work. I acknowledge the various funding sources of this Ph.D., including the National Science Foundation, Mazamas Organization, Central Oregon Geoscience Society, Thunderegg Rock Club, and Community Foundation for Southwest Washington. Additionally, thank you to the Volcano Science Center administration and Cascades Volcano Observatory staff who supported my year at CVO. I am especially grateful to my mentors who have helped me develop as a scientist. Thank you, Margi Rusmore, for instilling a sense of excitement, curiosity, and drive for perfection during my undergrad days at Occidental College. Thank you to my Ph.D. committee for helping me achieve this degree: Josef Dufek, Nathan Andersen, Katharine Cashman, James Watkins, and John Christian. Thank you, Kathy, for always happily looking at one more figure, providing fantastic research and career advice, and being my University of Oregon volleyball buddy. Thank you, Nathan, for sharing your expertise, pushing my writing to improve, and putting up with all my talking during the miles we’ve shared in the field. Finally, thank you, Joe, for teaching me the world of numerical modeling, supporting my efforts for community outreach, 8 and always answering my (many) questions. I appreciate every red mark on my drafts, celebration for the little victories, and questions posed to make me think deeper. Thank you! 9 TABLE OF CONTENTS Chapter Page I. CHAPTER 1: INTRODUCTION ............................................................................. 17 1.1 Goals of Modern Volcanology ........................................................................ 17 1.2 Three Sisters Volcanic Complex ..................................................................... 18 1.2.1 Indigenous Perspectives ......................................................................... 18 1.2.2 Eruptive History ...................................................................................... 18 1.2.3 Modern Unrest ........................................................................................ 19 1.2.4 Unanswered Questions ............................................................................ 20 1.3 Geochronology and Geochemistry ................................................................. 22 1.4 Geophysics ...................................................................................................... 23 1.5 An Interdisciplinary Approach ....................................................................... 24 II. CHAPTER 2: ZIRCON CONSTRAINTS ON THE ERUPTIVE SEQUENCE AND MAGMA EVOLUTION OF RHYOLITES AT SOUTH SISTER VOLCANO, OR ......................................................................................................... 25 2.1 Introduction ...................................................................................................... 25 2.2 Geologic Background ...................................................................................... 27 2.3 Samples and Methods ..................................................................................... 29 2.4 Results .............................................................................................................. 32 2.4.1 Zircon 230Th-238U Dates .......................................................................... 32 2.4.2 Compositional Trends ............................................................................. 34 2.4.3 Zircon Temperatures .............................................................................. 34 2.5 Discussion ........................................................................................................ 36 10 2.5.1 Revisited Eruption Chronology and Eruption Rate ................................ 36 2.5.1.1 Updated Eruption Chronology ....................................................... 36 2.5.1.2 Implications of Updated Chronology ............................................. 37 2.5.2 Trace Element Fingerprints of Magma Reservoir Structure .................. 39 2.5.3 Integrated Perspective of the Late Pleistocene Magma System ............. 42 2.5.3.1 Longevity and Structure ................................................................ 42 2.5.3.2 Regional Connections to the Magma System ................................ 43 2.5.4 Comparing the Late Pleistocene and Holocene Magma Systems .......... 43 2.5.5 Holocene Implications ............................................................................ 45 2.6 Conclusions ..................................................................................................... 46 2.7 Bridge .............................................................................................................. 47 III. MAGMA DYNAMICS AND ERUPTION MODELING LINKED TO THE MODERN BOUGUER GRAVITY ANOMALY REVEALS MAGMATIC EVOLUTION AND STRUCTURE OF THE THREE SISTERS VOLCANIC COMPLEX, OR, USA ................................................................................................ 48 3.1 Introduction ...................................................................................................... 48 3.2 Geologic Background ..................................................................................... 49 3.3 Methods ........................................................................................................... 54 3.3.1 Bouguer Gravity Survey ........................................................................ 54 3.3.2 Magma Dynamics Model ....................................................................... 54 3.3.3 Eruption Modeling and Comparison to Three Sisters Eruptive History 56 3.3.4 Background Crust and Intrusions Compositions ................................... 56 3.3.5 Modeling Space Explored ...................................................................... 60 11 3.3.6 Forward Geophysical Model .................................................................. 61 3.3.6.1 Gravity .......................................................................................... 62 3.3.6.2 Conductivity .................................................................................. 63 3.3.6.2 Seismic Velocities ......................................................................... 63 3.4 Results .............................................................................................................. 64 3.4.1 Example Realization .............................................................................. 64 3.4.2 Model Sensitivity for Forward Gravity Models ..................................... 67 3.4.3 Plutonic Models ..................................................................................... 68 3.4.3.1 Gravity .......................................................................................... 68 3.4.4 Effusive Models ..................................................................................... 69 3.4.4.1 Gravity ........................................................................................... 69 3.4.4.2 Eruptive History ............................................................................ 69 3.4.5 Explosive Models ................................................................................... 70 3.4.5.1 Gravity .......................................................................................... 70 3.4.5.2 Eruptive History ............................................................................ 71 3.5 Discussion ....................................................................................................... 72 3.5.1 Good Fits – Eruptive Volumes .............................................................. 72 3.5.2 Best Fits – Eruptive Volumes ................................................................ 74 3.5.3 Structure of the Magmatic System ......................................................... 76 3.5.4 Magmatic Flux ....................................................................................... 78 3.5.5 Past, Present, and Future Gravity Anomaly Patterns ............................. 82 3.5.6 Seismic and Resistivity Models ............................................................. 84 3.5.7 Comparison to Modern Inflation ........................................................... 85 12 3.6 Conclusions ........................................................................................................... 85 3.7 Bridge .................................................................................................................... 86 IV. MATCHING MODELED ZIRCONS TO NATURAL SAMPLES REVEALING MAGMA STORAGE AT SOUTH SISTER VOLCANO .................. 87 4.1 Introduction ...................................................................................................... 87 4.2 Methods............................................................................................................ 89 4.2.1 Magma Dynamics Model ....................................................................... 89 4.2.2 Background Crust .................................................................................. 91 4.2.3 Zircon Growth Model ............................................................................ 92 4.2.4 Modeling Space Explored ...................................................................... 94 4.2.5 South Sister Zircon Data ........................................................................ 95 4.3 Results ............................................................................................................. 96 4.3.1 Zircon Ages ............................................................................................ 96 4.3.2 Zircon Sizes ........................................................................................... 99 4.3.3 Zircon Trace Elements ........................................................................... 101 4.3.4 Zircon Saturation Temperature .............................................................. 103 4.3.5 Zircon Survivability ............................................................................... 103 4.4 Discussion ....................................................................................................... 104 4.4.1 Zircon Modeling Matches South Sister Age Data ................................. 104 4.4.2 Zircon Sizes and Model Assumptions ................................................... 105 4.4.3 Modeled Zircon Trace Element Compositions Differ from South Sister Data .............................................................................................. 106 4.4.4 Evolution and Structure of the Magmatic System ................................. 108 4.5 Conclusions ...................................................................................................... 110 13 V. CHAPTER 5: CONCLUSIONS ............................................................................ 111 APPENDICES ............................................................................................................. A. SUPPLEMENTAL MATERIAL FOR CHAPTER 2 ....................................... 114 B. SUPPLEMENTAL MATERIAL FOR CHAPTER 3 ....................................... 133 C. SUPPLEMENTAL MATERIAL FOR CHAPTER 4 ....................................... 158 REFERENCES CITED ................................................................................................ 163 14 LIST OF FIGURES Figure Page 1. Map of the Cascade Volcanic Arc ........................................................................ 17 2. Map showing the age and SiO2 trends of eruptions in the Three Sisters ............... 21 3. Simplified geologic map of the South Sister rhyolites .......................................... 27 4. Rank order plot of zircon crystallization ages ...................................................... 33 5. Probability density functions of the zircon populations ....................................... 33 6. Variation of selected zircon trace elements and ratios .......................................... 35 7. Comparison of South Sister zircon crystallization ages to zircon trace element ratios ....................................................................................................................... 36 8. Age and volume plot of South Sister for all erupted lava types ........................... 38 9. Comparison of the trace element composition and ratios for the main Holocene rhyolite lava, the rrm mini dome, and the late Pleistocene rhyolites .................... 44 10. Residual Bouguer gravity anomaly data (mGal) ................................................... 51 11. Density versus depth profile for the modeled background crust .......................... 59 12. Schematic setup for the intrusions of magma into the modeling space ................ 60 13. Example magma dynamics model outputs ........................................................... 65 14. Forward gravity model outputs ............................................................................. 66 15. Eruptive volumes during the four eruptive episodes ............................................ 67 16. Eruption volume misfit (km3) compared to the gravity misfit (mGal) ................. 70 17. Gravity and geology results for the 40 realizations that fit at least one modeled eruptive volume episodes ........................................................................ 73 18. Gravity and geology results for the 17 realizations that fit two modeled eruptive volume episodes ...................................................................................... 75 19. Modeled average magma flux and modeled total eruptive volume ...................... 81 15 20. Schematic of the magma dynamics and zircon growth model ............................. 90 21. Major variables set for each of the magma dynamics realizations ...................... 95 22. Age data for South Sister zircons and for mafic model realizations ..................... 97 23. Age data for South Sister zircons and for intermediate model realizations .......... 98 24. Size data for South Sister zircons and for mafic model realizations .................... 99 25. Size data for South Sister zircons and for intermediate model realizations ......... 100 26. Ti ppm data for South Sister zircons and for mafic model realizations ................ 101 27. Ti ppm data for South Sister zircons and for intermediate model realizations ..... 102 16 LIST OF TABLES Table Page 1. Table of the compositions for the modeled background crust and magmatic simulations ............................................................................................................ 58 17 CHAPTER I CHAPTER 1: INTRODUCTION 1.1 Goals of Modern Volcanology The US contains 161 potentially active volcanoes, which have experienced over 120 eruptions since 1980. Eighteen of these volcanoes are classified as “very high threat” by the U.S. Geological Survey (USGS) (Ewert et al., 2018). One of the challenges posed to volcanologists is “to document and understand the repose, unrest, precursors, and timing of eruptions during the entire life cycle of volcanoes” to help mitigate volcanic hazards (Committee on Improving Understanding of Volcanic Eruptions et al., 2017). In just ten years, from 1980 to 1990, 620,000 people were directly impacted by volcanic eruptions, and by 2000, at least 500 million people lived within a volcanic risk zone (Tilling and Lipman, 1993; Chester et al., 2000). With potential eruption hazards of pyroclastic density currents, emission of gas, ash plumes, volcanic floods, and lava flows, investigating past volcanic eruptions to inform future eruption hazards is essential. This includes studying eruptive products (e.g., lava compositions, tephra dispersal, dating methods), magmatic storage conditions of the past (e.g., pressure, depth, temperature, composition, volatile content, longevity), and geophysical imaging of the modern system (e.g., magnetotellurics, seismic velocities, gravity anomalies) to inform volcanic hazard and risk assessments. Inherently, all methods ask and answer distinct questions, but powerful results follow combining methods for a more holistic understanding of a volcano. This dissertation applies interdisciplinary approaches to a “very high threat” volcanic complex: the Three Sisters (Figure 1). 40°N 42°N 44°N 46°N 48°N 126°W 124°W 122°W 120°W USGS Threat Level Very High High Low USGS Threat Level Very High High Low Figure 1: Map of the Cascade Volcanic Arc with volcano symbology representing the USGS threat level following Ewert et al. (2018). 18 1.2 Three Sisters Volcanic Complex 1.2.1 Indigenous Perspectives Multiple indigenous groups resided in the central Oregon region, including the Confederated Tribes of Siletz Indians, Confederated Tribes of Warm Springs, Kalapuya, and Confederated Tribes of Grand Ronde (Native Land Digital, 2023). The Confederated Tribes of the Warm Springs have passed down stories and traditions about the peaks of their area—tales of mountains who were once people, fighting over jealousy or carrying food on their journeys. One story archived by Clark (1953) from The Confederated Tribes of the Warm Springs reveals the character of the Three Sisters, without needing to sample any rocks or take any measurements: Klah Klahnee, the Three Sisters, was once the biggest and highest mountain of all; it could be seen for many miles. One time, the earth shook for days, and the mountain boiled inside. It boiled over, and hot rocks came out of the top of it. Flames and smoke rose high in the air. Red-hot stones were thrown out in every direction. Many villages and many Indians were buried by the rocks. When the mountain became quiet again, most of it was gone. Only three points were left. That is why it is called Klah Klahnee, for it means “three points.” You can still see the black rocks all around the base of the three mountains. This vivid story records the past eruptive hazards posed in the region and discusses how volcanic eruptions directly impacted the Warm Spring tribe. I acknowledge the people of this region and appreciate their vivid description of the land I study through this dissertation. While I will refer to these peaks as the Three Sisters, I appreciate that the original name, Klah Klahnee, or “three peaks,” is an even more descriptive and meaningful name for the volcanoes. 1.2.2 Eruptive History There has long been interest in mapping and describing the geology of the Three Sisters region. This begins with the first geological overview of the area (Hodge, 1925), followed by a more volcanology-based study (Williams, 1944) that was built upon for a more holistic study of the central Oregon region (i.e., Taylor, 1990; Sherrod et al., 2004). The Three Sisters have been the topic of multiple dissertations (such as this one) (i.e., Wozniak, 1982; Clark, 1983; Hill, 1991; Webster, 1992; Schmidt, 2005). Notably, a U.S. Geological Survey group led a mapping and dating study of the region, laying the base for future, more detailed studies of the volcanic complex (Fierstein et al., 2011; Hildreth et al., 2012; Calvert et al., 2018). 19 Despite their proximity, each of the three volcanoes in the Three Sisters shows distinct characteristics. In only 20 km running north-south in the central Oregon Cascades lie the major stratovolcanoes North Sister, Middle Sister, and South Sister. The oldest and northernmost of the three main stratovolcanoes is North Sister. North Sister is built of uniform basaltic andesite from 130 – 50 ka (Schmidt and Grunder, 2009). The uniform composition through time is formed by repeated mantle-derived injections of low-potassium tholeiitic basalt, causing low-degree melting of crustal gabbro (Schmidt and Grunder, 2011). Middle Sister is next in space and time, with largely mafic and intermediate eruptions, but also a single rhyolite building the edifice from 49 – 14 ka (Calvert et al., 2018). Finally, South Sister initiated 39 ka, building up the edifice with largely alternating rhyolite, dacite, and andesite until 22 ka (Fierstein et al., 2011; Dechert et al., 2024). After a long pause, South Sister erupted two series of rhyolites around 2 ka on the east and west flanks of the volcano (Scott, 1987). Given the rarity of rhyolite in the high Cascades, numerous studies have postulated mechanisms to produce the more felsic materials with no consensus yet, such as partial melting (Waters et al., 2021), fractional crystallization (Brophy and Dreher, 2000; Parker et al., 2023), melt contributions from dehydration melting of the crust (Hill, 1991; Fierstein et al., 2011), and rejuvenation of residual felsic magmas (Calvert et al., 2018). While such previous work sets the stage at the Three Sisters, many questions remain. 1.2.3 Modern Unrest In addition to studies of the history of the Three Sisters, efforts have been made to probe the modern-day magmatic system. Deformation to the west of South Sister was first identified by Wicks (2002) suggesting that the uplift was caused by magma intrusion or pressurization of a geothermal system at 6.5 km depth with a volume increase of 0.023 km3. Following studies used InSAR, non-continuous GPS, and leveling methods to explore the deformation that started between 1995 and 1996 (Dzurisin et al., 2006, 2009; Riddick and Schmidt, 2011; Rodríguez- Molina et al., 2021; Lisowski et al., 2021). These studies all have variations in the depths, volumes, and source hypotheses around the inflation. However, the most updated time series suggests that the total uplift of 30 cm (by 2020) is sourced from 5.9 km below the surface and is caused by either 1) persistent magma input causing the surrounding crust to respond elastically to the source-volume increase or 2) an initial injection of magma before 2020 after which the surrounding viscoelastic shell has relaxed (Lisowski et al., 2021). Utilizing another method of geophysical monitoring, Zurek et al. (2012) conducted a microgravity survey over the deforming 20 area from 2002 to 2009, but did not find evidence of any mass changes, suggesting the ongoing deformation was in response to a past intrusion. To more completely image the hydrothermal or magmatic system at the Three Sisters, Le Mével et al. (submitted) conducted a Bouguer gravity survey, identifying a hydrothermal system of 5 km3 down to 3 km depth, but found no evidence of an eruptible magma body. Instead, they hypothesize that there is a felsic mush with 15-55% partial melt and that the mid-1990s intrusion intersected the base of the mush magmatic system. Besides geodetic monitoring, spring geochemistry in the Three Sisters also indicates a magmatic source. In the deformation area, low-temperature springs display magmatic carbon, mantle-like helium isotope ratios, and magmatic heat (Evans et al., 2004; Ingebritsen et al., 1994). However, given a chloride anomaly in the area that predates the deformation, Evans et al. (2004) suggest that these geochemical anomalies originate from an intrusion that pre-dates the mid-1990s intrusion. 1.2.4 Unanswered Questions Even with the fantastic previous research on the Three Sisters volcanic complex, I identify three main groups of unanswered questions: the relationship between the three main volcanoes, the Cascades rhyolite anomaly, and the current magmatic system regarding forecasting future eruptions. Not only does the Three Sisters Reach of the Cascades arc have one of the highest vent densities with 466 Quaternary volcanoes, three major stratovolcanoes, plus Broken Top, all lie in a 20 km arc segment, which is an anomaly in the arc (Hildreth, 2007). However, each volcano is distinct from the others. Open questions include: What controls the SiO2 progression to become more felsic to the south (Figure 2a)? Why do we see an age progressing from oldest in the north to youngest in the south, but also see the older Broken Top (300-150 ka) to the east (Figure 2b)? What causes the shift to high eruptive productivity from 50-20 ka, producing coeval Middle and South Sister? Finally, does a single, connected magmatic system form all three volcanoes, or does each system have a physically distinct source? Through this dissertation, I will address the eruptive sequence of South Sister in Chapter 2, creating a clearer picture of the age progression 21 of the system. In Chapter 3, I will address the connectivity of a single or multiple distributed intrusions sourcing the eruptions. Outside the rear arc, the Three Sisters and Lassen are the only locations with true high- silica rhyolite in the Cascades (Hildreth, 2007). In fact, the scarcity of rhyolites in the front arc settings is true worldwide. Middle Sister has one rhyolite, South Sister erupted 13, and just to the east of the Three Sisters lies a “silicic highland” of the Tumalo Volcanic Field with numerous rhyolites (Sherrod et al., 2004). As Hildreth (2007) says, “there are many roads to rhyolite,” and those roads have been hypothesized at South Sister (Hill, 1991; Brophy and Dreher, 2000; Waters et al., 2021; Parker et al., 2023). However, the root cause of why rhyolites can form here and not elsewhere along the arc is vastly missing. Some studies will suggest that central Oregon is the leading edge of the High Lava Plains rhyolites associated with the presence of the Yellowstone Plume (MacLeod et al., 1975; Jordan et al., 2004; Fierstein et al., 2011). Still, there seems to be no effort to connect that statement to the region's geochemistry or tectonics. As such, I believe the root cause of rhyolites in the Three Sisters region is a significant, outstanding question. Finally, I believe we have a vague understanding of the current magmatic system beneath the Three Sisters. The geodetic evidence suggests a magmatic intrusion at 6 km depth, but neither a microgravity nor a Bouguer gravity survey has been able to successfully image such a magma body (Zurek et al., 2012; Le Mével et al., submitted). The null hypothesis of Le Mével et Figure 2: Map showing the age and SiO2 trends of eruptions in the Three Sisters. a) SiO2 contents of eruptions from the Three Sisters. Note the more felsic eruptions to the south. b) Ages of eruptions from the Three Sisters. Black dots denote ages >50 ka. 22 al. (submitted) is an excellent step to understanding the system, highlighting that there is no perceptible, eruptible magma of large volumes in the shallow crust. Additionally, in Chapter 3, I show that even while the Three Sisters are in the highest eruptive state, the possible gravity anomaly due to the magma system is not very different from the modern signal, suggesting a magma system with a density close to the background crust. If this is true, and there is a magma body between 3-6 km with 15-55% partial melt, what sort of recharge event would create an eruptible magma body? More importantly, how would we monitor for a recharge event large enough to remobilize such a magma body (i.e., Cooper and Kent, 2014)? Such questions and answers directly impact the setup of a monitoring network and the interpretation of future unrest. 1.3 Geochronology and Geochemistry One method regularly utilized to probe the past of volcanic eruptions is a geochronology perspective, or the clocks of volcanoes. In short, we can use the decay of radioactive isotopes in minerals to measure time. The first date calculation through this method is attributed to Bertram Boltwood in 1907, who used the U/Pb series to date a uranium ore sample (White, 2015). Decay systems commonly used in geologic settings include the K-Ar-Ca, Rb-Sr, Sm-Nd, Lu-Hf, and U- Th-Pb systems (White, 2015). A common dating method in volcanology is the 40Ar/39Ar method, and in unaltered volcanic rocks, the date is interpreted as the eruption age (Schaen et al., 2021a). This is the primary dating method for recording the eruptive history of the Three Sisters (Fierstein et al., 2011; Calvert et al., 2018). However, a single date from a mineral is not necessarily helpful in discerning the evolution of or processes within the magmatic system. This led to the term petrochronology, or the study connecting time (i.e., ages) with rock and mineral-forming processes and physical conditions (Engi et al., 2017). Zircon is a prime mineral for petrochronology studies, with time from the U-Th-Pb series, temperature from the Ti-in-zircon calculation, and magmatic processes recorded in the crystal's trace elements and rare earth elements (REE) (Watson and Harrison, 1983; Finch and Hanchar, 2003; Hoskin and Schaltegger, 2003; Engi et al., 2017). Zircon studies have focused on large-volume systems (Reid et al., 1997; Simon and Reid, 2005; Rubin et al., 2017) and small-volume systems (Claiborne et al., 2010; Stelten and Cooper, 2012), uncovering the thermal, chemical, spatial, and temporal evolution of magmatic systems (Cooper, 2015). The petrochronologic perspective has revealed magmatic processes such as magma stalling and rejuvenation at Mount St Helens (i.e., Claiborne et al., 2010), independent evolution patterns rather than a long-lived magma body at 23 Long Valley caldera (i.e., Simon and Reid, 2005), near-solidus crystal storage with only rapid heating events at the Okataina Volcanic Center (i.e., Rubin et al., 2017), and heterogeneous hot and cold magma storage zones of the magmatic system at Laguna del Maule (i.e., Andersen et al., 2019). I utilize this advantageous mineral to explore South Sister volcano’s eruptive history and magmatic evolution in Chapter 2. 1.4 Geophysics While looking to the past of volcanoes through geochronology is useful, another method focuses on probing the magmatic system of the present day: geophysics. Rather than using past eruptive products to understand a volcano, geophysical methods utilize various wave-matter or field-matter methods to image physical differences in the subsurface, identifying magma or hydrothermal fluids beneath active volcanoes. These methods typically use the inverse method, which uses the data collected to infer the subsurface (Blakely, 1995). However, this generally leads to non-unique solutions of the best Earth system model (Magee et al., 2018). However, when used in parallel, these methods are informative about modern magmatic systems. Handy methods include gravity surveys, seismic imaging, and electromagnetic methods. Gravity surveys measure the relative gravitational field at a volcano, and relate the gravity anomaly to the distribution of mass in the subsurface, finding anomalously low-density regions that can be interpreted as melt or fluids (Seigel, 1995). Such a survey identified low-density material ranging from 8-30 km depth in the Campania Active Volcanic Area (i.e., Campi Flegrei, Vesuvius), suggesting long-term magmatic injection and crystallization (Fedi et al., 2018). Seismic imaging techniques measure seismic wave speeds to identify changes in wave arrival times, as S-waves cannot travel through fluids, which also slow down P-waves (Gudmundsson, 2012). A recent effort to image the magmatic systems along the entire Cascade Arc analyzed receiver functions to identify P-to-S wave scattering, finding partial melt under Lassen Peak, Crater Lake, Newberry Volcano, Mount Hood, Mount St Helens, and Mount Rainier (Pang et al., 2025). Electromagnetic techniques measure the resistivity or conductivity variations of the crust to identify higher conductivity melt bodies (Stanley et al., 1990; Magee et al., 2018). In the Altiplano-Puna volcanic complex of the central Andes, a magnetotelluric survey imaged low resistivity melts at 15 km depth and shallow melts, highlighting recent intrusions (Comeau et al., 2015). Geophysical studies of modern magmatic systems provide essential data to support future monitoring efforts. The volume, depth, and melt fraction of an active magmatic system impact 24 how monitoring networks are built and control the types and magnitude of signals that unrest could produce. I compare the data from a gravity survey of the Three Sisters to simulated gravity anomalies to understand magmatic evolution in Chapter 3. 1.5 An Interdisciplinary Approach Volcano science is inherently interdisciplinary. One of the main ways to continue strengthening volcanic science is through more multidisciplinary approaches, “integrating diverse types of observations from geophysics, geology, geochemistry, geodynamics, and remote sensing” (Committee on Improving Understanding of Volcanic Eruptions et al., 2017). Successful programs moving towards such endeavors include the Geodynamic Processes at Rifting and Subducting Margins (GeoPRISIMS) Program and the Subduction Zones in four Dimensions (SZ4D), which can help support projects such as imaging Magma Under St. Helens (iMUSH) to uncover the architecture of the Mount St. Helens magmatic system. The utility of such studies is shown at Kīlauea volcano, Hawaiʻi, where olivine diffusion timescales and seismic swarms align, connecting real-time monitoring data with timescales of magmatic processes (Lynn et al., 2024). Additionally, an effort to combine petrologic, geophysical, stratigraphic, and modeling approaches displays a fuller picture of the 79 CE Plinian eruption of Vesuvius, Italy, impacts on the surrounding communities, and how such an eruption would impact society today (Doronzo et al., 2022). I believe multidisciplinary approaches are essential for a complete understanding of volcanic processes and specific volcanic edifices. As such, this dissertation considers geochronology and geophysical data, examining the magmatic system at the Three Sisters volcanic complex through an interdisciplinary lens. 25 CHAPTER II CHAPTER 2: ZIRCON CONSTRAINTS ON THE ERUPTIVE SEQUENCE AND MAGMA EVOLUTION OF RHYOLITES AT SOUTH SISTER VOLCANO, OREGON From Dechert, A. E., Andersen, N. L., Dufek, J., & Jilly, C. E. (2024). Zircon constraints on the eruptive sequence and magma evolution of rhyolites at South Sister volcano, Oregon. Geochemistry, Geophysics, Geosystems, 25, e2024GC011680. https://doi.org/10.1029/2024GC011680 2.1 Introduction Globally, small volume but frequent explosive volcanic eruptions impact populations and aviation due to eruptive plumes, ashfall, lava flows, and lahars. Therefore, continuing interdisciplinary research at active volcanoes will help inform hazard assessments and mitigation strategies for future eruptions (e.g., Ewert et al., 2018). Geologic and petrologic approaches offer complimentary views of past eruptions and their implications for future behavior. On one hand, geologic studies yield information about the frequency, magnitude, and style of past activity to assess plausible future eruption scenarios. On the other, petrologic methods provide understanding of how the magmatic system has evolved and may offer clues about which past eruptions are more likely to be representative of the modern system and future behavior. The mineral zircon provides coupled temporal and compositional information that can inform both the timing of past eruptions and the evolution of the underlying magmatic system. The utility of zircon is due to its durability, high compatibility with elements useful for geochronology and petrology (e.g., U, Th, and the rare earth elements), and a well-characterized temperature dependence of saturation (Watson and Harrison, 1983; Finch and Hanchar, 2003; Hoskin and Schaltegger, 2003). Zircon petrochronology has been applied to large eruptive volume systems (Reid et al., 1997; Simon and Reid, 2005; Rubin et al., 2017) as well as small to medium eruptive volume systems (e.g., Claiborne et al., 2010) to distinguish chemical, thermal, spatial, and temporal characteristics of melt evolution and processes (Cooper, 2015). In the Cascade Range, zircon studies have informed models of magma evolution and storage at Mount St. Helens (e.g., Claiborne et al., 2010), South Sister (e.g., Stelten & Cooper, 2012), Crater Lake (e.g., Bacon & Lowenstern, 2005), and Lassen Peak (e.g., Klemetti & Clynne, 2014). However, the dearth of rhyolite and lower U in these arc magmas has, in part, hindered broader application of zircon geochronology in the Cascade Range (Hildreth, 2007). The previous zircon study at 26 South Sister focused only on the most recent eruptions (Stelten and Cooper, 2012); however, the rhyolites that erupted in the late Pleistocene remain little investigated. The U.S. Geological Survey has classified the Three Sisters volcanic complex as a “very high threat potential” to the surrounding communities, such as Bend (pop. 100,000), only 50 km away (Ewert et al., 2018). Past eruptions included hazards such as pyroclastic density currents, lava flows, and ash plumes (Scott et al., 2001). South Sister volcano has erupted approximately 12 km3 of volcanic rocks over its lifetime, including 13 rhyolitic eruptions ranging from 51 to 2 ka (Fierstein et al., 2011; Hildreth et al., 2012). Frequent eruption of rhyolite and its modest eruptive volume makes South Sister an ideal system to investigate the origin and pre-eruptive storage of rhyolite in smaller volume magmatic system. In addition to the recent eruptions, South Sister has experienced ground inflation since the mid-1990s focused 5 km west of the summit and only ~3 km from the late Holocene Rhyolite of Rock Mesa (rrm) vent (Lisowski et al., 2021). InSAR, continuous GPS, leveling, and modeling studies proposed multiple mechanisms of this uplift including the continuing intrusion of basaltic magma (hydraulic), a delayed response of visco-elastic crust to a past magma intrusion (viscoelastic), or the pressurizing of a geothermal system (poroelastic) (Wicks, 2002; Dzurisin et al., 2009; Riddick and Schmidt, 2011). These studies suggest the source of deformation is approximately 5-7 km deep. Additionally, a study of springs near South Sister found He and C isotope ratios and elevated heat flow consistent with contributions from mafic magma (Evans et al., 2004). As such, the working hypothesis suggests an intrusion of basaltic magma 5-7 km below South Sister. This intrusion could lead to a future eruption at South Sister, and understanding how the past magmatic system compares to the modern system can help inform eruption and hazard forecasting. We present 230Th-238U crystallization ages and trace element compositions of zircon crystals from rhyolites at South Sister volcano. Results of this study suggest (a) South Sister hosted a long-lived magmatic system primed for zircon crystallization from 100 to 20 ka, (b) the earliest rhyolite eruptions are more than 10 kyr younger than past studies suggested, (c) South Sister was characterized by a high eruption rate from 40-30 ka, and (d) the Pleistocene and Holocene rhyolites are derived from chemically distinct sources. 27 2.2 Geologic Background The Three Sisters volcanic complex comprises three stratovolcanoes (North Sister, Middle Sister, and South Sister) located in the central Oregon Cascade Range, USA (Figure 3). The regional tectonics of this reach of the Cascade Range is complicated by the interaction between the warm Cascadia subduction zone, the past influence of the Yellowstone plume, the edge of the Basin and Range extension, and may contribute to the unusual abundance of high- SiO2 volcanism for this region of the Cascade Range (Obrebski et al., 2010; McCrory et al., 2012; Cheng et al., 2017). The oldest and northernmost of the three, North Sister, erupted basaltic andesite with major cone construction from approximately 120-45 ka (Schmidt and Grunder, 2009; Hildreth et al., 2012). Growth of Middle Sister occurred from approximately 49- 15 ka with compositions ranging from basaltic andesite to rhyolite (Calvert et al., 2018). Building of the South Sister edifice alternated between rhyolite and intermediate products from ~51-22 ka, largely coevally with Middle Sister (Fierstein et al., 2011). The growth of all the Three Sisters stratovolcanoes was preceded by the growth of Broken Top stratovolcano, located just 5 km southeast of South Sister. Broken Top erupted 7-10 km3 of basaltic andesite to rhyodacite from approximately 300-150 ka (Hildreth, 2007). In total, the Three Sisters reach of the Cascade Range extends 90 km N-S in Central Oregon, and has at least 466 Quaternary volcanoes (Hildreth, 2007). Postglacial mafic eruptions are more common in this area than rdc interior (SC) rrm interior (SC) rrm interior (minidome) rct surface rgl interior rgl surface rdh surface rse surface rmc surface rsw interior rsw surface Summits ESRI World Hillshade Rhyolite of Devil's Chain (rdc) Rhyolite of Rock Mesa (rrm) Summits ESRI World Hillshade Rhyolite of South Sister climbers trail (rct) Rhyolite of Separation Creek (rsc) Rhyolite of Kaleetan Butte (rkb) Rhyolite of Green Lakes (rgl) Rhyolite of Park Creek (rpc) Rhyolite of Prouty Glacier (rpg) Rhyolite of Devils Hill (rdh) Rhyolite southest of Lewis Glacier (rse) Rhyolite of South Fork (rsf) Rhyolite of Mesa Creek (rmc) Rhyolite southwest of Lewis Glacier (rsw) Summits ESRI World Hillshade Zircon Samples South Sister Rhyolites Holocene Late Pleistocene a b Figure 3: a) Simplified geologic map of the South Sister rhyolites (modified from Hildreth et al., 2012) with locations of samples discussed in this paper. Brown triangles show South Sister, Middle Sister, North Sister, and Broken Top summits. Cyan and magenta triangles denote the samples of Stelten and Cooper (2012). Squares symbolize zircon interior analysis and circles symbolize zircon surface analyses for samples from this study. Colors signify samples from the same eruptive unit. See Dechert (2024) or supplemental data for sample coordinates. b) Inset shows the location of the Three Sisters in the Oregon Cascade Range, other notable Cascade Range volcanoes, and major cities. 28 elsewhere in the Cascade Range. These postglacial volcanic centers include the Sand Mountain volcanic field, Blue Lake volcano, LeConte Crater, Fish Lake lava from Nash Crater, Lost Lake cones, Four In One Cone, Yapoah Crater, Collier Cone, and Belknap Crater (Sherrod et al., 2004). Slightly farther south lies Mount Bachelor, which is the largest of approximately 50 N-S aligned basalt to basaltic andesite vents that erupted from 18-8 ka (Gardner, 1994). The eruption chronology at Three Sisters has largely been established by 40Ar/39Ar dates (Fierstein et al., 2011; Hildreth et al., 2012; Calvert et al., 2018). However, Fierstein et al. (2011) used the Taylor Creek sanidine standard age of 27.87 Ma while Calvert et al. (2018) used a more recently determined, slightly older age for the Taylor Creek sanidine standard of 28.345 Ma (Fleck et al., 2019). To ensure that all the 40Ar/39Ar dates are reported on the same basis and are directly comparable, we have recalculated the 40Ar/39Ar dates using the ArAR software of Mercer & Hodges (Mercer and Hodges, 2016) to an updated Taylor Creek age of 28.447 Ma (equivalent to a Fish Canyon Sanidine age of 28.201 Ma; (Fleck et al., 2019; Schaen et al., 2021a)) and decay constants of Min et al. (2000). For the remainder of this paper, we will report 40Ar/39Ar ages for South Sister based on this recalculation (Figure A1). South Sister’s eruptive history highlights a range of compositions and one of the rare instances of rhyolite in the Cascade Range (Hildreth, 2007). The oldest exposed lava at South Sister is rhyolitic at 51.4 ± 9.7 ka and eruptions alternated between rhyolitic, dacitic, and andesitic through ~37 ka (Fierstein et al., 2011). Then, a broad cone composed of andesitic to dacitic lavas grew from 37 to 30 ka. Activity became concentrated at the summit cone by ~27 ka culminating with the final and most mafic South Sister eruption at ~22 ka, which formed the present-day summit crater (Fierstein et al., 2011). Finally, two rhyolitic eruptions occurred at approximately 2 ka on the SW and E flanks of the volcano (Scott, 1987). Previous studies have hypothesized that the mechanism for rhyolite production at South Sister could include a) partial melting from different crustal compositions, b) fractional crystallization from mafic precursor magma, or c) the Holocene rhyolites could reflect the rejuvenation of residual rhyolitic magmas related to the late Pleistocene eruptions (Hill, 1991; Brophy and Dreher, 2000; Fierstein et al., 2011; Calvert et al., 2018; Waters et al., 2021; Parker et al., 2023). Hypotheses to explain the compositional gaps in the suite of Three Sisters eruptive products include fractional crystallization (Brophy & Dreher, 2000; Dufek & Bachmann, 2010) and partial melting (Waters et al., 2021); however, these gaps are not as dramatic at South Sister when considering the more 29 comprehensive data set of Hildreth et al. (2012). These hypotheses highlight the importance of further study as the mechanisms of silicic magma genesis at South Sister and the connection between the late Pleistocene and Holocene rhyolite eruptive episodes remain an open question. Zircon and plagioclase geochronology and trace element chemistry of the two Holocene rhyolite eruptions suggests two distinct magma sources. Rhyolite of Rock Mesa (rrm) erupted at 2,150 ± 150 14C yr B.P. (preferred 14C age of Scott, 1987) with a main coulee, pumice-fall deposits, and satellite domes on the SW flank of South Sister (Scott, 1987; Hildreth et al., 2012). Rhyolite of Devils Chain (RDC) erupted at 2,000 ± 200 14C yr B.P. in a series of seven lava flows and domes aligned N-S on the southeast flank of South Sister (Scott, 1987; Hildreth et al., 2012). Here, we note different naming of units between sources but consider them superseded by Hildreth et al. (2012). The zircon 230Th-238U crystal ages are predominantly 80-20 ka for both Holocene eruptions with a small number of grains in secular equilibrium (>350 ka) (Stelten & Cooper, 2012). The zircon crystals with resolvable dates are interpreted to be antecrystic, inherited from the longer-lived system (Stelten and Cooper, 2012). Conversely, bulk plagioclase 230Th-226Ra crystal ages trend closer to eruption age, 6.8-2.3 ka for rrm and 9.6-4.0 ka for rdc (Stelten and Cooper, 2012). The Ti-in-zircon geothermometer indicates zircon crystallization ranges from 640-837°C (Stelten and Cooper, 2012). The distinct trace element ratio of Hf/Nb of the rrm and rdc zircon crystals suggest these eruptions were sourced from two chemically and physically distinct reservoirs that coexisted temporally in the broader, shared magmatic system (Stelten and Cooper, 2012). While this work focused on the Holocene rhyolites, no such work has included the late Pleistocene South Sister rhyolites. 2.3 Samples and Methods Samples collected for this work include six late Pleistocene and one Holocene rhyolite from South Sister (Figure 3). The units span the southern spatial extent of the rhyolite eruptions and the known eruption age range. Zircon data were collected for the Rhyolite southwest of Lewis Glacier (rsw), Rhyolite of Mesa Creek (rmc), Rhyolite of southeast Lewis Glacier (rse), Rhyolite of Devils Hill (rdh), Rhyolite of Green Lakes (rgl), and Rhyolite of South Sister climbers trail (rct). These units erupted between 51.4 and 24 ka based on 40Ar/39Ar ages reported by Fierstein et al. (2011), Hildreth et al. (2012), and Calvert et al. (2018). Ages are based on incremental heating of plagioclase separates (rsw, rmc, and rct) or groundmass (rse, rdh, rgl) separates. Of the Pleistocene rhyolites, rse, rgl, rdh, and rdc are on the southeastern flank of the 30 volcano, which situates them near Broken Top. Our work also considers the Holocene rrm and rdc rhyolites erupted at 2,150 ± 150 14C yr B.P. and 2,000 ± 200 14C yr B.P. These eruption ages are based on 14C dates of wood fragments within rrm pumice fall and a soil directly underlying rdc tephra (Scott, 1987). We collected a sample of ejecta associated with the rrm satellite minidomes and compare our results to those of Stelten & Cooper (2012) samples from the main rrm coulee and two of the southern rdc domes (Figure 3). Zircon separates were produced by Zirchron LLC using standard magnetic and density separation techniques, including crushing, electro pulse disaggregation, water table, heavy liquids, and Frantz. Analysis of 238U-230Th dates and trace element compositions were conducted on the SHRIMP-RG (reverse geometry) ion microprobe co-operated by U.S. Geological Survey and Stanford University in the SUMAC facility at Stanford University during multiple sessions in 2023. Zircons were pressed into indium mounts and left unpolished. For these surface analyses, trace element data were collected followed by 230Th-238U analysis on the same spot. We also attempted to measure zircon interiors by mounting in epoxy and polishing the mount. A cathodoluminescence image of the polished interiors displays size, shape, and textures of the crystals analyzed (Figure A2). The average c-axis length is 87 µm and the average a-axis length is 49 µm. Crystals are euhedral and the most common shape is equant with additional shapes being tabular and prismatic. Approximately 30% of the crystals exhibit sector zoning, another 30% have patchy zoning, with oscillatory, homogeneous, diffuse, and irregular patterns for the rest. One sample, rrm minidome, does include resorbed cores on 6 of the 22 crystals analyzed. Unfortunately, the small size of most of the zircon crystals frequently resulted in the beam impinging on the epoxy during 230Th-238U analysis of crystal interiors yielding a significant carbide interference and unusable data. The U-Th analytical procedure and data reduction procedures follow methods developed by Williams (1997) and Ireland & Williams (2003). An O2- primary beam with accelerating voltage of 10 kV was used to sputter secondary ions from the sample surface with a 15-16 nA primary beam current focused to a ~37 µm spot. Seven masses were measured, including 90Zr216O, 238U+, 232Th12C+, 230Th16O+, background measured 0.050 AMU above the 230Th16O+ peak, 232Th16O+, and 238U16O+. Data were collected over 8 scans per spot, for a total run time of ~30 minutes, collected by magnet peak-jumping on an electron multiplier. Zircon U 31 concentration data were standardized against the well-characterized MAD-559 zircon standard (3940 ppm U; (Coble et al., 2018)). For trace element analyses, an O2- primary beam with accelerating voltage of 10 kV was used to sputter secondary ions from the sample surface with a 1.5 nA primary beam current focused to a ~16µm spot. The acquisition included analysis of 29 masses, all measured by peak jumping on an electron multiplier: 28Si16O+, 45Sc+, 48Ti+, 49Ti+, 56Fe+, 89Y+, 93Nb+, 92Zr1H+, 96Zr+, 139La+, 140Ce+, 141Pr+, 146Nd+, 147Sm+, 153Eu+, 155Gd+, 165Ho+, 159Tb16O+, 162Dy16O+, 166Er16O+, 169Tm16O+, 172Yb16O+, 175Lu16O+, 90Zr216O+, 180Hf16O+, 206Pb+, 207Pb+, 232Th16O, and 238U16O+. Count times ranged from 2–18 seconds per mass to optimize counting statistics for each isotope, for a total run time of ~13 minutes per spot. See Dechert et al. (2024) Text S1 for additional details of the SHRIMP-RG analytical methods. Zircon ages are calculated using the decay constant of (Cheng et al., 2013) and the method of Boehnke et al. (2016). The Boehnke et al. (2016) method calculates a date based on a mode equilibrium melt for each zircon spot based on the 238U/232Th measured in the zircon, an assumed U/Th partition coefficient ratio between the zircon and the melt, and assumes that melt variance from the equiline is ±15% (1s). We calculated this mean U/Th partition coefficient to be 5.4 ± 0.8 from the South Sister zircons and sample whole rock compositions (Dechert et al., 2024 Table S4). The dates reported by Stelten & Cooper (2012) for the Holocene rhyolites were calculated using a two-point model isochron approach (e.g., Reid et al., 1997) using the glass Th isotope ratios of the host rhyolites and the older decay constant of Cheng et al. (2000). To ensure our dates are comparable to those of Stelten & Cooper (2012) we first recalculated their model isochrons using the newer decay constant of Cheng et al. (2013) and find only 5% average difference with the updated decay constant after removing one crystal in secular equilibrium. We then compared these dates to those calculated from their data using the Boehnke et al. (2016) model and employed for our measurements and find a 38% difference for rdc (59% absolute difference within uncertainty) and 30% for rrm (78% absolute difference within uncertainty). Stelten & Cooper (2012) found the Holocene rhyolite whole rock compositions to be in 2% Th- excess. This contrasts to the average of secular equilibrium ±15% that is assumed in the Boehnke et al. (2016) model and is a potential source of bias. We tested the magnitude of this bias by calculating model isochron dates using a whole rock Th isotope composition based on the whole 32 rock U/Th ratio and assumptions of secular equilibrium with 2% Th excess. We find that assuming 2% Th-excess makes an average difference of 2.2% in the calculated dates, which is less than the analytical uncertainty. In the absence of whole-rock or glass Th isotope compositions for the Pleistocene rhyolites, we prefer the method of Boehnke et al. (2016) and find the assumption that the host melt is in secular equilibrium on average is simpler and no less reasonable than assuming the Th-excess measured for the Holocene rhyolites was homogeneous throughout the magma system and unchanging during the interval of zircon crystallization. All zircon dates discussed for the rest of this paper were calculated by the Boehnke et al. (2016) model without Th excess. Additionally, to also ensure that the trace element measurements from Stelten & Cooper (2012) are comparable to the late Pleistocene data, we recalculated the Holocene data with preferred concentrations of Coble et al. (2018) that intercalibrate MAD-1 (standard of Stelten & Cooper (2012)) and MAD-559 (standard of this work) zircon trace element standards. 2.4 Results 2.4.1 Zircon 230Th-238U Dates Individual 230Th-238U ages are model ages calculated using the Boehnke et al. (2016) model and represent zircon crystallization ages. Our data include 18 analyses in secular equilibrium, and thus unresolvable dates >350 ka. The 228 resolvable zircon dates range from 282.5-3.4 ka, but 199 (87%) are between 100-20 ka (Figure 4). There are three dates from rgl interiors that are nominally resolvable from 344.4-310.9 ka but are within error of secular equilibrium and therefore potentially older. The zircon dates for units rsw, rmc, rse, rdh, and rct all yield probability density maxima around 40 ka (Figure 5). Units rrm (combined) and rdc show a broad distribution with no distinct peaks. The unpolished surface analyses (rsw, rmc, rse, rdh, rgl, rct) yield a narrower range of younger dates compared to the crystal interiors (rdc, rrm coulee, rrm minidome, rgl, rsw). This is especially true for rgl interiors; however, we do not have a complete dataset of interior and surface analyses for the other units. Unit rgl shows two populations of zircon ages, one between 100-20 ka, probability density maxima of ~40 ka, like the other Pleistocene rhyolites, and one >150 ka, probability density maxima of ~250 ka (Figure 5). The coherent >150 ka population was not found in the rgl surface dates, nor the surface dates of the other rhyolites. While we only have four interior analyses for rsw, the zircon crystal date range of 70.2-27.6 ka matches that of the larger surface dataset, 96.7-25.4 ka. 33 In addition to individual zircon dates, we calculate weighted mean dates for the youngest distinguishable group of dates. We calculate this weighted mean age by including progressively older zircon crystal dates until the mean square weighted deviation (MSWD) exceeded the critical value. Given the spread of ages and the size of the dataset, we are not able to calculate a zircon weighed mean age for rdc (Figure A4). The zircon weighted mean age for the rrm main flow, rrm minidome, rct, and rgl interiors are all older than the associated 40Ar/39Ar or 14C date. The Figure 4: Rank order plot of zircon crystallization ages. Data of Stelten and Cooper (2012) symbolized as triangles and labeled SC. Zircon surface analyses are circles and interior analyses are squares. Black horizontal lines with gray fields show the associated 14C date (rdc and rrm) or 40Ar/39Ar date (all others) with 1s uncertainty (Fierstein et al., 2011; Scott, 1987). Red horizontal lines with salmon fields show the zircon weighted mean date with 1s uncertainty of the youngest group of zircons that yield an acceptable mean square weighted deviation. Note that the zircon weighted mean age for rmc and rsw are younger than the 40Ar/39Ar date at 1s . Geologic unit abbreviations are as in Figure 3. 5 20 40 80 100 25015060 Model Age (ka) Figure 5: Probability density functions of the zircon populations by unit based on model isochron slope calcaulted by the method of Boehnke et al. (2016). Surface and interior data are combined for each unit. Geologic unit abbreviations are as in Figure 3. 34 zircon weighted mean date for rgl surfaces, rdh, and rse are equivalent to the associated 40Ar/39Ar date. However, the zircon weighted mean date for rmc and rsw surface and rsw interior are younger than the associated 40Ar/39Ar dates at the 1σ level (Figure A4). The 40Ar/39Ar eruption date for rmc is 47.4 ± 8.2 ka and the zircon weighted mean date is 31.5 ± 2.1 ka (n = 30, MSWD = 1.02). The 40Ar/39Ar date for rsw is 51.4 ± 9.7 ka and the zircon weighted mean date for surfaces is 39.1 ± 2.4 ka (n = 30, MSWD = 0.89). For both rmc surfaces and rsw surfaces, the 40Ar/39Ar eruption date and the zircon weighted mean age are differentiated at 1s confidence level, such that the propagated uncertainty of the difference is less than the difference, but not at the 2s confidence level. 2.4.2 Composition Trends Nearly all the zircon crystals have U contents lower than 500 ppm (Figure S5); only 18 have higher concentrations up to 2,610 ppm. Ti concentrations cluster between 3 and 10 ppm, some reach up into low to mid 20s ppm Ti (Figure 6). Eu/Eu* is generally low, ranging from 0.02 to 0.34 (Figure 6). Outside of a few outliers, the Hf concentrations ranges between approximately 6,000 to 11,000 ppm. Zircon Hf concentrations correlate with some trace elements. We see decreasing Ti concentrations and Eu/Eu* with increasing Hf concentrations (Figure 4, Figure A5). However, elements such as Y, Sc, Nb, and ratios such as U/Th are not correlated with Hf concentrations, producing sub-vertical arrays (Figure A5). These elements and ratios that do not correlate with Hf do show some correlation with each other. For example, higher ratios of U/Th correlate with higher ratios of Hf/Nb and Yb/Gd. There is no apparent correlation between zircon age and composition for any of the Pleistocene rhyolites (Figure 7). 2.4.3 Zircon Temperatures Crystallization temperatures are calculated using zircon Ti concentrations and the Ti-in- zircon geothermometer calibration of (Ferry and Watson, 2007). This calibration requires activities of SiO2 and TiO2. However, determining these values is not straightforward for the South Sister rhyolites. We consider a range of these two activities, how those ranges impact the Ti-in-zircon temperature, and how phases in the South Sister system would impact these activities (see Dechert et al., 2024 Supplemental Text S1). The most reasonable assumptions are an activity of 1 for SiO2 and 0.7 for TiO2, and zircon temperatures presented here use these activities. The temperature range for all zircons spans from 647-855°C (Figure A6). The 35 Pleistocene eruptions all show a similar average Ti-in-zircon temperature from 724-768°C (Figure A6). Of these Pleistocene units, rgl yields the highest average while rmc and rdh show the coolest temperatures. For the Holocene rhyolites, rdc has the highest average temperature. Interestingly, rrm collected from the main coulee has a slightly higher average temperature than the sample we collected from the minidome ejecta (Figure 3). There is no correlation or trends between temperature and zircon age. In addition to calculating the Ti-in-zircon temperature, we can use the whole rock compositions to determine the temperature at which the melt would have reached zircon saturation (Boehnke et al., 2013). The range of temperatures that these melts reached saturation is 743-781°C (Figure A6). We also considered the calibration of (Watson and Harrison, 1983); 6000 7000 8000 9000 10000 11000 12000 13000 14000 Hf (ppm) 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Eu /E u* rdc interior (SC) rrm interior (SC) rrm interior (minidome) rct surface rgl interior rgl surface rdh surface rse surface rmc surface rsw interior rsw surface 1.0 1.5 2.0 2.5 3.0 3.5 U/Th 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Eu /E u* 6000 7000 8000 9000 10000 11000 12000 13000 14000 Hf (ppm) 5 10 15 20 25 Ti ( pp m ) 1.0 1.5 2.0 2.5 3.0 3.5 U/Th 0 20 40 60 80 100 Y/ Sc 6000 7000 8000 9000 10000 11000 12000 13000 14000 Hf (ppm) 1.0 1.5 2.0 2.5 3.0 3.5 U /T h 1.0 1.5 2.0 2.5 3.0 3.5 U/Th 10 20 30 40 Yb /G d 0 20 40 60 80 100 Y/Sc 0 10 20 30 N b (p pm ) 1.0 1.5 2.0 2.5 3.0 3.5 U/Th 0 10 20 30 N b (p pm ) a b c d e f g h zrc + plag xtl down-temperature zrc xtl zrc xtl titn xtl amph xtl zrc xtl amph + titn xtl zrc xtl zrc + plag xtl Figure 6: Variation of selected zircon trace elements and ratios. One anomalously high Hf/Nb ratio and Y/Sc ratio plots off scale. zrc = zircon, amph = amphibole, titn = titanite, xtl = crystalization. Geologic unit abbreviations are as in Figure 3. 36 however, the range of temperatures calculated is higher (791-824°C) and are not consistent with the Ti-in-zircon temperature onset as the saturation temperatures using Boehnke et al. (2013). We do not calculate saturation temperatures using Crisp & Berry (2022) because their experiments only go down to 750°C, which is not a low enough temperature for this system, and the model includes H2O values, which we do not have a constraint for. 2.5 Discussion 2.5.1 Revisited Eruption Chronology and Eruption Rate 2.5.1.1 Updated Eruption Chronology In addition to utilizing zircon for information about the magmatic system, zircon surface dates can represent an eruption age if zircon crystallization continued to the time of eruption, and growth was rapid enough to produce a thick enough rim that could be analyzed by ion probe (~4 µm for 230Th-230U dating). Examples include Lava Creek Tuff in Yellowstone National Park (Matthews et al., 2015) and Wilson Creek Formation of Mono Lake, California (Vazquez and Lidzbarski, 2012). Whereas a majority of the Pleistocene zircon dates are older than or indistinguishable from the associated 40Ar/39Ar date, in accordance with the typical assumptions of 40Ar/39Ar and 238U-230Th systematics, the weighted mean age of the youngest zircons from rsw and rmc are significantly younger than the corresponding 40Ar/39Ar date at the 1s level, regardless of the standard used. The 40Ar/39Ar dates for rsw and rmc are plateau ages of bulk plagioclase separates. The low K contents and low radiogenic yields make the ages susceptible to excess Ar—i.e., 40Ar 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Eu /E u* rdc interior (SC) rrm interior (SC) rrm interior (minidome) rct surface rgl interior rgl surface rdh surface rse surface rmc surface rsw interior rsw surface 0 50 100 150 200 250 300 350 Model Age (ka) 0 10 20 30 N b (p pm ) a b Figure 7: Comparison of South Sister zircon crystallization ages to zircon trace element ratios. There are no clear trends between zircon crystal age and compositions for these, or other, trace elements. Geologic unit abbreviations are as in Figure 3. 37 that is not attributable to in situ decay or atmospheric Ar (Kelley, 2002). The inverse isochron intercepts of the incremental heating experiments for rsw and rmc are within uncertainty of the 40Ar/36Ar ratio of atmospheric Ar; however, particularly the larger uncertainty of the isochron intercept for rsw may not be sufficiently sensitive to the small quantities of excess 40Ar needed to produce a spurious date due to the low 40Ar * yields of the plagioclase separates. Moreover, we note that the incremental heating spectra for rmc is somewhat “saddle shaped” (i.e., the ages of intermediate heating steps are consistently lower than the lower and high temperature heating steps), which is often interpreted as reflecting the presence of excess Ar (Kelley, 2002). It is also possible that inherited plagioclase crystals contributed to the older dates; however, we do not find consistent textural evidence that the plagioclase was inherited, and 230Th-226Ra ages for the two Holocene rhyolites suggest the plagioclase record near-eruption growth (Stelten & Cooper, 2012). While the 40Ar/39Ar plagioclase ages and 238U-230Th ages overlap at the 95% confidence level, we propose that the more precise, and younger, zircon surface weighted mean dates of 31.5 ± 2.1 and 39.1 ± 2.4 ka for rmc and rsw, respectively, represent the better determinations of the eruption age for these units. Zircon crystallization occurs in the melt body, and thus cannot occur after eruption, therefore, the zircon age must be a more accurate representation of eruption age. These younger eruption ages are consistent with field relationships. Because rws and rmc are the oldest known eruptions of South Sister, this revision of their eruption ages inverts the relative timing of the initiation of Pleistocene volcanism at South Sister and Middle Sister and compresses the overall period of cone building at South Sister (Figure 8). 2.5.1.2 Implications of Updated Eruption Chronology Early work at Three Sisters suggested that North Sister is the oldest of the three, Middle Sister started erupting next, and after Middle Sister was no longer erupting, South Sister’s eruptive sequence began (Wozniak, 1982; Clark, 1983; Hildreth et al., 2012). However, new mapping and dating by Fierstien et al. (2011), Hildreth et al. (2012), and Calvert et al. (2018) found the oldest Middle Sister lava eruption age of 48.5 ± 10 ka (mnf) and the oldest South Sister lava of rsw at 51.4 ± 9.7 ka were of similar age. They concluded that North Sister is still oldest, but that South Sister was the next edifice to be active, followed by Middle Sister as the last edifice to have its first eruption, noting that these dates are not distinguishable, and that activity at South and Middle Sisters largely spanned a similar period in the late Pleistocene. Our 38 revised eruption age for rsw, 39.1 ± 2.4 ka, now postdates the earliest known eruption from Middle Sister (Figure 8) but does not change the broader point that the Pleistocene South Sister and Middle Sister eruptions were largely contemporaneous. Products of earlier eruptions from both edifices that would further refine this timeline may have been removed by glaciation and/or buried, but our new dates indicate that Middle Sister initiated before South Sister. The revised rsw and rmc eruption ages fall in the general pause of Middle Sister activity ranging from approximately 37 to 27 ka (Fierstein et al., 2011). While Middle Sister and South Sister have overlapping eruptive periods, this shift in the South Sister eruption timeline suggests a higher productivity of silicic eruptions at South Sister while Middle Sister was eruptively quiet. The new, younger eruption ages of rsw and rmc condense the eruption sequence such that nine rhyolites erupted from 40-30 ka, rather than just seven, increasing the eruptive rate of the system significantly during this timeframe, and all the Pleistocene rhyolites erupted within 15.1 Figure 8: Age (ka) and volume (km3) plot of South Sister for all erupted lava types. The y-axis shows the cumulative South Sister lava exposed (km3) (Fierstein et al., 2011). The previously adopted timeline (Calvert et al., 2018; Fierstein et al., 2011; Hildreth et al., 2012) suggests South Sister erupted before Middle Sister initiated (gray box, 1s) and the first South Sister eruptions was 51.39 ka (red diamonds). Our zircon dates suggest that volcanism at South Sister began erupting after Middle Sister and South Sister experienced a higher rate of rhyolite eruptive productivity from 40-30 ka than previously recognized (black squares). 39 kyr rather than 27.4 kyr. To evaluate the changes in eruptive rate at South Sister due to the revised eruptive sequence, we consider the cone building phase of the late Pleistocene that concluded at ~ 22 ka (i.e., excluding the Holocene eruptions). Based on the 40Ar/39Ar dates, the eruptive rate for all lava types is 0.39 km3/kyr and for rhyolite lava is 0.17 km3/kyr over the cone building stage. However, with the updated eruption ages, the eruptive rate for all lava types is 0.67 km3/kyr and for rhyolite lava is 0.30 km3/kyr over the cone building stage, almost doubling the eruptive rate during this time interval. Calculated for all eruptions over the lifetime of the volcano, the eruptive rate of South Sister is 0.4 km3/kyr, which is a moderate rate compared to other well-characterized Cascade Range volcanoes: Mount Baker is 0.1-0.2 km3/kyr, Mount St. Helens is 2 km3/kyr, Mount Adams is 0.4-0.6 km3/kyr, and Mount Shasta is 0.75 km3/kyr (Hildreth, 2007). 2.5.2 Trace Element Fingerprints of Magma Reservoir Structure Correlations among zircon trace element compositions are commonly interpreted to reflect the compositional evolution of the host melt resulting from the crystallization or destabilization of zircon or other co-existing phases. The Hf composition is a common index of zircon crystallization, whereby an increase in Hf suggests progressive zircon crystallization (Hoskin & Schaltegger, 2003; Claiborne et al., 2006). Likewise, the Eu anomaly is used as a differentiation index for plagioclase, where a decrease in the Eu/Eu* ratio tracks a melt crystallizing plagioclase, and thus the fractional crystallization of the melt more generally. The partitioning of Ti in a zircon is temperature sensitive and as such, the concentration of Ti in zircon is used to understand the temperature of the system (Ferry and Watson, 2007). Additionally, U has a higher partition coefficient than the still compatible Th, and with high partition coefficients for the heavy rare earth elements (HREE) in zircon, progressive zircon crystallization should also decrease the U/Th ratio and the abundance of HREE relative to MREE and LREE (Hoskin and Schaltegger, 2003; Claiborne et al., 2018). The decrease in Ti and Eu/Eu* with increasing Hf concentrations for the Pleistocene rhyolite’s zircon data is consistent with down-temperature co-crystallization of zircon and plagioclase (Figure 6a, 6c). In contrast to Ti and Eu/Eu*, we do not see a correlation between Hf and the U/Th or rare earth elements (REE) ratios that would be expected if zircon crystallization was the principal process affecting these ratios (Figure 6e). Moreover, we do not see a zircon crystallization trend in the generally flat U/Th versus Eu/Eu* trend (Figure 6b). However, the 40 U/Th ratio and REE ratios such as Yb/Gd are correlated to each other, but not with the indication of zircon crystallization or general crystallization, suggesting that another process exerts greater leverage on the evolution of these elements. Processes that could also influence these concentrations outside of zircon crystallization include the crystallization or destabilization of phases such as amphibole, clinopyroxene, titanite, and/or Fe-Ti oxides. Broadly speaking, clinopyroxene and Fe-Ti oxides do not incorporate U, Th, or the REE in high enough concentrations to have a large impact on their evolution. Importantly we note that KDU/KDTh < 1 and KDMREE/KDHREE > 1 in titanite and amphibole while KDU/KDTh > 1 and KDMREE/KDHREE < 1 in zircon, so crystallization of these phases drives the U/Th and Yb/Gd ratios of the residual melt in opposite directions. The evolution produced by U/Th versus Yb/Gd could be caused by either zircon or titanite/amphibole crystallization. However, the lack of correlation with Hf, the flat U/Th versus Eu/Eu* trend, and the negative correlation between U/Th versus Y/Sc (which is the opposite trend if zircon crystallization were controlling the evolution of these elements (Bachmann et al., 2005)) suggests that zircon is not leveraging control of the evolution of these elements. Additionally, we see higher Nb concentrations and Y/Sc ratios at low U/Th ratios that indicate zircon fractionation, but this is not consistent with the whole rock compositional evolution (Hildreth et al., 2012). These data trends rule out zircon crystallization as the main influence on these trace elements so instead, we propose that the trends in U/Th and the REE are consistent with the involvement of amphibole and titanite (Figure 6f). While the trace element compositions of the surfaces are largely indistinguishable, some zircon interiors are marked by higher concentrations of Nb, Ti, and Sc. Most Nb concentrations are <5 ppm, Ti concentrations are ~2-10 ppm, and Sc concentrations are ~50-100 ppm (Figure 6, Figure A5). However, the zircon interiors for rgl and rsw reach higher concentrations of 15, 26, and 150 ppm for Nb, Ti, and Sc, respectively. Comparing Y/Sc versus Nb shows two distinct trends: one with lower Nb values and one with higher Nb values (Figure 6g). Nb can be moderately incompatible or compatible in zircon and amphibole, but is always compatible in titanite (Berlo et al., 2004; Bachmann et al., 2005; Tiepolo et al., 2007). Due to these differences in the affinity for Nb, we interpret the higher Nb trend to be due to titanite destabilization, and the lower Nb trend to be caused by amphibole destabilization (Bachmann et al., 2005). The high- Nb titanite signature is found more often in zircon interiors. 41 As such, we suggest that the destabilization of amphibole and titanite in response to magma recharge deeper in the magmatic system is the principal driver of the U, Th, and REE ratios thereby decoupling them from zircon and plagioclase co-crystallization. Very rare amphibole is present in the Holocene lavas, and its cryptic involvement in the Pleistocene lavas is inferred from trace element trends and is a common feature of intermediate to silicic arc magmatic systems more generally (Davidson et al., 2007). While titanite has not been identified in the South Sister lavas, cryptic titanite destabilization in the magma could still be contributing to these zircon compositions. Titanite is a common late-stage accessory mineral in felsic calc- alkaline plutons, thus, its presence in the cooler, more crystallized mush would be not surprising (Kohn, 2017). Various vocabulary and schematics have been postulated to describe the anatomy of the magmatic system beneath a volcano informed by both geophysical and geochemical evidence (e.g., Annen, 2009; Bachmann & Bergantz, 2004; Bacon & Druitt, 1988; Charlier et al., 2005; Cooper & Kent, 2014; Hildreth, 1981, 2004; Huber et al., 2009; Marsh, 1981). Here, we propose a framework involving a mush zone and melt-rich holding chambers. We consider a mush zone located in the upper crust where silicic melts focus and crystalize to a high-crystal content, the crystallinity, longevity, and temperature of which can vary (Hildreth, 2004; Dufek and Bachmann, 2010; Kent and Cooper, 2018). The melt-rich holding chambers represent the melt- segregation from the mush that can more readily be erupted (Bachmann and Bergantz, 2004; Hildreth, 2004). It has been suggested that there is magma distributed throughout the upper crust beneath South Sister (e.g., Brophy & Dreher, 2000; Hill, 1991) as well as in the lower crust (e.g., Hildreth, 2007; Hildreth et al., 2012). The presence of rare amphibole in the Holocene rhyolites, but not in the late Pleistocene rhyolites, suggests some level of interconnectivity between a slightly deeper melt body and the shallow, more evolved mush that is tapped for rhyolitic eruption. The common mush zone fingerprints the indistinguishable trace element trends across the Pleistocene eruptions, and interconnectivity of the magma system provides transport of amphibole and titanite into the crystal mush (Cashman et al., 2017). Periodic intrusions needed to maintain melt in the upper crust (Dufek and Bergantz, 2005; Karakas and Dufek, 2015) may periodically destabilize amphibole and titanite as is recorded by the trace element trends. The resulting melts are then assimilated into the larger magma reservoir. The locus of amphibole and titanite destabilization need not be coincident in space or time with the zircon crystallization, but 42 apparently was volumetrically significant to dominate the REE and U/Th budget of the magma system. 2.5.3 Integrated Perspective of the Late Pleistocene Magma System 2.5.3.1 Longevity and Structure In addition to the trace element fingerprints of the structure of the magma system, we can utilize the zircon dates to consider the longevity of the Pleistocene magma system. The continuity of zircon crystal dates for all Pleistocene rhyolites falling from 100-20 ka suggests a long-lived mush zone of silicic magmas and zircon crystallization. We suggest that this mush zone allowed for the destabilization of amphibole and titanite, catalyzed by new arrivals of more mafic magma from depth of the interconnected magmatic system and melt-rich holding zones of more readily eruptible material. These holding chambers could be distinct, one getting tapped for each Pleistocene eruption, or more interconnected and tapped over the Pleistocene rhyolite eruptions. The clustering of zircon surface dates around the eruption ages suggests that the magma body remained at zircon saturation composition and temperatures up to eruption, or long enough to grow the ~4 µm zircon depth to be analyzed, suggesting hundreds to thousands of years of storage based on zircon growth and dissolution kinetics (Watson, 1996). This common mush zone is also supported by the lack of temperature, temporal, or spatial trends in the zircon crystal ages and compositions. We can distinguish neither individual Pleistocene rhyolites nor geographic clusters of eruptions, for example, those east versus west of the summit, based on their zircon trace element compositions. While we do see a range of Eu/Eu* and Yb/Gd ratios, there is no systematic increase or decrease through time to suggest major melt body processes or evolution in the system. Additionally, using the Ti-in-zircon thermometer, we see no thermal correlation with zircon model age. The shortened interval of rhyolite eruptions in the Pleistocene suggested in this work and a lack of temporal or spatial trace element trends supports a single reservoir that was repeatedly tapped over the eruption sequence. Conversely, multiple physically distinct reservoirs could have been tapped and emptied via a single eruption. The zircon data cannot distinguish between a melt-rich or a mushy extraction source. No matter the process of eruption, the mush zone was likely a dynamic system that was continually replenished via the interconnected system but remained constant chemically through the Pleistocene rhyolite eruption sequence. 43 2.5.3.2 Regional Connections to the Magma System In addition to the proposed magma system beneath South Sister, there is evidence of xenocrystic zircons in one eruptive unit from a nearby volcano. While most of the zircon crystal populations range from 100-20 ka, Rhyolite of Green Lakes (rgl) includes an older subset of zircon ages ranging from 150 ka to secular equilibrium (>350 ka). The extent of the exposed rgl is on the eastern flank of South Sister with proximity to Broken Top volcano (Figure 1). While other lavas, such as rse and rdc, are also on the eastern flank, rgl has the most complete dataset including both surface and interior zircon analyses. Broken Top eruptive activity lasted from 300-150 ka and compositionally ranges from basalt to rhyodacite (Hildreth, 2007). The correlation between the older rgl zircon age population and the Broken Top eruptions suggests that either: (a) rgl inherited xenocrysts from the solidified Broken Top magma system or (b) the Broken Top magma system remained active for more than 100 ky following the final eruptions from Broken Top and also fed the rgl eruption. Scenario (a) is favored due to the bimodal age distribution of interior zircon model ages (Figures 4 and 5), whereas scenario (b) would have produced a more continuous age distribution from secular equilibrium through to 20 ka. Currently, Broken Top is an understudied volcano. Further mapping, geochronology, and petrologic investigation at Broken Top would allow an improved understanding of its magmatic system and support comparison to that of South Sister. Other volcanic centers in the area such as Todd Lake volcano, Tam McArthur Rim, and the Tumalo volcanic complex could also be the source of these older zircon crystals; however, the overlapping ages and proximity suggest Broken Top is the more likely candidate (Hill, 1991; Hildreth, 2007; Klemetti et al., 2023). 2.5.4 Comparing the Late Pleistocene and Holocene Magma Systems In addition to the Pleistocene dataset collected for this project, we can include comparison of the Holocene dataset of Stelten & Cooper (2012) to infer the state of the magmatic system over time. As discussed above, the Pleistocene rhyolites largely cannot be distinguished from one other based on zircon compositions, except for the higher Nb interiors of rgl and rsw. However, zircons from the Holocene rhyolites show a distinction in compositions when compared to those from the Pleistocene rhyolites (Figure A7). The highest Nb concentrations occur in the Holocene zircons and the Holocene zircons have predominately lower Eu/Eu* (Figure 9). 44 This distinction of trace element concentrations between the two main Holocene lavas and the Pleistocene lavas suggests that the mush zone and melt-rich bodies that fed these two eruption sequences are distinct. The overlap in zircon crystal dates from all units implies that the mush zones that fed these eruptions were largely coeval in the crust and likely relatively restricted in spatial area of the system. As such, it is likely that the melt source for the two main Holocene lavas was present during the Pleistocene eruptive sequence but was not tapped for an eruption. These could be two physically distinct melt-rich bodies (e.g., Stelten & Cooper, 2012) or a single body that was tapped for both Holocene eruptions. Additionally, compared to the late Pleistocene zircon crystals, the Holocene samples tend to have higher Nb compositions, reflecting likely more titanate dissolution (Figure 9c). The erupted Holocene magmas may have been derived from a slightly more fractionated source than the late Pleistocene magmas. Such constraints could also suggest that the late 6000 7000 8000 9000 10000 11000 12000 13000 14000 Hf (ppm) 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Eu /E u* rdc interior (SC) rrm interior (SC) rrm interior (minidome) rct surface rgl interior rgl surface rdh surface rse surface rmc surface rsw interior rsw surface 1.0 1.5 2.0 2.5 3.0 3.5 U/Th 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Eu /E u* 10 20 30 40 50 60 70 80 90 Y/Sc 0 10 20 30 N b (p pm ) a b c Figure 9: Comparison of the trace element compositions and ratios for the main Holocene rhyolite lava (data from Stelten and Cooper, 2012) (red), the rrm mini dome (pink), and the Pleistocene rhyolites (blue). The zircons from the Holocene and Pleistocene lavas can be distinguished based on their trace element compositions, except for the rrm minidome (pink squares), which displays some chemical affinity to both groups, see text. Geologic unit abbreviations are as in Figure 3. 45 Pleistocene rhyolites reflect a deeper part of the magmatic system relative to the Holocene rhyolites. While the compositional distinction between the Pleistocene eruptions and the two main Holocene eruptions is clear, the rrm minidome complicates the story. The rrm minidome is located just east of the main rrm coulee vent and shares compositional trends with both the main Holocene and late Pleistocene groups of lavas. In U/Th, Nb, and U, the rrm minidome compositions matches that of the Holocene group (Figure A7). In Eu/Eu* and Ti, the rrm minidome has similar values to the Pleistocene rhyolite zircons. However, there is also considerable overlap between the rrm minidome with both the late Pleistocene and Holocene groups in other compositions. While we do not have enough data to make a concrete conclusion, causation for these trends include cooling along the edge of the rrm reservoir that erupted as the rrm minidome or the interaction of the rrm reservoir with a remnant of the Pleistocene system. 2.5.5 Holocene Implications Not only does this work help us understand the evolution and structure of the past rhyolitic magmatic system at South Sister, but it can give us clues to how a modern rhyolitic magma system could evolve. Specifically, we highlight how the Pleistocene rhyolite eruptive sequence could be repeated in the Holocene. The Pleistocene rhyolite magmas all originated from the common South Sister magmatic system, evolved in a long-lived partially crystallized mush, and erupted from melt-rich region(s) for each eruption. The updated South Sister eruption chronology suggests a higher rhyolitic eruptive rate of 0.30 km3/kyr during the 39-22 ka cone building stage. Notably, the recurrence interval of rhyolitic eruptions from 40-30 ka is just over a thousand years, during which time nine of the 11 Pleistocene rhyolites erupted. The inflation at South Sister since the mid-1990s and the spring geochemistry suggest that magmatism continues beneath South Sister (Evans et al., 2004; Lisowski et al., 2021). Therefore, it must be considered if the Holocene eruptions are the beginning of a cluster of rhyolite eruptions similar to the Pleistocene episode. Continuing and additional efforts, such as additional monitoring and geophysical imaging to identify and characterize an active reservoir that could feed multiple future eruptions is currently beneath South Sister, would help explain if the late Pleistocene episode is a relevant guide to future volcanism. If so, this would suggest the eruptive volume, sequences, hazards, and hazard zones from the 40-30 ka eruptions would 46 become an important reference for mitigation and monitoring of future volcanic hazards at South Sister. 2.6 Conclusions 1. The two oldest South Sister rhyolites (rsw and rmc) are 12.3 and 15.9 kyr younger than previously thought, making South Sister the most recent of the Three Sisters to initiate eruptive activity. The resulting shorter duration of the edifice building stage at South Sister almost doubles the eruptive rate during this stage. 2. Zircon compositions record magma evolution by down-temperature plagioclase and zircon co-crystallization. However, they also record the influence of destabilized amphibole and titanite in the magma system. The destabilization of these phases exerted compositional leverage on U, Th, and REE composition of the melt and decoupled these elements from the down-temperature zircon and plagioclase crystallization recorded by Hf, Ti, and Eu/Eu*. 3. The Pleistocene rhyolites were derived from a shared crystal mush that imprinted the compositional similarities of all lavas. This mush remained at zircon saturation composition and temperature for hundreds to thousands of years for zircon rim crystallization prior to eruption of the melt-rich region(s). 4. The main flows and domes of the Holocene rhyolites have zircons with similar ages but distinctly lower Eu/Eu* ratios that zircons from the Pleistocene rhyolites, suggesting they crystallized from melts that underwent more plagioclase fractionation than the Pleistocene rhyolites. As such, these Holocene eruptions come from a separate reservoir or mush body than the Pleistocene eruptions. 5. The Holocene rhyolite eruptive sequence is distinct from the Pleistocene rhyolite eruptive sequence and thus could follow a unique eruptive future. Conversely, the Pleistocene rhyolite eruptive sequence could be repeated in the Holocene. If volcanic history repeats itself, we can expect a rhyolitic eruption recurrence interval of just over a thousand years and an eruptive rate of 0.30 km3/kyr. The Pleistocene rhyolite eruptive sequence is an important reference for evaluating future hazards at South Sister while acknowledging South Sister could be in a new, rhyolitic eruption sequence. 47 2.7 Bridge In Chapter 2, I utilized a petrochronology lens to explore rhyolites at South Sister volcano. This not only revealed updated eruption ages for two units, making South Sister the most