The plumbing systems of arc volcanoes comprise a vertically extensive network of melt reservoirs (Cashman et al., 2017). A step change from melt-controlled to crystal-controlled rheology naturally divides these reservoirs into (melt-rich) magma chambers and (melt-poor) mush zones (Sparks et al., 2019). It is widely accepted that magma differentiation occurs primarily in the deep crust (Annen et al., 2006), while shallow melt storage conditions control the style and volume of eruptions (Cashman et al., 2013;Popa et al., 2021). The question of whether a shallow reservoir is a magma chamber or a mush zone is critical for assessing eruption likelihood and dynamics. Presence of mobile melt (magma) allows for a rapid eruption (Caricchi et al., 2021;Seropian et al., 2021), while pre-eruptive crystallinity affects buoyancy-driven outgassing (Parmigiani et al., 2017). Just as important for hazard assessment is the reservoir geometry. The spatial distribution of melt can be used to estimate the magma recharge rate and the likely volume of future eruptions.

While petrology gives invaluable insights into thermal evolution and ascent history of magma (e.g., Druitt et al., 2012), constraining reservoir geometry and current storage conditions under active volcanoes requires geophysical imaging. Seismic tomography has provided strong evidence for partial melting beneath many active volcanoes and has helped constrain magma storage depths and volumes (e.g., Heath et al., 2015;Kiser et al., 2018;McVey et al., 2020;Paulatto et al., 2012;Ward et al., 2014). The seismological evidence is broadly consistent with low melt fraction storage (Lees, 2007), but the inherent resolution limits (Malcolm & Trampert, 2011;Paulatto et al., 2022) lead to significant over-and under-estimation of a reservoir's volume and melt fraction, respectively. Although the evidence for crystal mush reservoirs seems compelling, the picture is far from complete. It is unclear whether magma can be rapidly extracted from crystal mush and directly erupted or whether it must first accumulate in high melt fraction magma chambers or lenses (Sparks et al., 2019). The question of how eruptible magma can be extracted from mush reservoirs is also hotly debated (Holness et al., 2019). Persistent melt lenses may have been overlooked if they are smaller than the resolution limits of the geophysical methods used to date. Such "covert" melt lenses have been drilled at ∼2 km b.s.l. at several volcanoes (Rooyakkers et al., 2021). Although their prevalence is unknown, they may pose a significant hazard, affecting the availability of eruptible magma.

Traditional seismic imaging methods are unsuitable for addressing these questions since they lack the resolution required to detect melt lenses. We show that full-waveform inversion (FWI) can overcome this limitation and can provide robust constraints on melt fraction. FWI is well-suited for detecting crustal magma (Morgan et al., 2013), but requires dense sources and/or receivers. Focusing on a submarine volcano has the advantage of dense marine seismic acquisition. Our method can achieve a resolution of a few hundred metres, enough to detect small melt lenses that have so far remained beyond reach (Morgan et al., 2013).

Here, we use full-waveform seismic imaging to study Kolumbo, an active submarine volcano located in the Santorini volcanic field (SVF), along the Hellenic Arc (Figure 1a). Kolumbo is the most active center of seismic (Andinisari et al., 2021;Bohnhoff et al., 2006;Dimitriadis et al., 2009;Schmid et al., 2022) and high-temperature hydrothermal activity within the SVF (Carey et al., 2013;Hübscher et al., 2006;Kilias et al., 2013;Rizzo et al., 2016;Sigurdsson et al., 2006). It last erupted in 1650 AD in a highly explosive, tsunamigenic event (Ulvrova et al., 2016) which produced a few km 3 of hydrous, biotite-rich rhyolitic lava (Cantner et al., 2014;Konstantinou, 2020). Despite the high heat and gas fluxes (Rizzo et al., 2019) which suggest an elevated magmatic activity, little is known about the underlying magmatic system. Local-earthquake tomography showed a 5%-reduction in S-wave velocity at 5-7 km b.s.l. (Dimitriadis et al., 2010). Its interpretation as a magma chamber was utilized in petrological (Cantner et al., 2014) and numerical studies (Konstantinou, 2020). Active-source travel time tomography does not detect a significant anomaly in the uppermost 6 km (McVey et al., 2020), with just a reduction of 3% in P-wave velocity at 2-4 km b.s.l. directly below the edifice (Hooft et al., 2019).

