# Crustal thickness and magma storage beneath the Ecuadorian arc

Clinton D. Koch, Jonathan Delph, Susan L. Beck, Colton Lynner, Mario Ruiz, Stephen Hernandez, Pablo Samaniego, Anne Meltzer, Patricia Mothes, Silvana Hidalgo

► To cite this version:

Clinton D. Koch, Jonathan Delph, Susan L. Beck, Colton Lynner, Mario Ruiz, et al.. Crustal thickness and magma storage beneath the Ecuadorian arc. Journal of South American Earth Sciences, 2021, 110, 10.1016/j.jsames.2021.103331. insu-03708917

HAL Id: insu-03708917

https://insu.hal.science/insu-03708917v1

Submitted on 2 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 **Crustal Thickness and Magma Storage Beneath the Ecuadorian Arc**

2 Clinton D. Koch¹, clinton.d.koch@gmail.com

3 Jonathan Delph², jdelph@purdue.edu

4 Susan L Beck¹, slbeck@email.arizona.edu

5 Colton Lynner³, clynner@udel.edu

6 Mario Ruiz⁴, mruiz@igepn.edu.ec

7 Stephen Hernandez⁴, hrndz.stphn@gmail.com

8 Pablo Samaniego⁵, pablo.samaniego@ird.fr

9 Anne Meltzer⁶, ameltzer@lehigh.edu

10 Patricia Mothes⁴, pmothes@igepn.edu.ec

11 Silvana Hidalgo⁴, shidalgo@igepn.edu.ec

12

13 ¹Department of Geosciences, University of Arizona, Tucson, AZ

14 ²Department of Earth, Atmospheric, and Planetary Sciences, Purdue University, West Lafayette, IN

15 ³Department of Earth Sciences, University of Delaware, Newark, DE

16 ⁴Instituto Geofísico, Escuela Politécnica Nacional, Quito, Ecuador

17 ⁵Laboratoire Magmas et Volcans, IRD-CNRS-OPGC, Université Clermont Auvergne, Clermont-

18 Ferrand, France

19 ⁶Department of Earth and Environmental Sciences, Lehigh University, Bethlehem, PA

20 Corresponding Author: Clinton D. Koch, Clinton.d.koch@gmail.com, (605)209-6162

21
22 **Abstract**

The Northern Andes of Ecuador contain some of the most active volcanic systems in the Andes and extend over a broad region from the Western Cordillera to the Subandean Zone. While it is known that the arc straddles a range of basement compositions, from accreted mafic oceanic terranes in the west to silicic continental terranes in the east, the details of the crustal structure beneath the arc is unclear despite being critical for understanding magmatic and tectonic processes in this portion of the Andes. To gain insight into these processes, we create two 3D models of crustal and upper mantle seismic properties throughout the region. The first highlights the discontinuity structure using receiver functions, which allows for the recovery of crustal thickness beneath the Ecuadorian Andes. We observe a range from ~50 to 65 km under the high elevations, with thicker crust beneath the lower elevation Western Cordillera compared to the higher elevation Eastern Cordillera. This can largely be explained by density variations within the crust that are consistent with observed terranes at the surface, implying these terranes extend to depth. The second model combines our receiver functions with Rayleigh wave dispersion data from ambient noise measurements in a joint inversion to construct a 3-D shear wave velocity model. This model shows several mid-crustal (5-20 km below sea-level) low velocity zones beneath Ecuadorian arc volcanoes that contain a maximum of ~14% melt. These low velocity zones likely represent zones of long-term magma storage in predominantly crystalline reservoirs, consistent with "mush zones". Furthermore, the depth of the inferred reservoirs below several of the volcanic centers (e.g., Chiles-Cerro Negro and Tungurahua) are in broad agreement with previous geobarometry and geodetic modeling. Our results provide new observations of possible long-term magma reservoirs below other less-studied volcanic systems in the Ecuadorian arc as well, and further contributes to a mounting number of observations indicating long-term magma storage at low melt percentages in the mid-crust beneath active arc systems.
46 **Introduction**

47
The South American Andes hosts the longest continental volcanic arc and some of the
48
thickest crust in the world. The Andes have some of the thickest crust in the world and
49
understanding crustal thickness variations beneath the Andes has been the subject of many studies
50
(e.g., Yuan et al., 2000; Tassara and Echaurren, 2012; Poveda et al., 2015; Ryan et al., 2016; Condori
51
et al., 2017; Bishop et al., 2017; Rodriguez and Russo, 2020 and others). The crust-mantle
52
transition, where crustal rocks transition to peridotite, is often assumed to be synonymous with the
53
Mohorovičić discontinuity (the “Moho”), which is defined as the boundary where P-wave seismic
54
velocities increase to ~8.0 km/s (~4.5 km/s for S-waves; Christensen and Mooney, 1995). The
55
expected Moho signal beneath active orogenic belts can be complicated by magmatic underplating,
56
basal accretion/relamination, and delamination, thereby making the location and character of
57
crust-mantle transition difficult to constrain (Gilbert et al., 2006; Muntener and Ulmer, 2006;
58
Hacker et al., 2011, Frassetto et al., 2011; Beck et al., 2015). Constraining crustal thickness is
59
important for understanding many orogenic and magmatic processes.

60
Constraints on the structure of crustal magmatic systems is critical for understanding
61
where magma is formed and stored, and provides insight into the plausible pathways it takes as it
62
propagates to the surface. Seismological data can be used to image the in-situ architecture of
63
magnetic plumbing systems (e.g. Lees, 2007). Combining seismic images with other observations,
64
including those from petrologic, geochemical, and geodetic studies allows us to further constrain
65
some of these parameters (Cashman et al., 2017; Chaussard and Amelung, 2014). The presence of
66
low-velocity zones in seismic images beneath several volcanic systems, including several in the
67
Central Volcanic Zone of the Andes, have been associated with regions of partial melt and storage of
68
magma (Ward et al., 2014; Huang et al., 2015; Kiser et al., 2016; Delph et al., 2017; Ward et al.,
69
2017; Schmandt et al., 2019). Often geophysical models image crustal melt stored in the middle and
70
lower crust with melt fractions below the eruptible limit.
Despite a broad knowledge of crustal thickness beneath the Andes, local variations of crustal thickness within the Ecuadorian arc are poorly resolved. The Ecuadorian Andes are a ~150-200 km wide mountain belt in the Northern Andes with elevations reaching >6000 m and some of the most active volcanic systems in the Andes. Seismic data recorded by permanent instruments deployed by the Instituto Geofísico at the Escuela Politécnica Nacional (IG-EPN) along with other preexisting stations in Ecuador offer an opportunity to better resolve the crustal thickness, volcanism, and tectonics in this region. Previous studies in the region have estimated crustal thickness using gravity (Feininger and Seguin, 1983; Tamay et al., 2018) or local seismicity (Vaca et al., 2019). Additionally, seismic imaging has shown broad variations in the seismic velocities of the crust (Araujo 2016; Lynner et al., 2020). Here, we present the first receiver function-derived crustal thickness map beneath the Ecuadorian Andes. To better constrain arc processes, we combine receiver functions with ambient noise dispersion data to construct a 3D shear-wave velocity model that illuminates crustal properties. Our crustal thickness results show that the region is largely in isostatic equilibrium across all tectonic regimes. Additionally, we image several low shear velocity zones in the upper ~20 km of the crust beneath many active volcanic systems in both the Western and Eastern Cordilleras.

## Regional Geologic Setting

The Ecuadorian orogenic system can be divided into six major geologic provinces which are roughly oriented parallel the trench. From west to east these are: the forearc, Western Cordillera, Inter-Andean valley, Eastern Cordillera, Subandean Zone, and the Oriente basin (Fig. 1). The forearc of Ecuador consists of Mesozoic sedimentary basins separated by localized outcrops of mafic material. The mafic basement consists of accreted oceanic terranes that are thought to underlie much of the forearc region (Jaillard et al, 1995; Luzieux et al., 2006, Koch et al., 2020). The Western Cordillera basement consists of similar mafic lithologies, interlayered with oceanic sediments, which are thought to have accreted during a single event during the Late Cretaceous – Paleogene
96 (Jaillard et al., 2009; Vallejo et al., 2019), although multiple accretions have also been proposed
97 (Jaillard et al., 1995; Spikings et al., 2005). The Inter-Andean valley is a topographical depression
98 hosting thick Oligocene - Pleistocene volcanoclastic, fluvial, and lacustrine strata (Lavenu et al.,
99 1995). The Eastern Cordillera basement consists of Paleozoic through Jurassic metamorphosed
100 sedimentary and igneous rocks, Mesozoic granitoids, and metasedimentary rocks overlain by
101 Cenozoic volcanic deposits (Pratt et al., 2005; Chiaradia et al., 2009; Spikings et al., 2015). The
102 transition between the mafic basement beneath the Western Cordillera and the felsic basement
103 Eastern Cordillera is obscured by the Inter-Andean Valley. The Subandean Zone is a transitional
104 region including thin- and thick-skinned thrust belts. The Oriente basin consists of Mesozoic and
105 Cenozoic sedimentary sequences overlaying the Pre-Cambrian Guyanese craton (Aspden and
106 Litherland, 1992). Active arc volcanism is pervasive from the Western Cordillera to the Subandean
107 Zone and the landscape is covered with large stratovolcanoes (Hall et al., 2008).

108 *Crustal Thickness*

109 Previous crustal thickness estimates in Ecuador have come primarily from gravity
110 observations (Feininger and Seguin, 1983; Tamay et al., 2018), the joint determination of
111 earthquake locations and crustal structure (Vaca et al., 2019), and local tomography (Araujo, 2016).
112 Condori et al. (2017) and Poveda et al., (2015) used receiver functions and H-k stacking in Northern
113 Peru and Colombia respectively. These studies each extend crustal thickness interpretations into
114 southern and northern Ecuador, however, they had few constraints within most of the Ecuadorian
115 Andes. All of these studies have suggested crustal thicknesses ranging from ~45-75 km beneath the
116 cordilleras. Although the broad characteristics of the different crustal thickness models are similar,
117 with crustal thickness increasing beneath the high elevations of the Andes and thinning towards the
118 forearc and Subandean Zone, they exhibit variations in the detailed crustal structure of the region.
119 *Volcanism*

