MAGMA ACCUMULATION BENEATH SANTORINI VOLCANO FROM P-WAVE TOMOGRAPHY by BRENNAH GENE MCVEY A THESIS Presented to the Department of Earth Sciences and the Graduate School of the University of Oregon in partial fulfillment of the requirements for the degree of Master of Science June 2019 THESIS APPROVAL PAGE Student: Brennah Gene McVey Title: Magma Accumulation beneath Santorini Volcano from P-wave Tomography This thesis has been accepted and approved in partial fulfillment of the requirements for the Master of Science degree in the Department of Earth Sciences by: Emilie Hooft Chairperson Amanda Thomas Member Thomas Giachetti Member and Janet Woodruff-Borden Vice Provost and Dean of the Graduate School Original approval signatures are on file with the University of Oregon Graduate School. Degree awarded June 2019 ii © 2019 Brennah Gene McVey iii THESIS ABSTRACT Brennah Gene McVey Master of Science Department of Earth Sciences June 2019 Title: Magma Accumulation beneath Santorini Volcano from P-wave Tomography Despite multidisciplinary evidence for crustal magma accumulation below Santorini Volcano, the structure of the shallow magmatic system remains elusive. We use tomographic inversions of P-wave, active-source seismic data to identify a pronounced low-velocity anomaly, 2.8-5.4 km below the northern caldera basin, which causes seismic attenuation and ray-bending. After extensive synthetic testing, we constrain the shape, volume, and melt-content of shallow, inter-caldera magma storage. Low-velocities are elongated ~15 km NE-SW, attributed to a thin mush region (6-12% melt) confined by tectono-magmatic lineaments, which extends from an inter-caldera magma chamber (15- 25% melt). The main magma chamber is consistent with depth estimates of pre-eruptive storage and a recent inflation episode. We suggest magma accumulation under Santorini is controlled by crustal extension and local edifice loads. Our results provide the first structural constraints of Santorini’s upper crustal magma system during inter-caldera periods. This thesis includes previously unpublished co-authored material. iv CURRICULUM VITAE NAME OF AUTHOR: Brennah Gene McVey GRADUATE AND UNDERGRADUATE SCHOOLS ATTENDED: University of Oregon, Eugene, OR Colorado School of Mines, Golden, CO DEGREES AWARDED: Master of Science, Earth Sciences, 2019, University of Oregon Bachelor of Science, Geophysical Engineering, 2017, Colorado School of Mines AREAS OF SPECIAL INTEREST: Seismology Geophysics Natural Hazards PROFESSIONAL EXPERIENCE: Graduate Employee, University of Oregon, 2017-2019 Intern, United States Geological Survey, 2016-2017 Intern, Johnson & Prewitt Engineering, 2014 GRANTS, AWARDS, AND HONORS: Graduate Recruitment Award, University of Oregon, 2017 First Year Fellowship, University of Oregon, 2017 President’s Scholarship, Colorado School of Mines, 2013-2017 Magna Cum Laude, Colorado School of Mines, 2017 v PUBLICATIONS: Allstadt, K. E., McVey, B. G., & Malone, S. D. (2017). Seismogenic landslides, debris flows, and outburst floods in the western United States and Canada from 1977 to 2017. U.S. Geological Survey data release, https://doi.org/10.5066/F7251H3W. vi ACKNOWLEDGMENTS I wish to thank my wonderful advisor Emilie Hooft and my committee members Amanda Thomas and Thomas Giachetti, as well as Doug Toomey and Benjamin Heath for insightful discussion of the results. This work was funded by NSF grant OCE - #1459794 awarded to Emilie Hooft and Doug Toomey. vii TABLE OF CONTENTS Chapter Page I. OVERVIEW ................................................................................................. 1 II. MAGMA ACCUMULATION BENEATH SANTORINI VOLCANO FROM P-WAVE TOMOGRAPHY ............................................................. 2 APPENDIX ...................................................................................................... 18 REFERENCES CITED .................................................................................... 32 viii LIST OF FIGURES Figure Page 1. PROTEUS experiment geometry ................................................................ 6 2. Caldera-crossing seismic data ..................................................................... 8 3. P-wave velocity model ............................................................................. 10 4. Conceptual model of shallow magmatic structure ..................................... 15 S1. All raypaths added for this study .............................................................. 23 S2. P-wave velocity anomalies ....................................................................... 24 S3. Physical property calculation curves ......................................................... 25 S4. Synthetic test of Caldera Low Velocity (CLV) ......................................... 26 S5. Synthetic test of thin CLV ........................................................................ 27 S6. Synthetic test of deep CLV ....................................................................... 28 S7. Synthetic test of Linear Low Velocity (LLV) ........................................... 29 S8. Synthetic test of preferred model .............................................................. 30 S9. Past eruption centers controlled by tectono-magmatic lineaments ............. 31 ix CHAPTER I OVERVIEW The body of this thesis is original research prepared for submission to Geology with Dr. Emilie Hooft, Benjamin Heath, Dr. Doug Toomey, Dr. Michele Paulatto, Dr. Paraskevi Nomikou, Dr. Costas Papazachos, Dr. Joanna Morgan, and Dr. Michael Warner as coauthors. All authors discussed the results and their implications and contributed editorial support. E. H. and D. T. supervised the project. B. H. provided starting tomographic models and seismic data picks for the shallow velocity structure of the study region (1-3 km depth). M.P. provided python codes for physical property calculations. D. T. provided tomographic analysis software. E. H., D. T., J. M., C. P., P. N., and M. W. designed, funded, permitted, and executed the data collection. The remaining work presented is my own, including data picking and tomographic analysis to constrain the velocity structure at depths detailed in this manuscript (3-6 km), physical property calculations, writing the manuscript and drafting all figures. 1 CHAPTER II MAGMA ACCUMULATION BENEATH SANTORINI VOLCANO FROM P-WAVE TOMOGRAPHY Material in this chapter is being prepared for publication with Dr. Emilie Hooft, Benjamin Heath, Dr. Doug Toomey, Dr. Michele Paulatto, Dr. Paraskevi Nomikou, Dr. Costas Papazachos, Dr. Joanna Morgan, and Dr. Michael Warner. I. Summary Shallow magma chambers regulate volcanic eruptions and represent the topmost accumulation of melt within a transcrustal magmatic system. After caldera collapse, the formation and eruptibility of a small crustal chamber depends on melt injections from below. Santorini Volcano is a caldera in the Hellenic subduction zone, characterized by alternating effusive, shield-building eruptions and large, explosive eruptions sometimes culminating in caldera collapse. It is active on both geologic and historic timescales with well-documented eruptions, including a Plinian eruption (~3.6 ka) and several effusive events (most recently in 1950). Despite many geodetic and petrologic studies of Santorini, the structure of the magmatic system remains largely unknown. Using a 3D velocity model derived from active-source seismic data collected during the PROTEUS experiment and extensive synthetic testing, we constrain the shape, volume, and melt- content of shallow, inter-caldera magma storage. We image a magma body under the northern caldera basin, where a pronounced low-velocity anomaly from 2.8-5.4 km depth 2 causes seismic attenuation and ray-bending. The body has a minimum average melt- content of 10% and is consistent with depth estimates of pre-eruptive storage and a recent inflation episode. Synthetic tests show that the melt region likely extends below the resolution of our model and could contain up to 100% melt. Low-velocities are elongated ~15 km NE-SW, perpendicular to the maximum extensional stress direction from regional tectonism. This is attributed to a thin mush region (6-12% melt) confined by tectono-magmatic lineaments, which extends from an inter-caldera magma chamber (15- 25% melt). We suggest magma accumulation under Santorini is controlled by crustal extension and local edifice stresses. Our results provide the first seismic constraints on Santorini’s upper crustal magma system and its properties during inter-caldera periods. II. Introduction and Background Caldera forming eruptions require the pressurization and storage of large magma volumes in the shallow crust, sometimes thousands of cubic kilometers. After an explosive eruption, the shallow magmatic system reforms via small repeat injections of melt from the lower crust (Degruyter et al., 2016). These injections will accumulate where crustal stresses are least compressive due to topographic loads and regional tectonic extension (Corbi et al., 2015). Small shallow magma reservoirs govern inter- caldera activity but may only be temporary features that freeze after eruption since thermal conditions in the upper crust are not conducive to host small molten bodies. Conditions would be more favorable if the magma reservoir sat within a long-lasting, transcrustal mush system that maintains increased temperatures throughout the crust (Cashman et al., 2017). We use P-wave tomography to identify the volume, geometry, 3 and minimum melt-content of shallow magma storage below Santorini to constrain the anatomy of upper crustal magmatic systems during inter-caldera periods. Santorini is the most active volcano in the Hellenic island arc, sourced by the subduction of African oceanic lithosphere beneath the Aegean Sea plate (Fig. 1). For ~360 ka it has maintained a distinctive cycle of Plinian eruptions sometimes culminating in caldera collapse and interplinian periods of frequent, effusive shield-building events (Druitt et al., 1999). The Late Bronze Age (LBA aka Minoan) eruption was the last caldera-forming eruption, ~3600 years ago, of ~40-80 km3 dense rock equivalent of magma (Druitt, 2014). Subsequent effusive activity has formed the Kameni shield islands within Santorini caldera, including the most recent eruption of ~10-6 km3 in 1950 (Nomikou et al., 2014). Since the volcano is active on both geologic and historic time scales, it provides ideal conditions to image shallow inter-caldera storage that has quickly reformed since the last caldera collapse and is currently sourcing small eruptions. Santorini sits within a 40 km wide, N50E striking rift zone created by crustal extension of the Aegean Sea Plate (Nomikou et al. 2013). The rift zone is complemented by two volcanic lineaments that cross the caldera: the Kolumbo Line and Kameni Line (Fig. 1). Previously, the Kolumbo Line was defined by the alignment of Pleistocene vents within Santorini, and Kolumbo, the nearest active submarine cone (Druitt et al., 1999). The Kameni Line was defined by the alignment of post-LBA vents within the caldera and a small fault in the eastern caldera wall (Nomikou et al., 2014; Pyle and Elliott, 2006). For this study, we use a revised geometry of the lineaments, tomographically defined and discussed by Heath et al. (2019). 4 From 2011 to 2012, Santorini experienced an unrest period marked by increased H2 and CO2 emissions, earthquake swarms, and 10 cm of vertical displacement on the Kameni islands (Parks et al., 2015; Rizzo, 2015; Papadimitriou, 2015). Spherical Mogi source models of ground deformation measured by InSAR and GPS indicate 0.01-0.02 km3 of inflation north of the Kameni islands between 4 and 5 km depth (Parks et al., 2015), which is attributed to an intrusion of melt into a shallow magma chamber. The depth estimate of the intrusion agrees with a melt inclusion study of Kameni island 729 AD eruption products indicating a pre-eruptive storage depth of ~4 km (Druitt et al. 2016). Previous tomographic results of the upper 3 km of Santorini caldera suggest magma pooling could be concentrated north of the Kameni islands due to edifice stresses (Hooft et al., 2019). Based on past erupted volumes and a long-term average recharge rate of ~10-3 km2/yr, Degruyter et al. (2016) provide a best fit model of the current shallow chamber with 50% melt and a volume on the order of 10 km3. Despite multidisciplinary evidence of magma storage under Santorini, no geophysical studies have had sufficient resolution to characterize the existence of shallow magma accumulation, until now. Unlike most tomographic studies of arc volcanoes, our dense survey geometry provides the resolution required to constrain the shape and size of melt storage in the shallow crust, with constraints on melt-content and interactions between Santorini’s magmatic system and local tectonic lineaments. The unparalleled ray coverage also improves the understanding of the limitations and capabilities of seismic tomography when imaging low-velocity magmatic systems. Our results combine a high- resolution seismic experiment with ample geologic supporting evidence to understand the nature of shallow reservoirs at active caldera-forming arc volcanoes. 5 Figure 1: PROTEUS experiment geometry. Bathymetry (Hooft et al., 2017) and seismic data collected during PROTEUS experiment (Heath et al., 2019). Kameni Line and Kolumbo Line defined by Heath et al., (2019). Black square shows subset of model discussed in this study. Model is defined by x and y axis shown, rotated 22.5° from north. Inset map shows regional geography and volcanic arc. III. Data and Methods Santorini is semi-submerged, so exceptional ray coverage can be achieved through marine-land seismometers and airgun sources within and surrounding the caldera. In 2015, the PROTEUS experiment collected high-density seismic data utilizing 90 ocean bottom seismometers, 65 seismic land stations, and over 14,000 controlled- sound marine sources from the R/V Marcus Langseth (Fig. 1). We perform 3D tomography of over 230,000 Pg travel-times and present a subset of the full experiment in a 15 x 15 x 6 km velocity model of the caldera and surrounding area (Fig. 1). We build on the results of Heath et al. (2019), who provide a complete description of the 6 experiment, the seismic data, and the tomographic inversion, which we summarize and augment below. For this study, we added ~30,000 Pg first arrivals on shot-receiver combinations that probe depths >3 km below the caldera. Rays that intersect the northern caldera basin show signs of attenuation (i.e. decreased amplitudes, delayed arrivals, and lower signal- to-noise ratio (SNR)). To target the shallow magmatic system, we picked data on azimuthal swaths through the caldera to maintain lateral continuity between rays approaching the caldera and those passing through the magmatic system (Fig. 2). These data are filtered from 3-12 or 4-12 Hz and picked on the peak negative motion, with a median error visually assigned of 20 msec. The picking errors for the previous dataset (~200,000 arrivals) range from 5 msec to 30 msec with a median of 10 msec (Heath et al., 2019). We employed an iterative technique similar to Heath et al. (2019) to recognize mispicked arrivals. Data with high SNR are picked first and the generally exceptional data quality drives the velocity model. Subsequent picks on data with low SNR are compared to predicted travel-times through the model to identify and revise possible cycle-skips. We duplicate some of Heath et al.’s (2019) picks using the chosen lowpass filter and assigned pick phase, on stations with high SNR, and calculate a static offset to merge the datasets of travel-time picks. We utilize the same tomographic method (Toomey et al., 1994), smoothing and penalty parameters, and 3D starting model as Heath et al. (2019). The final velocity model utilized all ~230,000 arrivals, has horizontal and vertical node spacing of 200 m and 100 m, respectively, and was generated after 5 iterations with an RMS decrease from 107 msec for iteration 1 to 15.5 msec for iteration 5. 7 Figure 2: Caldera-crossing seismic data. (A) Azimuthal seismic data recorded at one station, flattened on predicted arrivals through the final velocity model. Data is a stack of hydrophone and vertical seismometer, filtered from 4-12 Hz. Travel-time picks shown in blue and red. Red travel-time picks are on attenuated data. (B) Raypaths of the same seismic data shown, probing the initial model, colored to match travel-time picks. Grey raypaths have no travel-time pick. Initial ray coverage is evenly distributed with no bending. (C) The same raypaths, now probing the final model. Those rays associated with attenuated data (red), bend around a volume in the northern caldera basin. Ray-bending and seismic attenuation provide direct evidence for a low-velocity volume north of the Kameni islands (caldera low-velocity, CLV). Raypaths can be calculated for each iteration of the inversion (Toomey et al., 1994). Rays travelling through the starting model have uniformly distributed paths between source and receiver with no bending (Fig. 2, Fig. S1; See Appendix for all supplemental figures). As the inversion progresses through each iteration, the same rays bend around a low-velocity region in the northern caldera basin to take the fastest path from source to receiver. These rays are consistently associated with the most attenuated data (Fig. 2, Fig. S1). No rays probe the CLV after 5 iterations, so we are limited to a minimum interpretation of the 8 physical properties within the volume. We run dozens of synthetic tests to further constrain the geometry, depth extent, and melt-content of the CLV and use geologically motivated synthetic models to probe the magmatic structure beneath Santorini and to test the resolving power of the tomographic results. IV. Results Below the northern caldera basin there are two distinct low-velocity bodies separated by a faster pinch-out (Fig. 3, Fig. S2). Hooft et al. (2019) describe the upper (<3 km depth) velocity structure and attribute the upper low-velocity body to a high- porosity column created by the collapse of a limited area of the caldera floor. The azimuthal data contributed for this study deepened model coverage to 6 km depth, revealing a second low-velocity volume (caldera low-velocity, CLV) between 2.8 and 5.4 km, directly below Hooft et al.’s (2019) high-porosity column (Fig. 3, Fig. S2). The porous column and the CLV are likely separated by a competent layer of solidified intrusives. 9 Figure 3: P-wave velocity model. (A) Depth slice at 3.4 km depth with lines showing cross section locations. Bounds of synthetic magma chamber shown as white oval over the CLV (Fig. S6). Kolumbo shown as white star in black circle. (B) Cross-section along profile A-A’. Bounds of synthetic magma chamber outlined in white, surrounding the CLV (Fig. S6). (C) Cross-section along profile B-B’, through the linear low-velocity (LLV). Bounds of synthetic magma chamber outlined in white (Fig. S6). All panels use the same color scale and contour interval of 2 km/s. Whited-out portions of the image show inadequate initial ray coverage. 10 The CLV is associated with a maximum P-wave velocity anomaly of -21% at 3.4 km depth (Fig. S2). This anomaly is significantly larger than those associated with magma storage at other arc volcanoes: Kiser et al. (2018) recover an anomaly of at most -9% attributed to magma storage 4 to 6 km below Mt. Saint Helens, Heath et al. (2015) find a maximum anomaly of -10% indicating storage between 3 and 5 km depth at Newberry Volcano, and Paulatto et al. (2012) report an anomaly of up to -12% below Monserrat with magma storage between 5.5 and 7.5 km. Surrounding the CLV to the NW and SE are regions of high velocity, similar to observations at Newberry, Mt. Saint Helens, and Monserrat (Heath et al., 2015; Kiser et al., 2018; Paulatto et al., 2012). These fast regions are likely solidified intrusives, which focus rays around the CLV, a common phenomenon when imaging magmatic systems (Lees, 2007). We also consider the vertical extent of these tomographically identified magmatic features. Below the CLV, relatively low-velocities continue below the resolution of our model. At Mt. Saint Helens, low-velocities form a nearly continuous feature from 3.5 to 14 km depth (Kiser et al., 2018) and the low-velocity volume under Monserrat extends below the resolution of the tomographic model (Paulatto et al., 2012). These observations may provide initial evidence of vertically extensive magmatic systems at arc volcanoes. We consider temperature and pressure conditions, magma composition, and pore space geometry to find the melt fraction required by the tomographically recovered CLV, (Text S1; See Appendix for all supplemental text; Berryman, 1980; Dunn et al., 2000). After correcting for a background geotherm of volcanic arcs, the temperature required to explain the anomaly exceeds 700°C, the solidus of andesite, indicating partial melt is required. We use the self-consistent effective medium method (Berryman, 1980) and 11 consider two end-member pore aspect ratios (0.05 and 0.2) to show the CLV requires an average of 10% melt from 2.8 to 5.4 km depth with a maximum of 8-17% at 3.4 km (Fig. S3). These melt fractions are minimum bounds since the tomographic inversion under- recovers the low velocities beneath Santorini. We use a suite of synthetic tests to explore magmatic structures that are consistent with tomographic inversion of the observed travel times (Text S2). These show that a fully molten body would be tomographically resolved to have the same minimum melt- content as the result (Fig. S4). Tests of depth extent of the low velocity body show that we can exclude a melt accumulation that is thinner than what is recovered (Fig. S5). Instead, melt must be present at least to 5.4 km depth. In fact, the tomographic recovery of a velocity anomaly of the same depth extent as the observed CLV requires that melt extend beyond the bottom of the recovered model (Fig. S6, white line on Fig. 3). Melt likely extends past 6 km depth but decreasing ray coverage, smoothing and penalty constraints limit the depth of the observed CLV to 5.4 km. The low-velocity region under Santorini is elongated ~15 km NE-SW towards Kolumbo seamount and perpendicular to regional extension direction (linear low- velocity, LLV). Within this feature, we recover two localized low-velocity concentrations with similar magnitudes to the CLV. However, synthetics show that these localized anomalies are an artifact of the tomographic inversion; complex raypaths interacting with a faster layer above the LLV are focusing low-velocities at these locations (Fig. S7). Synthetics do show that a linear extension of low-velocities out of the caldera is real and requires at least a small melt fraction. 12 The final synthetic model reflects the magmatic structure that we infer is consistent with the observed data and their tomographic signatures (Fig. S8). It includes a vertically extensive melt region beneath the northern caldera of Santorini (15-25% melt to at least 8 km depth) and a shallow NE-SW linear elongation of mush (6-12% melt to 5 km depth) northeast of the caldera (Fig. S8). This model is our preferred interpretation (Fig. 4) and the tomographic model recovered through this structure reproduces the main features of the observed model well; however, tomography is inherently nonunique so other interpretations may also fit the recovered model. The results support melt accumulation below the caldera with a higher melt-content and a deeper extent than the linear low-velocity (LLV). This is consistent with geologic observations and with pronounced attenuation and ray-bending through the caldera. V. Interpretation The caldera low-velocity (CLV) under Santorini’s northern caldera basin is centered at the approximate depth and location of the 2011-2012 inflation source (Parks et al., 2015; Fig. 4). Its depth is also consistent with pre-eruptive storage of Kameni island dacites (Druitt et al. 2016). Rarely do seismic velocity models, geodetic measurements, petrologic constraints, and geologic observations all converge on the location of magma storage in the upper crust. We infer this feature is the location of inter- caldera magma storage sourcing the Kameni island eruptions since 46 AD and likely since the LBA caldera-forming eruption ~3600 years ago. We estimate that the volume of this chamber is ~40 km3; an ellipsoid with 3 axes: the depth extent required by melt calculation curves, and the width and length of the 13 synthetic magma body (Fig. S6, cross-sections A and B, respectively). Using an average minimum melt-content of 10%, at least 4 km3 of melt currently resides below the northern caldera basin. This includes the volume of 0.01-0.02 km3 intruded in 2011-2012 (Parks et al., 2015) and is similar to the total volume of the Kameni islands (3.2 km3; Nomikou et al., 2014). This minimum estimate of 4 km3 of pure melt accumulated into a compact body in just 3,000 years indicates the shallow magmatic system is long-lived and quick to reform geologically. Thermo-mechanical models of the current state of inter-caldera magma storage estimate a volume of 10 km3 with 50% melt (Townsend et al., 2019; Degruyter et al., 2016) - within the range of possible volumes based on the seismic velocity model. The imaged body is consistent with the record of frequent, low volume eruptions that are expected from a small inter-caldera magma chamber (Townsend et al., 2019). The LLV from Santorini towards Kolumbo indicates a strong relationship between tectonism and magmatism. It is bounded by parallel, tomographically defined tectono-magmatic lineaments: the Kolumbo Line to the NW and the Kameni Line to the SE (Heath et al., 2019; Fig. S9). These lineaments are sub-parallel to active NE-SW trending regional faults and regional volcanic chains (Heath et al., 2019). All eruptive vents of the last 180 kyr are centered on, or between, these lineaments, as well as all dikes observed in the caldera walls (Fig. S9; Hooft et al., 2019; Druitt et al. 1999; Browning et al., 2015). Earlier earthquake tomography studies also recovered a linear low-velocity region elongated in this direction, but from 5 to 7 km depth, and attributed it to a tectonic zone which could provide pathways for ascending magma (Dimitriadis et al., 2010). In spite of the geometric alignment, we hesitate to attribute this feature to a 14 shallow magmatic connection between Santorini and Kolumbo since (1) it terminates short of, and slightly to the northwest of, Kolumbo, and (2) eruptive products from the two volcanoes indicate independent crustal differentiation (Klaver et al., 2016). Instead, this feature could be tomographic evidence of a series of dikes that have ruptured from the shallow magma chamber, consistent with the NE strike of all dikes observed in the caldera walls (Browning et al., 2015). It is also notable that recent volcanism is consistently located on the southern boundary of the feature, including the Kameni islands and Kolumbo. Figure 4: Conceptual model of shallow magmatic structure, from DV cross-sections (Fig. S2) with 2x vertical exaggeration. Pre-eruptive storage depth for Kameni island dacites, 4 km, shown as teal line (Druitt et al. 2016), and Mogi source estimates of recent inflation shown as yellow star with error bars (Parks et al., 2015). Melt accumulation sits below a high-porosity column from caldera collapse, separated by a competent layer (Hooft et al., 2019). A shallow mush region extends from the main chamber, perpendicular to regional extension, shown here in pink based on a dv contour of -10% past 2.4 km depth (Fig. S2). The mush region terminates short of and to the NW of Kolumbo. 15 Melt accumulation below Santorini is controlled by regional extension and topographic loads. The main magma body is located in the center of the caldera, directly below a low-density cylinder, where unloading rotates the least compressive stress to vertical (Corbi et al., 2015). The combination of topographic caldera rim loads, and localized unloading promotes post-caldera resurgence under the northern caldera basin. This is consistent with a recent influx of melt from the lower crust that is stored in the main magma body. While these localized top-down controls concentrate melt under the caldera, maximum regional horizontal stresses create a shallow NE-SW trending mush zone that extends from the main magma chamber. Shallow magma is confined by tectono-magmatic lineaments, created by crustal extension of the continental plate. These lineaments provide pathways for melt to migrate to the surface, as suggested by the location of past eruption centers (Fig. S9). Future shield-building eruptions will likely produce vents along the Kameni line, on the southeast boundary of the main magma body. VI. Conclusion During a Plinian eruption, large volumes of melt are evacuated from the shallow crust creating empty space that is often filled by caldera collapse. The shallow magmatic system must then reform through injections of melt from the lower crust, like the 2011- 2012 inflation episode. A long-lived, transcrustal mush system permits shallow magma to accumulate on geologically short timescales. At Santorini, melt injections pool in a small, crustal chamber from 2.8 to at least 5.4 km depth (Fig. 4). The magmatic structure within this ~40 km3 body is unknown; it could contain stacked sills, mush, or molten rock. 16 However, synthetic tests, melt estimates, ray bending, seismic attenuation, and geologic evidence all conclude that melt is most concentrated under the northern caldera basin. Magma accumulation continues from the main chamber, with decreased melt-content, for ~15 km NE (Fig. 4). This mush zone is bounded by tectono-magmatic lineaments, perpendicular to regional tectonic extension. The location of shallow magma pooling persists through the inter-caldera period (20-30 ka) fueling small, frequent shield-forming eruptions. Its eruptibility depends on the volume and rate of melt injections from the lower crust. This feature likely prevails through multiple eruption cycles since (1) extensional fracturing of the crust persistently focuses magmatism along tectono-magmatic lineaments and (2) edifice stresses generated by both low-density collapse and caldera topography repeatedly focuses inter-caldera pooling. These seismic constraints of upper crustal magma accumulation can be used to inform future studies of inter-caldera activity at Santorini. 17 APPENDIX SUPPLEMENTAL TEXT AND FIGURES Text S1: Physical Property Calculations We interpret the seismic velocities observed in the tomographic model through physical property calculations, including temperature, pressure, porosity, and the material filling pores. We start with the 1D velocity profiles in Fig. S3: a reference Vp profile from a metamorphic horst and a caldera Vp profile. We correct the reference velocity profile for temperature and pressure conditions accounting for both anharmonic and anelastic effects with no frequency dependence (Karato, 1993), described by Hooft et al., (2019)). We use a background geotherm of 40°C/km for the upper 0.2 km, where seawater mixing with sediments keeps the temperature cool. Then from 0.2-1.2 km depth, the geotherm increases to 160°C/km. This steep gradient is expected above the magmatic system. For the remainder of the profile we revert back to the background geotherm of 40°C /km, based on temperature properties of a volcanic arc (Hooft et al., 2019). We subtract the temperature and pressure corrected reference profile from the caldera Vp profile, to find the Vp anomaly under the caldera that can be attributed to excess temperature, porosity, and pore fluid. We calculate the temperatures that would be required to generate the Vp anomaly (Dunn et al., 2000) and find that they exceed the solidus of andesite, 700°C. Andesite has a range of solidus estimates (700-900°C), but we choose the minimum because the CLV 18 under Santorini contains mostly dacite. We correct the Vp anomaly for temperatures up to 700°C. The remaining anomaly is attributed to fluid filled pores. We use a self-consistent effective medium method where fluid inclusions are randomly oriented spheroids with a specified pore-aspect ratio (Berryman, 1980). At each depth, the velocity of the fluid filled pores considers the temperature/pressure at that depth. We use two different pore aspect ratios for the upper 2.4 km, 0.05 (thin cracks in fractured igneous rock) and 0.5 (more equant pores like tuffs and ignimbrites) filled with water (Hooft et al., 2019). There is no evidence (heat flux, inclusions, etc.) that melt would be stored this shallow. For depths greater than 2.8 km, we use aspect ratios of 0.05 and 0.2, filled with andesitic melt. We believe a competent layer separates water filled pores and melt storage, which is why we do not attribute any physical properties between 2.4 and 2.8 km depth. Text S2: Synthetic Modelling We use simple synthetic models to test possible magmatic features that could generate the low-velocity features in the observed velocity model. Using the tomographic image alone, we are limited to a minimum interpretation of the true subsurface structure. This is especially true when imaging low-velocity features that are typically under- recovered because of ray-bending and smoothing/penalty constraints. To go beyond a minimum-recover tomographic result, we use simple synthetic tests of various geologically motivated models of the subsurface structure. We generate a suite of synthetic models to explore tomographic recovery and the resolution of various subsurface features. To test the tomographic recovery of these 19 synthetic models, we generate travel times through the model (only for the rays used in the inversion). Then we follow the same inversion method used to generate our observed results, replacing the travel-times we picked with travel-times through the synthetic model. The synthetic results, or recovered model, can be compared to the observed model which is based on the true data. The Caldera Low-Velocity (CLV) • We first test thin magmatic bodies to see if the CLV could be generated by smoothing of a thinner, high melt-content volume. We test a molten sill (V = 2.3 km/s), from 2.4-3.4 km depth within the CLV, removing melt below the sill and from the shallow linear low-velocity (LLV) by replacing the velocities with average velocities (Fig. S5). This result, and others like it, do not reproduce the CLV amplitude of the observed model, or the depth extent. Instead, all depths of the CLV must contain melt. It is not an artifact of smoothing. This is supported by direct evidence of attenuation below the caldera at all depths of the model. • We then test different velocities to describe the CLV, knowing that it must be as deep as the observed model. We replace the CLV, from 2.4-5.4 km depth, with the velocity of molten andesite (V=2.3 km/s; Fig. S4), with a mush (15-25%, V=4 km/s), and a stack of molten layer and mush. All these tests equally match the observed model. Any reduction in velocity within the CLV, greater than that needed to exclude rays from the CLV volume, cannot be recovered. Any low- velocity variations within the CLV cannot be tomographically imaged. These tests 20 also show that a completely molten CLV is not recovered to have the same depth extent as the observed CLV, likely because rays are travelling below the body. • We test larger melt storage regions to represent the CLV. We test a deep, molten CLV and replace the LLV with average background velocities (Fig. S6). We find that the amplitude and depth extent of the observed CLV are best reproduced by synthetic models with melt below the resolution of the seismic model. We also test a mush (15-25%) instead of a molten vertical extension of the CLV, which also reproduce the observed CLV. The CLV in the observed model appears to terminate at 5.4 km depth because of smoothing to the faster background model. Ray coverage decreases with depth so low-velocities cannot be maintained towards the bottom of the model space. Our synthetics show that magma storage likely exceeds 6 km depth. The Linear Low-Velocity (LLV) • Our synthetic test of the deep, molten caldera low-velocity (CLV) with no LLV (Fig. S6) shows no sign of low-velocity smearing from the CLV to the NE. Therefore, the low-velocity trend must be real. • We replace the low-velocity elongation with a continuous body of molten andesite from 2.4-5 km depth (Fig. S7), keeping the molten CLV from Fig. S4. The recovered model shows the same two localized low-velocity concentrations found in the observed model. We test the same geometry within a 1D starting model and try an LLV with a mush region (6-12%, V=5 km/s). All results show the same two low-velocity concentrations. We conclude ray coverage and 21 complex raypaths are creating localized low-velocity bodies within the LLV. Raypaths are complicated by the fast region that sits above the LLV. • We test a deep LLV that extends beyond the resolution of the velocity model. A deeper LLV is recovered, which is not consistent with the observed model. The CLV is the only region that requires deeper melt accumulation. Preferred Model • Using all the evidence we have derived from these synthetic models, we create a preferred model (Fig. S8). We use this preferred model to support our interpretation of the tomographic results (Fig. 4). It includes a deep melt accumulation, 2.4-8 km below the northern caldera basin, of 15-25% melt. A shallow, linear elongation of mush (6-12%) from 2.4-5 km depth extends from the CLV out of the caldera to the NE. The recovered model does not exactly match the observed model; however, it does explain the main features found in the tomographic results. 22 Figure S1: All raypaths for travel-times picked for this study. (A) All rays probing the initial model, in grey, have even distribution and no bending. (B) All rays probing the final model bend around a volume in the northern caldera basin, correlated with the location of recent inflation. Figures are clipped to 6 km depth (where resolution is lost in the model). Stations used for this study in yellow. 23 Figure S2: P-wave velocity anomalies relative to average velocity, in percent. (A) Depth slice at 3.4 km depth with lines showing cross section locations. Bounds of synthetic CLV (main magma chamber) shown as black oval (Fig. S6). Kolumbo shown as white star in black circle. (B) Cross-section along profile A-A’. Bounds of synthetic CLV (main magma chamber) outlined in black (Fig. S6). (C) Cross-section along profile B-B’, through the LLV. Bounds of synthetic CLV (main magma chamber) outlined in black (Fig. S6). All panels use the same color scale and contour interval of 2%. Whited-out portions of the image show inadequate initial ray coverage. 24 Figure S3: Physical property calculation curves. (A) P-wave velocity model at 3.4 km depth with locations of 1D velocity profiles. Black box shows area used to average a background velocity profile for a metamorphic horst. Star corresponds to 1D profile of caldera low-velocity body. (B) Velocity profiles, colored to match locations on map. Green dashed line is the temperature corrected background profile, including temperatures up to 700C. Area between green dashed curve and red curve is explained by melt/water filled pores. (C) Background geotherm in black. Temperature required to generate the velocity anomaly under the caldera (observed – background) in green. Temperature curve, clipped to the solidus of andesite, in dashed green. (D) Physical property curves labeled with corresponding pore aspect ratio. From 0-2.4 km depth, pores are filled with water, from 4.8-6 km pores are filled with andesitic melt. From 2.4-2.8 km there is a competent layer. 25 Figure S4: Synthetic test of molten Caldera Low Velocity (CLV). (A) Synthetic model of CLV from 2.4–5.4 km depth. The LLV and depths below the molten CLV are replaced with average velocities. Travel-times are generated through this model. (B) Tomographic recovery of the synthetic model. (C) Observed model, same as Fig. 3, for comparison. All figures use the same depth slice at 3.4 km, color scale, cross-section location, and a contour interval of 0.2 km/s. Recovered model nearly reproduces the CLV in size and magnitude. 26 Figure S5: Synthetic test of molten sill. (A) Synthetic model of thin CLV, 1 km thick from 2.4-3.4 km depth. The LLV and depths below the molten sill are replaced with average velocities. Travel-times are generated through this model. (B) Tomographic recovery of the synthetic model. (C) Observed model, same as Fig. 3, for comparison. All figures use the same depth slice at 3.4 km, color scale, cross-section location, and a contour interval of 0.2 km/s. Recovered model does not reproduce the CLV in magnitude or in size. Melt must extend below 4 km depth to match the observed model. 27 Figure S6: Synthetic test of deep molten CLV. (A) Synthetic model of deep CLV to 8 km depth. The LLV is replaced with average velocities. Travel-times are generated through this model. (B) Tomographic recovery of the synthetic model. (C) Observed model, same as Fig. 3, for comparison. All figures use the same depth slice at 3.4 km, color scale, cross-section location, and a contour interval of 0.2 km/s. Recovered model best reproduces the CLV in magnitude and in size. Melt likely extends beyond the resolution of the velocity model (6 km depth) to match the observed model. 28 Figure S7: Synthetic test of molten Linear Low Velocity (LLV). (A) Synthetic model of LLV from 2.4–5 km depth and CLV from Fig. S4. Travel-times are generated through this model. (B) Tomographic recovery of the synthetic model. (C) Observed model, same as Fig. 3, for comparison. All figures use the same depth slice at 3.4 km, color scale, cross-section location, and a contour interval of 0.2 km/s. Recovered model shows the same two localized concentrations of low-velocity within the LLV. While the LLV is real, these localized bodies are an artifact of the tomographic inversion. 29 Figure S8: Synthetic test of preferred model. (A) Synthetic model of CLV with 15-25% melt (V=4 km/s to 8 km depth) and LLV with 6-12% melt (V=5 km/s to 5 km depth). Travel-times are generated through this model. (B) Tomographic recovery of the synthetic model. (C) Observed model, same as Fig. 3, for comparison. All figures use the same depth slice at 3.4 km, color scale, cross-section location, and a contour interval of 0.2 km/s. Recovered model reproduces the magnitude and size of the CLV and the LLV, serving as an interpretation of the observed model (Fig. 4). 30 Figure S9: Past eruptions controlled by tectono-magmatic lineaments. Bathymetry (Hooft et al., 2017) with Kolumbo and Kameni lineaments defined by Heath et al. (2019). These lineaments bound the CLV and LLV, shown in purple as a velocity contour of 5.3 km/s, and are perpendicular to extension direction. All dikes observed in the caldera walls (Browning et al., 2015) and all vents for the last 180 kyr (Hooft et al., 2019; Druitt et al., 1999) fall on or within these lineaments, suggesting tectonic extension is a dominant control on crustal magmatism. All dikes have a strike ~NE-SW. 31 REFERENCES CITED Berryman, J. G., 1980, Long wavelength propagation in composite elastic media I. Spherical inclusions and II. Ellipsoidal inclusions: Journal Acoustical Society of America, v. 68, p. 1809–1819. https://doi.org/10.1121/1.385172. Browning, J., Drymoni, K., and Gudmundsson, A., 2015, Forecasting magma-chamber rupture at Santorini volcano, Greece: Scientific Reports, v. 5, n. 15785. https://doi.org/10.1038/srep15785. Cashman, K.V., Sparks, R. S. J., and Blundy, J. D., 2017, Vertically extensive and unstable magmatic systems: a unified view of igneous processes: Science, v. 355 (6331), https://doi.org/10.1126/science.aag3055. Corbi, F., Rivalta, E., Pinel, V., Maccaferri, F., Bagnardi, M., and Acocella, V., 2015, How caldera collapse shapes the shallow emplacement and transfer of magma inactive volcanoes: Earth and Planetary Science Letters, v. 431, p. 287-293, https://doi.org/10.1016/j.epsl.2015.09.028. Degruyter, W., Huber, C., Bachmann, O., Cooper, K. M., and Kent, A. J. R., 2016, Magma reservoir response to transient recharge events: The case of Santorini volcano (Greece): Geology, v. 44, p. 23-26. Dimitriadis, I., Papazachos, C., Panagiotopoulos, D., Hatzidimitriou, P., Bohnhoff, M., Rische, M., and Meier, T., 2010, P and S velocity structures of the Santorini- Coloumbo volcanic system (Aegean Sea, Greece) obtained by non-linear inversion of travel times and its tectonic implications: Journal of Volcanology and Geothermal Research, v. 195(1), p. 13–30, https://doi.org/10.1016/j.jvolgeores.2010.05.013. Druitt, T. H., Mercier, M., Florentin, L., Deloule, E., Cluzel, N., Flaherty, T., Médard, E., and Cadoux, A., 2016, Magma storage and extraction associated with Plinian and Interplinian activity at Santorini caldera (Greece): Journal of Petrology, v. 57, p. 461–494, https://doi.org/10.1093/petrology/egw015. Druitt, T. H., Edwards, L., Mellors, R. M., Pyle, D. M., Sparks, R. S. J., and Lanphere, M., 1999, Santorini Volcano Geological Society Memoir (19th ed.): Geological Society Memoir. Druitt, T. H., 2014, New insights into the initiation and venting of the Bronze-Age eruption of Santorini (Greece), from component analysis: Bulletin of Volcanology, v. 76, p. 794, https://doi.org/10.1007/s00445-014-0794-x. 32 Dunn, R. A., Toomey, D. R., and Solomon, S. C., 2000, Three-dimensional seismic structure and physical properties of the crust and shallow mantle beneath the East Pacific Rise at 9°30'N: Journal of Geophysical Research, v. 105(B10), p. 23537- 23555, https://doi.org/10.1029/2000JB900210. Heath, B. A., Hooft, E. E. E., Toomey, D. R., and Bezada, M. J., 2015, Imaging the magmatic system of Newberry Volcano using joint active source and teleseismic tomography: Geochemitry, Geophysics, Geosystems, v.16, p. 4433–4448, https://doi.org/10.1002/2015GC006129. Heath, B. A., Hooft, E. E. E., Toomey, D. R., Papazachos, C. B., Nomikou, P., Paulatto, M., Morgan, J. V., and Warner, M. R., 2019, Tectonism and its relation to magmatism around Santorini volcano from upper crustal P-wave velocity: Journal of Geophysical Research. Hooft, E. E. E., Nomikou, P., Toomey, D. R., Lampridou, D., Getz, C., Christopoulou, M. E., O’Hara, D., Arnoux, G. M., Bodmer, M., Gray, M., Heath, B.A., and VanderBeek, B. P., 2017, Backarc tectonism, volcanism, and mass wasting shape seafloor morphology in the Santorini-Christiana-Amorgos region of the Hellenic Volcanic Arc: Tectonophysics, v. 712–713, p. 396–414, https://doi.org/10.1016/j.tecto.2017.06.005. Hooft, E. E. E., Heath, B. A., Toomey, D. R., Paulatto, M., Papazachos, C. B., Nomikou, P., Morgan, J. V., and Warner, M. R., 2019, Seismic imaging of Santorini: Subsurface constraints on caldera collapse and present-day magma recharge: Earth and Planetary Science Letters, v. 514, p. 48-61. https://doi.org/10.1016/j.epsl.2019.02.033. Karato, S., 1993, Importance of anelasticity in the interpretation of seismic tomography: Geophysical Research Letters, v. 20(5), p. 1623–1626. Kiser, E., Levander, A., Zelt, C., Schmandt, B., and Hansen, S., 2018, Focusing of melt near the top of the Mount St. Helens (USA) magma reservoir and its relationship to major volcanic eruptions: Geology, v. 46(9). p. 775-778. https://doi.org/10.1130/G45140.1. Klaver, M., Carey, S., Nomikou, P., Smet, I., Godelitsas, A., and Vroon, P., 2016, A distinct source and differentiation history for Kolumbo submarine volcano, Santorini volcanic field, Aegean arc: Geochemistry, Geophysics, Geosystems. https://doi.org/10.1002/2016GC006398. Lees, J. M, 2007, Seismic tomography of magmatic systems: Journal of Volcanology and Geothermal Research, v. 167, p. 37-56. 33 Nomikou, P., Papanikolaou, D., Alexandri, M., Sakellariou, D., and Rousakis, G., 2013, Submarine volcanoes along the aegean volcanic arc: Tectonophysics, v. 597–598, p. 123–146. https://doi.org/10.1016/j.tecto.2012.10.001. Nomikou, P., Parks, M.M., Papanikolaou, D., Pyle, D.M., Mather, T.A., Carey, S., Watts, A.B., Paulatto, M., Kalnins, M.L., Livanos, I., Bejelou, K., Simou, E., and Perros, I., 2014, The emergence and growth of a submarine volcano: the Kameni islands, Santorini (Greece): GeoResJ, v.1–2, p. 8–18. https://doi.org/10.1016/j.grj.2014.02.002. Papadimitriou, P., Kapetanidis, V., Karakonstantis, A., Kaviris, G., Voulgaris, N., and Makropoulos, K., 2015, The Santorini Volcanic Complex: a detailed multi- parameter seismological approach with emphasis on the 2011–2012 unrest pe- riod: Journal of Geodynamics, v. 85, p. 32–57, https://doi.org/10.1016/j.jog.2014.12.004. Parks, M. M., Moore, J., and Papanikolaou, X., 2015, From quiescence to unrest: 20 years of satellite geodetic measurements at Santorini volcano, Greece: Journal of Geophysical Research, v. 120. https://doi.org/10.1002/(ISSN)2169-9356. Paulatto, M., C. Annen, T. J. Henstock, E. Kiddle, T. A. Minshull, R. S. J. Sparks, and B. Voight, 2012, Magma chamber properties from integrated seismic tomography and thermal modeling at Montserrat: Geochemistry Geophysics Geosystems, v. 13, Q01014, doi:10.1029/2011GC003892. Pyle, D. M., and Elliott, J. R., 2006, Quantitative morphology, recent evolution, and fu- ture activity of the Kameni Islands volcano, Santorini, Greece: Geosphere, v. 2,253, https://doi.org/10.1130/GES00028.1. Rizzo, A. L., Barberi, F., Carapezza, M. L., Piazza, A. Di., Francalanci, L., Sortino, F., and D’Alessandro, W., 2015, New mafic magma refilling a quiescent volcano: Evidence from He-Ne-Ar isotopes during the 2011–2012 unrest at Santorini, Greece: Geochemistry Geophysics Geosystems, v. 16, p. 789-814, doi:10.1002/ 2014GC005653. Toomey, D. R., Solomon, S. C., and Purdy, G. M., 1994, Tomographic imaging of the shal-low crustal structure of the East Pacific Rise at 9◦30′ N: Journal of Geophysical Research, v. 99, p. 24135–24157, https://doi.org/10.1029/94JB01942. Townsend, M., Huber, C., Degruyter, W., and Bachmann, O., 2019, Magma chamber growth during intercaldera periods: Insights from thermos mechanical modeling with applications to Laguna del Maule, Campi Flegrei, Santorini, and Aso: Geochemistry, Geophysics, Geosystems, v. 19, https://doi.org/10.1029/2018GC008103 34