Numerically Calculated 3 D Space-Weighting Functions to Image Crustal Volcanic Structures Using Diffuse Coda Waves

Seismic coda measurements retrieve parameters linked to the physical characteristics of rock volumes illuminated by high frequency scattered waves. Space weighting functions (SWF) and kernels are different tools that model the spatial sensitivity of coda envelopes to scattering and absorption anomalies in these rock matrices, allowing coda-wave attenuation (Qcoda) imaging. This note clarifies the difference between SWF and sensitivity kernels developed for coda wave imaging. It extends the SWF previously developed in 2D to the third dimension by using radiative transfer and the diffusion equation, based on the assumption that variations of Qcoda depend solely on variations of the extinction length. When applied to active data (Deception Island, Antarctica), 3D SWF images strongly resemble 2D images, making this 3D extension redundant. On the other hand, diffusion does not efficiently model coda waveforms when using earthquake datasets spanning depths between 0 and 20 km, such as at Mount St. Helens volcano. In this setting, scattering attenuation and absorption suffer tradeoffs and cannot be separated by fitting a single seismogram energy envelope for SWF imaging. We propose that an approximate analytical 3D SWF, similar in shape to the common coda kernels used in literature, can still be used in a space weighted back-projection approach. While Qcoda is not a physical parameter of the propagation medium, its spatially-dependent modeling allows improved reconstruction of crustal-scale tectonic and geological features. It is even more efficient as a velocity independent imaging tool for magma and fluid storage when applied to deep volcanism.