120 Throughout the Western Cordillera, Inter-Andean Valley, Eastern Cordillera, and Subandean
121 Zone, the active continental volcanic arc is characterized by mostly calc-alkaline basaltic andesite to
122 dacitic volcanism north of ~2°S (Fig. 1; Hall et al., 2008; Ancellin et al., 2017; Bablon et al., 2020).
123 South of ~2°S, the volcanic arc is inactive due to flat-slab subduction of the Nazca Plate (Nur and
124 Ben-Avraham., 1983). This transition separates the Northern Andean and the Central Andean
125 Volcanic Zones. Within the active arc, volcanism exhibits across-arc geochemical variations that
126 suggest a decreasing slab input into the mantle wedge away from the trench (Barragan et al., 1998;
127 Bourdon et al., 2003; Bryant et al., 2006; Chiaradia et al., 2009; Hidalgo et al., 2012; Ancellin et al.,
128 2017). This decrease in the amount of slab-derived fluids and/or melts is concomitant with a
129 decrease in the amount of mantle melting away from the trench. More recently, noticeable along-
130 arc geochemical variations have been described for the frontal arc volcanoes (mostly those
131 constructed on the Western Cordillera), and have been interpreted as due to variations in slab and
132 mantle processes (Ancellin et al., 2017, Narvaez et al., 2018) or crustal processes (Chiaradia et al.,
133 2020).

134 In the Western Cordillera and Inter-Andean Valley, arc volcanism constructed on the mafic
135 basement is referred to as the frontal arc and includes active volcanoes such as Chiles- Cerro Negro,
136 Cotocachi-Cuicocha, and Pichincha, among others (Fig. 1). Where the volcanism erupts through the
137 more felsic basement of the Inter-Andean Valley and the Eastern Cordillera, the arc is referred to as
138 the main arc. Geochemical proxies for crustal contamination show an increase both west to east and
139 north to south (Ancellin et al., 2017). Back arc volcanism is observed locally in the Subandean zone
140 and is characterized by alkaline magmatic suites enriched in most incompatible elements,
141 suggesting the lowest degree of mantle melting within the Ecuadorean Andes (Ancellin et al., 2017).
142 **Data and Methods**

143 Data for this study primarily come from the IG-EPN permanent seismic network (Alvarado
144 et al., 2018) with one station from the Global Seismic Network (Albuquerque Seismological
145 Laboratory (ASL)/USGS, 1988). In total, we calculate receiver functions from 66 broadband stations
146 located mainly in the northern Ecuadorian Andes (Fig. 1). Dispersion data from Lynner et al. (2020)
147 used in this study use additional forearc stations from the IG-EPN and a temporary network in the
148 forearc region that was deployed following the 2016 Pedernales earthquake (Meltzer et al., 2019).
149 Many of the IG-EPN stations are clustered around volcanoes as they are mainly used for monitoring
150 purposes, and this high density of stations leads to high-fidelity results around several arc
151 volcanoes.

152 *Receiver Functions*

153 P-wave receiver functions (RFs) are a technique that isolates P to S conversions from
154 teleseismic P-waves that result from velocity contrasts beneath a seismic station (Langston, 1979;
155 Ammon, 1991). In this study, we calculate P-wave radial RFs using data from earthquakes with
156 magnitude ≥ 5.9 and epicentral distances between 30° and 90°. We window the data to 20 seconds
157 before the P-wave and 100 seconds after. These time series are band-pass filtered between 0.05
158 and 2 Hz and rotated to a ZRT coordinate system to isolate P-to-S conversions to the radial
159 component. Before the construction of the RFs, the seismic data is visually inspected for a clear P-
160 wave arrival on the vertical and radial components with minimal P-wave energy on the tangential
161 component. RFs are then calculated using a time-domain iterative deconvolution with a Gaussian
162 parameter of 2.8 (Ligorria and Ammon, 1999). For further quality control, the individual RFs for
163 each station are manually inspected, and data with strong positive arrivals at ~0 seconds, radial fits
164 ≥ 65%, and minimal amounts of ringing are retained for further analysis. In total, 4139 RFs were
165 further analyzed out of 8341 calculated.
166 *Adaptive Common Conversion Point Stacking*

167 Common conversion point stacking creates a 3D model of receiver function discontinuity
168 structure throughout a study area (Dueker and Sheehan, 1997). Common conversion point stacking
169 maps the amplitudes of each RF into bins along its theoretical ray path and stacks the amplitude
170 where bins overlap, resulting in a 3-D grid of RF amplitudes. Here we use the adaptive common
171 conversion point (ACCP) stacking method that allows the bins to increase in size until a minimum
172 number of RFs pass through the bin (Delph et al., 2015, 2017). In this study, we use bins spaced at
173 0.1° with a minimum size of 0.2° and minimum hit count of 5 RFs and allow bins to dilate to a
174 maximum size of 0.4° if less than 5 RFs are contained in a bin. Due to the bin dimensions, the
175 smallest anomaly in our models will appear at least 22 km wide, with anomaly edges within 16.5
176 km of their true location due to lateral smoothing (or 0.75x the bin width). Each bin is 0.5 km thick,
177 and the model extends to a depth of 100 km. To migrate from time to depth we use a 1-D velocity
178 model averaged from the arc region from ambient noise tomography (Figure S1; Lynner et al.,
179 2020) and a Vp/Vs ratio of 1.8. This Vp/Vs ratio is justified via through the H-k stacking on receiver
180 functions near the arc (Zhu and Kinamori, 2000; see supplemental material and Figure S2), and is
181 consistent with the ~1.65 – 1.90 range of Vp/Vs ratios in the Colombian and Northern Peru
182 cordilleras (Poveda et al., 2015; Condori et al., 2017). We note that local tomography has suggested
183 lower Vp/Vs values, closer to an average of 1.75 (Araujo, 2016). A Vp/Vs of 1.75 shift arrivals slightly
184 deeper compared to assuming a 1.8 Vp/Vs. For example, a conversion arriving at 2 seconds would
185 be ~1 km deeper and at ~7 seconds would be ~3 km deeper than what we observe in our ACCP
186 stacks.

187 *Joint Inversion of RFs and Ambient Noise Dispersion*

188 We utilize the RFs from this study combined with ambient noise-derived Rayleigh wave
189 dispersion data from Lynner et al. (2020) to construct a shear-wave velocity model of the region.
190 Rayleigh wave dispersion data can recover absolute velocities of vertically-polarized shear waves,
but the broad sensitivity kernels of Rayleigh waves result in weak sensitivity to velocity discontinuities. Conversely, RFs provide strong constraints on velocity changes associated with seismic discontinuities but are not sensitive to absolute velocities. By combining these highly complementary data types, we mitigate the shortcomings of each dataset as well as the dependence on assumptions about velocity structure (Julia et al., 2000). The resulting shear-wave velocity model is sensitive to both vertical and lateral absolute velocities and their variations. The dispersion periods used in this study have robust sensitivities between ~5 km – 50 km, with decreasing resolution below 50 km (Lynner et al., 2020). Outside of these depths, the sensitivity to absolute shear-wave velocities will be limited and structure will be primarily controlled by RF arrivals.

Lynner et al. (2020) used intermediate and broadband seismic data from stations deployed in the arc and forearc (AN stations; Fig. 1) to obtain dispersion data. This dataset included stations from the IG-EPN network as well as a temporary network that was deployed following the 2016 Pedernales earthquake (Alvarado et al., 2018; Meltzer et al., 2019). Dispersion curves were obtained after cross-correlation of day-long records and inversion of interstation phase velocity measurements on a 0.1° x 0.1° grid. We use dispersion curves from each grid point in the Lynner et al. (2020) model in our joint inversion (examples of dispersion curves are shown in Figure S3). Additional details of the processing and inversion of the ambient noise data can be found in Lynner et al. (2020, and references therein).

We extract RFs from the ACCP stacked model at the same grid spacing as the dispersion curves so that they can be directly paired in the joint inversion (Delph et al., 2017; Koch et al., 2020). We follow the methodology described in Delph et al. (2017) in extracting the RFs which migrates the ACCP volume back to time using the same velocity model and applies a filter to minimize changes that result from the discretization of the model into bins of constant thickness and the variable amount of data present in each bin.
The joint inversion of the dispersion curves and ACCP-derived RFs utilizes a non-linear least squares method that begins by assuming an initial velocity model, in this case, a 4.44 km/s half-space model discretized into 1 km thick layers, as this is the theoretical resolution of the receiver functions we use in this study. Assuming a simple starting velocity model ensures that emerging structures are the result of our data and not imposed by assumptions in the starting model (such as a 2-D or 3-D model of velocity or crustal thickness), which can bias the resulting model if the *a priori* constraints are inaccurate. The predicted RF and dispersion data are calculated using this initial velocity model and compared with the observed data (RF and dispersion curve) at each grid point. The misfit between the observed and predicted is calculated and used to iteratively update the velocity model. This process is repeated for 40 iterations or until the change in the model between iterations is negligible (<0.05%). To ensure the starting model is not strongly biasing our results, we also test starting models with 3.6 km/s and 4.8 km/s (Figures S4 and S5). These models show minimal changes compared to our final velocity model in the upper 50 km. The weighting between RFs and dispersion data is also an important parameter. We choose a 70%/30% weighting of the RFs and dispersion data, respectively, as this produces the highest average fit between the two datasets (Figure S6) and is consistent with previous studies (Delph et al., 2017; Koch et al., 2020).

## 233 Results

### 234 Receiver Functions