We apply 3D acoustic isotropic FWI (Morgan et al., 2013;Tarantola, 1984) to a wide-angle, multi-azimuth seismic data set (Hooft et al., 2017) consisting of early arriving waves generated by densely spaced airgun shots and recorded by ocean-bottom seismometers (Figure 1b). The analyzed data set is a subset of the data collected in the PROTEUS experiment (Hooft et al., 2017) contained in a 17 × 18 km wide box (Figure 1b), comprising 16 ocean bottom seismometers (OBSs) and 1,505 air-gun shots fired at ∼150 m intervals along the profile lines.

The local Cartesian coordinate system was obtained by rotating a local E-N coordinate system centered at 25.3971°E 36.4042°N by 25.5 counterclockwise. Traces where direct water waves interfere with the crustal refractions (shot-receiver offset less than 5 km) were excluded from the data set, followed by manual removal of bad-quality traces, resulting in an average of 977 traces per OBS. Each trace was 5-s long, sampled at 2.5 ms. The data were filtered with a minimum-phase, frequency-domain, band-pass Butterworth filter of corner frequencies 2-3-4.5-6.5 Hz followed by muting before, and 1 s after, first breaks calculated from the starting model (see below).

We used an explicit, time-domain, least-squares-optimized, finite-difference solver of the acoustic, isotropic, 3D wave equation implemented in the FULLWAVE3D package (Warner et al., 2013) with effective accuracy of sixth order in space and fourth order in time. As the number of air-gun shots was much greater than the number of OBSs, source-receiver reciprocity was used by labeling sources as receivers and receivers as sources.

Our starting P-wave velocity model was derived by first-arrival travel-time tomography (Heath et al., 2019) and discretized at 50-m interval in all dimensions. The density was calculated at each iteration from that model using an empirical relation (Gardner et al., 1974), except for a constant value in the water. The seabed was modeled implicitly as a velocity contrast. The sea surface was simulated as a flat, planar free surface located one node above the top of the model grid. Absorbing boundaries (Cerjan et al., 1985) of 2.5-km width were used elsewhere. A single-source wavelet assigned to all shots was extracted from the data by deghosting stacked, high-quality, near-offset traces using an initial-guess wavelet and a Wiener filter.

A least squares misfit between synthetic and field data, both normalized by their respective norms, was minimized with a steepest-descent algorithm preconditioned by an approximation of the diagonal Hessian of the misfit function (Warner et al., 2013). Both the misfit gradient and the step-length were calculated using the adjoint method (Lions, 1972;Tarantola, 1984). The regularization involved smoothing the misfit gradient by a Gaussian filter over two (local wavelengths) horizontally and one vertically. Both the source wavelet and field data were low-pass filtered at a cut-off frequency increasing by 0.5 Hz from 3 (the lowest inverted frequency) to 6.5 Hz across eight blocks of 20 iterations, totaling 160 iterations.

The gradient-based FWI algorithms based on the least squares misfit, such as the one used here, require a starting model that is good enough to match the majority of observed data to a half of the lowest inverted frequency (Hole et al., 2005). When this condition is not met, the prevailing "cycle-skipped" data may drive the inversion toward a wrong local minimum (Virieux & Operto, 2009). This means that travel-time data inverted for the starting model to be used in FWI needs to be obtained with a well-sampled seismic survey. The data behind our starting model was collected with the same dense geometry that was used for FWI. The high quality of the starting model was confirmed by reproducing most of the observed waveforms in the offset-range of interest down to a half the cycle at the lowest inverted frequency. This frequency was chosen to be 3 Hz based on the trade-off between the poor S/N ratio at lower frequencies and the worse accuracy of the model at higher frequencies (>4 Hz).