Introduction
Seismic attenuation imaging performed using coda waves provides novel information about tectonic structures and fluid content at crustal [1,2], regional [3], and local [4] scales.The attenuation coefficient is proportional to the sum of the inverse intrinsic (Q −1 i ) and scattering (Q −1 s ) quality factors.A separate estimate of scattering attenuation and absorption is crucial for understanding seismic wave propagation in highly-heterogeneous volcanic environments [5], or when targeting areas having different tectonic and scattering properties at crustal and lithospheric scales [6][7][8].A scattering ellipsoid has been adopted for decades by scientists to map the sensitivity of coda waves to Earth's heterogeneities, and map scattering attenuation and absorption in space [9].More recently, 2D and 3D coda sensitivity kernels based on multiple scattering propagation have been proposed to separate Q −1 i and Q −1 s [8,10] and to invert for attenuation in the subsurface at different scales and depths [2,11,12].These sensitivity kernels define the source parameters observed at a station as a space-weighted average of attenuation characteristics of the sampled medium, where the weights are defined via integral equations [10,12].Their application has led to absorption mapping at lithospheric scales [2] and are considered important for the evaluation of the effective sensitivity in ambient noise imaging [12].
The space-weighting functions (SWF) discussed in this note are designed to be applied in the back-projection (or regionalization) method to retrieve the attenuation parameters in space [9,13].In this case, Q −1 i and Q −1 s estimated for a single source-receiver couple characterize the entire space volume, weighted by SWF values between zero and one.The SWF are designed with a Monte Carlo simulation of the multiple scattering process, following the method of Yoshimoto [14].Each SWF value, associated with a point in space for a single station observation, is proportional to the probability that at that point the attenuation value is equal to the single sstation observation.At a point in space, we thus have as many probabilities as observations.The average of all the observed values weighted by the SWFs provides the value of the attenuation at the point.These SWFs have been expressively designed to map scattering attenuation and absorption in volcanoes using a diffusion model and active sources [4,15,16].In the resulting models, the high attenuation contrasts are often related to magma/fluid storage under volcanoes and ongoing volcano dynamics.
For a full discussion of the practice of attenuation mapping by weighted back-projection in volcanoes, the reader is refered to Del Pezzo et al. [17].These authors obtain SWF for Q −1 i and Q −1 s .The two parameters can be rewritten using the associated parameters; either the seismic albedo (B 0 ), the inverse extinction length (Le −1 ), or the intrinsic (η i ) and scattering (η s ) attenuation coefficients (Throughout this paper the syntactic rules used by the Wolfram Mathematica software for the use of parentheses is used: square brackets indicate the argument of a function; curly brackets indicate the elements of a matrix; round brackets indicate an algebraic grouping.): With an SWF, the spatial Q −1 i and Q −1 s are obtained using the following equations: where K 2D i and K 2D s are the intrinsic and scattering SWF, Q −1 ik and Q −1 sk represent the estimates calculated from the fit of experimental energy envelopes with the diffusion model, and k spans the available energy envelopes.The uncertainties of the estimates of Q −1 ik and Q −1 sk can be propagated in Equations ( 2) and (3) to estimate the variances of with the assumption of small covariance and null uncertainty in the determination of the weighting functions.
Del Pezzo et al. [17] additionally obtained that, in the case of a uniform half space and for diffusive propagation, the following function well approximates the numerically calculated SWF for both absorption and scattering attenuation: In Equation ( 4), D is the source receiver distance, x and y are the spatial coordinates, x s and y s the source coordinates, and x r and y r the receiver coordinates.The function fits the numerically calculated SWF reasonably well in the case of a short lapse time (around 15 s), a highly diffusive media, and δ x = δ y = 0.2.These parameters represent the spatial aperture of the weighting function.The two numerically evaluated SWF have approximately the same shape once the level of heterogeneity increases (i.e., when the scattering processes approach the diffusion regime).This is contrary to what happens for lower heterogeneity [18,19] and as a result is valid only for volcanoes and the active data geometry.
The spatial patterns described by the SWF depict the contribution of each cell to the coda formation and is thus proportional to the sensitivity kernels, for scattering and intrinsic dissipation.Equation ( 4) is indeed equal to that proposed at crustal scales for absorption mapping, however only at late lapse times [2,10,11].The sensitivity is maximum at the source and receiver stations, remains high across the area contouring the seismic ray, then decreases at a distance controlled by the extinction length.This similarity in shape goes even further, as the spatial pattern of the function is identical to the depth-dependent diffusive sensitivity kernels in 3D defined by Obermann et al. [12].The difference is in that the kernels do not assume a depth-dependent velocity structure, an approximation that is unfulfilled for shallow volcanic sources, but rather a constant velocity in half-space approximation.The analytical solution of Equation ( 4) is thus an approximate analytical equation for mapping Q coda , similar in shape and meaning to those developed to map absorption.The equation was re-framed as a forward problem in a 5 km deep volcanic medium [20] to map coda attenuation at the Campi Flegrei caldera.The results of the inversion show the increased illumination provided by this technique and additionally demonstrated important correlations of the coda attenuation anomalies with deformation sources at the volcano.
The present note investigates how effective the SWF are in illuminating multi-scale volcanism in 3D.It is divided into three parts: 1. Equation ( 4) is extended to the third dimension, maintaining the assumptions of shallow source and receiver in a diffusive Earth medium with no depth dependency-this is the case for the analysis of active seismic shots in volcanoes.2. We propose and discuss an SWF for mapping Q coda , calculated for a deep source in a non-diffusive medium and discuss its limits.3. We check the reliability and limits of the new approaches by applying 3D SWFs to published seismic databases.We use precalculated attenuation measurements for single source station paths from active data recorded at the Deception Island volcano (Antarctica) [21] and volcano tectonic earthquakes at the Mount St. Helens volcano (USA) [22].
In the Appendix, we report the main tests which were carried out in developing the applications.Test images are compared with previous tomography results obtained in the same areas using different seismic attributes, and show consistent features.

Diffusive Earth Media
We extended the numerical simulations described above to the third (depth) dimension, introducing the z-axis and keeping the half space approximation.For the assumption of no anomalous relevant depth dependency, we rely on the results of Reference [1].The weighting function remains symmetrical around the axis connecting the source to the receiver, in analogy with simulations using radiative transfer theory [10] and spectral elements methods [12].This symmetry allows one to evaluate the 3D SWF analytically for source (a shot) and receiver both placed at the surface.In the case of a uniform half space, the function: approximates the numerically calculated SWF in 3D to the first order (Figure 1).This analytical approximation is valid for the same range of Q i and Q s values and the same lapse time (15 s) used in Del Pezzo et al. [17].This approximated space weighting function is actually a "kernel" function.Different from the other diffusive kernels, it is valid solely for diffusive fields, short seismograms, and surface sources, like those recorded from shots fired in volcanoes for tomography purposes [21].