We present results of the P-wave radial RFs for stations within the arc (Fig. 2). The longest-running seismic station in Ecuador, OTAV, is from the Global Seismic Network and is located on the eastern edge of the Western Cordillera near the Inter-Andean valley in Northern Ecuador (Fig. 1). As shown in Fig. 2, a strong negative conversion ($Ps_1$) representing a decrease in velocity with depth, is present at ~3.5 seconds on the RFs. At 5.5 and 6.5 seconds, two positive conversions can also be seen ($Ps_2$ and $Pms$, Fig. 2). These represent increases in velocity with depth. $Ps_1$ and $Ps_2$ are
241 manually demarcated, where Pms is the predicted result based on a crustal thickness of 55±0.2 km
242 below the station and Vp/Vs ratio of 1.83±0.02 from the H-k results. Poveda et al. (2015) also apply
243 H-k stacking to this station and found the Moho conversion to be at ~50±2 km depth with a Vp/Vs
244 ratio =1.77±0.07, suggesting our Pms likely results from the same discontinuity. Another mid-
245 crustal structure is suggested by Ps2, and H-k stacking indicates that this arrival is associated with a
246 maximum at ~38 km below the station with an average Vp/Vs ratio of ~1.81 above this
247 discontinuity (Fig. 2). The predicted multiples for the Ps2 are showing in Figure S2 along with the H-
248 k stacking results and do not appear to interfere with the Pms multiples.

249 Station ANTG is located 13 km west of Antisana volcano in the main arc (Fig. 1). Data quality
250 at ANTG is somewhat representative of other stations located near volcanoes, although there is
251 significant interstation variability. A strong negative conversion (Ps₁) is observed at ~2.5 seconds
252 that is consistent across most ray parameters. At ~8 seconds, a strong conversion which we
253 interpret as the Moho-converted phase (Pms), is continuous across all ray parameters. H-k stacking
254 at ANTG shows the Pms conversion time corresponds to a depth of ~55 ± 1 km below the station
255 with an average Vp/Vs value of 1.85 ± 0.02. Several other positive and negative peaks are present
256 but are not continuous across ray parameters. These other peaks may be a result of reverberations
257 from shallow structures, anisotropy, or noise resulting from a low number of RFs sampling specific
258 ray parameters being stacked.

259 A third example of RFs is provided from station PIAT which is located in the Eastern
260 Cordillera at ~1°S. RFs at PIAT show a strong positive conversion at ~8 seconds (Fig. 2). H-k
261 stacking finds this conversion to correspond to a depth of 55 ± 1 km below the station and a Vp/Vs
262 of 1.78 ± 0.02. Similar to ANTG, strong positive and negative conversions are observed irregularly
263 at a variety of ray parameters.
264 *Adaptive Common Conversion Point Stacks*

265 Arc-perpendicular and arc-parallel cross sections through the ACCP stacks are shown in Fig. 3. These sections have been corrected for elevation, so all depths are indicated relative to sea level unless otherwise stated. In the upper ~20 km of the ACCP stacks, the most prominent features are several strong negative conversions indicating decreases in velocity with depth (i.e., the top of a low-velocity zone). These conversions, labeled Ps₁ in Fig. 3, are prominent beneath several active volcanic centers, including Chiles-Cerro Negro, Cayambe, Antisana, and Tungurahua. In all sections, Ps₁ represents a decrease in velocities at relatively shallow depths but does not imply continuity throughout the study area. In A-A', the peak of Ps₁ is present at ~10 km beneath the Chiles-Cerro Negro region and appears shallower farther eastward at ~5 km beneath Soche Volcano (Fig. 3). The conversion beneath Chiles-Cerro Negro can also be seen in E-E', which follows the Western Cordillera, deepening to ~18 km. In B-B' the Ps₁ conversion is deeper, at ~20 km beneath Cayambe Volcano, and is not as prominent beneath the Western Cordillera. A strong positive conversion, Ps₂, parallels Ps₁ at ~28 km. In C-C' and D-D', the Ps₁ conversion can be seen at ~8 km beneath Antisana and ~6 km beneath Tungurahua respectively.

279 A strong semi-continuous positive conversion (Pms) that we interpret as the Moho is observed in all sections between depths of ~45 km and 65 km (thick solid and dashed black line in Fig. 3). This conversion is deepest in the north and beneath the Western Cordillera (A-A' and E-E'). There are several areas where this conversion weakens or appears to separate into multiple conversions at different depths (dashed segments in Fig. 3). Examples of this can be seen in B-B', C-C', and E-E'. These results suggest a complicated velocity structure that indicates incremental increases in velocity between the crust and mantle. Multiples of conversions in the upper 20 km of the crust may also interfere with a primary Moho conversion. Fig. 3 shows the location of multiples from shallow discontinuities along each section (blue and red lines). For example, along the southern half of section E-E', a multiple from a shallow discontinuity at ~10 km depth arrives at
approximately the same time as the Pms (Fig. 3). However, the amplitude and continuity of this arrival to the north and south lead us to interpret that the signal at this depth is primarily due to a primary conversion. Multiples from primary conversions deeper than ~20 km arrive beneath the Pms conversion and hence do not interfere with our interpretation of the Moho.

## Joint Inversion

Cross sections through the upper 75 km of our shear-wave velocity model are shown in Fig. 4. Although the full model extends to 100 km, the dispersion data lose sensitivity below ~50 km so deeper features are primarily controlled by the RFs, implying that the absolute velocities are weakly constrained. The locations of cross-sections are dependent on station distribution, coverage of the ACCP stacks, and areas where ambient noise dispersion data has good ray coverage (see Fig. 2 in Lynner et al., 2020). As in the ACCP stacks, depths are given relative to sea-level.

The northernmost profile, A-A', parallels the Ecuador-Colombian border extending from the Western Cordillera into the Eastern Cordillera. Two regions of slow shear velocities are observed in A-A'. We use the 3.2 km/s contour to estimate the size of the low-velocity zones (LVZs), as this velocity is generally lower than what would be expected from crystalline crustal material (Christensen, 1996). The first LVZ is located beneath Chiles-Cerro Negro volcanoes, where the LVZ is ~8 km thick and extends from the western edge of the profile to ~35 km east. A second region of slow velocities is observed in A-A' beneath the region of Soche volcano at depths ranging from ~6-13 km. This LVZ extends from the eastern edge of the profile ~42 km northwest. Beneath these two LVZs, the velocities rapidly increase, reaching 3.5 km/s at ~18-20 km, and remain in the 3.5-4.1 km/s range down to the Moho. One exception to this is beneath the LVZ in the Chiles-Cerro Negro region where the velocities reach > 4.4 km/s at a depth of ~40 km (A-A'). These anomalously fast velocities at ~40 km depth are likely an artifact resulting from the inversion attempting to fit multiples in the RFs that result from shallower structures (Fig. 3).
Farther south, at ~0.1°N, cross-section B-B' transects the Western Cordillera (from 313 Pichincha and Pululahua volcanoes) and extends into the Eastern Cordillera (Cayambe volcano; 314 Figs. 1 and 4). This profile shows slower shear velocities in the upper ~15 km across the profile. 315 These slow velocities extend deeper beneath the Eastern Cordillera reaching ~24 km depth. Below 316 the LVZ and above the Pms conversion, velocities range from 3.5 – 4.1 km/s and appear faster in the 317 Western Cordillera than in the Eastern Cordillera. To explore this further, we averaged the 318 velocities from the shear-wave velocity model in the Western Cordillera and the Eastern Cordillera 319 north of ~2°S (Fig. 5). The average shear velocities indicate that the Western Cordillera generally 320 has a slower shallow crust (<10 km) and faster middle and lower crust (>10 km) than the Eastern 321 Cordillera (Fig. 5A). 322

Cross-section C-C' trends roughly SW-NE from the Inter-Andean Valley into the Eastern 323 Cordillera and transects Cotopaxi and Antisana volcanoes. Crustal velocities are slower in this 324 region, only reaching >3.5 km/s below ~40 km depth in the SW and ~30 km depth towards the NE. 325 LVZs are present beneath Antisana at depths between ~8-15 km, with a shallower LVZ around 5-10 326 km. As in A-A', multiples beneath these LVZs may impact the joint inversion, and the fast velocity 327 feature ~30 km beneath Antisana and adjacent slower velocities observed at the NE end of the 328 profile may not be robust. 329

The southernmost shear velocity profile D-D' extends from the Western Cordillera into the 330 Eastern Cordillera and crosses the region of Tungurahua (Fig. 4). We observe a LVZ beneath 331 Tungurahua at depths between 7 and 12 km that covers a lateral extent of ~40 km. Below this LVZ, 332 slow velocities extend down to around ~22 km in depth. Multiples may be interfering with the rest 333 of the profile deeper than ~40 km but given the expected low amplitudes of the multiples, we 334 suggest that we can still identify the Moho (Fig. 3). Seismic velocities in the upper 25 km 335 surrounding the primary LVZ are ~3.5-3.8 km/s (Fig. 4). Depths between ~25-40 km on the eastern 336 side of the LVZ exhibit faster velocities (3.95 – 4.25 km/s). The fast anomaly may be influenced by 337
multiples in the region (Fig. 3). A similar increase in seismic velocity is not observed on the western side of the LVZ, where shear-wave velocities are at intermediate values (~3.5-3.8 km/s) down to ~50 km. At ~55 km, beneath the Pms conversion, a region of relatively slow velocities is observed.

Given the spatial correspondence to the LVZ to these apparent upper mantle low velocities, it is likely that this is an artifact due to RF multiples, which is frequently observed in volcanic areas (Bianchi et al., 2015).

Cross-section E-E' parallels the Western Cordillera following the volcanic arc from the Ecuador-Colombian border to the southern extent of the Northern Volcanic Zone at ~1.5°S (Fig. 4). This section shows clear along-strike variation in crustal velocities. In the north, the upper ~20 km show two prominent LVZs. The first is described previously in section A-A'. This LVZ extends 25 km southward from the northern edge of the profile E-E' at approximately the same depths as in A-A'. A second LVZ is observed further south, centered ~9 km NE of Cuicocha. South of Cuicocha, slow velocities (3.05-3.5 km/s) are pervasive throughout most of the section in the upper 18 km. These slow velocities continue to shallow to ~85 km southward until the edge of the model. The deeper crust appears generally faster along the northern ~120 km of the profile.