The main result of this study is a P-wave velocity model of the subsurface at and around Kolumbo. The model is constrained to a 17 × 18 km horizontal domain (Figure 1b) and stretches down to 4 km below sea level, sampled regularly at a 50-m interval. Synthetic waveforms calculated through this model match the observed data much better than for the starting model and the inversion converged to a waveform-misfit minimum (Figure 2, Figures S6-S17 in Supporting Information S1).

We analyzed phase spectra (wrapped between -π and π radians) of the early synthetic and observed arrivals in order to estimate the S/N ratio of the field data and the adequacy of the input and forward modeling (Shah et al., 2012). Phase-residual plots of early synthetic and observed arrivals were used to assess the quality of the starting model, in particular, to detect problems with data cycle-skipping indicated by abrupt changes from -π to π (or vice versa) between neighboring stations (Figure 2).

To assess the effective resolving power of the inversion, in particular, to quantify the expected underestimation of velocity anomalies, we performed a suite of spike tests (Rawlinson & Spakman, 2016). We used the same strategy and input as for the field-data inversion. The only difference was replacing the observed data with acoustic synthetics calculated through an ensemble of 12 models constructed by perturbing the starting model with a suite of spherically symmetric 3D Gaussian anomalies centered roughly at the observed low-velocity anomaly (13.25; 8.75; 2.75 km), with full widths at half-maximum of 0.2, 0.8, and 1.4 km and peak amplitudes equal to 1%, 5%, 10% and 50% of the background model. Only anomalies stronger than 10% were detectable. The best recovery (80% of the true peak amplitude) was obtained for the 50% anomaly of 0.8-km width (Figure S5 in Supporting Information S1). Anomalies smaller and/or weaker were not recovered. Larger (1.4 km) but weaker (10%) anomalies yielded worse recovery than the best case, while the 50% anomaly of 1.4-km width led to artifacts characteristic of cycle-skipping.

To detect potential outliers in the OBS data and quantify the uncertainty of the final model we used jack-knife resampling (e.g., Rawlinson et al., 2014). We performed 16 field-data inversions, every time excluding a different OBS from the data set. The standard deviation of the models recovered after 160 iterations (Figure S4 in Supporting Information S1) was taken as a proxy for uncertainty.

We used the self-consistent scheme (Berryman, 1980) to calculate the effective bulk and shear moduli of a two-phase, partially molten granite with a composition that is the average of seven rock samples from the most recent Kolumbo eruption (supplement of Klaver et al., 2016). For the solid phase, we used v P = (6.25 ± 0.13) km/s, v S = (3.67 ± 0.12) km/s, and ρ = (2660 ± 20) kg/m 3 from the laboratory measurements at 200 MPa and 30°C of more than 100 granite samples (Christensen & Stanley, 2003). We corrected the velocities (Christensen, 1979) for lithostatic pressure (66 MPa) and a solidus temperature of 700°C (Boettcher & Wyllie, 1968), obtaining v P = 6.20 km/s and v S = 3.58 km/s. Thus, the temperature effect is a second-order factor and partial melting is required to explain the observed velocity anomaly. For the fluid phase, we use ρ = 2,110 kg/m 3 and bulk modulus K = 8.08 GPa from the thermodynamic model based on the Birch-Murnaghan equation (Ueki & Iwamori, 2016).

The assumption underlying the SCS calculation is spheroidal-shaped pores. We use 0.05-0.5 as the range of plausible aspect ratios (Takei, 2002), due to a lack of experimental constraints. We assumed the low-frequency, equilibrated limit (Biot, 1956) starting with air-filled pores followed by fluid substitution according to Gassman's equation (Gassmann, 1951). Sensitivity analysis showed that varying input parameters other than pore aspect ratio and fluid bulk modulus (which in turn is sensitive to water content) has a negligible effect on the final result.