Deep Sources (Natural Events) and Non-Diffusive Fields
In the case of deep earthquakes, the assumptions made in calculating the approximate SWF given by Equation ( 4) are invalid, and a multiple scattering regime better models coda wave propagation.We thus adopt the Paasschens [23] approximation of the energy transport equation solution in three dimensions to describe the seismogram energy envelope: where: and δ and H are the Dirac delta and the Heaviside step functions, respectively.Here, W 0 is the source energy and v is the seismic velocity.Fitting Equation ( 6) to the experimental energy envelopes, the single path separate estimates of B 0 and Le −1 are possible in principle; however, in 2D, a severe tradeoff affects the two parameters as discussed in Del Pezzo et al. [17].An alternative is the use of a simplified formula, which estimates Le −1 and B 0 by fitting the data to a first order approximation of the energy transport model equation, as given by Zeng et al. [24]: With such a fit, the severe tradeoffs disappear.Equation ( 7) is equivalent to the single scattering model developed by Sato [25] and is valid for low heterogeneity and short lapse times.In this case, intrinsic attenuation controls Le −1 , since η s is small.The physical meanings of the retrieved B 0 and Le −1 become controversial when the energy envelopes are recorded in media with high heterogeneity and modeled with Equation (7).In this case, the fit function is based on improper assumptions and Le −1 is proportional to the widely measured Q coda , the coda quality factor [25] used to map different tectonic settings at the crustal scale [1].The downside is that Q coda is not a physical parameter of the propagation medium; however, the Le −1 (or Q coda ) spatiak distribution can still depict attenuation properties, and the corresponding SWF, K coda , can be calculated.For this task, we use the hypothesis of Pacheco and Snieder [26], setting B 0 to an average value and Le coda .Hence: where is the point with coordinates {x, y, z}, T is the lapse time, τ is the integration variable (time).The integral can be numerically calculated.In Figure A1, we show the contour plot of Q coda as a function of Q −1 i and Q −1 s .For low scattering attenuation, Q −1 coda is independent of Q −1 s and similar to Q −1 i (see Appendix, Figure A1, left panel).An increase of scattering (right panel) increases the tradeoff.In Figure 2, we reproduce the SWF calculated using Equation (8).

Application Examples
The final Q coda image as a function of the spatial coordinates in a 3D space is thus obtained with a back-projection analogous to that used in Equations ( 2) and (3): where K k coda is the weighting function for the k-th source-receiver couple, Q −1 coda,k is the k-th Q coda estimate, and k spans over the available source-receiver couples.To avoid confusion with respect to the definition of source station kernels we remind the reader that: 1.The values of K 3D coda,k [x, y, z] express the probability that the Q −1 coda estimated at a station is equal to the one measured at [x, y, z]. 2. Equation ( 9) is to be used exclusively for back-projection.

The kernel K 3D
num [x, y, z, x r , y r , x s, y s ] in Equation ( 5) can still be used in an inversion for the space-dependent parameters, if the underlying hypotheses are fulfilled.

Deception Island Volcano-Diffusive Approximation
The Deception Island volcano (Antarctica) is an extraordinary natural laboratory, characterized by a horseshoe shape which permits one to design seismic active field surveys characterized by elaborate source and receiver geometries.To test the 3D SWF discussed in this note, we used data from the seismic experiment TOMO-DEC (Tomography at Deception) [21] publicly available from the Australian Antarctic Data Center (AADC) repository.The same data set was used by Prudencio et al. [16], who obtained a 2D attenuation image of this island using a simplified, Gaussian SWF (data and final models also are available from the AADC repository).Del Pezzo et al. [17] improved this image using the 2D weighting function of Equation ( 4), applied to data filtered in several frequency bands centered from 4 to 20 Hz.The present test is carried out using data filtered in the 4 Hz band, where the highest attenuation contrasts were previously observed.Using the 3D SWF of Equation ( 5), we show the attenuation coefficient space distribution calculated at depths of 2 and 4 km, with a horizontal grid of 4 km (Figure 3).The two panels are similar, with high absorption affecting the eastern and southwestern parts of the island.Since the SWF are practically null at 6 km, no images can be calculated below this depth.s (f) models from Del Pezzo et al. [17] are redrawn for comparison using the same color scale.We use the same distribution of sources and receivers as in Prudencio et al. [15].

