Volume 10 - 2022 | https://doi.org/10.3389/feart.2022.970131
Advances in seismic imaging of magma and crystal mush
- 1Department of Earth Science and Engineering, Imperial College London, London, United Kingdom
- 2Department of Earth Sciences, University of Oregon, Eugene, OR, United States
- 3National Tsunami Warning Center, Palmer, AK, United States
Seismic imaging methods have provided detailed three-dimensional constraints on the physical properties of magmatic systems leading to invaluable insight into the storage, differentiation and dynamics of magma. These constraints have been crucial to the development of our modern understanding of magmatic systems. However, there are still outstanding knowledge gaps resulting from the challenges inherent in seismic imaging of volcanoes. These challenges stem from the complex physics of wave propagation across highly heterogeneous low-velocity anomalies associated with magma reservoirs. Ray-based seismic imaging methods such as travel-time and surface-wave tomography lead to under-recovery of such velocity anomalies and to under-estimation of melt fractions. This review aims to help the volcanologist to fully utilize the insights gained from seismic imaging and account for the resolution limits. We summarize the advantages and limitations of the most common imaging methods and propose best practices for their implementation and the quantitative interpretation of low-velocity anomalies. We constructed and analysed a database of 277 seismic imaging studies at 78 arc, hotspot and continental rift volcanoes. Each study is accompanied by information about the seismic source, part of the wavefield used, imaging method, any detected low-velocity zones, and estimated melt fraction. Thirty nine studies attempted to estimate melt fractions at 22 different volcanoes. Only five studies have found evidence of melt storage at melt fractions above the critical porosity that separates crystal mush from mobile magma. The median reported melt fraction is 13% suggesting that magma storage is dominated by low-melt fraction crystal mush. However, due to the limits of seismic resolution, the seismological evidence does not rule out the presence of small (<10 km3) and medium-sized (<100 km3) high-melt fraction magma chambers at many of the studied volcanoes. The combination of multiple tomographic imaging methods and the wider adoption of methods that use more of the seismic wavefield than the first arriving travel-times, promise to overcome some of the limitations of seismic tomography and provide more reliable constraints on melt fractions. Wider adoption of these new methods and advances in data collection are needed to enable a revolution in imaging magma reservoirs.
Our understanding of magma differentiation and storage in the crust at arc and hotspot volcanoes is undergoing a paradigm shift. An increasing number of observations are at odds with the classical magma chamber model, which postulates that evolved magma is brewed by fractional crystallization, assimilation of host-rock, and crystal settling in high melt fraction magma chambers that are several cubic kilometres in volume (DePaolo, 1981; Marsh, 1989). According to the magma chamber model, the separation of crystals from melt in the magma chamber is responsible for the diversity in the chemical characteristics of magmas and, to some extent, the eruptive behaviour of volcanoes. However, seismic tomography studies have been unable to detect large melt-rich magma chambers and are generally consistent with low melt fractions (Lees, 2007). In addition, petrological studies of crystal residence and differentiation timescales indicate that most magma chambers are ephemeral (Turner and Costa, 2007; Druitt et al., 2012) and that magmatic systems are likely dominated volumetrically and over long time periods by crystal mush reservoirs (Cooper, 2017).
A similar paradigm shift occurred a few decades ago for mid-ocean ridge magmatic systems following the collection of detailed geophysical images (Solomon and Toomey, 1992; Sinton and Detrick, 1992; Hussenoeder et al., 1996). Those results showed that crustal magma systems at fast-spreading ridges are complex, layered, and mush-dominated and do not conform to predictions for a steady-state, melt-filled axial magma chamber that, based on thermal models and geological observations in ophiolites, were envisioned to be up to 10-km-wide and several kilometres in thickness (Cann 1974; Casey & Karson 1981; Pallister & Hopson 1981).
At fast-spreading ridges such as the East Pacific Rise, the melt distribution inferred from seismic refraction and reflection studies indicates that magma accumulates at several levels in the magmatic system. A narrow (∼1-km-wide), axially discontinuous, and by inference ephemeral, magma lens is imaged 1–2 km beneath the seabed (Detrick et al., 1987; Kent et al., 1990; Hooft et al., 1997; Singh et al., 1998). The axial magma lens overlies a wider (5–7 km) lower-crustal region dominated by crystal mush (e.g., Vera et al., 1990) that broadens significantly at, and below, the Moho transition zone (Dunn et al., 2000; Toomey et al., 2007). The melt fraction in the axial magma lens is spatially variable (Xu et al., 2014) and partially molten magma lenses may also exist off-axis (e.g., Aghaei et al., 2017). Seismic studies at the intermediate-spreading Juan de Fuca Ridge reinforced the general view of a composite magmatic system comprised of small, ephemeral melt lenses (Canales et al., 2006; Carbotte et al., 2006; Van Ark et al., 2007) that caps a crystal mush zone of variable width in the crust (Arnoux et al., 2019) and topmost mantle (VanderBeek et al., 2016). Similar studies at slow-spreading ridges revealed a much more discontinuous and temporally variable magmatic system. Seismic results from the Mid-Atlantic ridge indicate that mantle melts are focused toward the segment centre at mantle depths (e.g., Tolstoy et al., 1993; Hooft et al., 2000; Dunn et al., 2005). In the shallow crust, low-velocity bodies of limited size indicate that melt is redistributed along axis by lateral dike injection to supply seafloor eruptions along the ridge segment (Magde et al., 2000).
The new paradigm at arc and hotspot volcanoes, the mush model, is described in detail by Sparks et al. (2019). The mush model does not entirely eliminate the concept of a magma chamber. Instead, it recognizes that magmatic systems may include both regions of mush and magma chambers. Magma differentiation is inferred to happen primarily in the middle and lower crust in mushy hot zones, but magma chambers can play an essential role in the final stages of magma evolution leading up to eruption (Sparks and Cashman, 2017). A range of melt storage conditions, spanning from low to high melt fraction storage, may exist in systems of different geological settings and maturity (Giordano and Caricchi, 2022).
Seismological and other geophysical studies are often cited as a key body of evidence supporting mush-dominated systems (Cashman and Giordano, 2014; Sparks et al., 2019) but the primary literature is usually cautious about interpreting the presence of melt. A seminal review paper on seismic imaging of magmatic systems by Lees (2007) concluded that “there is a general lack of strong evidence for large regions of pure melt below most volcanoes”. We show that Lees’ conclusions are still valid in 2022 and that more recent experiments have not uncovered large magma chambers.
A lack of evidence does not necessarily equate to evidence of absence. Lees’ review emphasized that the evidence for low melt fraction storage is not entirely conclusive because seismic tomography models come with two important caveats. Firstly, the limited sensitivity of seismic waves and inherent smoothing of many inversion methods lead almost universally to the under-recovery of low-velocity anomalies (Malcolm and Trampert, 2011). Secondly, rock physics models that tie seismic velocities to melt fraction are fraught with large uncertainties and trade-offs, meaning that most published melt fraction estimates are “speculative at best” (quote from Lees, 2007). These caveats are still extremely relevant today, despite the increased quality of seismological data and the development of advanced inversion methods. We discuss seismic resolution and sensitivity in Seismic tomography in brief and the rock physics of partially molten rocks in Seismological signatures of mush and magma.
The mush model has far-reaching implications for a wide range of phenomena, from volcanic hazards to the search for precious metal deposits. Therefore, it is of great importance that its foundations are established on robust observations. To fully understand the seismological evidence for the mush model one needs an understanding of the different imaging methods, seismic wavefield sensitivity and model resolution, and the rock physics of partially molten rocks. There are several barriers which make the results of seismic tomography inaccessible to the non-specialist: i) The last couple of decades have seen a proliferation of tomography methods using different parts of the seismic signal, each with its advantages and drawbacks; ii) There is no standard tomography code used by the community; codes are often proprietary and are used by small groups; iii) The results of seismic tomography depend as much on the inversion method as on accurately predicting the travel-times and on the quality of the input data; iv) Resolution limits and uncertainties are often unclear and/or not properly accounted for.
The principal aims of this review are twofold: 1) to help the volcanologist and other non-specialists to understand and interpret seismic tomography, including its limitations; and 2) to provide a template for the volcano tomography practitioner to enhance the impact of their results and make them more accessible to the volcanological community. We focus in particular on one aspect that previous reviews have only touched on briefly: how can melt fractions be estimated from seismological data? We also summarize the main contributions of volcano tomography to the field of volcanology, particularly concerning the development of the mush model. Our review is accompanied by an online database of volcano tomography studies, which will be updated with new results as they become available.
2 Seismic tomography in brief
Seismic tomography is a family of subsurface imaging techniques that have the objective to determine the seismic properties and structure of a geological target. The key feature of seismic tomography is the reconstruction of the seismic properties of the subsurface (such as the P-wave velocity,
The most common type of seismic tomography is arrival-time tomography, and it can be used to illustrate the general underlying principles. The starting point is to identify known seismic phases and measure their time of arrival at the recording station. If one knows the origin time, the average speed of the seismic waves along the assumed propagation path can be determined by dividing the distance by the travel-time. If enough measurements are available to cover the target from multiple angles and azimuths, one can in principle reconstruct the three-dimensional distribution of seismic velocity using one of the well-established algorithms such as the regularized least-squares method or back projection (e.g., Toomey et al., 1994; Zelt and Barton, 1998). In practice, seismic tomography must be carried out iteratively due to the linear approximation that has to be applied to tackle real-world seismic problems that are almost always non-linear, particularly for volcanoes. The main source of non-linearity is the fact that the ray trajectories are dependent on the unknown velocity distribution. If earthquakes are being used, the source parameters are also unknown and must be determined; this results in additional tradeoffs between hypocentral parameters and velocity structure.
The well-established workflow of geophysical inversion starts with an initial guess of the distribution of seismic velocities (the starting model). Synthetic travel-times are calculated for the starting model. This is called solving the forward problem. For volcanoes, where multipathing, diffraction, and attenuation are common, a critical step in forward modelling is correctly identifying if a phase is primary or secondary and then solving the forward problem for the correct ray geometry. The synthetic travel-times are then compared to the observed travel-times (the data). The differences between the data and synthetics are called residuals. A correction to the starting model (called the model update) that minimizes, or at least reduces, the residuals are then calculated. This is the inversion or optimization step. New synthetic data are then computed and compared again to the observed data to find a new model update. The inversion continues until the residuals are reduced to a level comparable to the measurement uncertainty or until the model update becomes too small to be meaningful, i.e., convergence is achieved. The residuals have two sources of error, one is observational (uncertainty in the time assigned to an arrival) and one is numerical (error in time predicted by forward modelling) and it is not uncommon that the latter rivals the former. In practice, many tomographic methods are incapable of fitting data to the measurement uncertainty due to these inaccuracies and to trade-offs between model parameters (e.g., seismic velocity and hypocentral parameters).
Although the basic principle of seismic tomography is relatively straightforward for simple velocity structures, its application to field data from volcanoes is nothing but. For volcanoes, the primary challenge is the physics of wave propagation around and through pronounced low-velocity anomalies, a topic discussed below. Additionally, one common misconception is that one can find the “true” solution to the seismic inverse problem. Because the solution to the seismic inversion of real data is always non-unique, there are an infinite number of models that can fit the data equally well. The goal is to constrain a family or families of models that reduce the residuals. The tomographic method can also be used in a hypothesis-testing mode (e.g., Hammond and Toomey, 2003) to both find classes of models that are consistent with the data and to eliminate classes of models that are inconsistent with observations. Ideally, we should assign an uncertainty to each of these models; however, the inherent non-linearity makes this a challenging task. Because of nonuniqueness, additional constraints are introduced to stabilize the inversion, reduce the range of possible models, or test null hypotheses (e.g., Hammond and Toomey, 2003). These constraints fall under two main categories: i) a priori constraints i.e. any independent information about the velocity distribution, e.g., information about the possible range of seismic velocities expected; ii) regularization constraints that prevent the output model from being too complex. Regularization effectively limits the impact of ill-constrained perturbations that are both large and uncertain on the output model.
The main challenges in seismic tomography are:
1. Designing field strategies that best illuminate the chosen target
2. Finding accurate and computationally efficient ways to solve the forward problem and calculate synthetic data
3. Finding efficient strategies to solve the inverse problem and minimize the residuals
4. Limiting the effect of noisy data on the results
5. Estimating the uncertainties and resolution
The development of different solutions to these challenges has led to an array of different seismic tomography algorithms. Many seismic imaging methods have been applied to volcanic targets. The main differentiating features are the type of seismic source and the part of the seismic wavefield used. Above we have described travel-time tomography, which uses the travel-times of body waves. Other popular methods use surface waves, and more advanced methods use the full wavefield. Below we briefly review some of the most popular seismic tomography approaches. This list will be used to categorize the recent literature presented in Overview of recent volcano tomography results.
2.1 Active source travel-time tomography
One of the most reliable methods for seismic imaging relies on generating artificial seismic signals at known times and locations and recording their arrival times at an array of seismic receivers. Offshore experiments typically use air guns (Figure 1A) while experiments on land use controlled explosions (Figure 1B). The recorded time series are analyzed to identify P-wave and sometimes the S-wave arrivals and measure their travel-times. The advantage of active source tomography is that the origin times and locations of the sources are known. The geometry of the experiment can be tightly controlled to provide the best possible target illumination and achieve the desired resolution. Field campaigns can be relatively short and typically last a few weeks or months. The main disadvantages are high costs and the limitation of having all sources and receivers at or close to the surface. Active source tomography can provide high resolution isotropic and anisotropic
FIGURE 1. Schematic representation of the most common seismic imaging methods. (A) Marine active source travel-time tomography (marine ASTT). (B) Land active source travel-time tomography (land ASTT). (C) Local earthquake tomography (LET). (D and E) Teleseismic tomography (TST). (D) illustrates the far field propagation from the source to the boundary of the local domain. (E) Illustrates the ray-paths within the local domain. (F) Teleseismic receiver functions imaging (TRFI).
2.2 Local earthquake tomography
Local earthquake tomography uses a similar principle to active source tomography but employs earthquakes as sources (Figure 1C). The earthquakes must be located within the target region, meaning that this technique is only suitable for tectonically active regions. Active volcanoes are very often seismically active, and many are monitored with continuously-recording seismometer arrays. Therefore, volcanoes are prime targets for the application of LET and indeed there are many published LET studies in the last 2 decades (as evidenced in our database). Since the origin time and location of the sources are unknown, LET suffers from greater non-linearity than ASTT. This limitation usually result in lower resolution and a greater incidence of artifacts in the output models compared to ASTT.
In addition to the issue of non-linearity, the main disadvantage of LET compared to ASTT is that the location and timing of the sources cannot be controlled. Illumination of the target may not be optimal. Moreover, if the seismicity rate is low, instruments may need to be deployed in the field for long periods of time to achieve sufficient coverage. On the other hand, the main advantages of LET are that it is cheaper to carry out, it can more easily provide joint constraints on
2.3 Teleseismic tomography
Teleseismic (literally “far earthquake”) body wave tomography relies on recording the body waves from distant earthquakes on a local or regional network of seismic stations to image the structure beneath the network (Figures 1D,E). This technique can be applied almost anywhere in the world with the deployment of a dense network of seismometers. TST usually requires broad-band instruments (that are sensitive to low-frequency signals) to achieve a good signal-to-noise ratio because high-frequency seismic waves are more rapidly attenuated when propagating over long distances and are more easily swamped by noise in teleseismic data.
TST can image deeper than most other methods but provides lower resolution because of the low-frequency content and long propagation paths of teleseismic signals (which result in broad Fresnel zones, see Wavefront healing). For this reason, it is most successful in imaging large magmatic systems associated with large calderas (e.g., Steck et al., 1998; Bai et al., 2020) or continental rifts (e.g. Kounoudis et al., 2021). For a more complete comparison of LET and TST see Thurber (2003). A detailed discussion of TST applications to imaging volcanoes can be found in Iyer and Hirahara (1993). Popular algorithms for TST include the algorithm of VanDecar et al. (1995), the FMTT code by Rawlinson et al. (2006), and TomoLab (Hammond and Toomey, 2003).
2.4 Surface wave tomography and ambient noise surface wave tomography
Surface wave tomography uses a different part of the wavefield to image the subsurface. Surface waves are seismic waves that propagate close to the Earth’s surface. Most often SWT relies on Rayleigh waves, a type of surface wave that causes ground oscillations along a vertical plane parallel to the direction of propagation. They arrive later than P and S waves in seismic recordings, but their amplitude is typically stronger. The velocity of surface waves is related to
SWT can be carried out using local earthquakes (earthquakes with hypocenter inside the area of the recording network) or distant earthquakes. Surface wave dispersion curves can also be generated by cross-correlation of ambient noise signals (Shapiro et al., 2005). Ambient noise surface wave tomography (ANSWT) has the advantage of not requiring any sources at all and it can therefore be applied at locations where there are few or no earthquake sources.
The resolution of SWT depends on the frequency range considered. Higher frequencies lead to better resolution but cannot sample deep targets. The characteristic depth of penetration is proportional to the wavelength of the signal. It follows that SWT resolution decays with increasing depth and the resolution length is about half of the depth.
SWT and ANSWT have been applied to several magmatic targets including Aso (Huang et al., 2018), Yellowstone (Jiang et al., 2018) and Campi Flegrei (Guidarelli et al., 2002). Popular SWT algorithms include the algorithm of Barmin et al. (2001), the FMST code by Rawlinson (Bodin et al., 2012a), the ASWMS code by Jin and Gaherty (2015), and the two-plane wave approximation technique (Forsyth and Li, 2005; Yang and Forsyth, 2006).
2.5 Full waveform inversion
Full waveform inversion is based on modelling and matching the seismic recordings wiggle for wiggle instead of just measuring and fitting travel-times (e.g., Pratt, 1999). FWI relies on the ability to accurately model the wavefield to generate synthetic data. As in travel-time tomography, data and synthetics are compared and the residuals are minimized. In this case, however, the residuals are the difference between the recorded and modelled waveforms within a pre-determined time window. The main advantages of FWI are improved resolution and in some cases the ability to constrain multiple physical parameters including seismic attenuation and anisotropy. The main disadvantages are high computational costs and tight requirements on the survey parameters. The recording network or the source distribution needs to be closely spaced (of the order of the dominant wavelength), therefore FWI has seen the widest adoption in active source studies (ASFWI). Typically, FWI starts by inverting low-frequency data to determine the long-wavelength velocity distribution and then proceeds by gradually introducing higher frequencies. The resolution of FWI depends on the frequency range considered and the resolution length can be as low as half the wavelength (Virieux and Operto, 2009). However, inverting higher and higher frequencies requires finer modelling grids and denser survey geometry leading to rapidly increasing compute times.
Monetary and computational costs have in the past limited the wide adoption of FWI, but several recent successful applications have shown that FWI is feasible today for a wide range of problems (Morgan et al., 2013). In magma imaging, FWI has been applied to active source data (Arnulf et al., 2014; Morgan et al., 2016) and to surface wave Green’s functions derived from ambient noise cross-correlation (Flinders and Shen, 2017). Finite-frequency tomography which may be thought of as an intermediate step between travel-time and full-waveform inversion (Tape et al., 2007), has not yet been widely applied to volcanic targets (e. g., Koulakov & Shapiro, 2015) and will not be discussed here.
2.6 Other imaging methods
Several other seismic imaging methods that do not technically fall within the category of seismic tomography have been applied to magmatic systems. All these methods seek to identify sharp boundaries at magma/rock interfaces and are complementary to the information obtained from tomography.
Multi-channel seismic (MCS) reflection imaging is an active source seismic imaging technique that uses reflected seismic waves to build an image of the subsurface. It is best known for its widespread application in hydrocarbon exploration. MCS can provide the highest resolution constrain of any of the listed techniques but struggles to image deep targets and does not on its own provide reliable constraints on the seismic velocities. It is best suited for imaging shallow sub-horizontal targets. It has been successfully applied to image axial magma lenses and other magmatic sills beneath mid-ocean ridges (e.g., Hooft et al., 2000; Carbotte et al., 2020) and to image a melt layer at Campi Flegrei caldera (Zollo et al., 2008). Reflection imaging is best carried out with man-made seismic sources but can also be attempted using earthquake sources and a dense recording array (Byerly et al., 2010).
MCS data offer one of the most robust ways of determining the seismic properties of a low-velocity layer. The analysis of amplitudes vs. offset (AVO analysis) or amplitude vs. angle of incidence (AVA analysis) of the seismic waves reflected at the interface at the top of the layer can be used to accurately determine the seismic velocities of the layer (e.g., Castagna and Backus, 1993). This method has been used to distinguish between fully molten and mushy layers at mid-ocean ridges (Singh et al., 1998; Canales et al., 2006; Xu et al., 2014).
Teleseismic receiver function imaging (TRFI) uses signals from distant earthquakes and their reflections and conversions created by interfaces beneath the recording array (Figures 1D,F) (e.g., Burdick and Langston, 1977; Zhu and Kanamori, 2000). It is conceptually similar to MCS (Ryberg and Weber, 2000) but is better suited to imaging larger scale and deeper targets. It is most sensitive to horizontal layering and is widely used to determine Moho depth variations and the shape of subducting slabs (e.g., MacKenzie et al., 2010). Receiver functions are most often used to generate a direct image of subsurface discontinuities (e.g., Cornwell et al., 2010), but can also be inverted to obtain a velocity model (Lodge and Helffrich, 2009). When applied to volcano imaging TRFI can provide robust constraints on the deeper parts of magmatic systems (e.g., Hammond et al., 2020). In addition, reflection Green’s functions can be obtained from autocorrelations of seismic signals and/or of ambient noise to identify sharp boundaries associated with the top of magma systems (Heath et al., 2018).
An alternative approach for imaging volcanoes is attenuation tomography (AT). Attenuation tomography uses information contained in the Fourier amplitude spectrum of seismic phases to reconstruct a model of the spatial distribution of the attenuation parameter Q (e.g., Lees and Lindley, 1994). AT can use body waves from active source or earthquake sources, earthquake coda waves or ambient seismic noise. Magmatically active areas are typically strongly attenuative because high temperatures and melting can greatly increase viscous dissipation (e.g., Mavko, 1980; Hammond and Humphreys, 2000). Different AT approaches have been developed and have been applied to active volcanoes, including at Newberry caldera (Zucca and Evans, 1992), Mount St. Helens (De Siena et al., 2014) and Tenerife (Prudencio et al., 2015). We do not discuss attenuation tomography in detail in this paper as we think that it would be best treated in a dedicated stand-alone manuscript.
Some of the methods discussed above can be used in tandem to achieve better imaging than would be possible with one single method. For example, surface waves are often used in conjunction with receiver functions, taking advantage of their complementary sensitivity which leads to more robust images (Bodin et al., 2012b). When performed separately, SWT and TRFI suffer from low resolution and a trade-off between the depth of the discontinuity and the velocity of the overburden (Ammon et al., 1990), respectively. Joint inversion, unlike TRFI, constrains absolute shear velocities, and provides better vertical resolution than SWT (e.g., Chrapkiewicz et al., 2020). LET, ASTT and TST are also often used together to achieve better target illumination and more robust convergence (e.g., Battaglia et al., 2008; Heath et al., 2015; Díaz-Moreno et al., 2018). Seismic methods can also be coupled with other geophysical data, e. g gravity (Paulatto et al., 2019). Joint inversion can help reduce the null space of the solution (narrow the range of acceptable models).
2.7 Understanding seismic resolution
Seismic tomography studies should be accompanied by comprehensive resolution and uncertainty analysis so that the output model can be interpreted with confidence. In addition, since volcanoes are strongly heterogenous, extensive forward-inverse modeling should be conducted to understand the sensitivity of the results to starting models and to pronounced low-velocity anomalies. Estimating seismic resolution and model sensitivity is not straightforward and can be even more time consuming than carrying out the inversion in the first place.
Seismic resolution is a matrix that describes the linear averaging of the model parameters that results from the distribution of data coverage (e.g., Tarantola, 2005; Rawlinson and Spakman, 2016). For a model defined by N parameters, the resolution matrix is an NxN square, symmetric matrix. Each row of the resolution matrix corresponds to the linear averaging kernel for one model parameter. For example, if we estimate
For small inverse problems, e.g., 2D travel-time inversions, the resolution matrix can be directly calculated. However, for most modern large and/or 3D problems this is impractical. This is the case particularly for FWI, which has a much higher computational cost than travel-time tomography (Fichtner and Trampert, 2011). Therefore, in most cases, the calculation of the full resolution matrix is replaced with a sensitivity analysis.
Checkerboard recovery tests are the most common type of sensitivity analysis. These consist of carrying out a simulated inversion experiment using synthetic data to test the recovery of an alternating pattern of positive and negative anomalies given the geometry of sources and receivers of the experiment (Zelt, 1998; Rawlinson et al., 2014). By performing multiple tests with anomalies of different wavelengths and observing where the anomalies can be recovered and how they are smoothed and smeared it is possible to estimate how the resolution varies throughout the model. Checkerboard tests for volcanoes, however, can be misleading since the amplitude of the anomalies tested (nominally ± a few percent) is much smaller than the heterogeneity typical of volcanoes. As a result, the non-linearity for the checkerboard test is much less than for a geologically plausible model.
Less common, but often more informative are spike tests, which use isolated anomalies rather than a pattern of alternating anomalies and test the recovery of a single model parameter. In practice, they are a way of estimating one row of the resolution matrix in isolation. Rawlinson and Spakman (2016) argue that spike tests should be preferred to checkerboard tests, because the ray coverage of the checkerboard test may be significantly different from the ray coverage of the experiment. Checkerboard tests lead to the overlap of recovered patterns, and they can give a false impression of the resolving power (Rawlinson and Spakman, 2016).
Another important concept, closely related to the resolution matrix, is the resolution length, i.e., the minimum separation between two real features that can be distinguished as separate in the output model. From a checkerboard test, one can estimate the resolution length as the smallest anomaly that can be correctly distinguished and recovered, under the assumption that the magnitude of anomalies in the checkerboard and the actual medium are similar.
Particularly useful for volcanoes are forward-inverse tests to study the recovery of low-velocity anomalies. Typically these are geologically or thermodynamically informed models that test the melt fraction and spatial geometry of the magmatic system (e.g., Beachly et al., 2012; Paulatto et al., 2012; McVey et al., 2020). Such tests probe the family of models that fit a given data set and are crucial for expanding our understanding beyond a discussion of the minimal tomographically recovered structure. These studies also allow authors to estimate the amount of under-recovery in their studies, a typical problem with wavefront healing around low-velocity bodies (discussed below). A sensitivity analysis that tests a suite of hypotheses can be more informative than either checkerboards or spike tests.
For the purpose of discussing the reliability of tomographic images, we define several terms (Figure 2).
•Image fidelity is defined by the size of the image pixels (voxels in the case of a 3D model) or the spacing of grid nodes (Figure 2B). The image fidelity or grid resolution is specified by the user and it affects both the forward and inverse problems. Image fidelity is sometimes determined by stability and accuracy requirements of the forward problem or considerations of computational cost. For example, stable wave propagation modelling requires a grid spacing that is fine enough so that the wavefront does not travel more than about half of the length of a model cell over each timestep (Press et al., 1990). Similarly, finite difference Eikonal solvers or shortest path methods require a relatively fine grid to produce accurate travel-times. Image fidelity also trades off with model resolution such that improvements in image fidelity decrease the resolution of individual model parameters as measured by the resolution matrix. Jackson (1979) argues that low image fidelity can give rise to non-random errors, or bias in the output model.
•Physical sensitivity describes the sensitivity of the seismic waves to the velocity structure that they travel through or in other words the extent to which the waves “feel” any velocity anomalies (Figure 2C). The sensitivity of seismic waves can be described by a sensitivity kernel and is sometimes included in the inversion regularization. Different parts of the wavefield have different sensitivity to certain structures. For example, transmitted waves are sensitive to the seismic velocity along the travel path while reflected waves are sensitive to discontinuities. Higher-frequency waves have better sensitivity to small scale structures. In principle, one could increase the frequency of investigation to improve the physical resolution, but in practice, this can only go so far. The main reason is that higher frequency signals are attenuated more strongly than low-frequency signals and are therefore swamped by noise at longer offsets.
•Data coverage describes the ability of the data to sample the target (Figure 2D). This includes the effect of the acquisition geometry, source distribution, density of sampling and azimuthal coverage of the target. Even within regions of good data coverage, the model resolution can vary spatially.
•Inversion regularization describes the model damping, smoothing, or other pre-defined relations that are introduced during inversion. Most seismic tomography studies use some form of smoothing or damping of the solution and some incorporate physical sensitivity via kernels. Spatial smoothing constraints or sensitivity kernels give rise to similar effects by imposing relations on volumetric averaging of the velocity model (Tape et al., 2007). Some form of regularization is always necessary to stabilize the inversion and rule out models that have large and poorly constrained perturbations.
•Model resolution is the combined effects of the image fidelity, data coverage, and inversion regularization (Figure 2E). Since many tomographic algorithms do not include sensitivity kernels in their inversions, the effects of physical sensitivity may not be explicitly part of model resolution. In these cases, one must use prior understanding of the physical sensitivity to interpret the results. In parts of the model with excellent data coverage, the model resolution can theoretically approach the physical sensitivity. However, as discussed below, the physics of seismic wave propagation around and through pronounced low-velocity anomalies means that excellent data coverage is not often attainable in a volcanic setting.
FIGURE 2. The effect of different quantities on model recovery. (A) Idealized true velocity model of an active volcano. (B) Representations of the true model with different image fidelities. The numbers on the top left corners indicate number of pixels on each side of the image. (C) The true model as “felt” by the seismic waves due to the limits of the physical sensitivity. (D) The data coverage. (E) The model resolution is affected by image fidelity, the data coverage, and the inversion regularization.
2.8 Resolution of travel-time tomography and surface wave tomography
To better understand seismic resolution, we use travel-time tomography as an example once again. Most travel-time tomography studies rely on the ray approximation i.e., on the assumption that the travel-time depends only on the seismic velocities along a ray. This is equivalent to an infinite frequency approximation. In practice, field data contain a limited range of frequencies and travel-times are picked on band-pass filtered data. In this case, scattering and diffraction effects become important, leading to inaccuracies in the reconstructed model (Williamson, 1991; Williamson and Worthington, 1993). The ray approximation leads to incorrect travel-times and wrong model updates when the size of the velocity anomaly is comparable to or smaller than the seismic wavelength. Scattering can be particularly pronounced in volcanic settings, where the Earth’s structure is highly heterogeneous (e.g., De Siena et al., 2014).
The physical sensitivity of travel-time tomography is approximated by the Fresnel radius or the radius of the first Fresnel zone. The first Fresnel zone is the area from which energy arriving at the receiver has a phase differing no more than half a period, and so interferes constructively with the direct arrival. For refracted ray paths the Fresnel radius can be approximated as
The resolution of traditional SWT methods is also limited by the Fresnel radius since ray theory is used to invert the phase velocity maps for the 3D
The overall effect of the physical sensitivity and model resolution is that the recovered image is blurred, and anomalies are diminished in amplitude compared to the true model (Dahlen, 2004). Anomalies with a size comparable to the physical sensitivity cannot be fully recovered. The sensitivity and resolution typically worsen with depth since longer offsets or lower frequencies are needed to sample deeper. The resolution length is not a clear-cut limit and smaller features can still be detectable, although the anomaly amplitude will be adversely affected (Sheriff and Geldart, 1995; Spetzler and Snieder, 2004).
2.9 Wavefront healing
A more intuitive way to understand why anomalies comparable to or smaller than the seismic wavelength cannot be recovered is to look at wavefront healing, which is the result of scattering (Malcolm and Trampert, 2011). Seismic waves crossing a low-velocity anomaly are slowed down creating a delay in the wavefront behind the anomaly (Figure 3). As the wavefront propagates further, scattering of the waves around the anomaly causes the shadow of the anomaly to “heal”. If the anomaly is similar in size to the wavelength of the signal and the propagation distance is significantly larger than the wavelength, wavefront healing can effectively hide the anomaly, i.e., the anomaly results in a negligible observed travel-time delay (Figure 3H). However, even when an anomaly is fully hidden from travel-times by wavefront healing, it affects the amplitude of arrivals and generates secondary waves (Figure 3) which can be exploited by forward waveform modelling (e.g., Beachly et al., 2012), explicit ray tracing of secondary arrivals, or more advanced methods like full-waveform inversion. There is an asymmetry in the evolution of the wavefield around positive vs. negative travel-time anomalies. Negative delays will initially grow, or ‘unheal’, but then heal at a rate faster than positive delays (Nolet and Dahlen, 2000). This means that slow anomalies are less resolvable than fast anomalies of the same size using travel-times.
FIGURE 3. Illustration of the wavefront healing effect. (A–C) Wavefield snapshots at 1.0 s, 1.5 s, and 2.0 s after the source origin time. A square low-velocity anomaly with side equal to 0.5 times the wavelength is located between the source (star) and the receivers (P1, and P2). The strength of the anomaly is −50% of the background velocity. (D–F) Wavefield snapshots at 1.2 s, 1.6 s, and 2.0 s after the source origin time with a square low-velocity anomaly with side equal to 1.5 times the wavelength. The first arrival wavefront is highlighted with a dashed line. (G) Time-series recorded at points P1 with the ×0.5 wavelength anomaly (green) and without the anomaly (black). (H) Same as (G) at point P2. (I) Same as (G) but for the ×1.5 wavelength anomaly. (J) Same as (I) at point P2. Wavefields were calculated using an acoustic finite difference method (Warner et al., 2013).
FWI resolution is not limited by the Fresnel radius and can recover anomalies as small as half the wavelength (Virieux and Operto, 2009). This is because FWI uses not only the travel-time but also the amplitude of the waves and can use secondary phases like the diffractions and converted phases generated by velocity perturbations. The resolution length of FWI is controlled by the seismic frequency and by the geometry of the data acquisition (via the scattering angle). It increases with the frequency of the wavefield and decreases with the scattering angle (Eq. 30 of Virieux and Operto, 2009). For the scattering angle approaching zero (normal-incidence reflection), the resolution length approaches half of the seismic wavelength.
In the FWI study of Chrapkiewicz et al. (2022) the frequencies used are up to 6 Hz leading to a physical resolution limit of ∼0.4 km. In practice, the model regularization introduced to stabilize the inversion due to data sparsity means that the model resolution is about 1 km. The resolution length for seismic imaging of magma reservoirs in the upper crust is typically 2 km–10 km for travel-time tomography (ASTT and LET), 4 km–20 km for SWT, and 0.4 km–2 km for ASFWI. We will discuss the implication of resolution limits on melt estimates in Melt fractions in the crust.
3 Seismological signatures of mush and magma
Using an analogy from hydrocarbon exploration, we can think of composition, temperature, and melt fraction within a magma reservoir as reservoir parameters. These are the parameters of interest to the volcanologist to understand the dynamics of the system and interpret unrest. However, the primary outputs of seismic tomography are constraints on the elastic properties, like seismic velocities and anisotropy. Determining the reservoir parameters requires an understanding of the physical relationships that tie them to the measured elastic properties. These relationships are sometimes called constitutive relationships and can be derived either theoretically or experimentally. In this section, we focus on the constitutive relationships that tie melt fraction to the elastic moduli and seismic velocities. Before we tackle this topic, we need to introduce a few rock physics concepts.
3.1 Elastic moduli and seismic velocities
The elastic moduli are the parameters that describe how a material reacts to stress under the assumption of linear elasticity. The most intuitive to understand is the bulk modulus (
Stress is defined as force per unit area, while strain describes the deformation.
For a general anisotropic material, 81 elastic moduli can be defined. The elasticity tensor or stiffness tensor C is a matrix containing the 81 elastic moduli. Of the 81 components only 21 are independent because of the symmetry properties of the stress tensor, therefore 21 constants are needed to fully describe the elastic response of a general anisotropic material. For isotropic materials, the number of independent constants is reduced to two. Popular choices are the Lamé parameters
The seismic velocities are related to the components of the stiffness tensor and to density. In a solid isotropic material two types of seismic body waves can propagate, P-waves and S-waves, and these are characterized by the P-wave velocity (
The seismic waves used in crustal-scale seismic tomography have wavelengths ranging from a few hundred meters to several kilometres. As seen in Seismic tomography in brief, the sensitivity of the seismic waves in the best-case scenario is about half of the signal wavelength. The Earth contains heterogeneity at all length scales, and magma reservoirs are no exception. Seismic waves passing through a partially molten region are insensitive to small-scale heterogeneities and will “feel” or react to the average or effective seismic velocities and the related effective elastic moduli.
3.2 Effective elastic properties of partially molten rocks
The effective elastic moduli of partially molten rocks depend on the elastic moduli of the components (crystals/melt/gases), their respective abundances, and on their geometric arrangement and connectivity (Figure 4A). In general, the effective elastic moduli will be frequency-dependent as the flow of fluid phases causes viscous energy dissipation and frequency dispersion. The response may also be anisotropic if any of the constituents are anisotropic (as are most minerals) and/or if the microstructure has a preferential alignment. In this paper, we focus on isotropic elastic properties. For a discussion of the effect of magma on seismic anisotropy see e.g. Hammond et al. (2014).
FIGURE 4. (A) Classification of partially molten rocks on a scale from solid rock to pure melt adapted from Nur et al. (1998). The critical porosity is the melt fraction where the mechanical behavior transitions from matrix supported mush to fluid supported suspension. (B) Spheroidal shapes used to model the effect of melting on elastic properties. (C) Crystal mush can be represented as a distribution of spheroidal melt pockets embedded in a solid matrix. Magma can be represented as a suspension of solid particles in melt.
3.2.1 Bounds and averages
The simplest approach to estimating the physical properties of composite materials is to take the volumetric average of the properties of the components and ignore any geometric effects. The arithmetic mean of the elastic moduli is also called the Voigt bound (Voigt, 1928), while the geometric mean (mean of the reciprocals) is called the Reuss bound (Reuss, 1929). It can be proved that the Voigt and Reuss bounds are maximum and minimum bounds on the possible elastic properties of an isotropic mixture (Hill, 1952). The Reuss bound is particularly useful as it is exactly valid for an idealized dilute suspension of solid particles in a fluid. The Reuss bound is therefore a good approximation for calculating the elastic moduli of magma with suspended crystals. For solid rocks and solid/fluid composites, the mean of the Voigt and Reuss bounds, known as the Voigt-Reuss-Hill (VRH) average (Chung et al., 1963) is sometimes used as a first approximation in the absence of any information on the microstructure. The Hashin-Shtrikman (HS) bounds [Hashin & Shtrikman, 1963, eqs (4.1–4.4)] are a bit more complex than the Voigt and Reuss bounds and are the tightest bounds that can be derived for the elastic moduli of an isotropic random mixture (Figure 5). The upper HS bound is exactly valid for an idealized distribution of fluid spheres within a solid. The lower HS bound is equivalent to the Reuss bound.
FIGURE 5. Example constitutive relationships for partially molten rocks. (A) Bulk modulus (K) and shear modulus (G) as a function of melt fraction for a two-phase composite with spheroidal melt inclusions (aspect ratio = 0.2) calculate with the SCS (solid line) and the DEM (dashed line). (B) Elastic moduli as a function of melt fraction calculated with the SCS for different values of the inclusion aspect ratio; blue: 0.05, black: 0.2, red: 1.0. (C) Elastic moduli as a function of melt fraction calculated with the critical porosity approach of Nur et al. (1998). (D–F) Same as (A–C) but showing Vp and Vs. instead of the K and G. The grey lines mark the upper and lower Hashin-Shtrikman bounds (HS+ and HS-).
3.2.2 Effective medium theories
Progressing beyond the simple bounds approach requires some information on the geometrical shape and arrangement of the components. A popular approach is to assume that the material is composed of a matrix made of one material with one or more additional phases embedded in the matrix in the form of inclusions with a defined simple geometrical shape. For example, a magma can be represented by a melt “matrix” with suspended crystal grains represented as sub-spherical solid inclusions. Similarly, a crystal mush can be represented as a solid matrix with interconnected melt inclusions (Figure 4C). This approach was pioneered by Eshelby (1957), who calculated the effect on the elastic properties of an ellipsoidal inclusion (Figure 4B). Based on his work, a wide range of effective medium theories have been developed to calculate the properties of materials containing inclusions with different shapes including spheroids (Wu, 1966; Kuster and Toksöz, 1974), films (Walsh and Grosenbaugh, 1979), and cuspate tubes (Mavko, 1980). For more in-depth information on the elastic properties of composite materials see e.g. Watt et al. (1976) and Mavko et al. (2020).
The two most widely used effective medium approaches are the differential effective medium (DEM) approach and the self-consistent scheme (SCS). These are both implicit semi-analytical schemes that require integration or iteration to reach a solution. The DEM approach (Norris, 1985; Avellaneda, 1987) simulates a composite material and calculates the effective elastic moduli by gradually adding inclusions to a matrix. The method is solved by numerical integration and is suitable when the inclusions are sparse enough and do not form an interconnected network (Berge et al., 1993). The SCS approach (Budiansky, 1965; O’Connell and Budiansky, 1974; Berryman, 1980; Mavko, 1980; Schmeling, 1985) considers inclusions embedded in an effective matrix with properties to be determined. The properties of the medium are calculated by iteration, starting with an approximate initial guess (e.g the VRH average). Its advantage is that it implicitly accounts for interaction between inclusions and predicts the existence of a critical porosity (see discussion below).
3.2.3 Inclusion aspect ratios
Figures 5A,D shows the melt-velocity relationships calculated with the DEM and SCS approaches for spheroidal inclusions using typical solid and melt properties for a granite. The velocity reduction for a given melt fraction can be vastly different depending on the chosen inclusion geometry, described by the aspect ratio
3.2.4 Critical porosity
An alternative approach that avoids the problem of choosing an aspect ratio altogether relies on the observation that the mechanical behaviour of partially molten rocks is bimodal. At low melt fractions, below a critical porosity
3.3 Estimating melt fraction from seismic velocities
If the seismic velocities of the magma reservoir can be measured with sufficient accuracy one may attempt to estimate melt fraction. This is not as simple as converting the velocity with a simple scaling constant. Multiple inter-related factors need to be carefully considered. First, one needs to have an approximate idea of the bulk composition of the melt and of the host rock. The composition of the components will dictate their elastic moduli, densities, melting temperatures, and likely micro-geometry or critical porosity.
Estimating these parameters requires some geological knowledge of the volcanic system under study. The seismic properties of the host rock can be estimated by analyzing plutonic xenoliths found in the erupted lavas (Kiddle et al., 2010). Plutonic xenoliths are fragments of crystalline rocks carried to the surface during an eruption. The assumption is that they represent an approximation of the solid component of the magmatic system. Their seismic velocity can be determined by laboratory measurements or as the VRH average of the properties of the constitutive minerals. This type of analysis is only available for a few of the most well-studied volcanoes. More often, the average composition of the host rock is estimated based on the bulk composition of the eruptive products and any available general information about the tectonic setting. This might lead to a broad classification in terms of silica content (e.g., granitic or gabbroic). The seismic velocities of the host rock can then be determined from tabulated laboratory measurements (e.g., Christensen and Stanley, 2003) and then corrected for temperature and pressure. The composition of the melt component is easier to determine. For example, one may assume that the bulk composition of the melt in the magmatic system is similar to the composition of the quenched glass and groundmass in the eruptive products. From the major element composition of the melt, one can then estimate the elastic properties and density using for example the model of Ueki and Iwamori (2016), which was developed by fitting an extensive suite of experimental measurements.
With the composition of the components constrained, one may then make an educated guess of the melt geometry. The different approaches described in Effective elastic properties of partially molten rocks require different parameters and will require different strategies:
1. Critical porosity - Yu and Lee (2016) estimate that the critical porosity at a dacite volcano is ∼30%. This value can be used to tune a critical porosity curve or the SCS approach for volcanoes with similar composition. Van der Molen and Paterson (1979) obtain a value between 30% and 35% for granite. The latter value is used by Chu et al. (2010) to estimate melt fractions beneath Yellowstone caldera.
2. Dihedral angles and spheroidal models - If the mineral composition of the solid component is known at least approximately, dihedral angles can be estimated based for example on the compilation of Holness (2006). One can then use the model of Takei (2002) to relate the dihedral angles to an equivalent aspect ratio (or a range of aspect ratios) and then use a spheroidal model e.g. the SCS to calculate a constitutive relationship.
Having chosen and tuned an appropriate constitutive relationship, it can then be inverted, for example with a line search method. This is equivalent to finding the melt fraction where the constitutive curve matches the observed value of the elastic property considered. The predicted melt fraction should be accompanied by an uncertainty estimate. The largest uncertainty contribution typically arises from uncertainties in the measured seismic properties and in the melt geometry. As mentioned above, the microscopic melt geometry constrained from exhumed samples may not be entirely meaningful if a significant amount of melt is stored in macroscopic melt layers or dikes. More work is needed to determine the best way to deal with multi-scale melt distribution. In the meantime, the safest approach is to consider a range of melt geometries and propagate the resulting uncertainties (e.g., McVey et al., 2020). Other significant uncertainties may arise from the choice of host rock and melt seismic properties and from temperature effects.
4 Overview of recent volcano tomography results
Lees (2007) covered in some detail the most significant volcano tomography results published from 1993 to 2007. Since then, the number of published volcano tomography studies has more than doubled (Figure 6). We have collated over 270 volcano tomography studies published between 1975 and 2022. At least 78 volcanoes have been imaged with seismic methods (Figure 7). We have chosen not to include mid-ocean ridge magmatic systems in our database; for context we briefly summarize the mid-ocean ridge results in Introduction.
FIGURE 6. (A) Volcano imaging studies listed by publication year and grouped based on type of seismic source used. “Joint” means joint inversion of active source and earthquake data. “Other” includes mostly ambient noise studies. (B) Same as (A) but grouped based on imaging method used. ASTT, Active Source Travel-time Tomography; LET, Local Earthquake Tomography; TST, Teleseismic Tomography; TRFI, Teleseismic Receiver Function Inversion; ANSWT, Ambient Noise Surface Wave Tomography.
From each paper, we have extracted information about the recording array, sources, part of the wavefield used, imaging method, quantitative data about any observed low-velocity anomalies, and reported melt estimates. We have strived to be as complete as possible in our literature search, but a few studies may have been missed. The data are collected in an online database (https://github.com/kmch/VolcanoTomo) which is open to new submissions from authors and will be kept up to date with new results. We encourage readers to report any inaccurate or missing information in the database. We particularly welcome input from the authors of the cited studies.
Below we summarize some of the most significant volcano imaging results published before 2022. We focus particularly on those studies that reported clear evidence of melt and provided quantitative and well documented estimates of melt fraction (Table 1).
4.1 Santorini and Kolumbo
Santorini (Hellenic arc, Greece) was first imaged using LET by Dimitriadis et al. (2010). This study found low velocities beneath the northern part of the island extending in the direction of known tectonic lineaments, but the resolution afforded by the sparse dataset was not sufficient to estimate melt fractions. In 2015, a large, dense amphibious active source seismic study (the PROTEUS experiment) recorded ∼14,300 controlled sound sources from a large airgun array on 90 ocean-bottom bottom and 65 land seismometers covering an area of 120 × 60 km2 (Heath et al., 2019). Hooft et al. (2019) used ASTT to generate a model of the shallow subsurface, revealing evidence for a strong low-
Kolumbo is a small submarine volcanic cone, a part of the Santorini volcanic system. Dimitriadis et al. (2010) found evidence for a −3%
FIGURE 8. Examples of seismic tomography results. (A) P-wave velocity model of Kolumbo from active source FWI (adapted from Chrapkiewicz et al., 2022). (B) P-wave velocity model of Montserrat from joint inversion of active source P-wave travel-times and gravity data (Paulatto et al., 2019). (C) S-wave velocity model of Mount St. Helens from joint LET and ASTT (adapted from Ulberg et al., 2020). (D) S-wave velocity model of Laguna del Maule from SWT (adapted from Wespestad et al., 2019).
FWI has yet to be applied to the larger magmatic system directly beneath Santorini (McVey et al., 2020). For this target, the spatial extent, heterogeneity, and attenuation of the magmatic system significantly impact data quality. As such, it remains to be seen whether active-source FWI can image the larger magma bodies expected beneath volcanoes that are associated with significantly reduced signal-to-noise ratio data.
The Soufriere Hills Volcano (SHV) on Montserrat (Lesser Antilles) was surveyed in 2007 and imaged with ASTT. Shalev et al. (2010), published a
4.3 Mount St. Helens
Mount St. Helens (Washington, United States) was first imaged in detail using LET by Lees and Crosson (1989) and Lees (1992) who found a possible shallow magma reservoir, a high-velocity core, and a deeper low-velocity volume at 9 km depth. The shallow reservoir at 2 km–5 km depth is more evident in Lees (2007) who republished the model of Lees (1992) using higher quality colour figures. Waite and Moran (2009) carried out a LET study with a more extensive earthquake database and imaged two low-
Newberry volcano (Oregon, United States) was the subject of one of the earliest 3D volcano tomography studies (Achauer et al., 1988). They used ASTT and imaged a low-
4.5 Campi Flegrei
The Campi Flegrei caldera (Italy) has been imaged with LET (Aster & Meyer, 1988; Vanorio et al., 2005; Calò and Tramelli, 2018), ASTT (Zollo et al., 2003; Chiarabba and Moretti, 2006; Battaglia et al., 2008) and SWT (Guidarelli et al., 2002). At shallow depth, the system is characterized by low-velocity caldera fill surrounded by a high-velocity ring (e.g., Chiarabba and Moretti, 2006). At 4 km depth, there is evidence for a low-
The Toba caldera (Sumatra, Indonesia) was first imaged by Mastruyono et al. (2001) using LET. They imaged a low
4.7 Aso caldera
Aso caldera (Japan) was first imaged by Sudo and Kong (2001). Their study is an exemplary application of LET, containing detailed quality control and a quantitative estimate of melt fraction. They found evidence for a −26%
4.8 Uturuncu/Altiplano Puna
Uturuncu (Bolivia) is the location of some of the most spectacular volcano tomography images worldwide. An early TSRFI study found evidence for a low-
Volcán de Colima (Mexico) and the surrounding Colima Volcanic Complex were imaged using both LET (Ochoa-Chávez et al., 2016; Watkins et al., 2018; Sychev et al., 2019) and SWT (Spica et al., 2014; Escudero and Bandy, 2017; Spica et al., 2017). All these studies found evidence for low seismic velocities in the crust beneath Colima at depth of about 10 km–25 km, although depth ranges from each study differ. Particularly notable are the results of Spica et al. (2017), who imaged a low-velocity volume at 10 km–20 km depth.
4.10 Long Valley
Long Valley caldera (California, United States) has been the subject of many seismic studies already summarized in Lees (2007). Early studies found evidence for large, small-amplitude low-velocity anomalies. Dawson et al. (1990), Steck and Prothero (1994), and Weiland et al. (1995) used TST and found low
Katla (Iceland) was first imaged in 1994 by Gudmundsson et al. (1994). They collected active source data along a 2D profile crossing the caldera and the glacier that covers it. They used iterative forward modelling of P-waves to build a model of
4.12 Mt. Paektu/Changbaishan
Mt. Paektu (North Korea), also known as Changbaishan has been imaged by Kyong-Song et al. (2016) using receiver functions. Their study revealed a negative polarity arrival and high
The first attempts at imaging Yellowstone caldera (Wyoming, United States) were carried out in the 1970s and early 1980s by Iyer and others (e.g. Iyer, 1979; Iyer et al., 1981). By analysing teleseismic travel-times using one of the earliest applications of TST they found evidence for low-velocity bodies within the crust and upper mantle beneath the caldera. A subsequent LET study by Benz and Smith (1984) also found a substantial low-velocity body with
4.14 Canary Islands
Several of the Canary Islands volcanoes (Spain) have been imaged with seismic methods. Their large volumes, long histories and complex geological structure make it hard to distinguish melt bodies. Krastel and Schmincke (2002) imaged Gran Canaria using ASTT but found no evidence for significant melt regions. García-Yeguas et al. (2012) used ASTT to image Tenerife. Their
4.15 Hawaii volcanoes
The active volcanoes on Hawaii’s Big Island (United States) have been the subject of many seismic investigations starting with the TST study of Ellsworth and Koyanagi (1977). Early studies found a high-velocity core beneath Mauna Loa (Okubo et al., 1997) and both low and high velocity regions beneath Kilauea (Thurber, 1984). Haslinger et al. (2001) used LET to image Kilauea and found evidence for low velocities beneath the East Rift Zone, interpreted as a possible magma reservoir with 10% melt fraction. Park et al. (2007, 2009) found evidence for low velocities in the upper mantle beneath Lohi Seamount. Lin et al. (2014) used LET to generate joint
5 Melt fractions in the crust
We found 39 studies that reported quantitative melt fraction estimates associated with 22 different volcanoes (Table 1). The median and mean reported melt fractions are 13% and 18% respectively. The distribution is highly skewed, with the majority of studies reporting melt fractions between 0% and 20% (Figure 9). This result seems to unequivocally support a model where melt storage is dominated by low-melt-fraction processes; however, the story is not so simple.
FIGURE 9. (A) Histogram of reported melt fraction. (B) Reported melt fraction versus depth of low velocity anomalies. (C) Reported melt fraction vs estimated anomaly volume.
Recalling the discussion of seismic resolution in Seismic tomography in brief, we note that at any of the studied volcanoes low-velocity anomalies smaller than the seismic resolution length may have gone undetected. Given that the resolution length is typically 2 km–10 km (Table 1), anomalies as large as several tens of cubic kilometres may have been missed by geophysical imaging even at some of the best-studied volcanoes. In addition, any anomaly comparable to the seismic resolution length would be strongly affected by blurring and amplitude reduction. Many of the reported low-velocity volumes are 4 km–10 km in diameter, a size comparable to the resolution length.
Most studies only take into account some of the sources of error and typically do not account for the under-recovery of low-velocity anomalies. This issue is often glossed over because the under-recovery cannot be easily corrected with a multiplicative constant. Beachly et al. (2012) showed how a −15%
Only five studies reported melt fractions above 35%. These are listed below, together with the estimated melt fractions and imaging methods.
1. Chrapkiewicz et al. (2022), Kolumbo volcano (Hellenic arc), full waveform inversion, 26%–53%
2. Seccia et al. (2011), Long Valley caldera, teleseismic receiver function imaging, 30%–60%
3. Gudmundsson et al. (1994), Katla, active source P-wave forward modelling, 50%
4. Beachly et al. (2012), Newberry volcano, waveform forward modelling, 10%–100%
5. Zollo et al. (2008), Campi Flegrei, amplitude-vs-offset analysis, 65%–90%
It is worth noting that none of these five studies used travel-time tomography and only Chrapkiewicz et al. (2022) used a true inverse method. Beachly et al. (2012) used waveform forward modelling of secondary arrivals, Gudmundsson et al. (1994) used travel-time forward modelling, and the other studies used reflected or converted phases that are particularly sensitive to velocity contrasts across a solid-liquid interface. At all the five above locations, traveltime tomography (either LET or ASTT) was applied by the same or other authors and found small or negligible amounts of melt. Higher melt fractions were detected only when alternative methods, that are more sensitive to low-velocity anomalies, were applied.
There are some unequivocal cases where large anomalies have been observed that are clearly significantly larger than the resolution limit. At Uturuncu, Ward et al. (2014) imaged the 500,000 km3 Altiplano Puna magma body. The seismic velocities within this low-velocity anomaly are consistent with melt fractions below 25%. The diameter of the body is about 200 km, about 4–5 times larger than the horizontal resolution length which is about 30 km–50 km. When this evidence is combined with electrical resistivity data from magneto-telluric inversion (Comeau et al., 2016), it points to the existence of a large and heterogeneous region of crystal mush with melt fraction >15%. Little is known about the internal structure of this crystal mush reservoir. The 15%–25% melt fraction estimate is an average over the scale of the resolution length. Future investigation should focus on determining the multi-scale distribution of melt within mush reservoirs. The distribution of melt fraction is tightly tied to permeability and therefore can have a large effect on the dynamics of melt transport and on the rheology of the mush.
The seismic tomography evidence collected so far does not rule out the presence of small (<10 km3) and medium-sized (<100 km3) magma chambers. These may be present either as isolated bodies or embedded within regions of mush. The only study that reports the presence of a large (>100 km3) reservoir with a melt fraction possibly above the critical porosity is the teleseismic receiver function study of Long Valley Caldera by Seccia et al. (2011). None of the other studies analyzed has uncovered large magma chambers and all are generally in line with the results summarized by Lees (2007). There is an absence of seismological evidence for the existence of large magma chambers, but is this the same as stating that large magma chambers do not exist? Proving the non-existence of something is a notoriously difficult ontological problem since the lack of evidence does not necessarily equate to evidence of absence. However, we can conclude that large magma chambers are at least geologically rare since none have been imaged beyond doubt and we have no reason to believe that the present day is special in any way.
Many authors have argued that large magma chambers must form quickly before an eruption (e.g., Parks et al., 2012; Cashman et al., 2017; Tavazzani et al., 2020). Different mechanisms have been proposed for how low crystallinity magma may be extracted from mush zones. Transfer of heat and volatiles from fresh magma injections may rapidly bring the crystallinity below the lock-up threshold and prime a reservoir for eruption (Weber et al., 2020; Zou and Ma, 2020). Mush remobilization initiated by injection of mafic magma has been suggested as the trigger of eruptions at Volcan de Colima (Hughes et al., 2021) and Lassen Volcanic field (Klemetti and Clynne, 2014). Unrest and eruptions have been linked to the destabilization of the mush at Montserrat (Christopher et al., 2015; Didonna et al., 2022). More work is needed to understand the structure and dynamics of mushy magmatic systems to more confidently interpret geophysical signals that may be associated with the destabilization of the mush.
Outstanding knowledge gaps that are crucial to understanding eruption precursors and magma dynamics and can be addressed by advances in geophysical imaging include. (Figure 10C):
1. The prevalence of volatile-rich and high melt fraction lenses (Rooyakkers et al., 2021)
2. The geometry of eruptive conduits (Cassidy et al., 2018)
3. The melt fraction, internal structure, and eruptive potential of magma reservoirs (Rosi et al., 2022)
4. The melt fraction of deep crustal magma storage regions and the nature of deep melt ascent pathways (Giordano and Caricchi, 2022)
FIGURE 10. (A) Schematic representation of the typical level of detail provided by traditional travel-time tomography. (B) Level of detail that can be achieved by combining the best seismic imaging methods currently available. (C) Modern conceptual model of magmatic systems. Deeper and higher resolution constraints are needed to address outstanding knowledge gaps including: the distribution of magmatic volatiles, the geometry of eruptive conduits, the presence of melt lenses, the internal structure of mush reservoirs, the nature of deep ascent pathways, the properties of deep crustal storage regions.
6 Recommendations for best practice
Our extensive literature review uncovered an impressive body of work on volcano imaging. We identified a few desirable features that in our opinion make a volcano tomography paper particularly useful for the volcanological community. We have distilled these features into a set of recommendations for best practice. These are intended to allow the non-expert reader to extract the most important information, like the depth, size, and strength of low-velocity anomalies, with relative ease. We recommend that authors:
1. Provide a complete description of the dataset used, including the number and type of stations and sources and the frequency bandwidth of the data. Include a map displaying the geometry of the sources and receivers.
2. Describe the inversion method and procedure in sufficient detail to allow replication of the study or at least to allow the expert reader to confidently assess the robustness of the results in comparison to other studies.
3. Carry out and clearly describe a resolution and sensitivity analysis. Provide a quantitative estimate of the resolution length, preferably using spike tests instead of or in addition to checkerboard tests. Provide synthetic studies that specifically test hypotheses about the ability to recover the spatial extent and magnitude of velocity anomalies that are interpreted.
4. Describe any significant and interesting velocity anomaly in the text, listing its depth, size, and amplitude. Provide a rationale to justify the reference velocity used to calculate the anomaly. List both the absolute velocity and the velocity perturbation at the centre of the anomaly. Discuss how the resolution limits and physical sensitivity affect the recovery of the observed anomaly. Refrain from interpreting anomalies that are not robustly constrained.
5. If a low-velocity anomaly is detected and the resolution is sufficient to warrant further investigation, attempt a quantitative estimate of the melt fraction. Report the method used, and all the assumptions made, including the thermal structure, the assumed composition and elastic properties of the melt and of the host rock, and the melt geometries considered. Report the estimated ranges of melt fraction and melt volume with associated uncertainties.
7 Final remarks
Volcano tomography has contributed greatly to the development of the modern paradigm for magma storage and differentiation. An increasing number of active volcanoes have been imaged with seismic methods and many studies have found evidence of clear low-velocity anomalies that can only be explained by the presence of partial melt. The observation of low-velocity anomalies at multiple levels in the crust supports the emerging concept of transcrustal magmatic systems (Sparks and Cashman, 2017). Most seismological observations seem to be consistent with magma storage dominated by crystal-rich mush; however, estimates of melt fraction are still highly speculative.
We argue that traditional travel-time tomography and surface wave tomography provide weak constraints on melt fractions, with a potential bias toward underestimating. They remain powerful and useful tools for understanding the three-dimensional structure of volcanoes and they can inform the location of magma storage. When applied with care and interpreted within the bounds imposed by resolution and sensitivity limits, they can provide broad constraints on possible magma volumes. To recover robust constraints on melt fractions and answer the many important outstanding questions concerning the nature of magma reservoirs other methods must be employed that use more information derived from the full seismic wavefield.
Determining melt fraction from seismic properties is conceptually very similar to estimating porosity and hydrocarbon saturation in hydrocarbon reservoirs. Therefore, the volcano seismologist has the opportunity to borrow from the array of methods developed by the geophysical exploration industry. In traditional hydrocarbon exploration, tomography is used routinely to recover the long-wavelength structure and multi-channel seismic imaging is used to determine the high-resolution tectonic structure. The most common approach to determine reservoir porosity and hydrocarbon saturation is AVO analysis (Smith and Gidlow, 1987; Castagna and Backus, 1993; Fawad et al., 2020). This method has been used successfully to constrain melt fractions beneath mid-ocean ridges (e.g., Singh et al., 1998). More recently, the exploration industry has embraced full-waveform inversion, which can simultaneously provide high-resolution structural constraints and robust reservoir parameters (Prieux et al., 2013; Warner et al., 2013; Zhang and Alkhalifah, 2020). Where these methods have been applied to magma imaging, they have been successful in revealing small melt bodies and constraining melt fractions more robustly than possible with traditional methods. (e.g., Xu et al., 2014; Chrapkiewicz et al., 2022).
Wider adoption of these techniques is needed to shed new light on magma reservoir properties. Achieving this objective will also require new strategies for collecting dense broad-band seismic data on active volcanoes. Higher-resolution seismic imaging and robust melt fraction constraints have the potential to transform seismic imaging into a diagnostic tool for volcanic hazards. Deeper and higher resolution imaging should be integrated into a multi-disciplinary strategy including information from geophysical monitoring, the geological record, petrology, and numerical modelling to determine the potential for future eruptions and the likely eruption style and size (Giordano and Caricchi, 2022).
MP assembled and analysed the database, generated the figures and wrote the manuscript. KC produced the wavefront healing analysis and created the first version of the GitHub repository and the included Jupyter notebooks. EH, BH, DT, and JM helped interpret the result and write the manuscript.
MP was supported by a NERC Independent Research Fellowship (grant number NE/R015708/1) and EH by the National Science Foundation (grant number OCE-2023338).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Abe, Y., Ohkura, T., Shibutani, T., Hirahara, K., and Kato, M. (2010). Crustal structure beneath Aso Caldera, Southwest Japan, as derived from receiver function analysis. J. Volcanol. Geotherm. Res. 195 (1), 1–12. doi:10.1016/j.jvolgeores.2010.05.011
Abe, Y., Ohkura, T., Shibutani, T., Hirahara, K., Yoshikawa, S., and Inoue, H. (2017). Low‐velocity zones in the crust beneath Aso caldera, Kyushu, Japan, derived from receiver function analyses. J. Geophys. Res. Solid Earth 122 (3), 2013–2033. doi:10.1002/2016jb013686
Achauer, U., Evans, J. R., and Stauber, D. A. (1988). High‐resolution seismic tomography of compressional wave velocity structure at Newberry Volcano, Oregon Cascade Range. J. Geophys. Res. 93 (B9), 10135–10147. doi:10.1029/jb093ib09p10135
Aghaei, O., Nedimović, M. R., Marjanović, M., Carbotte, S. M., Canales, J. P., Carton, H., et al. (2017). Constraints on melt content of off‐axis magma lenses at the East Pacific Rise from analysis of 3‐D seismic amplitude variation with angle of incidence. J. Geophys. Res. Solid Earth 122 (6), 4123–4142. doi:10.1002/2016jb013785
Annen, C., Paulatto, M., Sparks, R. S. J., Minshull, T. A., and Kiddle, E. J. (2014). Quantification of the intrusive magma fluxes during magma chamber growth at Soufrière Hills Volcano (Montserrat, Lesser Antilles). J. Petrology 55 (3), 529–548. doi:10.1093/petrology/egt075
Arnoux, G. M., Toomey, D. R., Hooft, E. E. E., and Wilcock, W. S. D. (2019). Seismic imaging and physical properties of the Endeavour segment: Evidence that skew between mantle and crustal magmatic systems governs spreading center processes. Geochem. Geophys. Geosyst. 20 (3), 1319–1339. doi:10.1029/2018GC007978
Arnulf, A. F., Harding, A. J., Singh, S. C., Kent, G. M., and Crawford, W. C. (2014). Nature of upper crust beneath the Lucky Strike volcano using elastic full waveform inversion of streamer data. Geophys. J. Int. 196 (3), 1471–1491. doi:10.1093/gji/ggt461
Aster, R. C., and Meyer, R. P. (1988). Three-dimensional velocity structure and hypocenter distribution in the Campi Flegrei caldera, Italy. Tectonophysics 149 (3-4), 195–218. doi:10.1016/0040-1951(88)90173-4
Bai, T., Thurber, C., Lanza, F., Singer, B. S., Bennington, N., Keranen, K., et al. (2020). Teleseismic tomography of the Laguna del Maule volcanic field in Chile. JGR. Solid Earth 125 (8), e2020JB019449. doi:10.1029/2020jb019449
Barmin, M. P., Ritzwoller, M. H., and Levshin, A. L. (2001). “A fast and reliable method for surface wave tomography,” in Monitoring the comprehensive nuclear-test-ban treaty: Surface waves (Basel: Birkhäuser), 1351–1375.
Bastow, I. D., Pilidou, S., Kendall, J. M., and Stuart, G. W. (2010). Melt‐induced seismic anisotropy and magma assisted rifting in Ethiopia: Evidence from surface waves. Geochem. Geophys. Geosyst. 11 (6). doi:10.1029/2010gc003036
Battaglia, J., Zollo, A., Virieux, J., and Iacono, D. D. (2008). Merging active and passive data sets in traveltime tomography: The case study of Campi Flegrei caldera (southern Italy). Geophys. Prospect. 56 (4), 555–573. doi:10.1111/j.1365-2478.2007.00687.x
Beachly, M. W., Hooft, E. E., Toomey, D. R., and Waite, G. P. (2012). Upper crustal structure of Newberry volcano from P‐wave tomography and finite difference waveform modeling. J. Geophys. Res. 117 (B10). doi:10.1029/2012jb009458
Bedrosian, P. A., Peacock, J. R., Bowles-Martinez, E., Schultz, A., and Hill, G. J. (2018). Crustal inheritance and a top-down control on arc magmatism at Mount St Helens. Nat. Geosci. 11 (11), 865–870. doi:10.1038/s41561-018-0217-2
Behr, Y., Townend, J., Bannister, S., and Savage, M. K. (2011). Crustal shear wave tomography of the Taupo Volcanic Zone, New Zealand, via ambient noise correlation between multiple three‐component networks. Geochem. Geophys. Geosyst. 12 (3). doi:10.1029/2010gc003385
Benz, H. M., and Smith, R. B. (1984). Simultaneous inversion for lateral velocity variations and hypocenters in the Yellowstone region using earthquake and refraction data. J. Geophys. Res. 89 (B2), 1208–1220. doi:10.1029/jb089ib02p01208
Bodin, T., Sambridge, M., Tkalčić, H., Arroucau, P., Gallagher, K., and Rawlinson, N. (2012b). Transdimensional inversion of receiver functions and surface wave dispersion. J. Geophys. Res. 117 (B2). doi:10.1029/2011jb008560
Burdick, L. J., and Langston, C. A. (1977). Modeling crustal structure through the use of converted phases in teleseismic body-wave forms. Bull. Seismol. Soc. Am. 67 (3), 677–691. doi:10.1785/bssa0670030677
Canales, J. P., Singh, S. C., Detrick, R. S., Carbotte, S. M., Harding, A., Kent, G. M., et al. (2006). Seismic evidence for variations in axial magma chamber properties along the southern Juan de Fuca Ridge. Earth Planet. Sci. Lett. 246 (3-4), 353–366. doi:10.1016/j.epsl.2006.04.032
Carbotte, S. M., Detrick, R. S., Harding, A., Canales, J. P., Babcock, J., Kent, G., et al. (2006). Rift topography linked to magmatism at the intermediate spreading Juan de Fuca Ridge. Geol. 34 (3), 209–212. doi:10.1130/G21969.1
Carbotte, S. M., Arnulf, A., Spiegelman, M., Lee, M., Harding, A., Kent, G., et al. (2020). Stacked sills forming a deep melt-mush feeder conduit beneath Axial Seamount. Geology 48 (7), 693–697. doi:10.1130/g47223.1
Caricchi, L., Burlini, L., and Ulmer, P. (2008). Propagation of P and S-waves in magmas with different crystal contents: Insights into the crystallinity of magmatic reservoirs. J. Volcanol. Geotherm. Res. 178 (4), 740–750. doi:10.1016/j.jvolgeores.2008.09.006
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 355 (6331), eaag3055. doi:10.1126/science.aag3055
Chambers, E. L., Harmon, N., Rychert, C. A., and Keir, D. (2021). Variations in melt emplacement beneath the northern East African Rift from radial anisotropy. Earth Planet. Sci. Lett. 573, 117150. doi:10.1016/j.epsl.2021.117150
Chiarabba, C., and Moretti, M. (2006). An insight into the unrest phenomena at the Campi Flegrei caldera from Vp and Vp/Vs tomography. Terra nova. 18 (6), 373–379. doi:10.1111/j.1365-3121.2006.00701.x
Chrapkiewicz, K., Wilde-Piórko, M., Polkowski, M., and Grad, M. (2020). Reliable workflow for inversion of seismic receiver function and surface wave dispersion data: A “13 BB star” case study. J. Seismol. 24, 101–120. doi:10.1007/s10950-019-09888-1
Chrapkiewicz, K., Paulatto, M., Heath, B., Hooft, E. E. E., Nomikou, P., Papazachos, C. B., et al. (2022). Magma chamber detected beneath an active arc volcano with high-resolution velocity images. pre-print. doi:10.31223/X5934R
Christopher, T. E., Blundy, J., Cashman, K., Cole, P., Edmonds, M., Smith, P. J., et al. (2015). Crustal‐scale degassing due to magma system destabilization and magma‐gas decoupling at S oufrière H ills V olcano, M ontserrat. Geochem. Geophys. Geosyst. 16 (9), 2797–2811. doi:10.1002/2015gc005791
Chung, D. H., Swica, J. J., and Crandall, W. B. (1963). Relation of single‐crystal elastic Constants to polycrystalline isotropic elastic moduli of MgO. J. Am. Ceram. Soc. 46 (9), 452–457. doi:10.1111/j.1151-2916.1963.tb11775.x
Comeau, M. J., Unsworth, M. J., and Cordell, D. (2016). New constraints on the magma distribution and composition beneath Volcán Uturuncu and the southern Bolivian Altiplano from magnetotelluric data. Geosphere 12 (5), 1391–1421. doi:10.1130/ges01277.1
Cornwell, D. G., Maguire, P. K. H., England, R. W., and Stuart, G. W. (2010). Imaging detailed crustal structure and magmatic intrusion across the Ethiopian Rift using a dense linear broadband array. Geochem. Geophys. Geosyst. 11 (1). doi:10.1029/2009gc002637
Costa, A., Caricchi, L., and Bagdassarov, N. (2009). A model for the rheology of particle‐bearing suspensions and partially molten rocks. Geochem. Geophys. Geosyst. 10 (3), 2138. doi:10.1029/2008gc002138
Dawson, P. B., Evans, J. R., and Iyer, H. M. (1990). Teleseismic tomography of the compressional wave velocity structure beneath the Long Valley region, California. J. Geophys. Res. 95 (B7), 11021–11050. doi:10.1029/jb095ib07p11021
De Siena, L., Thomas, C., Waite, G. P., Moran, S. C., and Klemme, S. (2014). Attenuation and scattering tomography of the deep plumbing system of Mount St. Helens. J. Geophys. Res. Solid Earth 119 (11), 8223–8238. doi:10.1002/2014JB011372
Detrick, R. S., Buhl, P., Vera, E., Mutter, J., Orcutt, J., Madsen, J., et al. (1987). Multi-channel seismic imaging of a crustal magma chamber along the East Pacific Rise. Nature 326 (6108), 35–41. doi:10.1038/326035a0
Díaz-Moreno, A., Barberi, G., Cocina, O., Koulakov, I., Scarfì, L., Zuccarello, L., et al. (2018). New insights on Mt. Etna’s crust and relationship with the regional tectonic framework from joint active and passive P-wave seismic tomography. Surv. Geophys. 39 (1), 57–97. doi:10.1007/s10712-017-9425-3
Didonna, R., Costa, F., Handley, H., Turner, S., and Barclay, J. (2022). Dynamics and timescales of mafic–silicic magma interactions at Soufrière Hills Volcano, Montserrat. Contrib. Mineral. Pet.177 (2), 28–17. doi:10.1007/s00410-022-01891-z
Dimitriadis, I., Papazachos, C., Panagiotopoulos, D., Hatzidimitriou, P., Bohnhoff, M., Rische, M., et al. (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. J. Volcanol. Geotherm. Res. 195 (1), 13–30. doi:10.1016/j.jvolgeores.2010.05.013
Druitt, T. H., Costa, F., Deloule, E., Dungan, M., and Scaillet, B. (2012). Decadal to monthly timescales of magma transfer and reservoir growth at a caldera volcano. Nature 482 (7383), 77–80. doi:10.1038/nature10706
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. J. Geophys. Res. 105 (B10), 23537–23555. doi:10.1029/2000jb900210
Dunn, R. A., Lekić, V., Detrick, R. S., and Toomey, D. R. (2005). Three‐dimensional seismic structure of the Mid‐Atlantic Ridge (35 N): Evidence for focused melt supply and lower crustal dike injection. J. Geophys. Res. 110 (B9). doi:10.1029/2004jb003473
Eshelby, J. D. (1957). The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R. Soc. Lond. Ser. A. Math. Phys. Sci. 241 (1226), 376–396. doi:10.1098/rspa.1958.0133
Farrell, J., Smith, R. B., Husen, S., and Diehl, T. (2014). Tomography from 26 years of seismicity revealing that the spatial extent of the Yellowstone crustal magma reservoir extends well beyond the Yellowstone caldera. Geophys. Res. Lett. 41 (9), 3068–3073. doi:10.1002/2014gl059588
Flinders, A. F., Shelly, D. R., Dawson, P. B., Hill, D. P., Tripoli, B., and Shen, Y. (2018). Seismic evidence for significant melt beneath the Long Valley caldera, California, USA. Geology 46 (9), 799–802. doi:10.1130/g45094.1
Forsyth, D. W., and Li, A. (2005). “Array analysis of two-dimensional variations in surface wave phase velocity and azimuthal anisotropy in the presence of multipathing interference,” in Seismic Earth: Array analysis of broadband seismograms. Geophysical monograph. Editors A. Levander, and G. Nolet (Washington, D.C: AGU), 157, 81–97. doi:10.1029/157GM0
Garapić, G., Faul, U. H., and Brisson, E. (2013). High‐resolution imaging of the melt distribution in partially molten upper mantle rocks: Evidence for wetted two‐grain boundaries. Geochem. Geophys. Geosyst. 14 (3), 556–566. doi:10.1029/2012gc004547
García‐Yeguas, A., Koulakov, I., Ibáñez, J. M., and Rietbrock, A. (2012). High resolution 3D P wave velocity structure beneath Tenerife Island (Canary Islands, Spain) based on tomographic inversion of active‐source data. J. Geophys. Res. 117 (B9), 8970. doi:10.1029/2011jb008970
García-Yeguas, A., Ibánez, J. M., Koulakov, I., Jakovlev, A., Romero-Ruiz, M. C., and Prudencio, J. (2014). Seismic tomography model reveals mantle magma sources of recent volcanic activity at El Hierro Island (Canary Islands, Spain). Geophys. J. Int. 199 (3), 1739–1750. doi:10.1093/gji/ggu339
Giordano, G., and Caricchi, L. (2022). Determining the state of activity of transcrustal magmatic systems and their volcanoes. Annu. Rev. Earth Planet. Sci. 50, 231–259. doi:10.1146/annurev-earth-032320-084733
Godey, S., Snieder, R., Villaseñor, A., and Benz, H. M. (2003). Surface wave tomography of north America and the caribbean using global and regional broad-band networks: Phase velocity maps and limitations of ray theory. Geophys. J. Int. 152 (3), 620–632. doi:10.1046/j.1365-246x.2003.01866.x
Greenfield, T., White, R. S., and Roecker, S. (2016). The magmatic plumbing system of the Askja central volcano, Iceland, as imaged by seismic tomography. J. Geophys. Res. Solid Earth 121 (10), 7211–7229. doi:10.1002/2016jb013163
Gudmundsson, O., Brandsdóttir, B., Menke, W., and Sigvaldason, G. E. (1994). The crustal magma chamber of the Katla volcano in south Iceland revealed by 2-D seismic undershooting. Geophys. J. Int. 119 (1), 277–296. doi:10.1111/j.1365-246x.1994.tb00928.x
Guidarelli, M., Sarao, A., and Panza, G. F. (2002). Surface wave tomography and seismic source studies at Campi Flegrei (Italy). Phys. Earth Planet. Interiors 134 (3-4), 157–173. doi:10.1016/s0031-9201(02)00154-1
Hammond, J. O., Kendall, J. M., Wookey, J., Stuart, G. W., Keir, D., and Ayele, A. (2014). Differentiating flow, melt, or fossil seismic anisotropy beneath Ethiopia. Geochem. Geophys. Geosyst. 15 (5), 1878–1894. doi:10.1002/2013gc005185
Hammond, J. O., Wu, J. P., Ri, K. S., Wei, W., Yu, J. N., and Oppenheimer, C. (2020). Distribution of partial melt beneath Changbaishan/Paektu Volcano, China/democratic people's Republic of Korea. Geochem. Geophys. Geosyst. 21 (1), e2019GC008461. doi:10.1029/2019gc008461
Hammond, W. C., and Toomey, D. R. (2003). Seismic velocity anisotropy and heterogeneity beneath the mantle electromagnetic and tomography experiment (MELT) region of the East pacific rise from analysis of P and S body waves. J. Geophys. Res. 108 (B4), 1789. doi:10.1029/2002jb001789
Haslinger, F., Thurber, C., Mandernach, M., and Okubo, P. (2001). Tomographic image of P‐velocity structure beneath Kilauea's East Rift Zone and South Flank: Seismic evidence for a deep magma body. Geophys. Res. Lett. 28 (2), 375–378. doi:10.1029/2000gl012018
Heath, B. A., Hooft, E. E., Toomey, D. R., and Bezada, M. J. (2015). Imaging the magmatic system of Newberry Volcano using joint active source and teleseismic tomography. Geochem. Geophys. Geosyst. 16 (12), 4433–4448. doi:10.1002/2015GC006129
Heath, B. A., Hooft, E. E., and Toomey, D. R. (2018). Autocorrelation of the seismic wavefield at Newberry Volcano: Reflections from the magmatic and geothermal systems. Geophys. Res. Lett. 45 (5), 2311–2318. doi:10.1002/2017GL076706
Heath, B. A., Hooft, E. E. E., Toomey, D. R., Papazachos, C. B., Nomikou, P., Paulatto, M., et al. (2019). Tectonism and its relation to magmatism around Santorini Volcano from upper crustal P wave velocity. J. Geophys. Res. Solid Earth 124 (10), 10610–10629. doi:10.1029/2019jb017699
Heath, B. A., Hooft, E. E. E., Toomey, D. R., Paulatto, M., Papazachos, C. B., Nomikou, P., et al. (2021). Relationship between active faulting/fracturing and magmatism around Santorini: Seismic anisotropy from an active source tomography experiment. JGR. Solid Earth 126 (8), e2021JB021898. doi:10.1029/2021jb021898
Hobro, J. W., Singh, S. C., and Minshull, T. A. (2003). Three-dimensional tomographic inversion of combined reflection and refraction seismic traveltime data. Geophys. J. Int. 152 (1), 79–93. doi:10.1046/j.1365-246x.2003.01822.x
Hooft, E. E., Detrick, R. S., and Kent, G. M. (1997). Seismic structure and indicators of magma budget along the southern East Pacific Rise. J. Geophys. Res. 102 (B12), 27319–27340. doi:10.1029/97jb02349
Hooft, E. E. E., Detrick, R. S., Toomey, D. R., Collins, J. A., and Lin, J. (2000). Crustal thickness and structure along three contrasting spreading segments of the Mid‐Atlantic Ridge, 33.5–35 N. J. Geophys. Res. 105 (B4), 8205–8226. doi:10.1029/1999jb900442
Hooft, E. E. E., Heath, B. A., Toomey, D. R., Paulatto, M., Papazachos, C. B., Nomikou, P., et al. (2019). Seismic imaging of Santorini: Subsurface constraints on caldera collapse and present-day magma recharge. Earth Planet. Sci. Lett. 514, 48–61. doi:10.1016/j.epsl.2019.02.033
Huang, H. H., Lin, F. C., Schmandt, B., Farrell, J., Smith, R. B., and Tsai, V. C. (2015). The Yellowstone magmatic system from the mantle plume to the upper crust. Science 348 (6236), 773–776. doi:10.1126/science.aaa5648
Huang, Y. C., Ohkura, T., Kagiyama, T., Yoshikawa, S., and Inoue, H. (2018). Shallow volcanic reservoirs and pathways beneath Aso caldera revealed using ambient seismic noise tomography. Earth Planets Space 70 (1), 169–221. doi:10.1186/s40623-018-0941-2
Hughes, G. E., Petrone, C. M., Downes, H., Varley, N. R., and Hammond, S. J. (2021). Mush remobilisation and mafic recharge: A study of the crystal cargo of the 2013–17 eruption at volcán de Colima, Mexico. J. Volcanol. Geotherm. Res. 416, 107296. doi:10.1016/j.jvolgeores.2021.107296
Husen, S., Smith, R. B., and Waite, G. P. (2004). Evidence for gas and magmatic sources beneath the Yellowstone volcanic field from seismic tomographic imaging. J. Volcanol. Geotherm. Res. 131 (3-4), 397–410. doi:10.1016/s0377-0273(03)00416-5
Hussenoeder, S. A., Collins, J. A., Kent, G. M., and Detrick, R. S. (1996). Seismic analysis of the axial magma chamber reflector along the southern East Pacific Rise from conventional reflection profiling. J. Geophy. Res.: Solid Earth 101 (B10), 22087–22105.
Indrastuti, N., Nugraha, A. D., McCausland, W. A., Hendrasto, M., Gunawan, H., Kusnandar, R., et al. (2019). 3-D seismic tomographic study of Sinabung volcano, northern sumatra, Indonesia, during the inter-eruptive period october 2010–july 2013. J. Volcanol. Geotherm. Res. 382, 197–209. doi:10.1016/j.jvolgeores.2019.03.001
Iyer, H. M., Evans, J. R., Zandt, G., Stewart, R. M., Coakley, J. M., and Roloff, J. N. (1981). A deep low-velocity body under the Yellowstone caldera, Wyoming: Delineation using teleseismic P-wave residuals and tectonic interpretation. GSA Bull. 92, 1471–1646. doi:10.1130/gsab-p2-92-1471
Jackson, M. D., Blundy, J., and Sparks, R. S. J. (2018). Chemical differentiation, cold storage and remobilization of magma in the Earth’s crust. Nature 564 (7736), 405–409. doi:10.1038/s41586-018-0746-2
Janiszewski, H. A., Wagner, L. S., and Roman, D. C. (2020). Aseismic mid-crustal magma reservoir at Cleveland Volcano imaged through novel receiver function analyses. Sci. Rep. 10 (1), 1780–1789. doi:10.1038/s41598-020-58589-0
Jaxybulatov, K., Shapiro, N. M., Koulakov, I., Mordret, A., Landès, M., and Sens-Schönfelder, C. (2014). A large magmatic sill complex beneath the Toba caldera. science 346 (6209), 617–619. doi:10.1126/science.1258582
Jeddi, Z., Tryggvason, A., and Gudmundsson, Ó. (2016). The Katla volcanic system imaged using local earthquakes recorded with a temporary seismic network. J. Geophys. Res. Solid Earth 121 (10), 7230–7251. doi:10.1002/2016jb013044
Jónsdóttir, K., Tryggvason, A., Roberts, R., Lund, B., Soosalu, H., and Böðvarsson, R. (2007). Habits of a glacier-covered volcano: Seismicity patterns and velocity structure of Katla volcano, Iceland. Ann. Glaciol. 45, 169–177. doi:10.3189/172756407782282499
Kiddle, E. J., Edwards, B. R., Loughlin, S. C., Petterson, M., Sparks, R. S. J., and Voight, B. (2010). Crustal structure beneath Montserrat, Lesser Antilles, constrained by xenoliths, seismic velocity structure and petrology. Geophys. Res. Lett. 37 (19), 42145. doi:10.1029/2009gl042145
Kiser, E., Palomeras, I., Levander, A., Zelt, C., Harder, S., Schmandt, B., et al. (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 (6), 411–414. doi:10.1130/g37591.1
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 46 (9), 775–778. doi:10.1130/G45140.1
Klemetti, E. W., and Clynne, M. A. (2014). Localized rejuvenation of a crystal mush recorded in zircon temporal and compositional variation at the Lassen Volcanic Center, Northern California. PloS one 9 (12), e113157. doi:10.1371/journal.pone.0113157
Korenaga, J., Holbrook, W. S., Kent, G. M., Kelemen, P. B., Detrick, R. S., Larsen, H. C., et al. (2000). Crustal structure of the southeast Greenland margin from joint refraction and reflection seismic tomography. J. Geophys. Res. 105 (B9), 21591–21614. doi:10.1029/2000jb900188
Koulakov, I., Kasatkina, E., Shapiro, N. M., Jaupart, C., Vasilevsky, A., El Khrepy, S., et al. (2016). The feeder system of the Toba supervolcano from the slab to the shallow reservoir. Nat. Commun. 7 (1). doi:10.1038/ncomms12228
Koulakov, I., Yudistira, T., Luehr, B. G., and Wandono, (2009). P, S velocity and VP/VS ratio beneath the Toba caldera complex (Northern Sumatra) from local earthquake tomography. Geophys. J. Int. 177 (3), 1121–1139. doi:10.1111/j.1365-246x.2009.04114.x
Kounoudis, R., Bastow, I. D., Ebinger, C. J., Ogden, C. S., Ayele, A., Bendick, R., et al. (2021). Body-wave tomographic imaging of the Turkana Depression: Implications for rift development and plume–lithosphere interactions. Geochem. Geophys. Geosyst. 22 (8), e2021GC009782. doi:10.1029/2021gc009782
Krastel, S., and Schmincke, H. U. (2002). Crustal structure of northern Gran Canaria, Canary Islands, deduced from active seismic tomography. J. Volcanol. Geotherm. Res. 115 (1-2), 153–177. doi:10.1016/s0377-0273(01)00313-4
Kukarina, E., West, M., Keyson, L. H., Koulakov, I., Tsibizov, L., and Smirnov, S. (2017). Focused magmatism beneath Uturuncu volcano, Bolivia: Insights from seismic tomography and deformation modeling. Geosphere 13 (6), 1855–1866. doi:10.1130/ges01403.1
Kyong-Song, R., Hammond, J. O., Chol-Nam, K., Hyok, K., Yong-Gun, Y., Gil-Jong, P., et al. (2016). Evidence for partial melt in the crust beneath Mt. Paektu (changbaishan), democratic people’s republic of Korea and China. Sci. Adv. 2 (4), e1501513. doi:10.1126/sciadv.1501513
Lees, J. M., and Crosson, R. S. (1989). Tomographic inversion for three‐dimensional velocity structure at Mount St. Helens using earthquake data. J. Geophys. Res. 94 (B5), 5716–5728. doi:10.1029/jb094ib05p05716
Lin, G., Shearer, P. M., Matoza, R. S., Okubo, P. G., and Amelung, F. (2014). Three‐dimensional seismic velocity structure of Mauna Loa and Kilauea volcanoes in Hawaii from local seismic tomography. J. Geophys. Res. Solid Earth 119 (5), 4377–4392. doi:10.1002/2013jb010820
Luehr, B. G., Koulakov, I., Rabbel, W., Zschau, J., Ratdomopurbo, A., Brotopuspito, K. S., et al. (2013). Fluid ascent and magma storage beneath Gunung Merapi revealed by multi-scale seismic imaging. J. Volcanol. Geotherm. Res. 261, 7–19. doi:10.1016/j.jvolgeores.2013.03.015
MacKenzie, L. S., Abers, G. A., Rondenay, S., and Fischer, K. M. (2010). Imaging a steeply dipping subducting slab in Southern Central America. Earth Planet. Sci. Lett. 296 (3-4), 459–468. doi:10.1016/j.epsl.2010.05.033
Magde, L. S., Barclay, A. H., Toomey, D. R., Detrick, R. S., and Collins, J. A. (2000). Crustal magma plumbing within a segment of the Mid-Atlantic Ridge, 35 N. Earth Planet. Sci. Lett. 175 (1-2), 55–67. doi:10.1016/s0012-821x(99)00281-2
Martí, J., Villaseñor, A., Geyer, A., López, C., and Tryggvason, A. (2017). Stress barriers controlling lateral migration of magma revealed by seismic tomography. Sci. Rep. 7 (1). doi:10.1038/srep40757
Martinez‐Arevalo, C., Musumeci, C., and Patane, D. (2009). Evidence of a partial melt zone beneath Stromboli volcano (Italy) from inversion of teleseismic receiver functions. Terra nova.21 (5), 386–392. doi:10.1111/j.1365-3121.2009.00894.x
Mastruyono, M., McCaffrey, R., Wark, D. A., Roecker, S. W., Ibrahim, G., Fauzi, , et al. (2001). Distribution of magma beneath the Toba caldera complex, north Sumatra, Indonesia, constrained by three‐dimensional P wave velocities, seismicity, and gravity data. Geochem. Geophys. Geosyst. 2 (4). doi:10.1029/2000gc000096
McVey, B. G., Hooft, E. E. E., Heath, B. A., Toomey, D. R., Paulatto, M., Morgan, J. V., et al. (2020). Magma accumulation beneath Santorini volcano, Greece, from P-wave tomography. Geology 48 (3), 231–235. doi:10.1130/g47127.1
Meléndez, A., Korenaga, J., Sallarès, V., Miniussi, A., and Ranero, C. R. (2015). TOMO3D: 3-D joint refraction and reflection traveltime tomography parallel code for active-source seismic data—synthetic test. Geophys. J. Int. 203 (1), 158–174. doi:10.1093/gji/ggv292
Miller, D. S., and Smith, R. B. (1999). P and S velocity structure of the Yellowstone volcanic field from local earthquake and controlled‐source tomography. J. Geophys. Res. 104 (B7), 15105–15121. doi:10.1029/1998jb900095
Morgan, J., Warner, M., Bell, R., Ashley, J., Barnes, D., Little, R., et al. (2013). Next-generation seismic experiments: Wide-angle, multi-azimuth, three-dimensional, full-waveform inversion. Geophys. J. Int. 195 (3), 1657–1678. doi:10.1093/gji/ggt345
Morgan, J., Warner, M., Arnoux, G., Hooft, E., Toomey, D., VanderBeek, B., et al. (2016). Next-generation seismic experiments–II: Wide-angle, multi-azimuth, 3-D, full-waveform inversion of sparse field data. Geophys. J. Int. 204 (2), 1342–1363. doi:10.1093/gji/ggv513
Natale, M., Nunziata, C., and Panza, G. F. (2005). verage shear wave velocity models of the crustal structure at Mt. Vesuvius. Phys. Earth Planet. Inter. 152 (1–2), 7–21. doi:10.1016/j.pepi.2005.03.011
Nurfiani, D., Wang, X., Gunawan, H., Triastuty, H., Hidayat, D., Wei, S. J., et al. (2021). Combining petrology and seismology to unravel the plumbing system of a typical arc volcano: An example from Marapi, West Sumatra, Indonesia. Geochem. Geophys. Geosyst. 22 (4), e2020GC009524. doi:10.1029/2020gc009524
Ochoa-Chávez, J. A., Escudero, C. R., Núñez-Cornú, F. J., and Bandy, W. L. (2016). P-wave velocity tomography from local earthquakes in Western Mexico. Pure Appl. Geophys. 173 (10), 3487–3511. doi:10.1007/s00024-015-1183-x
Ohlendorf, S. J., Thurber, C. H., Pesicek, J. D., and Prejean, S. G. (2014). Seismicity and seismic structure at okmok volcano, Alaska. J. Volcanol. Geotherm. Res. 278, 103–119. doi:10.1016/j.jvolgeores.2014.04.002
Okubo, P. G., Benz, H. M., and Chouet, B. A. (1997). Imaging the crustal magma sources beneath Mauna Loa and Kilauea volcanoes, Hawaii. Geol. 25 (10), 867–870. doi:10.1130/0091-7613(1997)025<0867:itcmsb>2.3.co;2
Pallister, J. S., and Hopson, C. A. (1981). Samail ophiolite plutonic suite: Field relations, phase variation, cryptic variation and layering, and a model of a spreading ridge magma chamber. J. Geophys. Res. 86 (B4), 2593–2644. doi:10.1029/jb086ib04p02593
Park, J., Morgan, J. K., Zelt, C. A., Okubo, P. G., Peters, L., and Benesh, N. (2007). Comparative velocity structure of active Hawaiian volcanoes from 3-D onshore–offshore seismic tomography. Earth Planet. Sci. Lett. 259 (3-4), 500–516. doi:10.1016/j.epsl.2007.05.008
Park, J., Morgan, J. K., Zelt, C. A., and Okubo, P. G. (2009). Volcano‐tectonic implications of 3‐D velocity structures derived from joint active and passive source tomography of the island of Hawaii. J. Geophys. Res. 114 (B9), B09301. doi:10.1029/2008jb005929
Parks, M. M., Biggs, J., England, P., Mather, T. A., Nomikou, P., Palamartchouk, K., et al. (2012). Evolution of Santorini Volcano dominated by episodic and rapid fluxes of melt from depth. Nat. Geosci. 5 (10), 749–754. doi:10.1038/ngeo1562
Paulatto, M., Annen, C., Henstock, T. J., Kiddle, E., Minshull, T. A., Sparks, R. S. J., et al. (2012). Magma chamber properties from integrated seismic tomography and thermal modeling at Montserrat. Geochem. Geophys. Geosyst. 13 (1). doi:10.1029/2011gc003892
Paulatto, M., Moorkamp, M., Hautmann, S., Hooft, E., Morgan, J. V., and Sparks, R. S. J. (2019). Vertically extensive magma reservoir revealed from joint inversion and quantitative interpretation of seismic and gravity data. J. Geophys. Res. Solid Earth 124 (11), 11170–11191. doi:10.1029/2019jb018476
Prieux, V., Brossier, R., Operto, S., and Virieux, J. (2013). Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from the valhall field. Part 2: Imaging compressive-wave and shear-wave velocities. Geophys. J. Int. 194 (3), 1665–1681. doi:10.1093/gji/ggt178
Prudencio, J., Ibáñez, J. M., Del Pezzo, E., Martí, J., García-Yeguas, A., and De Siena, L. (2015). 3D attenuation tomography of the volcanic island of Tenerife (Canary Islands). Surv. Geophys. 36 (5), 693–716. doi:10.1007/s10712-015-9333-3
Rosi, M., Acocella, V., Cioni, R., Bianco, F., Costa, A., De Martino, P., et al. (2022). Defining the pre-eruptive states of active volcanoes for improving eruption forecasting. Front. Earth Sci. (Lausanne). 10, 795700. doi:10.3389/feart.2022.795700
Schmeling, H. (1985). Numerical models on the influence of partial melt on elastic, anelastic and electric properties of rocks. Part I: Elasticity and anelasticity. Phys. earth Planet. interiors 41 (1), 34–57. doi:10.1016/0031-9201(85)90100-1
Schmitz, M., Heinsohn, W. D., and Schilling, F. R. (1997). Seismic, gravity and petrological evidence for partial melt beneath the thickened central Andean crust (21–23 S). Tectonophysics 270 (3-4), 313–326. doi:10.1016/s0040-1951(96)00217-x
Seccia, D., Chiarabba, C., De Gori, P., Bianchi, I., and Hill, D. P. (2011). Evidence for the contemporary magmatic system beneath Long Valley Caldera from local earthquake tomography and receiver function analysis. J. Geophys. Res. 116 (B12), B12314. doi:10.1029/2011jb008471
Shalev, E., Kenedi, C. L., Malin, P., Voight, V., Miller, V., Hidayat, D., et al. (2010). Three‐dimensional seismic velocity tomography of Montserrat from the SEA‐CALIPSO offshore/onshore experiment. Geophys. Res. Lett. 37 (19). doi:10.1029/2010gl042498
Singh, S. C., Kent, G. M., Collier, J. S., Harding, A. J., and Orcutt, J. A. (1998). Melt to mush variations in crustal magma properties along the ridge crest at the southern East Pacific Rise. Nature 394 (6696), 874–878. doi:10.1038/29740
Sparks, R. S. J., Annen, C., Blundy, J. D., Cashman, K. V., Rust, A. C., and Jackson, M. D. (2019). Formation and dynamics of magma reservoirs. Phil. Trans. R. Soc. A 377 (2139), 20180019. doi:10.1098/rsta.2018.0019
Spica, Z., Cruz-Atienza, V. M., Reyes-Alfaro, G., Legrand, D., and Iglesias, A. (2014). Crustal imaging of Western Michoacán and the Jalisco Block, Mexico, from ambient seismic noise. J. Volcanol. Geotherm. Res. 289, 193–201. doi:10.1016/j.jvolgeores.2014.11.005
Stachnik, J. C., Dueker, K., Schutt, D. L., and Yuan, H. (2008). Imaging Yellowstone plume‐lithosphere interactions from inversion of ballistic and diffusive Rayleigh wave dispersion and crustal thickness data. Geochem. Geophys. Geosyst. 9 (6). doi:10.1029/2008gc001992
Stankiewicz, J., Ryberg, T., Haberland, C., and Natawidjaja, D. (2010). Lake Toba volcano magma chamber imaged by ambient seismic noise tomography. Geophys. Res. Lett. 37 (17). doi:10.1029/2010gl044211
Stauber, D. A., Green, S. M., and Iyer, H. M. (1988). Three‐dimensional P velocity structure of the crust below Newberry Volcano, Oregon. J. Geophys. Res. 93 (B9), 10095–10107. doi:10.1029/jb093ib09p10095
Steck, L. K., and Prothero, W. A. (1994). Crustal structure beneath Long Valley caldera from modeling of teleseismic P wave polarizations and Ps converted waves. J. Geophys. Res. 99 (B4), 6881–6898. doi:10.1029/93jb03284
Steck, L. K., Thurber, C. H., Fehler, M. C., Lutter, W. J., Roberts, P. M., Baldridge, W. S., et al. (1998). Crust and upper mantle P wave velocity structure beneath valles caldera, new Mexico: Results from the jemez teleseismic tomography experiment. J. Geophys. Res. 103 (B10), 24301–24320. doi:10.1029/98jb00750
Sychev, I. V., Koulakov, I., Egorushkin, I., Zhuravlev, S., West, M., El Khrepy, S., et al. (2019). Fault‐associated magma conduits beneath Volcán de Colima revealed by seismic velocity and attenuation tomography studies. J. Geophys. Res. Solid Earth 124 (8), 8908–8923. doi:10.1029/2019jb017449
Tape, C., Liu, Q., and Tromp, J. (2007). Finite‐frequency tomography using adjoint methods—methodology and examples using membrane surface waves. Geophys. J. Int. 168 (3), 1105–1129. doi:10.1111/j.1365-246X.2006.03191.x
Tavazzani, L., Peres, S., Sinigoi, S., Demarchi, G., Economos, R. C., and Quick, J. E. (2020). Timescales and mechanisms of crystal-mush rejuvenation and melt extraction recorded in Permian plutonic and volcanic rocks of the Sesia magmatic system (Southern Alps, Italy). J. Petrology 61 (5), egaa049. doi:10.1093/petrology/egaa049
Tolstoy, M., Harding, A. J., and Orcutt, J. A. (1993). Crustal thickness on the Mid-Atlantic Ridge: Bull's-eye gravity anomalies and focused accretion. Science 262 (5134), 726–729. doi:10.1126/science.262.5134.726
Toomey, D. R., Solomon, S. C., and Purdy, G. M. (1994). Tomographic imaging of the shallow crustal structure of the East Pacific Rise at 9° 30′ N. J. Geophys. Res. 99 (B12), 24135–24157. doi:10.1029/94jb01942
Toomey, D. R., Jousselin, D., Dunn, R. A., Wilcock, W. S., and Detrick, R. S. (2007). Skew of mantle upwelling beneath the East Pacific Rise governs segmentation. Nature 446 (7134), 409–414. doi:10.1038/nature05679
Ulberg, C. W., Creager, K. C., Moran, S. C., Abers, G. A., Thelen, W. A., Levander, A., et al. (2020). Local source Vp and vs tomography in the Mount St. Helens region with the iMUSH broadband array. Geochem. Geophys. Geosyst. 21 (3), e2019GC008888. doi:10.1029/2019gc008888
Van Ark, E. M., Detrick, R. S., Canales, J. P., Carbotte, S. M., Harding, A. J., Kent, G. M., et al. (2007). Seismic structure of the endeavour segment, juan de Fuca ridge: Correlations with seismicity and hydrothermal activity. J. Geophys. Res. 112 (B2), B02401. doi:10.1029/2005JB004210
Van Avendonk, H. J., Harding, A. J., Orcutt, J. A., and Holbrook, W. S. (2001). Hybrid shortest path and ray bending method for traveltime and raypath calculations. Geophysics 66 (2), 648–653. doi:10.1190/1.1444955
VanDecar, J. C., James, D. E., and Assumpção, M. (1995). Seismic evidence for a fossil mantle plume beneath South America and implications for plate driving forces. Nature 378 (6552), 25–31. doi:10.1038/378025a0
Vanorio, T., Virieux, J., Capuano, P., and Russo, G. (2005). Three‐dimensional seismic tomography from P wave and S wave microearthquake travel times and rock physics characterization of the Campi Flegrei Caldera. J. Geophys. Res. 110 (B3), B03201. doi:10.1029/2004jb003102
Vera, E. E., Mutter, J. C., Buhl, P., Orcutt, J. A., Harding, A. J., Kappus, M. E., et al. (1990). The structure of 0‐to 0.2‐my‐old oceanic crust at 9° N on the East Pacific Rise from expanded spread profiles. J. Geophys. Res. 95 (B10), 15529–15556. doi:10.1029/jb095ib10p15529
Waite, G. P., and Moran, S. C. (2009). VP Structure of Mount St. Helens, Washington, USA, imaged with local earthquake tomography. J. Volcanol. Geotherm. Res. 182 (1-2), 113–122. doi:10.1016/j.jvolgeores.2009.02.009
Ward, K. M., Porter, R. C., Zandt, G., Beck, S. L., Wagner, L. S., Minaya, E., et al. (2013). Ambient noise tomography across the central andes. Geophys. J. Int. 194 (3), 1559–1573. doi:10.1093/gji/ggt166
Ward, K. M., Zandt, G., Beck, S. L., Christensen, D. H., and 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 Planet. Sci. Lett. 404, 43–53. doi:10.1016/j.epsl.2014.07.022
Weiland, C. M., Steck, L. K., Dawson, P. B., and Korneev, V. A. (1995). Nonlinear teleseismic tomography at Long Valley Caldera, using three‐dimensional minimum travel time ray tracing. J. Geophys. Res. 100 (B10), 20379–20390. doi:10.1029/95jb01147
Wespestad, C. E., Thurber, C. H., Andersen, N. L., Singer, B. S., Cardona, C., Zeng, X., et al. (2019). Magma reservoir below Laguna del Maule volcanic field, Chile, imaged with surface‐wave tomography. J. Geophys. Res. Solid Earth 124 (3), 2858–2872. doi:10.1029/2018jb016485
Xu, M., Pablo Canales, J., Carbotte, S. M., Carton, H., Nedimović, M. R., and Mutter, J. C. (2014). Variations in axial magma lens properties along the East Pacific Rise (9° 30′ N–10° 00′ N) from swath 3‐D seismic imaging and 1‐D waveform inversion. J. Geophys. Res. Solid Earth 119 (4), 2721–2744. doi:10.1002/2013jb010730
Yang, B., Lin, W., Hu, X., Fang, H., Qiu, G., and Wang, G. (2021). The magma system beneath Changbaishan-Tianchi Volcano, China/North Korea: Constraints from three-dimensional magnetotelluric imaging. J. Volcanol. Geotherm. Res. 419, 107385. doi:10.1016/j.jvolgeores.2021.107385
Yang, Y., and Forsyth, D. W. (2006). Regional tomographic inversion of the amplitude and phase of Rayleigh waves with 2-D sensitivity kernels. Geophys. J. Int. 166, 1148–1160. doi:10.1111/j.1365-246X.2006.02972.x
Yu, X., and Lee, C. T. A. (2016). Critical porosity of melt segregation during crustal melting: Constraints from zonation of peritectic garnets in a dacite volcano. Earth Planet. Sci. Lett. 449, 127–134. doi:10.1016/j.epsl.2016.05.025
Zelt, C. A., and Barton, P. J. (1998). Three‐dimensional seismic refraction tomography: A comparison of two methods applied to data from the faeroe basin. J. Geophys. Res. 103 (B4), 7187–7210. doi:10.1029/97jb03536
Zhang, Z. D., and Alkhalifah, T. (2020). High-resolution reservoir characterization using deep learning-aided elastic full-waveform inversion: The North Sea field data example. Geophysics 85 (4), WA137–WA146. doi:10.1190/geo2019-0340.1
Zollo, A., Judenherc, S., Auger, E., D'Auria, L., Virieux, J., Capuano, P., et al. (2003). Evidence for the buried rim of Campi Flegrei caldera from 3‐d active seismic imaging. Geophys. Res. Lett. 30 (19), 2002. doi:10.1029/2003gl018173
Zollo, A., Maercklin, N., Vassallo, M., Dello Iacono, D., Virieux, J., and Gasparini, P. (2008). Seismic reflections reveal a massive melt layer feeding Campi Flegrei caldera. Geophys. Res. Lett. 35 (12). doi:10.1029/2008gl034242
Zou, B., and Ma, C. (2020). Crystal mush rejuvenation induced by heat and water transfer: Evidence from amphibole analyses in the jialuhe composite pluton, East kunlun orogen, northern tibet plateau. Lithos 376, 105722. doi:10.1016/j.lithos.2020.105722
Keywords: volcanoes, crystal mush, magma storage, geophysics, seismic tomography, resolution, melt fraction, wavefront healing
Citation: Paulatto M, Hooft EEE, Chrapkiewicz K, Heath B, Toomey DR and Morgan JV (2022) Advances in seismic imaging of magma and crystal mush. Front. Earth Sci. 10:970131. doi: 10.3389/feart.2022.970131
Received: 15 June 2022; Accepted: 05 October 2022;
Published: 20 October 2022.
Edited by:Antonio Costa, National Institute of Geophysics and Volcanology, Italy
Reviewed by:Michael Marani, National Research Council (CNR), Italy
Simona Gabrielli, Istituto Nazionale di Geofisica e Vulcanologia (INGV), Italy
Copyright © 2022 Paulatto, Hooft, Chrapkiewicz, Heath, Toomey and Morgan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Michele Paulatto, email@example.com