The most striking feature of the final model is a prominent low-velocity zone (LVZ) beneath Kolumbo extending from 2.1 to at least 4 km b.s.l., with a diameter of 0.6 km for the 4 km/s contour and a minimum value of 3.4 km/s at 2.55 km b.s.l. (40% reduction relative to the starting model, see below). To highlight its geometry, we display it as a negative anomaly relative to the starting model, with the absolute values (for both the starting and final model) shown in the insets (Figure 3). Its shape is plume-like in the SW-NE cross-section (Y axis), thinner at the base where the wavelength is larger, and thus not produced by the directional smoothing which we used to regularize the inversion. The regularization leads to recovery of only 80% of the true amplitude as shown by synthetic tests (Figure S5 in Supporting Information S1). The maximum recovered anomaly was 1.6 km/s relative to the starting 5.0 km/s (Figure 3). Correcting for regularization gives approximately a 2.0 km/s anomaly (40% reduction), corresponding to a P-wave velocity of 3.0 km/s within the LVZ.

Several factors make this anomaly a robust feature required by the data. First, the synthetic tests show that a similar anomaly is recoverable by our method (Figure S5 in Supporting Information S1). Second, the anomaly, although much more smeared, was already present in the high-quality starting model (Figure S1 in Supporting Information S1) obtained by travel-time tomography (Heath et al., 2019), which has lower resolution but is more robust against artifacts compared to FWI (Virieux & Operto, 2009). Third, the largest enhancement of the anomaly amplitude, corresponding to the largest improvement of waveform fit, was achieved during the first, low-frequency iterations, indicating that the anomaly is not produced by over-fitting of the high-frequency waveforms. Fourth, the feature is robust against different inversion strategies and data subsets. In particular, jack-knife resampling shows that it is not caused by potential outliers (Figure S4 in Supporting Information S1). Although the recovered final model is likely to have errors due to data noise, inversion regularization, and the simplified physics (acoustic, isotropic, non-attenuative medium) used to calculate synthetic data, we argue that these errors are small and do not affect our conclusions (Text S2 in Supporting Information S1).

The observed low-velocity anomaly can be explained by the presence of either melt or hydrothermal fluids (magmatic gases and/or seawater). We favor the former possibility based on independent seismic reflection observations of the hydrothermal system that place the top of the hydrothermal reservoir at ∼0.6 km b.s.l. (Hübscher et al., 2006). Geochemical analysis of the hydrothermal fluids suggests that interaction of magmatic gases and seawater occurs at ∼1 km b.s.l. (Rizzo et al., 2019). The hydrothermal system is not visible in our model because the shallow structure is poorly constrained by the data (see Figures S4 andS5 in Supporting Information S1). Nonetheless, the weak low-velocity layer at 1.4-1.9 km b.s.l. may be evidence of a volatile-rich lens. Although a multi-level hydrothermal system extending to the base of the observed anomaly is conceivable, the persistent seismic unrest down to 16 km b.s.l. (Schmid et al., 2022) makes this possibility less likely, as explained below.

We can also rule out the hypothesis that the LVZ is caused by a purely thermal anomaly. The temperature required to explain such a strong velocity perturbation greatly exceeds the solidus temperature for granite, implying partial melting. We interpret the LVZ as a partially molten intrusion and we assume it to be granitic (rhyolitic) based on the bulk composition of lava from the 1650 AD eruption (Cantner et al., 2014).