3D SWF at Mount St. Helens Volcano-Non-Diffusive Media
Mount St. Helens volcano (US) is a central-cone stratovolcano, characterized by 0-7 km deep earthquakes (under the central cone) and lateral fault seismicity (down to 20 km).A 3D Q −1 c attenuation model of the area has been calculated using the SWF described in Equations ( 8) and ( 9) at Mount St. Helens, with a test passive dataset of 451 waveforms ( [27], available from the PANGAEA Data Centre).We use the single-path Q coda estimates obtained by De Siena et al. [3] at 6 Hz.In the Appendix, a sensitivity test, carried out to check the reliability of the method, is described.Equation (9) has been applied to a spatial grid with points separated by a distance of 4 km.In this way, we obtain the Q coda space values at 500 3D grid points.
The Q −1 coda 3D spatial distribution is plot for two horizontal slices, crossing the z-axis at depths of 0.5 km and 4 km (Figure 4, uppermost panels).The vertical slice (lower panel) intersects the surface along the white line drawn in the upper left panel.A sensitivity test using a hemispherical anomaly centered in the middle of the study area is described in the Appendix.The input test is only roughly reproduced: the small number of data available would correspond to an underdetermined inversion problem.This strongly reduces the sensitivity of the method to small anomalies; making any checkerboard test unsuccessful.coda spatial distribution has been interpolated before plotting the average inverse Q coda , < Q −1 coda >.All panels are drawn using Mathematica 10 TM .

Discussion
Figure 1 shows that, in areas of high heterogeneity (diffusion approximation) and for data shots fired at the surface, the sensitivity of the SWF method strongly reduces on increasing depth; since the SWF values strongly decrease with depth.Coda waves recorded from shots fired at the surface in a diffusive Earth media and recorded at short distances, as for the Deception Island case study, propagate mainly in the upper 3-4 km of the crust (Figure 3).The northern part of the island is associated with the crystalline basement and shows low attenuation, while high-attenuation bodies, spatially-correlated to high-velocity structures [28,29], characterize the southern part of the volcano.There is a consistent agreement between low/high coda attenuation and high/low velocity structure from the first scattering/absorption separations [16].The correlation between the SWF-dependent 2D models [17] and the 3D models indicates that coda attenuation estimates are stable using this dataset.Comparing the present 3D attenuation images with the total-Q images obtained by Prudencio et al. [30], which used direct-P coda-normalized waves (MuRAT code-De Siena et al. [22]), we observe a good match between the 3D intrinsic-Q and the total-Q distributions.The location of the main total high attenuation body retrieved by Prudencio et al. [30] spatially fits the main absorption anomaly.
To investigate greater depths, deeper sources (passive data) are necessary.In this case, the diffusion equation is inappropriate, as the Earth's heterogeneity strongly reduces with depth, and a theory based on multiple scattering is necessary.However, inverting a multiple scattering model for the energy envelopes associated with single source-receiver couples prevents the recovery of separate scattering and intrinsic tomography images due to the tradeoff between B 0 and Le −1 .The only possibility is thus the use of an approximate kernel to invert for a unique parameter, Le −1 , a quantity proportional to the widely measured Q coda parameter.In this case, we proposed to calculate the corresponding SWF by using the approach described by Pacheco and Snieder [26].Despite the controversial physical meaning of Q coda , images of the spatial variations of Q coda are still retrievable, like those recently described by Mayor et al. [2], which depict the attenuation structure of the Alps.
Following this approach, we calculated the 3D Q coda image of Mount St. Helens volcano (Figure 4).We compared it to a map of the 2D Q coda spatial distribution obtained by De Siena et al. [3].The authors used maps of late lapse time Q coda -assuming they measured absorption, and energy envelope peak delays, which are proportional to scattering-Q-to separate scattering attenuation from absorption.They back-project the single station Q coda value assuming that it is distributed on a strip connecting the source and ray, derived from pre-calculated 3D rays.
At both depths shown in Figure 4, the low inverse Q coda west of Mount St. Helens is a major feature, similar to that observed by De Siena et al. [3] (see their Figure 5, 6 Hz panel).Nevertheless, this area is a unique anomaly in our analysis, located west of the volcano at a depth of 4 km, and extends to the south at 500 m.A wide area inside this anomaly was not sampled in the previous study, as it assumed a back-projection of the single station Q coda along a strip.In the case of Mount St. Helens, many of the seismic sources are located at, or below, 8 km.The SWF theoretically produce an improved resolution at this depth range due to the wider illumination at near-source nodes.The example reported in the present paper is made with a limited amount of data.The images obtained for Mount St. Helens are thus defocused and need to be improved by using a latger data set.Despite this limitation, the use of SWF is promising in enlightening the space attenuation contrasts.We are confident that it may become a useful tool to complement tomography images achieved with different techniques, especially due to its independence of velocity tomography results.
Author Contributions: Edoardo Del Pezzo wrote the paper together with Francesca Bianco, Jesus Ibanez, and Luca De Siena.Edoardo del Pezzo, Simona Gabrielli, and Luca De Siena analyzed the data and calculated the 3D SWF derived models for all the examples reported.Angel de La Torre wrote the software reported in the Appendix implementing the SWF.