### Validation of low-velocity zones

We examine the LVZs in more detail to ensure that they are robust features. Rayleigh wave phase velocities are predicted to monotonically increase given an increasing velocity structure with depth. Thus, we inspect our dispersion curves in regions of these LVZs for a decrease in phase velocity with increasing periods, which necessitates a decrease in shear-wave velocity with depth assuming isotropic media. Figure S3 shows dispersion curves through the LVZs beneath Antisana (C-C') and Tungurahua (D-D'). Each of these dispersion curves show a decrease in velocity between 10 and 12 second periods, clearly indicating that a LVZ is present and required by the surface wave data. For the LVZs observed in B-B' beneath Cayambe, the dispersion curves do not show a clear decrease, however, the phase velocities are notably low, likely a result of the slow shear velocities
extending into the upper crust. Cross-sections through a Vs model created using only dispersion
data (i.e., no RF information included) also show a LVZ in the same region, suggesting that this is a
robust feature (see Figure S7). The LVZs in the region of Chiles-Cerro Negro and Soche are not
observed in the dispersion data or the dispersion-only Vs model, suggesting that these LVZs
primarily result from RF waveforms. This however does not imply that they are artifacts, as these
structures can be explained by low-velocity layers that are too thin to affect the broad sensitivities
of the dispersion data. Synthetic models of LVZs with varying thickness through a 50 km thick crust
are shown in Figure S8. These synthetics suggest that LVZs characterized by a Vs of 2.8 km/s with a
thickness less than ~9 km will not be clearly seen in dispersion data. Given that these are the
lowest velocities we observe, and the maximum observed thickness is ~10 km, it is not surprising
that a decrease in the dispersion curves is not observed.

These LVZs are also observed in our joint inversions that used a different starting model
(Supplemental Figures S4 and S5). The fact that the structure and velocities observed in these LVZs
are very similar independent of the starting model ensures that they are resolvable features and
not dependent on the starting model.

378 Discussion

379 *Crustal Thickness Beneath the Ecuador Arc*

380 The Pms conversion throughout our study region has a large amplitude, is continuous in
381 nature, and is consistent with previous studies (Feininger and Seguin, 1979; Poveda et al., 2015;
382 Araujo, 2016; Condori et al., 2017; Tamay et al., 2018; and Vaca et al., 2019). We, therefore,
383 interpret this Pms conversion as the crust-mantle boundary, where crustal rocks transition into the
384 peridotites that comprise the lithospheric mantle (Moho; Fig. 3). Due to low sensitivity to depths
385 greater than ~50 km, we do not interpret the nature of the mantle beneath the arc in this study.
A crustal thickness map derived from our ACCP stack is shown in Fig. 6. The depths of the conversion are tracked in the ACCP stack and interpolated across our study area. In regions lacking ACCP coverage, crustal thickness was calculated assuming the crust to be in isostatic equilibrium. For this calculation, we used a crustal density of 2.8 g/cm³ and a mantle density of 3.3 g/cm³. An optimal isostatic baseline of 36.5 km (corresponding to sea-level) was found by minimizing the difference between the observations and the isostatic calculations (Schmandt et al., 2015). These areas are shaded on the map in Fig. 6. We have good sampling of the Moho north of 2°S in the cordilleras but very little coverage between 2°S and 5°S due to sparse station coverage (Fig. 6). Where we do have reliable results from our ACCP stacks, the crust appears largely near isostatic equilibrium for a reasonable range of densities (see Fig. 7).

The thickest crust (~66 km) is observed in northern Ecuador in the regions near the Chiles-Cerro Negro and Soche volcanoes. Araujo (2016) found comparable thicknesses (65-70 km) in this region. In southern Colombia, Poveda et al., (2015) found a crustal thickness of ~56-58 km and Vp/Vs ratio of 1.74 from receiver function H-k stacking at station GCUF, ~50 km north of the Ecuador/Colombian border. In Ecuador, Vaca et al. (2019) estimated an overall thinner crust (averaging ~45 km beneath the arc), but still found some of the thickest crust in this northern region (~50-55 km). The average elevation in this region ranges from ~3200 - 3800 meters. If we apply a simple Airy isostatic calculation for crustal thickness between 62 and 67 km, this calculation suggests crustal densities of ~2.9 g/cm³. These densities are reasonable for the mafic terranes that are thought to make up the basement of the Western Cordillera (Jaillard et al., 2009; Vallejo et al., 2019).

At latitudes between ~0.25°S and ~2.0°S, the Western Cordillera exhibits crustal thicknesses between ~55 and 60 km while the Eastern Cordillera thicknesses are ~5-10 km thinner. Although the average elevation of the Western Cordillera is ~300 m lower than the Eastern Cordillera, this result is not entirely unexpected due to the presence of a denser, mafic, crust
411 beneath the Western Cordillera. Feininger and Seguin (1983) estimated a much thinner Western
412 Cordillera crust compared to the Eastern Cordillera crust using a gravity profile at ~2°S.

413 Conversely, Tamey et al. (2018) used gravity and estimated crustal thicknesses of up to 70 km
414 beneath both the Western and Eastern Cordilleras. Araujo (2016) found a slightly thicker Eastern
415 Cordillera ~1-2°S and a thicker Western Cordillera to the north. Vaca et al. (2019) estimated a
416 thicker Western Cordillera crust between ~1-2°S. South of 2°S stations are primarily located at
417 lower elevations on the edges of the cordilleras, and we observe correspondingly thinner crustal
418 thickness, ~41-48 km thick crust for stations on the edges of the cordilleras and between ~30 – 35
419 km thick crust for stations in the forearc and Oriente Basin.

420 Lateral variations in the Vp/Vs ratio or the velocity of the crust could explain some of the
421 variations in crustal thickness. Varying the Vp/Vs ratio used to migrate the CCP stack by ± 0.03
422 results in ±2 km shift for a 50 km conversion. The observations from average velocities in our joint
423 inversion (Fig. 5A) suggest an overall faster crust in the Western Cordillera versus the Eastern
424 Cordillera. Using a faster velocity model to migrate the RF data would increase the crustal
425 thickness, and a slower velocity model would decrease the crustal thickness. Thus, we find this
426 difference of thicker Western Cordillera crust to thinner Eastern Cordillera crust to be robust.

427 In multiple areas (dashed black line in Fig. 3) we observe a deep crustal conversion that lies
428 above our interpreted Moho (Pms) conversion (dashed black lines; Fig. 3). In these areas, we map
429 the crust-mantle transition as the deeper conversion based on its correspondence to where
430 velocities increase to >4.2 km/s in the joint inversion (Fig. 5B) although we note that the presence
431 of multiples may complicate the interpretation in some areas. The shallower conversion likely
432 represents the top of a distinct lower crust, possibly resulting from the tectonic or magmatic
433 underplating (Beck and Zandt, 2002; ANCORP Working Group, 2003; Muntener and Ulmer, 2006;
434 Frassetto et al., 2011). The velocities in these regions (~3.8-4.1 km/s) are consistent with a mafic
435 lower crust and appear to be distributed beneath several active volcanic systems.
437 Seismic wave speeds can be affected by many factors, such as composition, anisotropy,
438 temperature, and fluids. We approximate the expected isotropic velocities of the basement beneath
439 the arc using a program by Hacker et al. (2003), updated by Aber and Hackers (2016). For the
440 frontal arc, we assume a composition based on the mafic terranes that comprise the basement
441 (Table S1; Hughes and Pilatasig, 2002), whereas, for the main arc, where the basement is less well
442 constrained, we use the compositions of recent volcanic products (Table S1; Samaniego et al., 2011;
443 Nauret et al., 2018) which has been shown to approximately correspond with the bulk composition
444 of intrusive magmatism (Glazner et al., 2015). For the frontal arc (Western Cordillera), the
445 predicted basement velocity is ~3.75 km/s, and for the main arc (Eastern Cordillera) it is ~3.6
446 km/s. Observations of anisotropy of mid-crustal rocks in the Central Andes show maximum values
447 slightly higher than ~10%, with the strongest regions of anisotropy occurring beneath magmatic
448 centers (Lynner et al., 2018). This value also encompasses the anisotropic values from most
449 laboratory measurements in mid-lower crustal lithologies (Brownlee et al., 2017). 10% anisotropy
450 will reduce the isotropic velocity by ~5%, which would decrease our calculated values to ~3.55
451 km/s in the frontal arc and ~3.45 in the main arc, which we take as a lower estimate for the
452 expected velocities at depth in a melt-free system. As we observe shear velocities below ~3.2 km/s
453 beneath several active volcanic systems in our study area (Figure 4 and 8), these velocities indicate
454 the likely presence of melt in these systems. Previous studies have interpreted similar mid-crustal
455 LVZs below other volcanic systems as regions that contain a low percentage of melt (e.g., Ward et
456 al., 2014 and 2017; Kiser et al., 2016; Delph et al., 2017; Schmandt et al., 2019). These LVZs could
457 represent crystal-rich magmatic reservoirs (“mush” zones; Marsh, 1981) that act as long-term
458 magmatic accumulation and storage centers beneath their respective magmatic arc (Costa, 2008).
459 The tectonomagmatic environment of our study area indicates that this is the simplest explanation
460 for the LVZs within the Ecuadorian arc (Fig. 9).
The frontal arc LVZs underlie the potentially active volcanos Chiles-Cerro Negro, Yanaurcu-Chachimbiro, and Cotacachi-Cuicocha. Geodetic measurements near Chiles-Cerro Negro have suggested a source at depths greater than 9 km and less than 22.5 km below sea level (Ebmeier et al., 2016). These depths are consistent with our findings that suggest a reservoir between depths of ~7 – 15 km (Fig. 8A). There have been several instances of volcanic tremor swarms in the upper ~10 km in 2013 and 2014 at this system, (Ebmeier et al., 2016) and it is considered potentially active. South of ~0.1°N the upper ~20 km of the Western Cordillera shows relatively slow velocities that may indicate the presence of regional anisotropy or a very low percentage of diffuse partial melt in the upper crust. Petrologic studies have estimated amphibole crystallization at 5 km below the surface at Pichincha, likely corresponding to the upper part of the reservoir (Samaniego et al., 2010). Similar velocities are observed in the Altiplano low-velocity zone (ALVZ) in the Central Andes, which is characterized by a highly reflective mid-crust, negative conversions in RFs, and positive radial anisotropy (Yuan et al., 2000; Ward et al., 2016; Lynner et al., 2018). The ALVZ has been interpreted as a zone of ductile, perhaps flowing crust that may contain some partial melt. Given the active volcanism in the Western Cordillera and the history of accretion, it is possible that we are seeing a similar structure beneath the Western Cordillera.