We estimate the maximum melt fraction of the intrusion by comparing the maximum observed P-wave velocity anomaly of 3.0 km/s with values calculated for partially molten granite using the self-consistent effective medium approach (Berryman, 1980), assuming melt-filled spheroidal pores (see Methods). There is a strong trade-off between the pore geometry (aspect ratio) and melt fraction. Here we assume thin, oblate pores (melt films) of aspect ratio between 0.05 and 0.5 (Takei, 2002). Recalling the uncertainty in velocity of 0.3 km/s (from jack-knifing) yields melt fractions ranging from 0.26 to 0.53 (Figure 4). Assuming the velocity error to be normally distributed (with 𝐴𝐴 𝐴𝐴 = 0.1 km/s), the melt fraction is 𝐴𝐴 0.42 ± 0.04 .

The total volume of the low-velocity anomaly equals 6.2 km 3 , which corresponds to a sphere of 1.1 km radius. This is probably a lower bound, as the anomaly seems to extend below the bottom of the model. The total volume of melt stored in the reservoir is ( 𝐴𝐴 1.38 ± 0.32 ) km 3 , where the error bounds come from the aspect-ratio uncertainty. Note that these values are also affected by difficult-to-quantify errors due to the simplified physics used in FWI (see e.g., Marjanović et al., 2019) and should probably be used only in comparison to other acoustic inversions.

A recent seismicity study by Schmid et al. (2022) shows that the observed anomaly coincides with the tip of a cone-shaped region hosting enhanced seismic activity in the 2006-2007 period. The earthquakes were clustered in several swarms starting with the strongest events at 6-9 depth and migrating upwards, terminating at the base of the LVZ. The swarm propagation was faster than expected for hydrothermal gases diffusing from a point source or for active diking (Schmid et al., 2022). However, the speed and focal mechanisms associated with the swarms remain consistent with melt ascent along pre-existing, subvertical pathways through a rheologically strong layer at 4-9 km b.s.l. Schmid et al. (2022) interpret this seismic activity as a result of either transport (Hashin & Shtrikman, 1963). Colored lines correspond to different values of the pore aspect ratio listed in the key, with rhombi marking the critical porosity (Nur et al., 1998). The maximum observed anomaly with three error bounds from jackknifing is denoted in red.

of melt from larger depths to the observed reservoir, or formation of fractures enabling such transport which occurred after the recording period. These findings support our interpretation of a melt-filled reservoir below Kolumbo. Similar seismic activity has been recorded during several other time periods (Andinisari et al., 2021;Bohnhoff et al., 2006;Dimitriadis et al., 2009) suggests ongoing reservoir recharge by frequent injections of magma from greater depth.

The melt fraction of 𝐴𝐴 0.42 ± 0.04 at the center of the velocity anomaly may be overestimated by neglecting the 6% of volatiles reported for the 1650 AD eruptive products (Cantner et al., 2014). However, it is likely above the critical porosity (Figure 4) corresponding to the mush-magma transition (Sparks et al., 2019), implying rheology controlled by melt rather than crystalline framework. The observed >40% reduction of P-wave velocity relative to the regional average is larger than for anomalies obtained by travel time tomography at other volcanoes studied to date. For example, 10% for Newberry (Heath et al., 2015), Yellowstone (Huang et al., 2015) and Mt. St. Helens (Kiser et al., 2018), 17% for Montserrat (Paulatto et al., 2019), and 21% for Santorini (McVey et al., 2020). For this reason, despite the remaining uncertainty in melt fraction, we interpret the imaged reservoir as a small magma chamber rather than a mush zone, or a hybrid reservoir in which a small magma chamber is embedded in a larger mush zone.

Petrological (Cantner et al., 2014) and numerical models (Konstantinou, 2020) of the Kolumbo magmatic system have been based on the interpretation of the low S-wave velocity anomaly from local-earthquake tomography (Dimitriadis et al., 2010), which is located 5-7 km b.s.l., offset from the edifice. This interpretation is inconsistent with the seismicity data (Schmid et al., 2022) which suggests that the source of melt feeding the shallow chamber may be located deeper than 9 km b.s.l., beneath a rheologically strong, seismically active layer (Schmid et al., 2022). The composition of the ascending melt is at least partially mafic, as indicated by the enclaves present in the 1650 AD lavas (Cantner et al., 2014;Klaver et al., 2016). We propose that frequent ascent of mafic melt from deeper levels leads to heating and remobilization (Burgisser & Bergantz, 2011) of the shallow felsic melt reservoir, and consequently to formation of a magma chamber.