Appendix C. Numerical Integration of Equation (8)
The function to be integrated (Equation ( 8)) is a product of two functions, each one including a delta and a continuously decaying term, which here we call "coda".Hereafter we drop out the dependence of Equation ( 8) on B 0 , Le −1 and v leaving unaltered and t.Therefore, the integral K ss [ , t] can be decomposed into four integrals (i.e., delta•delta, delta•coda, coda•delta and coda•coda): each of them null for t < (t a + t b ) where t a and t b are respectively the time the perturbation reaches position from the source and the time from to the receiver.These integrals are defined as: where E 1 [r, t] refers to the wavefront (or delta) contribution, E 2 [r, t] refers to the coda contribution, [r a , t a ] refers to the source-impulsive response, and [r b , t b ] refers to the receiver-impulsive response.
Taking into account the sampling property of the Dirac delta function:

Figure 1 .
Figure 1.Plot of the 3D kernel function obtained using Equation (5).The source and receiver are set at [xs = 5 km, ys = 2 km] and [xr = 5 km, yr = 8 km], respectively.The color scale marks the isosurfaces.The kernel function is normalized to its value at [x = 5 km, y = 5 km].The vertical sections correspond to the white lines shown on the x-y plane.

Figure 2 .
Figure 2. Vertical section showing the 3D kernel function obtained using Equation (8).The color scale marks the isosurfaces.The kernel function is normalized to its value at [x = 5 km, z = −2.5 km].

Figure 3 .
Figure 3.The three-dimensional images obtained for Deception Island at 4 Hz are compared with the bi-dimensional images in Del Pezzo et al., 2016.Horizontal slices cut the inverse intrinsic, Q −1 i

Figure 4 .
Figure 4. Q −1 coda spatial distribution at Mount St. Helens.Horizontal slices are calculated at depths of 0.5 km and 4.0 km.The vertical section intersects the horizontal plane along the white line in the upper left panel.Topography isolines (only in the zone of Mount St. Helens) are superimposed.The discrete Q −1 coda spatial distribution has been interpolated before plotting the average inverse Q coda , < Q −1 coda >.All panels are drawn using Mathematica 10 TM .

Figure A3 .
Figure A3.Synthetic test at Mount St. Helens volcano.(Left) panels: test input, where the contrasts are expressed as a percentage with respect to the average, and correspond to Q = 50 (red) and Q = 500 (blue).(Right) panels: output, where the attenuation contrast is reproduced only in the center of the area.