Compared to the frontal arc, the main arc (mainly constrained to the Eastern Cordillera) shows several crystal mush zones beneath nearly all of the instrumented volcanic systems including Soche, Cayambe, Antisana, and Tungurahua (Figs 4, 8B, and 9). In contrast, we do not image a distinct mid-crustal mush zone beneath Cotopaxi (Fig. 8C), despite evidence for magma storage from geodetic modeling, deformational modeling, and geochemistry (Hickey et al., 2015; Morales Rivera et al., 2017; Martel et al., 2018; Saalfield et al., 2020). No estimates of depths of magma storage have been published for Soche, Cayambe, or Antisana volcanic centers, however, geochemical studies suggest an increased component of upper crustal contamination in the main arc compared to the frontal arc (Yin et al., 2020). Specifically, Ancellin et al., (2017) estimated the
percent upper- and lower-crustal (UC and LC) contribution to magmas throughout the arc, and
found the highest amounts of upper crustal contamination at Antisana (~1.8% UC, ~22% LC),
Cayambe (~1.2% UC and ~13% LC) and Tungurahua (~1.3% UC and ~27% LC) and lower amounts
of upper crustal contamination at Soche (~0.5% UC and ~9% LC), Cotopaxi (~0.8% UC and ~10%
LC), and Chiles-Cerro Negro (~0.2% UC and ~6% LC). This could be explained by replenishment of
mid-crustal magma reservoirs by rejuvenation events, potentially causing increased crustal
assimilation by extending the lifetime of the reservoir. It is difficult to know if this process is
ubiquitous across the arc, as these proxies for upper crustal contamination do not appear to show
systematic variations with the individual crystal mush zones that we interpret. On a broader scale,
however, Ancellin et al. (2017) found a general trend of increasing crustal contamination from west
to east and north to south throughout the arc, which is consistent with a general decrease in crustal
shear-wave velocities from north to south and west to east (Figures 4 and 5).

The higher shear velocities in the lower crust are consistent with a mafic lower crust
beneath Cayambe, Antisana, and Cotopaxi. The source and implications of a mafic lower crust in
these regions remain unclear. Beneath Cotopaxi, Garrison et al. (2011) suggested a fragment of the
Raspas terrane (a mafic terrane with oceanic plateau affinities) may be present in the lower crust.
Beneath Antisana and Cayambe some evidence of an oceanic crustal signature has been observed
(Samaniego et al., 2002; Bourdon et al., 2002; Samaniego et al., 2005, Garrison et al., 2011), which
suggests that the higher velocities we observe is an extension of this terrane. Alternately, magmatic
underplating of mafic melts fueling widespread assimilation and fractional crystallization may
deplete the lower crust of felsic minerals leaving behind a mafic lower crust (Müntener and Ulmer,
2018). While the geochemical signatures are consistent with a mafic lower crust, we cannot rule out
other processes, such as partial melt of a subducting slab as proposed by others (Bourdon et al.,
2002; Samaniego et al., 2002).
Tungurahua and Cotopaxi are well-studied volcanoes, and both have multiple lines of evidence for mid-crustal magma storage. For Tungurahua, one of the most active volcanic systems in Ecuador (Fig. 8B), a long-term storage zone was modeled at 11.5 km below sea level based on surface uplift measured from InSAR data (Champenois et al., 2014). Furthermore, phase equilibrium experiments on Tungurahua andesites suggest long-term storage of melt at depths 15-16 km beneath the edifice (~12 – 13 km below sea level) and a shorter-lived temporary storage between at 7.5 – 9.5 km below the summit (Andujar et al., 2017; Samaniego et al., 2011). Both of these are reasonably consistent with the reservoir we interpret at ~7-12 km (Figs. 4 and 8B). At Cotopaxi, magmatic storage regions have been suggested at depths ranging from 4-5 km below the summit (Hickey et al., 2015) and 3 – 17 km below sea level (Martel et al., 2018; Saalfield et al., 2019). Our results image a generally slow crust beneath Cotopaxi, however, no distinct LVZs are observed (Fig. 8C). The lack of such distinct LVZs may indicate that storage in these regions is shorter-lived or in smaller distributed volumes that are not detectable at our resolutions, or that Cotopaxi is laterally tapping into other magma storage systems, such as beneath Chalupas to the SW of Cotopaxi, which has been suggested by petrologic studies (Garrison et al., 2011).

The presence of magma reservoirs beneath the Ecuadorian arc is supported by petrologic, geodetic, and seismic observations. The seismic images we present add to the growing body of literature that images sizeable zones of mid-crustal magmatic storage. Figure 9 shows our preferred interpretation for several of the cross-sections presented in this study, highlighting the placement of these mid-crustal magma storage zones in the Ecuadorian arc. In the next section, we explore in more detail what are the plausible melt percentages within these LVZs.

## Estimating Melt Percentages

Multiple studies have used Gassmann's relations to estimate the percent of partial melt based on seismic velocities assuming equilibrated pore pressures (Chu et al., 2010; Ward et al., 2014; Kiser et al., 2018, and more). While useful, this method neglects possible anisotropy and
535 anelastic effects of high temperatures within magmatic systems. To approximate elastic properties
536 of the solid-melt system, we use petrologic constraints from basement rock and eruptions over the
537 past ~3000-4000 years at Pichincha and Tungurahua, located at the Western and Eastern
538 Cordilleras respectively (Table S1; Hughes and Pilatasig, 2002; Samaniego et al., 2010; Nauret et al.,
539 2018) using programs by Hacker et al., (2003; updated in Abers and Hacker, 2016) and Ukei and
540 Iwamori (2016). Different compositions of host rock and melt were tested and produced minor
541 variations in the resulting melt percent (Samaniego et al., 2010 and 2011). Pressure conditions are
542 approximated based on the depth of the LVZ (~250-350 MPa) and temperature conditions from
543 petrologic studies (Andujar et al., 2017; Samaniego et al., 2011).

544 As discussed in the previous section, the basement rock isotopic velocities in the frontal arc
545 (Western Cordillera) and main arc (Eastern Cordillera) are ~3.75 km/s and ~3.6 km/s respectively.
546 We then follow the approach of Chu et al., (2010) using the aforementioned assumptions to
547 calculate melt percentages for the velocities observed in our LVZs (Fig. 8D). To be conservative, we
548 also calculate melt percentages assuming we are sampling the slow axis of a ~10% anisotropic
549 medium with our methods. In this case, the slowest possible velocity for a melt free rock of the
550 appropriate composition is ~3.55 km/s for the Western Cordillera and ~3.45 km/s for the Eastern
551 Cordillera. In the frontal arc below Chiles-Cerro Negro, and the Cotacachi-Cuicocha volcanic
552 complex we observe are among the slowest in Ecuador, reaching lower than 2.8 km/s in the center
553 of the anomalies (Fig. 8), much slower than is expected for even the highly anisotropic case. This
554 represents up to a ~25% decrease in shear velocity relative to the isotropic velocity (~21% for the
555 anisotropic case) and corresponds to a maximum melt percentage in the LVZ of 14% (~12.5 for the
556 anisotropic case). In the main arc, Soche, Cayambe, Antisana, and Tungurahua all show their
557 slowest velocities ~2.9 km/s corresponding to maximum melt percentages of 13% (~11.5 for the
558 anisotropic case) (Fig. 8D).
The melt percentages we estimate further support the interpretation of the LVZs in the Ecuadorian arc as highly crystalline magmatic reservoirs consistent with crystal mush zones (Fig. 9). These are thought to be non-eruptible bodies due to their high crystallinity state, which renders the bulk system immobile (Marsh, 1981). However, they can allow for magmatic accumulation, geochemical hybridization, and melt transport to shallower depth within the mush zone through percolative flow (Jackson et al., 2018), which may eventually make its way to the surface. It is important to note that these crystal mush zones likely result from the long-term (tens of thousands to millions of years; Petrelli et al., 2020) amalgamation of magma within the crust of an arc and are likely characterized by amalgamated plutonic bodies rather than representing an individually emplaced plutonic body. Similar crystal mush zones have been interpreted beneath a variety of systems, both in other geologic systems (Kiser et al., 2016; Paulatto et al., 2019; Chen et al., 2020) and elsewhere in the Andes (Ward et al., 2017; Wespestad et al., 2019). Our findings provide further evidence that the mid crust acts as a zone of long-term magmatic storage and accumulation.

## Conclusions