We propose that the imaged shallow melt reservoir (Figure 5) represents the uppermost level of a multi-level system. Surprisingly, the imaged melt reservoir is not sill-like as it is elongated vertically, with a broader top. However, this shape may not be well constrained, as indicated by the relatively large uncertainty associated with this part of the model (Figure S4 in Supporting Information S1) and the smearing at the top of the recovered synthetic anomaly (Figure S5 in Supporting Information S1). If we were to treat the mushroom-like shape of the anomaly as a real feature, its broad top may represent a lens of brine (Hill et al., 2020) separated from magma ponding below the high-temperature hydrothermal system (Hübscher et al., 2015;Rizzo et al., 2016;Sigurdsson et al., 2006). In turn, the narrow base may reflect the predominant path of magma supplied through a rheologically strong layer (Schmid et al., 2022) that overlies a yet-to-be-imaged deeper mush zone. This may indicate a more complex multi-level system than numerical simulations which show small (hundreds of meters) magma chambers forming directly at the top of larger (kilometres) mush zones (Booth et al., 2019;Jackson et al., 2018).

The long-term evolution of a magma chamber depends on the time-scales of melt injection, cooling, and viscous relaxation of the surrounding crust (Degruyter & Huber, 2014). For a chamber to grow, mass injection must be faster than cooling, but slower than the viscous relaxation of elastic stress. Assuming that all of the observed melt has accumulated since the last eruption (1650 AD) and taking 2,600 kg/m 3 as the mafic melt density (Konstantinou, 2020), the average melt flux is equal to 4 × 10 3 km 3 yr -1 , or 330 kg/s. This is of the same order of magnitude as the estimate for nearby Santorini (Parks et al., 2012) and Montserrat (Annen et al., 2013). Calculations of eruption frequency of upper crustal silicic magma chambers using the approach of Degruyter and Huber (2014) suggest that this flux can result in tens of eruptions before mechanical lock-up, for chamber volumes of up to 20 km 3 and for fluxes as small as 50 kg/s (Figure 6).

The volume of mobile magma beneath Kolumbo is up to 0.02 km 3 (assuming a pore aspect-ratio of 0.05). This is twice the volume of magma injected in 2011-2012 beneath Santorini (Parks et al., 2012). The total volume of melt accumulated in the reservoir beneath Kolumbo (∼1.4 km 3 ) is close to the ∼2 km 3 of dry-rock equivalent (DRE) volume of the last eruption (Nomikou et al., 2014). For comparison, the DRE volume of the Mt. St. Helens 1980 eruption was 3 km 3 (Kiser et al., 2018). The eruption recurrence interval at Kolumbo is unknown, as only the last of five recognised eruptive phases has been dated with confidence (Nomikou et al., 2016). A sequence stratigraphic analysis indicates that all these phases occurred during the last 1.2 Ma, suggesting an average recurrence interval of ∼240 ka (Preine et al., 2021). A lower bound of 10 ka was proposed based on seismic reflection data and the tectonic history of the region (Hübscher et al., 2015). While this implies a considerably longer recurrence time than the interval since Kolumbo's last eruption, a sizable lens of mobile melt ∼2 km below the seafloor poses a serious threat of a highly explosive, tsunamigenic eruption in the near future. The high melt fraction implies that the reservoir could be readily mobilized by a larger influx of magma or an external trigger (Caricchi et al., 2021). The shallow reservoir depth lowers the critical overpressure needed to form and propagate a dyke to the eruptive vent (Jellinek & DePaolo, 2003). The relatively shallow water (500 m at the crater floor) is likely to enhance rather than suppress explosivity (Dürig et al., 2020). A tsunami and an eruptive column of tens of km with significant ashfall and extensive pumice rafts can be expected as the aftermath of the eruption (Kusky, 2022;Nomikou et al., 2014). Although it will probably be an order of magnitude less powerful than the recent Hunga Tonga-Hunga Ha'apai event, it may have a greater impact, considering CO 2 accumulated at the crater floor (Carey et al., 2013), the dense population on nearby Santorini, and the sea/air traffic in the Aegean.

