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 midocean 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-kmwide), 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) lowercrustal 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 traveltimes 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.
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, V P and the S-wave velocity, V S ) by solving a geophysical inverse problem using recordings of seismic waves propagating in the subsurface. Excellent reviews of the theory and methods behind it can be found e.g., in Iyer and Hirahara (1993) and Nolet (2008), while the methodological limitations are highlighted by Foulger et al. (2013). Earlier reviews that specifically focus on seismic tomography of magmatic systems can be found in Lees (2007) and Koulakov and Shapiro (2015).
The most common type of seismic tomography is arrivaltime 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 traveltime. 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 realworld 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 illconstrained 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
The development of different solutions to these challenges has led to an array of different seismic tomography algorithms. Many seismic imaging methods Values in bold were not reported directly in the cited manuscripts and were instead estimated by the authors of this review based on the published figures. Melt fractions reported in red exceed the critical porosity (assumed to be 35%). The column "Melt mtd" lists the methods used to estimate melt fraction from seismic velocities. Exp, experimental relationship (Caricchi et al., 2008); DA, Dehhidral angles (Takei, 2002); SCS, Self-consistent scheme (Berryman, 1980;Schmeling, 1985;Paulatto et al., 2012); HS, Hashin-Shtrikman bounds (Hashin and Shtrikman, 1963); CP, Critical porosity (Chu et al., 2010); ID, Iyer and Hirahara (1993); SCS/DEM, Taylor and Singh, 2002;OB77, O'Connell and Budiansky (1977); M80, Mavko (1980); HH, Hammond and Humphreys (2000); BS79, Bottinga and Steinmetz (1979); S97, Schmitz et al. (1997). The complete list of volcano imaging studies examined in this paper is available in the companion database (https://github.com/VolcanoTomo).
Frontiers in Earth Science frontiersin.org 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 traveltimes 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.
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 V P models, particularly for submarine volcanoes or island volcanoes where dense marine seismic acquisition can be employed; the method can also constrain compressional wave attenuation. Offshore experiments have the advantage of being able to use airguns, with many closely spaced and nearly identical sources, but they have difficulties in constraining V S since sources in the water column do not directly generate S-waves. ASTT has been applied to imaging Montserrat (Shalev et al., 2010;Paulatto et al., 2012), Santorini (McVey et al., 2020;Heath et al., 2021), Mount St. Helens (Kiser et al., 2018) and many other volcanoes (Table 1). Popular algorithms for ASTT include the FAST code by Zelt and Barton (1998); Jive3D by Hobro et al. (2003), TomoLab by Toomey et al. (1994), Tomo2D by Korenaga et al. (2000), Tomo3D by Meléndez et al. (2015), and the code of Van Avendonk et al. (2001).
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 V P and V S or V P /V S , and it can provide deeper constraints if deep earthquakes are present. LET is the most widely adopted volcano imaging technique and has been applied to Aso caldera (Sudo and Kong, 2001), Merapi (Luehr et al., 2013), Askja (Greenfield et al., 2016) and Sinabung (Indrastuti et al., 2019) among many others. Popular LET codes include the SimulPS family of algorithms (Thurber, 1983) and the LOTHOS code by Koulakov (2009).
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).
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 V S but is slower than both V P and V S . Lower frequency surface waves sample deeper, therefore deeper targets require lower frequency data. V S near the Earth's surface typically increases with depth therefore low-frequency surface waves, which are sensitive to deeper structure, propagate faster than high-frequency surface waves. This principle is at the core of SWT. SWT starts with determining surface wave dispersion curves (the frequency dependence of the surface wave phase velocity). These are used to construct horizontal 2D phase velocity maps at different frequencies. The phase velocity maps are then inverted to obtain a 3D model of V S .
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) andCampi 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).
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.
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 andEvans, 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).
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 forwardinverse 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 V P by measuring P-wave travel-times using rays travelling only in the z-direction, the resulting model will be smeared, and the resolution matrix will represent that smearing or vertical averaging of information.
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 nonlinearity 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 Frontiers in Earth Science frontiersin.org 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 underrecovery 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 nonrandom 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.
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 (J) Same as (I) at point P2. Wavefields were calculated using an acoustic finite difference method (Warner et al., 2013).
Frontiers in Earth Science frontiersin.org 11Fresnel 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 r Lλ √ , where L is the length of the propagation path and λ is the wavelength (Williamson, 1991). For a typical crustal-scale LET study like the one by Greenfield et al. (2016), where the path length is 10 km-20 km and the dominant frequency is 5 Hz-10 Hz, the Fresnel radius is estimated to be 2.5 km-5 km. This range matches the resolution length determined using checkerboard tests.
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 V S structure (e.g., Godey et al., 2003). SWT is often carried out using relatively low frequencies and therefore typically has relatively low resolution. For example, Huang et al. (2018) use Rayleigh wave phase velocity measurements at a maximum frequency of 1 Hz and 5-50 km long ray paths, leading to a resolution length of 3.5 km-11 km.
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).
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 fullwaveform 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.
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.
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.
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 (K), which is defined as the pressure increase resulting from a decrease in volume (compression). A material with a high bulk modulus will strongly resist compression (it has a low compressibility). Similarly, the shear modulus (G) describes the material's response to shear deformation. More generally, the response of a material to an external force is described by an elastic 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 λ and μ or the bulk and shear modulus. For a detailed description of the stiffness tensor see e.g. Mavko et al. (2020).
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, Frontiers in Earth Science frontiersin.org P-waves and S-waves, and these are characterized by the P-wave velocity (V P ) and the S-wave velocity (V S ):
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.
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).
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.
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 selfconsistent 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).
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 α of the spheroids (Figures 5B,E). It is impossible to recommend one simplified geometry that is optimal for all applications. In real rock samples, the geometry of the melt is related to the melt fraction, the composition of the melt and crystals, and the temperature of the reservoir. Garapić et al. (2013) used nanoscale imaging to show that low percentages of melt in the upper mantle are best described by a combination of cuspate tubes and low aspect ratio films (α 10 -3 -10 -4 ). Takei (1998) derived a set of theoretical constitutive relationships tying melt fraction and the dihedral angle to the elastic properties assuming thermodynamic equilibrium between the melt and crystal grains. These relationships were then compared to the SCS approach to derive equivalent aspect ratios for a wide range of dihedral angles (Takei, 2002). These results provide a rationale for selecting a particular aspect ratio when the approximate range of dihedral angles is known. Typical dihedral angle values of 20 °to 40 °correspond to equivalent aspect ratios of 0.1-0.15. The link to dihedral angles is extremely useful to narrow down the range of plausible aspect ratios. However, it relies on two important assumptions that are rarely valid in crustal systems. Firstly, the melt is rarely in thermal equilibrium with the surrounding crystals; mush reservoir models predict that the melt and crystals are efficiently segregated by porous flow (Jackson et al., 2018). Secondly, a significant proportion of melt is likely stored in mesoscopic pockets, lenses, or channels that have no relation to the grain-scale equilibrium geometry. For these reasons, even when the microgeometry of the melt is well constrained, it is necessary to consider a range of aspect ratios.
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 ϕ c , the medium behaves like a solid, with stresses supported by the solid matrix. Above the critical porosity, the solid matrix loses its cohesion and the medium behaves like a fluid suspension (Costa et al., 2009). The elastic properties above the critical porosity can be described by the Reuss bound (fluid suspension), while below the critical porosity they can be approximated by a modified Voigt bound (Figures 5C,F) (Nur et al., 1998). The advantage of this approach is that it does not require knowledge of the microgeometry of the melt and instead relies on estimating the critical porosity, which can be achieved by fitting experimental data (Chu et al., 2010). This approach results in constitutive curves that contain a kink at the critical porosity and are similar to those obtained with the SCS approach (Figure 5). By using a modified Voigt bound there is a risk of overestimating the seismic velocities at low melt fractions compared to the SCS (Figure 5F), so this approach is best suited to systems where the melt fraction is expected to be higher than about 10%.
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. (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 midocean ridge magmatic systems in our database; for context we briefly summarize the mid-ocean ridge results in Introduction.
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). Frontiers in Earth Science frontiersin.org
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 km 2 (Heath et al., 2019). Hooft et al. (2019) used ASTT to generate a model of the shallow subsurface, revealing evidence for a strong low-V P anomaly at 2 km depth. They interpreted this anomaly as fractured caldera fill saturated with hydrothermal fluids. McVey et al. (2020) used the same ASTT approach to extend the V P model to greater depth and found evidence for a deeper lowvelocity volume at 2.8 km-5 km depth. The anomaly can be explained by a 4%-13% melt fraction, but higher melt fractions cannot be ruled out because of the effect of wavefront healing.
Kolumbo is a small submarine volcanic cone, a part of the Santorini volcanic system. Dimitriadis et al. ( 2010) found evidence for a -3% V P anomaly and a -5% V S anomaly at 5 km-9 km depth using LET. McVey et al. (2020) confirmed the presence of a small V P anomaly at 3 km-4 km depth using ASTT. They estimated the melt fraction required to explain the data as <1%. The PROTEUS dataset was the first volcano study specifically designed for active source FWI and Chrapkiewicz et al. (2022) used this method to process the same dataset as McVey et al., 2020. The resulting high-resolution model revealed the presence of a strong but small low V P volume at 2.2 km-4.0 km depth beneath Kolumbo (Figure 8A). The value of V P at the centre of the -40% anomaly is 3 km/s. The inferred melt fraction is 44 ± 4% based on detailed calculations using the SCS approach. This result is particularly significant for two reasons: i) it is the clearest and most robust image of a highmelt fraction reservoir or melt lens obtained with seismic tomography, ii) it starkly highlights the limits of travel-time tomography. These results are best explained by a multi-level melt reservoir topped by a melt lens.
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 V P model based on travel-times from offshore airgun shots to land stations. The model extends to 6 km depth and revealed the shallow structure of the island. SHV and the two neighbouring and extinct volcanoes Centre Hills and Silver Hills all have high V P cores in the top 4 km. High-velocity cores are observed at many volcanoes worldwide (e.g., Lees, 1992) and are attributed to the presence of intrusive igneous rocks contrasting with the surrounding volcaniclastic apron. Paulatto et al. (2012) complemented the land station data of Shalev et al. (2010) with OBS data. Thanks to the longer offset OBS data that undershot the island, their model was able to constrain the deeper structure at up to 8 km depth. This deeper model revealed a large low-V P volume beneath SHV. The anomaly has a volume of at least 20 km 3 and a maximum anomaly amplitude of -12%. They used the SCS method to estimate a melt fraction of 3%-10%, but acknowledged that this represents a minimum melt fraction, due to wavefront healing and smoothing. They attempted to better constrain the melt fraction by integrating the tomography results with thermal modelling and concluded that the observed seismic data could be explained by a magma reservoir with a 35% melt fraction and a total volume of 13 km 3 . Annen et al. (2014) further expanded the thermal modelling to constrain the longterm magmatic flux. Paulatto et al. (2019) reprocessed the seismic dataset and jointly inverted the seismic travel-times together with land and marine gravity data (Figure 8B). They used V P /density cross-plots and rock physics calculations to re-estimate the melt fraction. They conclude that the melt fraction in the shallow reservoir is at least 3% but could be 28% or higher depending on the severity of the wavefront healing effect. Their best estimate for the melt fraction is 17%.
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-V P volumes at 2 km-3.5 km and 5.5 km-8 km depth. They estimate the melt fraction in the shallow melt reservoir to be at least 5%. In 2014 Mount St. Helens was the subject of the most complex field campaign ever conducted on a continental arc volcano (the iMush experiment). The deployment involved a total of about 6000 recording stations and the firing of 23 explosive shots. Several studies were published using the data from this experiment. Kiser et al. (2016) processed the active source data along two perpendicular 2D profiles and found generally higher-thanaverage V P beneath the volcano. They also found high V P /V S at 4 km-13 km depth and interpreted it to correspond to a magma reservoir. Kiser et al. (2018) used 3D ASTT and confirmed the presence of a low-V P anomaly at 3.4 km-14 km depth. The low-V P volume shows a complex structure with multiple lobes and layers and probably encompasses the two reservoirs postulated by Lees (1992) and Waite and Moran (2009). The melt fraction is estimated to be 10%-12% using the critical porosity method of Chu et al. (2010). Flinders and Shen (2017) imaged the crust at Mount St. Helens as part of a larger-scale surface wave waveform inversion. They confirmed the presence of a low V S anomaly of at least -15% corresponding to ~5% melt fraction. They also found low-velocity volumes beneath the neighbouring Mount Rainer and Mount Adams. The most complete model of the Mount St. Helens magmatic system comes from joint ASTT and LET by Ulberg et al. (2020), providing constraints on both V P and V S (Figure 8C). They interpret the observed low-V P /low-V S as due to 3% melt fraction. An independent constraint on melt fraction was obtained by Bedrosian et al. (2018) using MT inversion. The overall picture at Mount St. Helens is of a complex multi-level magmatic system dominated by low melt fractions.
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-V P anomaly of -10% at depth of 2 km-5 km but could not confidently attribute it to partial melt. At the same time, Stauber et al. (1988) used TST to image the same target to a greater depth. Their model extends deeper and resolves a high-velocity core beneath the volcano but does not have the resolution to detect any robust low-velocity anomalies. The seismic attenuation study of Zucca and Evans (1992) seemed to contradict the findings of Achauer et al. (1988) and concluded that no magma was present. More recently, Beachly et al. (2012) combined the 1980s data with a new active source 2D dataset and inverted them together using a more modern tomography code (Toomey et al., 1994) to obtain a better resolved 3D V P model. Their tomography model finds a weak low-V P anomaly, in agreement with Achauer et al. (1988). More importantly, they conducted a detailed analysis of wavefront healing, using two-dimensional finite-difference wave propagation modelling. Their work showed that the travel-time results and even waveform modelling results are non-unique, and that the data could be explained equally well by a low melt fraction mush or a melt sill, with possible melt fractions ranging from 11% to 100%. Heath et al. (2015) combined the active source data from Beachly et al. ( 2012) with teleseismic earthquake data. Their V P model once again confirms the presence of a low velocity volume at 3 km-5 km depth, requiring a minimum of about 10% partial melt. A sharp boundary at the top of this magmatic system was mapped with a coherent P-wave reflections at ~2.5 km depth using auto-correlations of the ambient and seismic wavefields (Heath et al., 2018).
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-V P /V S ratio volume, interpreted as a volatile-rich fractured rock (e.g Vanorio et al., 2005). These studies did not find evidence for significant amounts of magma in the top 5 km of the crust. The breakthrough in terms of magma imaging came from Zollo et al. (2008), who used Amplitude vs. Offset analysis (AVO), a method commonly used in the oil industry to determine hydrocarbon saturation. They identified seismic reflections in marine active source data collected in the Bay of Naples. The reflections are generated by a discontinuity at 7.5 km depth, and the AVO analysis indicates that the underlying layer has a V P of about 2,800 m/s and a V P /V S of 2.35. They then used the Hashin-Shtrikman bounds to calculate a melt fraction of 80-90%. This low-V P layer might represent a large tabular magma body at the top of the magma reservoir, perhaps similar to the high melt fraction layers predicted by Jackson et al. (2018) based on modelling of reactive melt flow in mush reservoirs.
The Toba caldera (Sumatra, Indonesia) was first imaged by Mastruyono et al. (2001) using LET. They imaged a low V P anomaly in the upper crust with a magnitude of up to -37%. They concluded that melt must be present but did not attempt to quantify the melt fraction. The TRFI study of Sakaguchi et al. (2006) supported this conclusion and placed the magma reservoir at 8 km-14 km depth. They also found increased V P /V S but again did not attempt to estimate melt fraction. Koulakov et al. (2009) generated a regional V P and V S model using LET and confirmed the presence of an extensive region of low velocities beneath Toba. Stankiewicz et al. (2010) and Jaxybulatov et al. (2014) came to the same conclusion using SWT. Jaxybulatov et al. (2014) found evidence for strong positive V S radial anisotropy in the lower crust, indicative of horizontal layering, similar to observations at Yellowstone and Long Valley caldera (Jiang et al., 2018). Koulakov et al. (2016) revisited LET at Toba and found evidence for a multi-level magmatic system comprising at least two separate large reservoirs located in the crust and in the upper mantle.
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% V P anomaly and a -31% V S anomaly at 6 km depth. They attempted to estimate melt fraction using the SCS method of Mavko (1980) and inferred a minimum melt fraction of 10%. Abe et al. (2017) imaged the deeper structure of the system using TRFI and found evidence for two connected low-V S volumes, one at 8 km-15 km depth and one at 15 km-23 km depth. They estimate melt fraction to be about 10%-14%. An ANSWT study by Huang et al. (2018) found more evidence for a low-velocity volume at about 6 km depth and also detected an additional low-velocity volume at 2 km depth. The evidence at Aso is of a large transcrustal magmatic system comprising multiple levels of magma storage.
Uturuncu (Bolivia) is the location of some of the most spectacular volcano tomography images worldwide. An early TSRFI study found evidence for a low-V S (<0.5 km/s) layer at ~19 km depth (Chmielowski et al., 1999), which they interpret as a sill-like magma body, and they named it the Altiplano-Puna magma body (APMB). Ward et al. (2013) carried out ambient noise SWT and confirmed the presence of a large low-V S layer. Ward et al. (2014) combined SWT with receiver functions to better define the extent and properties of the APMB. They concluded that the APMB extends from 4 to 25 km depth and has a radius of ~100 km, corresponding to a colossal volume of ~500.000 km 3 . They use a critical porosity approach (Chu et al., 2010) to estimate the melt fraction to be 4%-25% and interpret it as a region of layered crystal mush. Independent evidence from inversion of magnetotelluric (MT) data (Comeau et al., 2016) puts the depth of the APMB at 13 km-17 km and indicates that the melt fraction is above 15%. A shallow magma reservoir may also be present. Kukarina et al. (2017) found high V P /V S values at 6 km-80 km depth using LET. The evidence at Uturuncu clearly points to the presence of a vast magmatic system, likely dominated by low-melt fractions, but sizable high-meltfraction bodies may also be present within the mush region.
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. V S at the centre of the anomaly is <2.3 km/ s, a velocity decrease of at least 12%. This low-V S anomaly and the region below it are characterized by negative radial anisotropy, suggesting vertically aligned heterogeneity likely related to vertical faults or dikes associated with the extensional tectonic environment. This result is unusual compared to other volcanically active regions discussed in this paper (e.g., Yellowstone and Toba), where seismic studies have found positive radial anisotropy caused by horizontal layering. Negative radial anisotropy, however, has been observed in other rift settings and is thought to be associated with vertically aligned melt or cracks (e.g., Bastow et al., 2010;Hammond et al., 2014;Hammond and Kendall, 2016;Chambers et al., 2021).
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, smallamplitude low-velocity anomalies. Dawson et al. (1990), Steck andProthero (1994), andWeiland et al. (1995) used TST and found low V P anomalies at about 7 km-20 km depth with magnitude as high as -30%. Melt estimates range from 8% to 25%. Almost every conceivable imaging method has been applied to Long Valley since those early experiments. Perhaps the most successful results came from combined LET and TRFI by Seccia et al. (2011). Their LET study found a weak low-V P anomaly at 6 km-7 km depth, indicating 2%-3% melt. More significantly, their receiver function analysis suggests the presence of a much stronger low-V S anomaly at 7 km-11 km depth, requiring melt fractions of 30-60%. Further evidence came from a regional surface wave waveform inversion by Flinders et al. (2018). Their V S model includes a low-V S anomaly of at least -20% at 5 km-20 km depth, requiring a melt fraction of 23% ± 4%. The evidence at Long Valley is suggestive of a large magmatic system including both regions of crystal mush and some highmelt fraction layers.
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 V P and found evidence for a shallow (1 km-2 km depth) low-V P volume with V P of 3 km/s. They estimate a melt fraction of 50% using the lower Hashin-Shtrikman bound and a melt volume of 5 km 3 . This result should be taken with a grain of salt because the forward modelling approach used by Gudmundsson et al. (1994) can be highly non-unique. A more recent experiment using LET did not have good enough resolution to confirm these results (Jonsdottir et al., 2007). Jeddi et al. (2016) used a denser earthquake dataset and combined it with the active source data from Gudmundsson et al. (1994) to perform joint ASTT/ LET. They found velocities of 3.5 km/s-4.0 km/s beneath the caldera, suggesting a smaller but still significant amount of melt. A SWT study (Jeddi et al., 2017) could not find evidence of melt.
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 V P /V S in the crust indicating the presence of melt. The model was extended and refined by Hammond et al. (2020) who carried out TRFI and detected a low V S layer at 4 km-8 km depth with a V P /V S as high as 2.1. Independent evidence from 3D MT inversion revealed a multi-level system with 7%-30% melt (Yang et al., 2021).
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 V P as low as 4.9 km/s, a 21% reduction compared to a solidified granite. Several additional LET studies were carried out in the following decades including notably by Kissling (1988), Miller and Smith (1999), Husen et al. (2004), Farrell et al. (2014), and Huang et al. (2015). They all found evidence for low-V P anomalies of various magnitudes and volumes at depths of 5 km-15 km. Huang et al. (2015) combined TST and LET to extend the velocity model to the upper mantle, providing a complete view of the magmatic system. The emerging picture at Yellowstone is of a system comprising two large reservoirs: a rhyolitic reservoir in the upper crust and a basaltic reservoir in the lower crust. Stachnik et al. (2008) andSchutt et al. (2008) used SWT to image the system down to depths of 150 km and came to similar conclusions. In addition, they imaged low velocities in the upper mantle attributed to the mantle plume. Chu et al. (2010) used a joint inversion of surface waves and teleseismic receiver functions to better constrain the melt fractions in the crust. They concluded that the upper crustal reservoir has a melt fraction of about 32%. Jiang et al. (2018) analyzed vertically polarized and horizontally polarized surface waves to determine seismic anisotropy beneath the caldera. They found positive radial anisotropy at 5 km-18 km depth indicative of horizontal layering as would be expected in the presence of a large sill complex or a layered mush reservoir.
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 V P model is highly heterogeneous but no clear evidence for magma storage was found. Tenerife's complex structure was also imaged with attenuation tomography by Prudencio et al. (2015). Several high-attenuation bodies stand out in their model, but no clear magma body can be identified. El Hierro was imaged using LET by García-Yeguas et al. (2014) and later by Martí et al. (2017). Both studies found evidence for a low velocity volume at the base of the crust or in the upper mantle at depth of about 12 km-22 km, but neither attempted to estimate melt fraction. More recently, D'Auria et al. ( 2022) used LET to image the source of the recent eruption on La Palma. They imaged a low-velocity volume at 7 km-25 km depth an interpreted it as the active magma reservoir.
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. (2007Park et al. ( , 2009) ) found evidence for low velocities in the upper mantle beneath Lohi Seamount. Lin et al. (2014) used LET to generate joint V P and V P /V S models. They discuss several low and high velocity anomalies but do not find evidence for large magma accumulations. Kilauea has erupted multiple times in recent decades; therefore, the lack of significant magma bodies is surprising. It is likely that shallow magma reservoirs exist but are too small to be detected with the methods used so far. Larger magma bodies may exist in the deeper crust or upper mantle.
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.
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 lowvelocity 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% V P anomaly and a -55% V P anomaly can produce the same observed traveltime delay and are therefore indistinguishable with travel-time tomography. A regularized least-squares travel-time inversion that has been carefully tuned aims for the simplest model and thus would recover the smallest anomaly that is compatible with the data. Consequently, any derived melt fraction should be interpreted as a lower bound. Therefore, we can expect a skew towards low melt fractions. This conclusion is further complicated when the uncertainties in melt fraction calculations are considered. It follows that reported melt estimates are often highly inaccurate.
Only five studies reported melt fractions above 35%. These are listed below, together with the estimated melt fractions and imaging methods. 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 lowvelocity 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 km 3 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 km 3 ) and medium-sized (<100 km 3 ) 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 km 3 ) 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):
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.
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 highresolution 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).
Frontiers in Earth Science frontiersin.org
Paulatto et al. 10.3389/feart.2022.970131
Paulatto et al.
The prevalence of volatile-rich and high melt fraction lenses(Rooyakkers et al., 2021)
The geometry of eruptive conduits(Cassidy et al., 2018)
The melt fraction, internal structure, and eruptive potential of magma reservoirs(Rosi et al., 2022)
The melt fraction of deep crustal magma storage regions and the nature of deep melt ascent pathways(Giordano and Caricchi, 2022)