Our study examines the crustal structure of the Ecuadorian Arc using adaptive common conversion point receiver function stacks and a 3D shear-wave velocity model derived from the joint inversion of receiver functions and surface wave dispersion data. These new images show that crustal thicknesses range from ~50-65 km beneath this portion of the Andean arc, with the Western Cordillera being ~5-10 km thicker than the Eastern Cordillera. The Western Cordillera becomes faster at mid-to-lower crustal depths, in agreement with evidence for a mafic basement beneath the Western Cordillera versus a silicic basement in the Eastern Cordillera. This likely implies a variation in bulk crustal density as well, allowing the 5-10 km thicker crust of the Western Cordillera to be at a slightly lower elevation while remaining in isostatic balance. We also observe that the mid-crust beneath several volcanic systems in the Ecuadorian arc is characterized by low
583 velocity zones at ~5 - 20 km below sea level. Below the well-studied Chiles-Cerro Negro and
584 Tungurahua volcanic centers, the depth extent of these low velocity zones is consistent with
585 previous estimates of magma storage from petrologic and geodetic studies. We also image low
586 velocity zones beneath other volcanic centers that have not been studied in as much detail and may
587 merit further investigation. The lowest velocities we observe in the mid-crust reach ~2.8 km/s and
588 correspond to a maximum of ~14% melt, indicating that these low velocity zones likely represent
589 predominantly crystalline magma reservoirs consistent with mush zones. These mush zones likely
590 result from the long-term (tens of thousands to millions of years) accumulation and storage of
591 magmas in the mid-crust. Our results provide new insights into the crustal structure of the
592 Ecuadorean Andes and the magmatic architecture of this broad portion of the Andean arc.

593 Acknowledgements

594 We used data from the Instituto Geofísico at the Escuela Politécnica Nacional (IG-EPN)
595 seismic network and the Global Seismic Network. Open and available seismic data (Albuquerque
596 Seismological Laboratory (ASL)/USGS, 1988; Alvarado et al., 2018) were accessed via the Data
597 Management Center of the Incorporated Research Institutions for Seismology. Additional data was
598 obtained directly from the IG-EPN. We thank the dedicated IG-EPN staff that operate and maintain
599 the EC network. This work was supported by NSF EAR-1723065, NSF EAR-1723042, and NSF
600 EAR-1565475. We used the Generic Mapping Tools developed by Wessel and Smith (1998) to
601 make figures.

602 References

603 Abers, G.A., Hacker, B.R., 2016. A MATLAB toolbox and Excel workbook for calculating the densities,
604 seismic wave speeds, and major element composition of minerals and rocks at pressure and
605 temperature. Geochemistry, Geophysics, Geosystems 17, 616-624.
606 Albuquerque Seismological Laboratory (ASL)/USGS (1988). Global Seismograph Network (GSN - IRIS/USGS). https://doi.org/10.7914/SN/IU

608 Alvarado, A., Ruiz, M., Mothes, P., Yepes, H., Segovia, M., Vaca, M., Ramos, C., Enríquez, W., Ponce, G., Jarrín, P., 2018. Seismic, volcanic, and geodetic networks in Ecuador: Building capacity for monitoring and research. Seismological Research Letters 89, 432-439.

611 Ammon, C.J., 1991. The isolation of receiver effects from teleseismic P waveforms. Bulletin of the seismological Society of America 81, 2504-2510.

613 Ancellin, M.A., Samaniego, P., Vlastélic, I., Nauret, F., Gannoun, A., Hidalgo, S., 2017. Across-arc versus along-arc S r-Nd-Pb isotope variations in the E cuadorian volcanic arc. Geochemistry, Geophysics, Geosystems 18, 1163-1188.

616 Andújar, J., Martel, C., Pichavant, M., Samaniego, P., Scaillet, B., Molina, I., 2017. Structure of the plumbing system at Tungurahua volcano, Ecuador: insights from phase equilibrium experiments on July-August 2006 eruption products. Journal of Petrology 58, 1249-1278.

619 Araujo, S.n., 2016. Travel time tomography of the crust and the mantle beneath Ecuador from data of the national seismic network. Grenoble Alpes.

621 Aspden, J.A., Litherland, M., 1992. The geology and Mesozoic collisional history of the Cordillera Real, Ecuador. Tectonophysics 205, 187-204.

623 Bablon, M., Quidelleur, X., Samaniego, P., Le Pennec, J.-L., Santamaria, S., Liorzou, C., Hidalgo, S., Eschbach, B., 2020. Volcanic history reconstruction in northern Ecuador: insights for eruptive and erosion rates on the whole Ecuadorian arc. Bulletin of Volcanology 82, 11.

626 Barragan, R., Geist, D., Hall, M., Larson, P., Kurz, M., 1998. Subduction controls on the compositions of lavas from the Ecuadorian Andes. Earth and Planetary Science Letters 154, 153-166.
628 Beck, S., Zandt, G., Ward, K.M., Scire, A., 2015. Multiple styles and scales of lithospheric foundering
629 beneath the Puna Plateau, central Andes, Memoir of the Geological Society of America. Geological
630 Society of America, pp. 43-60.

631 Beck, S.L., Zandt, G., 2002. The nature of orogenic crust in the central Andes. Journal of Geophysical
632 Research: Solid Earth 107, ESE 7-1-ESE 7-16.

633 Bishop, B.T., Beck, S.L., Zandt, G., Wagner, L., Long, M., Antonijevic, S.K., Kumar, A., Tavera, H., 2017.
634 Causes and consequences of flat-slab subduction in southern Peru. Geosphere 13, 1392-1407.

635 Bourdon, E., Eissen, J.-P., Gutscher, M.-A., Monzier, M., Hall, M.L., Cotten, J., 2003. Magmatic response
636 to early aseismic ridge subduction: the Ecuadorian margin case (South America). Earth and
637 Planetary Science Letters 205, 123-138.

638 Bourdon, E., Eissen, J.-P., Monzier, M., Robin, C., Martin, H., Cotten, J., Hall, M.L., 2002. Adakite-like
639 lavas from Antisana Volcano (Ecuador): evidence for slab melt metasomatism beneath Andean
640 Northern Volcanic Zone. Journal of Petrology 43, 199-217.

641 Brownlee, S. J., Schulte-Pelkum, V., Raju, A., Mahan, K., Condit, C., & Orlandini, O. F. (2017).

642 Characteristics of deep crustal seismic anisotropy from a compilation of rock elasticity tensors and
643 their expression in receiver functions. Tectonics, 36(9), 1835-1857.

644 Bryant, J., Yogodzinski, G., Hall, M., Lewicki, J., Bailey, D., 2006. Geochemical constraints on the origin
645 of volcanic rocks from the Andean Northern Volcanic Zone, Ecuador. Journal of petrology 47, 1147-
646 1175.

647 Cashman, K.V., Sparks, R.S.J., Blundy, J.D., 2017. Vertically extensive and unstable magmatic
648 systems: a unified view of igneous processes. Science 355.
649 Champenois, J., Pinel, V., Baize, S., Audin, L., Jomard, H., Hooper, A., Alvarado, A., Yepes, H., 2014.

650 Large-scale inflation of Tungurahua volcano (Ecuador) revealed by Persistent Scatterers SAR

651 interferometry. Geophysical Research Letters 41, 5821-5828.

652 Chaussard, E., Amelung, F., 2014. Regional controls on magma ascent and storage in volcanic arcs.

653 Geochemistry, geophysics, geosystems 15, 1407-1418.

654 Chiaradia, M., 2009. Adakite-like magmas from fractional crystallization and melting-assimilation of

655 mafic lower crust (Eocene Macuchi arc, Western Cordillera, Ecuador). Chemical Geology 265, 468-

656 487.

657 Chiaradia, M., Müntener, O., Beate, B., 2020. Effects of aseismic ridge subduction on the

658 geochemistry of frontal arc magmas. Earth and Planetary Science Letters 531, 115984.

659 Christensen, N.I., 1996. Poisson's ratio and crustal seismology. Journal of Geophysical Research:

660 Solid Earth 101, 3139-3156.

661 Christensen, N.I., Mooney, W.D., 1995. Seismic velocity structure and composition of the continental

662 crust: A global view. Journal of Geophysical Research: Solid Earth 100, 9761-9788.

663 Chu, R., Helmberger, D.V., Sun, D., Jackson, J.M., Zhu, L., 2010. Mushy magma beneath Yellowstone.

664 Geophysical Research Letters 37.

665 Condori, C., França, G.S., Tavera, H.J., Albuquerque, D.F., Bishop, B.T., Beck, S.L., 2017. Crustal

666 structure of north Peru from analysis of teleseismic receiver functions. Journal of South American

667 Earth Sciences 76, 11-24.

668 Delph, J.R., Ward, K.M., Zandt, G., Ducea, M.N., Beck, S.L., 2017. Imaging a magma plumbing system

669 from MASH zone to magma reservoir. Earth and Planetary Science Letters 457, 313-324.
670 Delph, J.R., Zandt, G., Beck, S.L., 2015. A new approach to obtaining a 3D shear wave velocity model of the crust and upper mantle: An application to eastern Turkey. Tectonophysics 665, 92-100.

672 Dueker, K.G., Sheehan, A.F., 1997. Mantle discontinuity structure from midpoint stacks of converted P to S waves across the Yellowstone hotspot track. Journal of Geophysical Research: Solid Earth 102, 8313-8327.

675 Ebmeier, S.K., Elliott, J.R., Nocquet, J.-M., Biggs, J., Mothes, P., Jarrín, P., Yépez, M., Aguaiza, S., Lundgren, P., Samsonov, S.V., 2016. Shallow earthquake inhibits unrest near Chiles-Cerro Negro volcanoes, Ecuador-Colombian border. Earth and Planetary Science Letters 450, 283-291.

678 Feininger, T., Seguin, M., 1983. Simple Bouguer gravity anomaly field and the inferred crustal structure of continental Ecuador. Geology 11, 40-44.

680 Frassetto, A.M., Zandt, G., Gilbert, H., Owens, T.J., Jones, C.H., 2011. Structure of the Sierra Nevada from receiver functions and implications for lithospheric foundering. Geosphere 7, 898-921.

682 Garrison, J.M., Davidson, J.P., Hall, M., Mothes, P., 2011. Geochemistry and petrology of the most recent deposits from Cotopaxi Volcano, Northern Volcanic Zone, Ecuador. Journal of Petrology 52, 1641-1678.