Our study is the first application of active-source, seismic full-waveform imaging to an active volcano. High-resolution imaging combined with high-sensitivity earthquake monitoring (Schmid et al., 2022) provides unique evidence for magma transport to the shallow crust and a direct image of the resulting magma chamber beneath Kolumbo volcano. The current state of the reservoir indicates that an explosive eruption of high societal impact in the future is possible (though not imminent), thus we suggest establishing a permanent observatory involving continuous earthquake monitoring (Schmid et al., 2022) and seafloor geodesy.

Our image bridges the gap between seismic data and conceptual models of relatively small but high-melt fraction silicic bodies forming at the uppermost level of transcrustal mush systems. Remarkably, the chamber was almost entirely missed by dense travel-time tomography (Heath et al., 2019;McVey et al., 2020), which used the same (Hübscher et al., 2015) and geochemical data (Rizzo et al., 2019). The approximate earthquake locations are based on Schmid et al. (2022). The velocity anomaly is extracted from the final model along the DD' profile (Figure 3b). No vertical exaggeration is applied. active-source data set. Similar melt reservoirs would have gone undetected at other active volcanoes, raising the possibility that they may be more common than previously thought. Further thermomechanical modeling should give insights into the formation and longevity of such chambers.

Wider adoption of full-waveform imaging may usher a new era of volcano geophysics, increasing detectability of melt reservoirs by providing detailed images of elastic parameters, as well as more accurate background models for earthquake relocation and 3D seismic migration (Magee et al., 2018). With ever-increasing computing resources, elastic FWI (e.g., Marjanović et al., 2019) should soon be commonplace, allowing to study subaerial volcanoes as well. These targets will require modeling of irregular land surface either by using finite-element-like solvers (which require a challenging meshing step, e.g., Komatitsch & Tromp, 2002) or the augmented finite-difference method (e.g., Chrapkiewicz, 2021). Whether marine or land, elastic FWI will also require detailed a priori knowledge of the shallow subsurface and/ or suppressing some of the elastic effects present in the data (e.g., Agudo et al., 2020).

Increasingly tight constraints on the spatial distribution of melt and volume of eruptible magma in the shallow crust will substantially improve our understanding of magmatic systems and revolutionize volcanic hazard assessment. (Konstantinou, 2020), crust viscosity equal to log 10 = (20.5 ± 1.5) log 10 (Pa•s) (Konstantinou, 2020), and thermal diffusivity equal to (3.0 ± 0.5) 10 6 m 2 /s, where the uncertainties correspond to 3σ of the respective distributions. This thermal diffusivity is relatively high (Whittington et al., 2009) and was motivated by the cold upper crust and efficient heat removal via the hydrothermal system at Kolumbo. the National Science Foundation. The Geophysical Instrument Pool Potsdam provided 60 land seismometers. The Aristotle University of Thessaloniki contributed five land seismometers and the Greek military donated helicopter time for installations on the smaller islands. The data collection was funded by the National Science Foundation Grant OCE-1459794 to the University of Oregon. The data analysis was funded by Leverhulme Trust Grant RPG-2015-363 to Imperial College London. FS received funding from Deutsche Forschungsgemeinschaft, Grant 453685621. MP was funded by a NERC Independent Research Fellowship (NE/R015708/1).