685 Gilbert, H., Beck, S., Zandt, G., 2006. Lithospheric and upper mantle structure of central Chile and Argentina. Geophysical Journal International 165, 383-398.

687 Glazner, A.F., 2015. Magmatic life at low Reynolds number: REPLY. Geology 43, e359-e359.

688 Group, A.W., 2003. Seismic imaging of a convergent continental margin and plateau in the central Andes (Andean Continental Research Project 1996 (ANCORP'96)). J. geophys. Res 108, 2328.

690 Hacker, B.R., Kelemen, P.B., Behn, M.D., 2011. Differentiation of the continental crust by relamination. Earth and Planetary Science Letters 307, 501-516.
692 Hall, M.L., Samaniego, P., Le Pennec, J.-L., Johnson, J., 2008. Ecuadorian Andes volcanism: A review of
693 Late Pliocene to present activity. Journal of Volcanology and Geothermal Research 176, 1-6.

694 Hickey, J., Gottsmann, J., Mothes, P., 2015. Estimating volcanic deformation source parameters with
695 a finite element inversion: The 2001–2002 unrest at Cotopaxi volcano, Ecuador. Journal of
696 Geophysical Research: Solid Earth 120, 1473-1486.

697 Hidalgo, S., Gerbe, M.C., Martin, H., Samaniego, P., Bourdon, E., 2012. Role of crustal and slab
698 components in the Northern Volcanic Zone of the Andes (Ecuador) constrained by Sr-Nd-O
699 isotopes. Lithos 132, 180-192.

700 Huang, H.-H., Lin, F.-C., Schmandt, B., Farrell, J., Smith, R.B., Tsai, V.C., 2015. The Yellowstone
701 magmatic system from the mantle plume to the upper crust. Science 348, 773-776.

702 Hughes, R.A., Pilatasig, L.F., 2002. Cretaceous and Tertiary terrane accretion in the Cordillera
703 Occidental of the Andes of Ecuador. Tectonophysics 345, 29-48.

704 Jackson, M., Blundy, J., Sparks, R., 2018. Chemical differentiation, cold storage and remobilization of
705 magma in the Earth's crust. Nature 564, 405-409.

706 Jaillard, E., Lapierre, H., Ordonez, M., Alava, J.T., Amortegui, A., Vanmelle, J., 2009. Accreted oceanic
707 terranes in Ecuador: southern edge of the Caribbean Plate? Geological Society, London, Special
708 Publications 328, 469-485.

709 Jaillard, É., Ordoñez, M., Benitez, S., Berrones, G., Jiménez, N., Montenegro, G., Zambrano, I., 1995.

710 Basin development in an accretionary, oceanic-floored fore-arc setting: southern coastal Ecuador
711 during late Cretaceous-late Eocene time.

712 Julia, J., Ammon, C., Herrmann, R., Correig, A.M., 2000. Joint inversion of receiver function and
713 surface wave dispersion observations. Geophysical Journal International 143, 99-112.
714 Kiser, E., Levander, A., Zelt, C., Schmandt, B., 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 46, 775-778.

717 Kiser, E., Palomeras, I., Levander, A., Zelt, C., Harder, S., Schmandt, B., Hansen, S., Creager, K., Ulberg, C., 2016. Magma reservoirs from the upper crust to the Moho inferred from high-resolution Vp and Vs models beneath Mount St. Helens, Washington State, USA. Geology 44, 411-414.

720 Koch, C.D., Lynner, C., Delph, J., Beck, S.L., Meltzer, A., Font, Y., Soto-Cordero, L., Hoskins, M., Stachnik, J.C., Ruiz, M., 2020. Structure of the Ecuadorian forearc from the joint inversion of receiver functions and ambient noise surface waves. Geophysical Journal International.

723 Langston, C.A., 1979. Structure under Mount Rainier, Washington, inferred from teleseismic body waves. Journal of Geophysical Research: Solid Earth 84, 4749-4762.

725 Lavenu, A., Winter, T., Dávila, F., 1995. A Pliocene-Quaternary compressional basin in the Interandean Depression, central Ecuador. Geophysical Journal International 121, 279-300.

727 Lees, J.M., 2007. Seismic tomography of magmatic systems. Journal of Volcanology and Geothermal Research 167, 37-56.

729 Ligorria, J.P., Ammon, C.J., 1999. Iterative deconvolution and receiver-function estimation. Bulletin of the seismological Society of America 89, 1395-1400.

731 Luzieux, L., Heller, F., Spikings, R., Vallejo, C., Winkler, W., 2006. Origin and Cretaceous tectonic history of the coastal Ecuadorian forearc between 1 N and 3 S: Paleomagnetic, radiometric and fossil evidence. Earth and Planetary Science Letters 249, 400-414.

734 Lynner, C., Beck, S.L., Zandt, G., Porritt, R.W., Lin, F.C., Eilon, Z.C., 2018. Midcrustal deformation in the Central Andes constrained by radial anisotropy. Journal of Geophysical Research: Solid Earth 123, 4798-4813.
737 Lynner, C., Koch, C., Beck, S.L., Meltzer, A., Soto-Cordero, L., Hoskins, M.C., Stachnik, J.C., Ruiz, M.,

738 Alvarado, A., Charvis, P., 2020. Upper-plate structure in Ecuador coincident with the subduction of

739 the Carnegie Ridge and the southern extent of large mega-thrust earthquakes. Geophysical Journal

740 International 220, 1965-1977.

741 Marsh, B., 1981. On the crystallinity, probability of occurrence, and rheology of lava and magma.

742 Contributions to Mineralogy and Petrology 78, 85-98.

743 Martel, C., Andújar, J., Mothes, P., Scaillet, B., Pichavant, M., Molina, I., 2018. Storage conditions of the

744 mafic and silicic magmas at Cotopaxi, Ecuador. Journal of Volcanology and Geothermal Research

745 354, 74-86.

746 Meltzer, A., Beck, S., Ruiz, M., Hoskins, M., Soto-Cordero, L., Stachnik, J.C., Lynner, C., Porritt, R.,

747 Portner, D., Alvarado, A., 2019. The 2016 Mw 7.8 Pedernales, Ecuador, earthquake: rapid response

748 deployment. Seismological Research Letters 90, 1346-1354.

749 Morales Rivera, A. M., Amelung, F., Mothes, P., Hong, S. H., Nocquet, J. M., & Jarrin, P. (2017). Ground

750 deformation before the 2015 eruptions of Cotopaxi volcano detected by InSAR. Geophysical

751 Research Letters, 44(13), 6607-6615.

752 Müntener, O., Ulmer, P., 2006. Experimentally derived high-pressure cumulates from hydrous arc

753 magmas and consequences for the seismic velocity structure of lower arc crust. Geophysical

754 Research Letters 33.

755 Narvaez, D.F., Rose-Koga, E.F., Samaniego, P., Koga, K.T., Hidalgo, S., 2018. Constraining magma

756 sources using primitive olivine-hosted melt inclusions from Puñalica and Sangay volcanoes

757 (Ecuador). Contributions to Mineralogy and Petrology 173, 80.
758 Nauret, F., Samaniego, P., Ancellin, M.-A., Tournigand, P.-Y., Le Pennec, J.-L., Vlastélic, I., Gannoun, A.,
759 Hidalgo, S., Schiano, P., 2018. The genetic relationship between andesites and dacites at
760 Tungurahua volcano, Ecuador. Journal of Volcanology and Geothermal Research 349, 283-297.

761 Nur, A., Ben-Avraham, Z., 1983. Volcanic gaps due to oblique consumption of aseismic ridges.
762 Tectonophysics 99, 355-362.

763 Paulatto, M., Moorkamp, M., Hautmann, S., Hooft, E., Morgan, J.V., Sparks, R.S.J., 2019. Vertically
764 extensive magma reservoir revealed from joint inversion and quantitative interpretation of seismic
765 and gravity data. Journal of Geophysical Research: Solid Earth 124, 11170-11191.

766 Petrelli, M., Zellmer, G.F., 2020. Rates and timescales of magma transfer, storage, emplacement, and
767 eruption. Dynamic Magma Evolution, 1-41.

768 Poveda, E., Monsalve, G., Vargas, C.A., 2015. Receiver functions and crustal structure of the
769 northwestern Andean region, Colombia. Journal of Geophysical Research: Solid Earth 120, 2408-
770 2425.

771 Pratt, W.T., Duque, P., Ponce, M., 2005. An autochthonous geological model for the eastern Andes of
772 Ecuador. Tectonophysics 399, 251-278.

773 Rodriguez, E., Russo, R., 2020. Southern Chile crustal structure from teleseismic receiver functions:
774 Responses to ridge subduction and terrane assembly of Patagonia. Geosphere 16, 378-391.

775 Ryan, J., Beck, S., Zandt, G., Wagner, L., Minaya, E., Tavera, H., 2016. Central Andean crustal structure
776 from receiver function analysis. Tectonophysics 682, 120-133.

777 Saalfeld, M.A., Kelley, D.F., Panter, K.S., 2019. Insight on magma evolution and storage through the
778 recent eruptive history of Cotopaxi Volcano, Ecuador. Journal of South American Earth Sciences 93,
779 85-101.
780 Samaniego, P., Le Pennec, J.-L., Robin, C., Hidalgo, S., 2011. Petrological analysis of the pre-eruptive
781 magmatic process prior to the 2006 explosive eruptions at Tungurahua volcano (Ecuador). Journal
782 of Volcanology and Geothermal Research 199, 69-84.

783 Samaniego, P., Martin, H., Monzier, M., Robin, C., Fornari, M., Eissen, J.-P., Cotten, J., 2005. Temporal
784 evolution of magmatism in the Northern Volcanic Zone of the Andes: the geology and petrology of
785 Cayambe Volcanic Complex (Ecuador). Journal of petrology 46, 2225-2252.

786 Samaniego, P., Martin, H., Robin, C., Monzier, M., 2002. Transition from calc-alkalic to adakitic
787 magmatism at Cayambe volcano, Ecuador: insights into slab melts and mantle wedge interactions.
788 Geology 30, 967-970.

789 Samaniego, P., Robin, C., Chazot, G., Bourdon, E., Cotten, J., 2010. Evolving metasomatic agent in the
790 Northern Andean subduction zone, deduced from magma composition of the long-lived Pichincha
791 volcanic complex (Ecuador). Contributions to Mineralogy and Petrology 160, 239-260.

792 Schmandt, B., Jiang, C., Farrell, J., 2019. Seismic perspectives from the western US on magma
793 reservoirs underlying large silicic calderas. Journal of Volcanology and Geothermal Research 384,
794 158-178.

795 Schmandt, B., Lin, F.C., Karlstrom, K.E., 2015. Distinct crustal isostasy trends east and west of the
796 Rocky Mountain Front. Geophysical Research Letters 42, 10,290-210,298.

797 Spikings, R., Cochrane, R., Villagomez, D., Van der Lelij, R., Vallejo, C., Winkler, W., Beate, B., 2015.
798 The geological history of northwestern South America: from Pangaea to the early collision of the
799 Caribbean Large Igneous Province (290-75 Ma). Gondwana Research 27, 95-139.

800 Spikings, R.A., Winkler, W., Hughes, R., Handler, R., 2005. Thermochronology of allochthonous
801 terranes in Ecuador: Unravelling the accretionary and post-accretionary history of the Northern
802 Andes. Tectonophysics 399, 195-220.
803 Tamay, J., Galindo-Zaldívar, J., Martos, Y., Soto, J., 2018. Gravity and magnetic anomalies of ecuadorian margin: Implications in the deep structure of the subduction of Nazca Plate and Andes
804
805 Cordillera. Journal of South American Earth Sciences 85, 68-80.

806 Tassara, A., Echaurren, A., 2012. Anatomy of the Andean subduction zone: three-dimensional density model upgraded and compared against global-scale models. Geophysical Journal
807
808 International 189, 161-168.

809 Ueki, K., Iwamori, H., 2016. Density and seismic velocity of hydrous melts under crustal and upper mantle conditions. Geochemistry, Geophysics, Geosystems 17, 1799-1814.

811 Vaca, S., Vallée, M., Nocquet, J.-M., Alvarado, A., 2019. Active deformation in Ecuador enlightened by a new waveform-based catalog of earthquake focal mechanisms. Journal of South American Earth Sciences 93, 449-461.

814 Vallejo, C., Spikings, R.A., Horton, B.K., Luzieux, L., Romero, C., Winkler, W., Thomsen, T.B., 2019. Late Cretaceous to Miocene stratigraphy and provenance of the coastal forearc and Western Cordillera of Ecuador: Evidence for accretion of a single oceanic plateau fragment, Andean Tectonics. Elsevier, pp. 209-236.

818 Ward, K.M., Delph, J.R., Zandt, G., Beck, S.L., Ducea, M.N., 2017. Magmatic evolution of a Cordilleran flare-up and its role in the creation of silicic crust. Scientific reports 7, 1-8.

820 Ward, K.M., Zandt, G., Beck, S.L., Christensen, D.H., McFarlin, H., 2014. Seismic imaging of the magmatic underpinnings beneath the Altiplano-Puna volcanic complex from the joint inversion of surface wave dispersion and receiver functions. Earth and Planetary Science Letters 404, 43-53.

823 Wespestad, C.E., Thurber, C.H., Andersen, N.L., Singer, B.S., Cardona, C., Zeng, X., Bennington, N.L., Keranen, K., Peterson, D.E., Cordell, D., 2019. Magma reservoir below Laguna del Maule volcanic
825 field, Chile, imaged with surface-wave tomography. Journal of Geophysical Research: Solid Earth
826 124, 2858-2872.

827 Yin, S., Ma, C., Xu, J., 2020. Recycling of K-feldspar antecrysts in the Baishiya porphyritic
828 granodiorite, East Kunlun orogenic belt, northern Tibet Plateau: Implications for magma
differentiation in a crystal mush reservoir. Lithos, 105622.
829

830 Yuan, X., Sobolev, S.V., Kind, R., Oncken, O., Bock, G., Asch, G., Schurr, B., Graeber, F., Rudloff, A.,
831 Hanka, W., 2000. Subduction and collision processes in the Central Andes constrained by converted
832 seismic phases. Nature 408, 958-961.

833 Zhu, L., Kanamori, H., 2000. Moho depth variation in southern California from teleseismic receiver
834 functions. Journal of Geophysical Research: Solid Earth 105, 2969-2980.

835
836 **Figures**

<figure><img src="image_2.png" /></figure>

837

**Figure 1.** Topographic map showing the major tectonic features, seismic stations (inverted triangles) and Holocene volcanic centers (red triangles) within our study area. Seismic stations used for receiver functions only are shown in inverted black triangles. Seismic stations used by Lynner et al. (2020) for ambient noise dispersion data are shown as small pink triangles and seismic stations used for both are shown as pink triangles with a black outline. Inset map in upper left corner shows the location of the map. Tectonic provinces are labeled or colored in red (Western Cordillera), yellow (Inter-Andean Valley), green (Eastern Cordillera), and blue (Subandean Zone). Dark blue lines indicate the Calacali Puijili Fault Zone (CPFZ) and the Peltetec Fault (PF). The red box indicates the regions of maps shown in Figs. 3, 4, and 5.

847
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Examples of RFs plotted as a function of ray parameter from representative stations in Ecuador. Station OTAV is in the Western Cordillera (left). Station ANTG is near Antisana volcanic center in the Eastern Cordillera (middle). Station PIAT is in the Eastern Cordillera (right). RF ray parameter plots are stacked into bins of .0025 s/km. Black lines are manually picked P-to-S conversions. Red, cyan, and blue lines show the phases calculated using results from H-k stacking; Ps (red), PpPs (dashed cyan), and PsPs + PpSs (dashed blue).</figcaption></figure>
<figure><img src="image_4.png" /></figure>

**Figure 3.** Cross-sections through the ACCP receiver function stack. Parameters and the velocity model used to calculate the ACCP stacks are described in the text. Map on the left shows the location of the profiles with seismic stations (black inverted triangles) and volcanos (red triangles). A topographic profile is plotted above each section with stations and volcanos projected from ± 20 km on either side of the profile. Green lines show strong Ps conversions in the upper ~20 km and the red and blue lines are the predicted location of the multiples in the velocity model (the color of the line corresponds to the predicted amplitude of the multiple; red = positive, blue = negative). Thick black line indicates the interpreted Moho discontinuity in each profile. The black dashed line is a shallower discontinuity most likely in the lower crust.
<figure><img src="image_5.png" /></figure>

**Figure 4.** Cross-sections through the Joint Inversion (JI) shear-wave velocity model. Map on the left shows the location of the profiles with seismic stations (black inverted triangles) and volcanos (red triangles). Each cross-section shows the JI fit to the RF data (red square) and surface wave dispersion data (SW; blue square) for the final iteration of the model at that grid point. A topographic profile is plotted above each section with stations and volcanos projected from ± 20 km on either side of the profile. The black line and dashed black line are described in Fig. 3. The contour values are indicated in the color bar, with the 3.2 km/s contour in bold.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (A) Average shear velocity through the Western (blue) and Eastern (red) Cordillera from the joint inversion. Standard deviations for each cordillera are shown in the shaded region. The areas that are averaged for each cordillera are the outlines in Fig. 1 north of ~1.5°S. (B) Average velocity profiles through three different regions showing variability in the Moho and lower crust. The sharp Moho is taken from stations in the Western Cordillera between ~0.7°S and ~1.6°S.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Crustal thickness map derived from the ACCP stack. Shaded areas show regions that were calculated from topography assuming isostatic equilibrium (averaging the topography over 0.2°). Seismic stations are shown in black inverted triangles and volcanoes are shown in red triangles. The crustal thickness contour interval is 5 km.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Graph of crustal thickness plotted against the topography averaged over a 0.2° area centered at each seismic station. Lines indicate the predicted crustal thickness using Airy isostasy for different crustal densities (blue = 2.6 g/cm³; green = 2.8 g/cm³, black = 3.0 g/cm³), a mantle density of 3.3 g/cm³, and an isostatic baseline of 36.5 km (crustal thickness at sea-level) with the compensation depth at the Moho. Different colors indicate stations located in different tectonic provinces. The largest source of uncertainty is lateral variations in Vp/Vs which result in ±2 km at 50 km for variations of ±0.03 Vp/Vs.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** 1-D velocity models located at the bins containing (A) Chiles-Cerro Negro, (B)</figcaption></figure>

Tungurahua, and (C) Cotopaxi. Red shaded region indicates the upper and lower limits where we expect partial melt (velocities < 3.2 km/s) from our joint inversion model. Red stars and brackets represent the depths at which other studies have suggested magma storage. References are (1)

Ebmeier et al. (2016), (2) Andujar et al. (2017), (3) Champenois et al, (2014); (4) Samaniego et al. (2011), (5) Rivera et al. (2017); (6) Saalfield et al. (2020), (7) Martel et al. (2018), (8) Hickey et al.

(2015). (D) Partial melt as a function of the isotropic (solid line) and anisotropic (dashed line) shear velocity calculated for the compositions from the frontal (black) and main (green) arc (Table S1).

Red shaded region indicates increasing melt percentage as velocity decreases up to the lowest velocities we observe (2.8 km/s).
<figure><img src="image_10.png" /><figcaption>Figure 9. Interpreted cross-sections B-B', C-C', and E-E' through the joint inversion shear velocity model. Polygons indicating the accreted mafic terranes (black), mafic lower crust (blue), South American continental crust (green), slow velocity upper crust (orange), and regions of partial melt (crystal mush; red). We do not interpret below the crust-mantle boundary due to lack of resolution of the dispersion data.</figcaption></figure>