Investigating the Apparent Seismic Diffusivity of Near-Receiver Geology at Mount St. Helens Volcano, USA

We present an expanded approach of the diffusive approximation to map strongly scattering geological structures in volcanic environments using seismic coda intensities and a diffusive approximation. Seismic data from a remarkably consistent hydrothermal source of Long-Period (LP) earthquakes, that was active during the late 2004 portion of the 2004–2008 dome building eruption of Mount St. Helens Volcano, are used to obtain coefficient values for diffusion and attenuation, and describe the rate at which seismic energy radiates into the surrounding medium. The results are then spatially plotted as a function of near-receiver geology to generate maps of near-surface geological and geophysical features. They indicate that the diffusion coefficient is a marker of the near-receiver geology, while the attenuation coefficients are sensitive to deeper volcanic structures. As previously observed by other studies, two main scattering regimes affect the coda envelopes: a diffusive, multiple-scattering regime close to the volcanic edifice and a much weaker, single-to-multiple scattering regime at higher source-receiver offsets. Within the diffusive, multiple-scattering regime, the spatial variations of the diffusion coefficient are sufficiently robust to show the features of laterally-extended, coherent, shallow geological structures.


Introduction
Reflection (or scattering) at an acoustic impedance surface is a principal process by which seismic energy propagates through the subsurface.Interfaces characterised by higher impedance contrasts will reflect more energy.In the framework of ray theory, the seismic source, the path the ray takes through the subsurface to the receiver, and the near-surface geology close to the receiver all contribute to energy propagation.While the separation of source and path effects is an extensively studied problem, receiver (or site) effects are often underestimated [1]; however, these effects make a significant contribution to the waveform [2] and, especially in volcanic environments, they may introduce errors into models and interpretations [3,4].This study links these site effects with the near-receiver geological characteristics, with the aim of discovering if seismic coda wave sensitivities are able to define shallow geological structures.
Given sufficient lapse-time from an initial P-arrival, a coda wave may be considered as a sequence of arrivals of reflected energy from the numerous heterogeneities encountered from source to receiver [5,6].The amplitude of these arrivals were found to be very similar at different stations for the same event, Geosciences 2017, 7, 130; doi:10.3390/geosciences7040130www.mdpi.com/journal/geosciencesregardless of distance from the source.What does vary, however, is the length of the coda.This forms the basis of the S-wave coda normalisation method (e.g., Aki 1980 [6]) and a theoretical basis to perform scattering tomography.By normalising the amplitude spectra at a given lapse time, energy is assumed to be uniform around the source at this point, allowing heterogeneities to be mapped spatially using a single-scattering assumption.The key assumption is that, at shallow depths, smaller-scale structures dominate over larger-scale variations, so much so that each arrival seen in the coda is less a function of a single reflection, and more a result (superposition) of multiple scattering events [7].These shallow effects are especially relevant when studying volcanoes, such as Mount St. Helens (MSH).The geology of this stratovolcano is complicated since it is the result of thousands of years of eruptions and landslides.Such complex geology results in a highly inhomogeneous medium, which then causes significant multiple scattering events as seismic energy propagates through it.By treating the entirety of the volcanic structure as a region of strongly scattering material, it can be assumed that, due to the stochastic nature of the medium, any seismic energy passing through will behave in a diffusive manner [8,9].By modelling seismic arrivals as the forefront of a diffusive "cloud" of seismic energy, the decay of the P and S energy is heavily affected by the strongly scattering structure close to the seismic receiver.Wegler (2001 [8] and 2003 [9]) argued that traditional seismic travel time tomography was not good enough to model the expected complexity of a volcanic interior.Hence, they developed the concept of multiple scattering in order to image small-scale heterogeneities.Using the assumption that a volcano can be described as a diffusive regime, the authors considered the decay of the seismic coda to obtain coefficient values for the diffusion and attenuation of seismic energy using an active source experiment at Merapi and Vesuvius volcanoes.
The coefficient values obtained with a diffusive approximation describe the way in which energy is scattered around a medium; the way in which this energy is dispersed changes depending on different lithological properties, even inside the same rock matrix.Figure 1 shows a cartoon representation of scattering behaviour for two different types of media.In panel 1a, a loosely consolidated structure results in very strong scattering of the propagating energy.On the other hand, panel 1b, shows how a small contrast between the acoustical properties of different crystals leads to very weak scattering behaviour.These processes are, evidently, frequency dependent.
the coda.This forms the basis of the S-wave coda normalisation method (e.g., Aki 1980 [6]) and a theoretical basis to perform scattering tomography.By normalising the amplitude spectra at a given lapse time, energy is assumed to be uniform around the source at this point, allowing heterogeneities to be mapped spatially using a single-scattering assumption.The key assumption is that, at shallow depths, smaller-scale structures dominate over larger-scale variations, so much so that each arrival seen in the coda is less a function of a single reflection, and more a result (superposition) of multiple scattering events [7].
These shallow effects are especially relevant when studying volcanoes, such as Mount St. Helens (MSH).The geology of this stratovolcano is complicated since it is the result of thousands of years of eruptions and landslides.Such complex geology results in a highly inhomogeneous medium, which then causes significant multiple scattering events as seismic energy propagates through it.By treating the entirety of the volcanic structure as a region of strongly scattering material, it can be assumed that, due to the stochastic nature of the medium, any seismic energy passing through will behave in a diffusive manner [8,9].By modelling seismic arrivals as the forefront of a diffusive "cloud" of seismic energy, the decay of the P and S energy is heavily affected by the strongly scattering structure close to the seismic receiver.Wegler (2001 [8] and 2003 [9]) argued that traditional seismic travel time tomography was not good enough to model the expected complexity of a volcanic interior.Hence, they developed the concept of multiple scattering in order to image small-scale heterogeneities.Using the assumption that a volcano can be described as a diffusive regime, the authors considered the decay of the seismic coda to obtain coefficient values for the diffusion and attenuation of seismic energy using an active source experiment at Merapi and Vesuvius volcanoes.
The coefficient values obtained with a diffusive approximation describe the way in which energy is scattered around a medium; the way in which this energy is dispersed changes depending on different lithological properties, even inside the same rock matrix.Figure 1 shows a cartoon representation of scattering behaviour for two different types of media.In panel 1a, a loosely consolidated structure results in very strong scattering of the propagating energy.On the other hand, panel 1b, shows how a small contrast between the acoustical properties of different crystals leads to very weak scattering behaviour.These processes are, evidently, frequency dependent.Attenuation can be considered a conductive process, where energy is transferred via physical contact rather than through radiative processes [10].As the seismic wavefield passes through a medium, energy is lost due to conversion to kinetic energy, via elastic perturbations of structural features such as sedimentary grains, or to thermal energy, due to friction between adjacent particles.How quickly energy is attenuated can be seen as a function of how much individual grains can move.In geological terms, freshly deposited sandstone, which is poorly consolidated, will attenuate seismic energy quicker than a much more consolidated, or even crystalline, rock.The presence of fluids may additionally enhance attenuation due to the loss of shear energy.
Diffusion can be described as a result of the spaces between the grains [11].As a process it is similar to seismic reflection.Energy is reflected off the sides of grains due to high acoustic impedances.The more of these high acoustic impedance surfaces there are, the more the seismic energy is 'scattered' around the medium.Some energy is lost at each scatterer due to transmission, but to a far lower degree as compared with the attenuation coefficient.For a review of how seismic diffusion and attenuation can be measured and behave in rock samples we refer to Mavko et al. (2009) [11].
By treating the diffusion and attenuation coefficients together as descriptors of the rate at which energy is lost in the medium, we can describe the geology in terms of consolidation and ability to attenuate energy.Attenuation is conceptually a simple process: higher coefficient values mean higher rates of attenuation and vice versa.The diffusion coefficient meanwhile, is closer to the inverse of 'scattering'.A lower diffusion value shows that energy is not leaving the local medium-just being constantly scattered around locally.Both are difficult to measure in the field but have wide implications for geophysical imaging in volcanoes [12,13].
This assumption forms the basis of a new way of using the diffusive approximation to depict surface geology, as presented in this paper.By applying the method of Wegler (2003) [9], we process a unique seismic dataset from the late 2004 portion of the 2004-2008 dome building eruption of MSH.The subsequent results are then compared and contrasted to the known near-surface geological structure as documented in the literature.

Mount St. Helens Volcano, USA: Geological and Seismic Structure
Mount St. Helens Volcano (MSH) is a region of diverse and complex geology, and one of several active volcanoes located in the Cascade Volcanic Arc of North America.The volcanoes in this area, and particularly MSH, have been extensively studied due to the potential hazard they pose on the metropolitan areas nearby and the destructive potential of their explosive eruptions.Significant volcanism in the region started around 37 million years ago [14] due to the subduction of the Juan de Fuca plate beneath the North American Plate, in a region known today as the Cascadia Subduction Zone.This process has resulted in a complicated network of stratovolcanoes, linked by intrusive bodies and interacting with a diverse range of hydrothermal fluids and partial melt that reach the shallow crust [15,16].
Due to the complex geological structure of the volcanic edifice and plumbing systems at MSH, the two broad scattering regimes identified in Wegler and Lühr (2001) [8] can be observed at the volcano using passive seismicity.Firstly, type A events, located at some distance from the edifice (typically greater than 10 km) and in the frequency range of 5-15 Hz [12].Secondly, type B events, which originate from within the volcanic edifice itself, and are of lower frequency, generally between 1 Hz and 5 Hz [15,17].These differences are dependent on both the source dynamics and media properties close to the source [12,13,15]: whilst type A events show clear P-and S-wave arrivals, type B do not show a clear S arrival.Since S-waves do not propagate through fluids this suggests that they are present close to the source of type B events, that is, fluid and molten materials are responsible for their generation [15].Here, we focus on the near-receiver effects affecting coda waves, that is, those created by strong near-receiver scattering, observed by De Siena et al. (2014) [12].Using type B events, the authors observe high scattering and high attenuation anomalies up to 6 km depth, indicative of fluid-rich zones within the volcanic plumbing system during the last eruption of the volcano.Kiser et al. (2016) [16] confirmed these results using travel-time velocity tomography.
Far from the source, path effects are primarily manifested in P-wave velocity models as high or low velocity anomalies.Lees (1992) [18] applies P-wave tomography (e.g., Nolet 1987 [19]) to Mount.St. Helens data recorded before and after its main eruption (1980) using a non-linear inversion method of hypocentral relocation and velocity recovery.His purpose was to try to determine the extent of magmatic bodies beneath the volcano, i.e., regions of low velocity, which can be attributed to regions of partial melt occurring in patches.The results show a ~1.5 km deep high-velocity zone, beneath the crater floor interpreted as a volcanic plug, directly above a low-velocity zone-interpreted as the magma conduit.This conduit broadens into a "chamber", which extends down to around 7.5 km depth.Beneath this region, data coverage was poor, but a second chamber was inferred to exist between 9 and 16 km depth by the same author.
Building on these results and adding seismicity recorded around the 2004 eruption, Waite and Moran (2009) [20] further improve the resolution in the upper 10 km of the crust.They found two igneous intrusive bodies NE and NW of the edifice as well as a structural boundary that runs parallel to the St. Helens Seismic zone (SHZ).These structures had already been retrieved by Moran et al. (1999) [17] and interpreted as a significant NNW trending structural anomaly, which ties in with the regional plate tectonics of a converging plate margin.The main difference between the Lees (1992) [18] and Waite and Moran (2009) [20] models is the location of a high-velocity anomaly.between depths of 3.5 km and 6.5 km.Waite and Moran (2009) [20] showed that this high-velocity body is 1 km deeper and ~1 km further to the east with respect to that of Lees (1992) [18].Whilst this may seem trivial, it does bring up implications for the interpretation of the extent of magmatic material associated with the seismicity in this region.This difference could be explained by the differences in the techniques used by the different authors, which is largely a matter of advancements in geophysics made in the last 25 years.
Recent velocity imaging from active surveys shows the most up to date picture of the plumbing system of the volcano [15].By utilising active source data from the iMUSH experiment, a low-Vp columnar anomaly was observed beneath and to the south-east of an inferred magma sill extending to 40 km depth [16].At the boundary of this region, deep long period seismicity is attributed to the injection of magmatic fluids into the crust.High-density accumulates, which may have a magmatic origin, are also observed and are thought to play an important role in the transport of such fluids within the crust.The low-velocity vertical anomaly gently dips towards SE; a low-scattering, high-intrinsic scattering anomaly retrieved at ~10 km depth by scattering tomography [12]; an high-conductivity anomaly retrieved by means of magnetotelluric tomography and interpreted as an extended crustal magma sill [21]; and a deep low-velocity anomaly, obtained using full-waveform tomography and interpreted as the shallowest extension of a mantle wedge [22].
Even though the geology close to the receivers has the obvious advantage of a visible surface exposure, the effects of the uppermost edifice on a seismic trace are very poorly constrained.Figure 2 shows a simplified geological map, which divides the area into 2 regions.Firstly, a younger, loosely consolidated group of Pleistocene-age volcanic rocks and debris flows.Secondly, an older, much more consolidated region of Oligocene-Miocene deposits.Overlain onto this map is the shallow subsurface extent of the Spud Mountain and Spirit Lake Plutons, along with a third body to the south, which is thought to be part of the Spirit Lake Pluton [23].Also shown is the extent of the 1980s debris flows, including some of the youngest deposits in this region, which likely show the least amount of consolidation.Poppeliers (2015) [24] demonstrated how a changing consolidation close to receivers has a disproportionally strong effect on the final waveform.By considering the partitioning of strain energy between P and S waves, it was found that sites with a pronounced stratification (layering), in particular, showed a strong frequency dependence on the ratios between direct energy and coda values.

"Drumbeat" Seismicity
Long-period (LP) earthquakes play a key role for driving main eruptions at Mount St. Helens.The 18 May 1980 eruption was preceded by the occurrence of two swarms of low-frequency seismic events and high values of the harmonic tremor indicated the action of interior pressurization promoting the weakening of the volcanic edifice [26].
Repetitive and shallow LP earthquakes were extensively recorded during the 2004 MSH eruption and quickly became known as 'drumbeats'.Although initially attributed to the large amounts of dome building occurring on MSH at the time-structurally manifested as a vertical strikeslip system [27]-the suggested mechanism for generating these 'drumbeat' events was that of

"Drumbeat" Seismicity
Long-period (LP) earthquakes play a key role for driving main eruptions at Mount St. Helens.The 18 May 1980 eruption was preceded by the occurrence of two swarms of low-frequency seismic events and high values of the harmonic tremor indicated the action of interior pressurization promoting the weakening of the volcanic edifice [26].
Repetitive and shallow LP earthquakes were extensively recorded during the 2004 MSH eruption and quickly became known as 'drumbeats'.Although initially attributed to the large amounts of dome building occurring on MSH at the time-structurally manifested as a vertical strike-slip system [27]-the suggested mechanism for generating these 'drumbeat' events was that of volumetric changes within a fluid-filled crack/conduit system.Waite et al. (2008) [15] showed how the first motions seen on the observed seismic data are dilatational (and so cannot be part of a strike-slip system) and originate beneath the dome itself at around 2000 m elevation.The actual 'crack' itself can be thought of as a fluid-filled fracture between two old lava flows.The LP events are then likely to be generated when magmatic activity below causes the pressure to rise within the fracture; once this increases past a threshold dictated by the frictional forces present, fluid then escapes, generating a characteristic LP event (drumbeat) and also allowing the crack to close [27].The crack then refills with more fluid and the process repeats itself.Figure 3 is a schematic illustration of the geometry of the source location, modified from previous studies [15,28].
volumetric changes within a fluid-filled crack/conduit system.Waite et al. (2008) [15] showed how the first motions seen on the observed seismic data are dilatational (and so cannot be part of a strikeslip system) and originate beneath the dome itself at around 2000 m elevation.The actual 'crack' itself can be thought of as a fluid-filled fracture between two old lava flows.The LP events are then likely to be generated when magmatic activity below causes the pressure to rise within the fracture; once this increases past a threshold dictated by the frictional forces present, fluid then escapes, generating a characteristic LP event (drumbeat) and also allowing the crack to close [27].The crack then refills with more fluid and the process repeats itself.Figure 3 is a schematic illustration of the geometry of the source location, modified from previous studies [15,28].The timing and energy of this repeating source was remarkably consistent, and so it presents an opportunity to treat this dataset as an active source experiment.Utilising the Pacific Northwest Seismic Network, data recorded between 26 September 2004 and 29 December 2004 were selected at stations located as shown in Figure 4 (red triangles).Approximately 2400 vertical component traces were selected, based on the timings of seismic arrivals and signal to noise ratios.This was done to ensure that we were processing "drumbeat" arrivals and that sites with particularly noisy datasets (e.g., YEL on the summit crater) did not distort the final results.The timing and energy of this repeating source was remarkably consistent, and so it presents an opportunity to treat this dataset as an active source experiment.Utilising the Pacific Northwest Seismic Network, data recorded between 26 September 2004 and 29 December 2004 were selected at stations located as shown in Figure 4 (red triangles).Approximately 2400 vertical component traces were selected, based on the timings of seismic arrivals and signal to noise ratios.This was done to ensure that we were processing "drumbeat" arrivals and that sites with particularly noisy datasets (e.g., YEL on the summit crater) did not distort the final results.

Diffusion Model Theory
One of the principle methods currently used in modelling wave propagation through random media is radiative transfer theory (e.g., Sato and Fehler, 1998 [7]).This technique assumes a velocity structure described by a random medium, and characterised by its autocorrelation function and the mean value of the background velocity.The propagating seismic energy is divided into two alternative processes: the free propagation in an effectively homogenous background medium, and the scattering caused by point-like scatterers.
At the extreme limit of strong multiple scattering, Wegler (2003) demonstrated that the diffusion model (e.g., Dainty and Toksov (1981) [29]) can describe a homogeneous and isotropic random medium by taking an acoustic approximation of the coda in Equation ( 1).

Diffusion Model Theory
One of the principle methods currently used in modelling wave propagation through random media is radiative transfer theory (e.g., Sato and Fehler, 1998 [7]).This technique assumes a velocity structure described by a random medium, and characterised by its autocorrelation function and the mean value of the background velocity.The propagating seismic energy is divided into two alternative processes: the free propagation in an effectively homogenous background medium, and the scattering caused by point-like scatterers.
At the extreme limit of strong multiple scattering, Wegler (2003) demonstrated that the diffusion model (e.g., Dainty and Toksov (1981) [29]) can describe a homogeneous and isotropic random medium by taking an acoustic approximation of the coda in Equation (1).
where W(r,t) is the density of seismic energy and p is a dimensional constant.For body waves p = 3 and for surface waves p = 2.As surface waves have been observed to strongly affect the coda at low frequencies (<3 Hz) [30,31], the dominance of body waves is assumed.The diffusion coefficient is related to the physical medium by the transport mean free path l tr and seismic velocity v by: where the transport mean free path is given by Gusev and Abubakirov (1996) [27]: In the case of non-isotropic scattering as assumed by this study, the angular dependent scattering coefficient is described by g(ϑ), where ϑ is the scattering angle and g tr describes how energy loses memory of its initial direction of propagation.With the strong multiple scattering observed in random highly-heterogeneous media, it is the parameter g tr that primarily controls coda generation [32].The reciprocals of g 0 and g tr are calculated to connect these coefficients to the physical medium.These are also known as the scattering mean free path l 0 and the transport mean free path l tr respectively.By taking the reciprocal of g tr we can connect the results to the physical medium in the form of the transport mean free path to fit the data to the diffusion model in the time domain, the method of Wegler and Lühr (2001) [8] is used.Multiplying Equation ( 1) by the geometrical factor t p/2 and taking the logarithm results in: This equation demonstrates that, for a fixed distance r as a function of time t, parameter W depends linearly on three base functions: 1, t, and 1/t.To calculate the three associated unknowns, a 1 , a 2 and a 3 , a linear least squares inversion is applied and using the following equations the physical parameters of E 0 , b and d are then given by: To ensure consistent data processing, an automated routine was developed in MatLab (see supplementary material) to apply the diffusion approximation to the seismic arrivals.However, some differences do apply depending on the location of the recording station.At MSH, we assume two separate scattering regimes that are both dominated by body waves.Firstly, a strongly scattering region near to the volcano and secondly, outside this area, a much weaker scattering environment.Since we cannot know the true length of the raypath within the strongly scattering volcanic edifice, distances have been simplified into two groups.Those that are found at stations off the volcano (CDF, ELK, FL2, MTM and TDL are the station short names) and stations that are found on the volcanic edifice itself (EDM, HSR, JUN, SHW, STD, SOS and YEL are the station's short names).For use in Equation ( 4), ray path lengths (r) are initially set to 10 km for the first group and 4 km for the second, respectively.
To calculate seismic energy density W, t min is set as the P-arrival and t max (see Figure 5) is defined when the amplitude approaches zero after the initial onset.The width of the energy windows used to measure W is defined as t max minus t min divided by 5.By using a small number of windows, we better depict the near-receiver geological overprint.

Data Processing
The first step to processing the inversion results is to apply an averaging procedure to calculate a mean value for each seismic recording station.Results from each station are separated and NaN or infinite values are removed.A histogram count is then applied, where the number of bins is defined as the number of traces for that station divided by two, minus one.This is performed to minimise the effects of a small number of arrivals at certain stations distorting the results with outliers.A mean is then calculated from those values which are within the limits of the most populous bin.By applying this procedure, results are weighted towards to the most common values.If two or more bins have the same number of values, the averaging limits are defined as the minimum and maximum values of the outer bins.An example of histograms plots is shown in Figure 6.Table 1 is the output of this procedure where on-volcano stations have been marked in italic.

Data Processing
The first step to processing the inversion results is to apply an averaging procedure to calculate a mean value for each seismic recording station.Results from each station are separated and NaN or infinite values are removed.A histogram count is then applied, where the number of bins is defined as the number of traces for that station divided by two, minus one.This is performed to minimise the effects of a small number of arrivals at certain stations distorting the results with outliers.A mean is then calculated from those values which are within the limits of the most populous bin.By applying this procedure, results are weighted towards to the most common values.If two or more bins have the same number of values, the averaging limits are defined as the minimum and maximum values of the outer bins.An example of histograms plots is shown in Figure 6.Table 1 is the output of this procedure where on-volcano stations have been marked in italic.The same data are plotted in Figure 7 with associated error bars.Where the error is defined as the standard deviation of the values used for averaging.For both groups of stations, the diffusion error is observed to increase with frequency but a far larger range in values is seen for off-volcano stations, which can be attributed to the breakdown of the diffusive approximation at far offsets.Another important difference is that values for on-volcano stations do not greatly vary as frequency increase.Whilst increasing the offset and change into the off-volcano regime, there is an obvious trend for values to increase with frequency.Values for the attenuation coefficient behave similarly for both scattering regimes in that there is a negative trend as frequency increases.However, the error range for on-volcano stations is typically larger.
The same data are plotted in Figure 7 with associated error bars.Where the error is defined as the standard deviation of the values used for averaging.For both groups of stations, the diffusion error is observed to increase with frequency but a far larger range in values is seen for off-volcano stations, which can be attributed to the breakdown of the diffusive approximation at far offsets.Another important difference is that values for on-volcano stations do not greatly vary as frequency increase.Whilst increasing the offset and change into the off-volcano regime, there is an obvious trend for values to increase with frequency.Values for the attenuation coefficient behave similarly for both scattering regimes in that there is a negative trend as frequency increases.However, the error range for on-volcano stations is typically larger.To assess the level of confidence that can be ascribed to these results as well as the impact of the ray-path length, r, simplification, theoretical diffusion curves are compared with the observed seismic energy density measurements, which have been stacked for each station.Figure 8 shows an example of a good model fit at the two offsets.Depending on the frequency, the fit becomes poorer as errors in coefficient measurement and uncertainty in r estimation increase.Fit values at the central frequency of 9 Hz, which showed the best fit at all stations, in our analysis range are shown in Table 2. To assess the level of confidence that can be ascribed to these results as well as the impact of the ray-path length, r, simplification, theoretical diffusion curves are compared with the observed seismic energy density measurements, which have been stacked for each station.Figure 8 shows an example of a good model fit at the two offsets.Depending on the frequency, the fit becomes poorer as errors in coefficient measurement and uncertainty in r estimation increase.Fit values at the central frequency of 9 Hz, which showed the best fit at all stations, in our analysis range are shown in Table 2.   To quantify the level of fit between the theoretical data and the observed, a two-sample Kolmogorov-Smirnov test was performed to assess the likelihood that both datasets came from the same distribution.This test quantifies the distance between the empirical distribution function of the observed data and the cumulative distribution of the theoretical.The outputs is a value of probability that both datasets came from the same distribution.With testing, it was observed that source-receiver distance played a large role in determining the diffusion coefficient and so controlled the theoretical seismic energy density.A simple routine was devised with fixed coefficient parameters that iterated values for r until a sufficient level of probability was acquired.Values for estimated r began at 20 km and were reduced by 0.1 km for each iteration.If the estimated value dropped below half the simplified distance, i.e., 4 km or 10 km, the routine was stopped, the probability threshold was lowered and the process restarted.New values for r are then applied to the inversion procedure to generate new diffusion and attenuation coefficients.By doing this we improve the fit of measured data to the theoretical diffusion approximation and provide a better estimation of the true raypath lengths that would occur within this model.Table 3 displays the new fit values, including values for estimated r and the diffusion coefficient.Figure 9 is a plot comparing the distribution of values of diffusion for the simplified and estimated results.Whilst there is some change in values, at the limited resolution of this methodology and the uncertainty at far offsets, the distribution of data points is sufficient to say that the use of simplified values for r is adequate.To quantify the level of fit between the theoretical data and the observed, a two-sample Kolmogorov-Smirnov test was performed to assess the likelihood that both datasets came from the same distribution.This test quantifies the distance between the empirical distribution function of the observed data and the cumulative distribution of the theoretical.The outputs is a value of probability that both datasets came from the same distribution.With testing, it was observed that source-receiver distance played a large role in determining the diffusion coefficient and so controlled the theoretical seismic energy density.A simple routine was devised with fixed coefficient parameters that iterated values for r until a sufficient level of probability was acquired.Values for estimated r began at 20 km and were reduced by 0.1 km for each iteration.If the estimated value dropped below half the simplified distance, i.e., 4 km or 10 km, the routine was stopped, the probability threshold was lowered and the process restarted.New values for r are then applied to the inversion procedure to generate new diffusion and attenuation coefficients.By doing this we improve the fit of measured data to the theoretical diffusion approximation and provide a better estimation of the true raypath lengths that would occur within this model.Table 3 displays the new fit values, including values for estimated r and the diffusion coefficient.Figure 9 is a plot comparing the distribution of values of diffusion for the simplified and estimated results.Whilst there is some change in values, at the limited resolution of this methodology and the uncertainty at far offsets, the distribution of data points is sufficient to say that the use of simplified values for r is adequate.

Diffusion
At the short lapse times from the P-wave arrival we consider, Calvet and Margerin (2013) [33] demonstrate that coda measurements are most sensitive to the effects of scattering, with temporal decay of the seismogram at later lapse times being mostly dependent on absorption effects.This is not necessarily valid at 3 Hz, where the LP sources show a resonant coherent component [15].In strongly heterogeneous media and using volcano-tectonic events, Wegler (2003) [9], Gao and Zhang (2013) [34] and De Siena et al. (2014) [12] confirm that at short lapse times scattering may be affected by local topography and is frequency dependent, with lower frequencies being most sensitive to heterogeneity at the surface.De Siena et al. (2016) [13] go further, showing the potential of coda waves for imaging the 30-m-thick debris flows of the 1980 MSH eruption.
These observations provide the framework for Figure 10.Results for the diffusion coefficient at 3 Hz are plotted at each recording station and gridded using an inverse distance method in Voxler (see supplementary material for an example data plot).The maps show strong local minima/maxima but can be seen to correlate with known near-surface geological features.A low-diffusivity anomaly is observed around the cone; as the offset from the volcanic centre increases, diffusivity is seen to increase rapidly.Since the diffusion coefficient is derived from the mean free path, diffusivity is the inverse of scattering.The strongest scattering is observed at stations positioned on Pleistocene-Recent deposits.Slightly smaller amounts of scattering are found at stations EDM and STD, which are located on landslide material from the 1980's eruption.Station HSR, and to a degree YEL, is notable in that is understood to be close to glacial deposits but not close to recent landslide material [35].In media where ice and water are trapped, diffusivity increases consistently resulting in the lower amounts of scattering observed [36].Outside of the strongly scattering regime, in the outskirts of MSH, the diffusion coefficient shows values that would be expected for single-to-weakly scattering media.Stations FL2 and CDF have average values, with the highest values being located at sites close to known plutonic intrusions.The variance that is seen within the two regimes is strongly dependent on the source-receiver distance and the coherency of larger-scale structures close to the seismic receiver.By increasing the offset value r, the diffusion coefficient rapidly increases, and so the simplification introduced in the model may mean that differences can be attributed to the resonant source signature or an incorrect estimation of r.
On the other hand, diffusion coefficient results at 3 Hz suggest that such differences can also be attributed to large coherent reflecting surfaces dominating the principle scattering direction.This is the case close to plutonic structures, such as station ELK.Here, the strong acoustic impedance contrast between the intrusion and country rock reduces the apparent diffusivity, as the interaction with the interface overshadows the effects of smaller-scale geological structure [37].The same result can be observed at stations positioned on landslide material, where instead of an increased amount of scattering that would be expected for such chaotic deposits, there is a relative decrease.In our interpretation, this is due to the creation of landslide blocks during the 1980s eruption [13], resulting in several large coherent slip planes for seismic energy to be scattered at.
Plotting results using a local polynomial method in Figure 11 yields tomographic maps that better correlate with previously observed regional trends.At 3 Hz a scattering anomaly following a NE-SW trend is observed.It encompasses the volcanic edifice and stretches to the NE.Between 6 Hz and 9 Hz this feature is shifted down and to the east.From 12 Hz to 15 Hz the amount of scattering rapidly increases in amplitude and area, fully comprising the volcanic edifice and the north of the study area.AT 18 Hz the anomaly is reduced and is in a E-W orientation but is shifted to the right of the volcanic edifice as observed at the lower frequencies.Further confirmation of the invalidity of the diffusion approximation at far-offsets is found when calculating the average transport mean free path, = 3 / 0, for on and off volcano stations.
Maintaining the assumption that there are two separate scattering regimes, velocity 0 is defined as 2 km/s and 5 km/s respectively.Within the frequency range 3-18 Hz, a value of 250 m is obtained for on-volcano stations and 430 m for off volcano.250 m is comparable with typical values of between 100 m and 1000 m found at other volcanic sites [8,9,[38][39][40][41][42] where scattering effects have been observed to dominate over intrinsic attenuation.However, away from the strongly scattering regime close to a volcano, typical crustal values are in the range of 100 km and so 430 m for off-volcano stations is several orders of magnitude lower than it should be.Again, this can be attributed to the uncertainty in ray path length and unsuitability of the diffusion model at these far offsets.Further confirmation of the invalidity of the diffusion approximation at far-offsets is found when calculating the average transport mean free path, l tr = 3d/v0, for on and off volcano stations.Maintaining the assumption that there are two separate scattering regimes, velocity v0 is defined as 2 km/s and 5 km/s respectively.Within the frequency range 3-18 Hz, a value of l tr ≈ 250 m is obtained for on-volcano stations and l tr ≈ 430 m for off volcano.l tr ≈ 250 m is comparable with typical values of between 100 m and 1000 m found at other volcanic sites [8,9,[38][39][40][41][42] where scattering effects have been observed to dominate over intrinsic attenuation.However, away from the strongly scattering regime close to a volcano, typical crustal values are in the range of l tr ≈ 100 km and so l tr ≈ 430 m for off-volcano stations is several orders of magnitude lower than it should be.Again, this can be attributed to the uncertainty in ray path length and unsuitability of the diffusion model at these far offsets.Increasing depth is associated with increasing frequency, and so the strongly scattering anomaly observed between 12 Hz and 15 Hz is attributed to a rapidly broadening magma chamber beneath and slightly to the NE of the volcanic edifice.
The interpretation follows that of De Siena et al. (2016) [13], where the attenuation and scattering parameters observed at different frequencies are depth dependent, with increasing frequency sampling deeper structure.At shallow depths (<200 m) the low frequency coda is strongly affected by surface waves [30,31] and so at 3 Hz the scattering anomaly is associated with the scattering volcanic edifice but is also strongly affected (as are all measured frequencies) by the presence of high scattering anomalies to the NE and E of the volcano.From 6 Hz to 9 Hz features are associated with the SHZ reducing the amount of observed scattering, suggesting a strongly coherent boundary forming this feature at 2.9 km depth.The increased amount of scattering that is observed between 12 Hz and 15 Hz can be correlated with the interpretation of De Siena et al. (2014) [12] of a very strongly scattering magma chamber.In the study, the authors observe a rapid broadening of this chamber which then floors off after 4 km depth.At 18 Hz features switch to a E-W trend that is associated with deeper regional magmatic structure and feeding systems [21].The interpretation follows that of De Siena et al. (2016) [13], where the attenuation and scattering parameters observed at different frequencies are depth dependent, with increasing frequency sampling deeper structure.At shallow depths (<200 m) the low frequency coda is strongly affected by surface waves [30,31] and so at 3 Hz the scattering anomaly is associated with the scattering volcanic edifice but is also strongly affected (as are all measured frequencies) by the presence of high scattering anomalies to the NE and E of the volcano.From 6 Hz to 9 Hz features are associated with the SHZ reducing the amount of observed scattering, suggesting a strongly coherent boundary forming this feature at 2.9 km depth.The increased amount of scattering that is observed between 12 Hz and 15 Hz can be correlated with the interpretation of De Siena et al. (2014) [12] of a very strongly scattering magma chamber.In the study, the authors observe a rapid broadening of this chamber which then floors off after 4 km depth.At 18 Hz features switch to a E-W trend that is associated with deeper regional magmatic structure and feeding systems [21].

Attenuation
Unlike the diffusion coefficient, which is structurally and source controlled, the attenuation coefficient is more sensitive to media properties.Quantities like pore fluid and relative seismic velocity control attenuation by either directly attenuating a propagating ray or by increasing the length of time the seismic energy is within a medium [11].However, in the diffusive approximation, the large uncertainties associated to average mean values (e.g., Figure 7) as well as the onset of "negative attenuation", which is suggestive of resonant processes [35] both translate in high uncertainty in the properties of the sampled medium.Some correlations with known geological features can be seen at 9 Hz, in Figure 12.As observed by Waite and Moran (2009) [20], MSH sits atop a significant NW-SE anomaly, which is related to regional plate tectonics.Between depths of 10 km and 14 km, the boundary shows the contrast between the south-western high-velocity and low-attenuation Siletz Terrane and the north-eastern low-velocity, high attenuation Cascade Arc Crust [12,43].
Spatial correlations can also be seen at 6 Hz with the results from De Siena et al. ( 2016) [13].In the study, the authors observe high coda attenuation values to the NW and SE of the volcanic edifice at 18 Hz, which correlate with regions of high attenuation as seen in Figure 10.The authors show a region where waveforms are characterised by higher S-wave peak delays, which stretches across the volcanic edifice in the NE-SW direction.High peak-delays are indicative of strong forward scattering and correlate with the lowest values of attenuation in our results; in our interpretation, seismic energy is quickly scattered away from these regions before being attenuated by local heterogeneities.

Attenuation
Unlike the diffusion coefficient, which is structurally and source controlled, the attenuation coefficient is more sensitive to media properties.Quantities like pore fluid and relative seismic velocity control attenuation by either directly attenuating a propagating ray or by increasing the length of time the seismic energy is within a medium [11].However, in the diffusive approximation, the large uncertainties associated to average mean values (e.g., Figure 7) as well as the onset of "negative attenuation", which is suggestive of resonant processes [35] both translate in high uncertainty in the properties of the sampled medium.Some correlations with known geological features can be seen at 9 Hz, in Figure 12.As observed by Waite and Moran (2009) [20], MSH sits atop a significant NW-SE anomaly, which is related to regional plate tectonics.Between depths of 10 km and 14 km, the boundary shows the contrast between the south-western high-velocity and lowattenuation Siletz Terrane and the north-eastern low-velocity, high attenuation Cascade Arc Crust [12,43].
Spatial correlations can also be seen at 6 Hz with the results from De Siena et al. (2016) [13].In the study, the authors observe high coda attenuation values to the NW and SE of the volcanic edifice at 18 Hz, which correlate with regions of high attenuation as seen in Figure 10.The authors show a region where waveforms are characterised by higher S-wave peak delays, which stretches across the volcanic edifice in the NE-SW direction.High peak-delays are indicative of strong forward scattering and correlate with the lowest values of attenuation in our results; in our interpretation, seismic energy is quickly scattered away from these regions before being attenuated by local heterogeneities.As before with the diffusion coefficient, gridding attenuation results with a local polynomial means that we lose the local scale sensitivity as assumed by the method, but better correlation can be observed with regional trends.In Figure 13, High attenuation areas correlate well with areas of low rise-times and Qc −1 at 6 Hz and 9 Hz, with high-rise times being associated with regions of very low As before with the diffusion coefficient, gridding attenuation results with a local polynomial means that we lose the local scale sensitivity as assumed by the method, but better correlation can be observed with regional trends.In Figure 13, High attenuation areas correlate well with areas of low rise-times and Qc −1 at 6 Hz and 9 Hz, with high-rise times being associated with regions of very low attenuation as shown in De Siena et al. (2016) [13].With increasing frequency, the rotation and shift of tectonic trend is observed to go from a strong NW-SE orientation to E-W with areas of high attenuation having a more westwards expression in the NE of the study area.3 Hz is very notable in that unlike the diffusion coefficient, it does not correlate well with the surface geology but instead follows the regional trend of the SHZ separating the Siletz and Cascade Arc terranes.Again, this can be attributed to surface waves having a strong effect on the seismic coda at these low frequencies, resulting in a combination of shallow and deep structure information being present in the results.attenuation as shown in De Siena et al. (2016) [13].With increasing frequency, the rotation and shift of tectonic trend is observed to go from a strong NW-SE orientation to E-W with areas of high attenuation having a more westwards expression in the NE of the study area.3 Hz is very notable in that unlike the diffusion coefficient, it does not correlate well with the surface geology but instead follows the regional trend of the SHZ separating the Siletz and Cascade Arc terranes.Again, this can be attributed to surface waves having a strong effect on the seismic coda at these low frequencies, resulting in a combination of shallow and deep structure information being present in the results.Attenuation (1/s) mapping using a local polynomial method.As with Figure 12, results correlate with peak delay times but show a strong dependence on large scale regional trends.

Conclusions
We use the diffusive approximation of the dispersion of seismic energy within strongly scattering media and data from a highly repetitive, natural source to measure station-dependent diffusion and attenuation coefficients at Mount St. Helens volcano, USA.Both factors are then spatially mapped as a property of the near-receiver medium.The results show important correlation with the local geology and improve our interpretation of up to 5 km deep volcanic structures.
At 3 Hz, values for the diffusion coefficient correlate with known near-surface geological structures.A weakly diffusive (strongly scattering) volcanic centre is observed close to the source,  12, results correlate with peak delay times but show a strong dependence on large scale regional trends.

Conclusions
We use the diffusive approximation of the dispersion of seismic energy within strongly scattering media and data from a highly repetitive, natural source to measure station-dependent diffusion and attenuation coefficients at Mount St. Helens volcano, USA.Both factors are then spatially mapped as a property of the near-receiver medium.The results show important correlation with the local geology and improve our interpretation of up to 5 km deep volcanic structures.

Figure 1 .
Figure 1.Representative scattering behaviour of seismic energy due to differences in medium structure.(a) Loosely consolidated media exhibit strongly multiple-to-diffusive scattering ray behavior; (b) Highly consolidated/crystalline media typically show weakly-to-single scattering characteristics.

Figure 1 .
Figure 1.Representative scattering behaviour of seismic energy due to differences in medium structure.(a) Loosely consolidated media exhibit strongly multiple-to-diffusive scattering ray behavior; (b) Highly consolidated/crystalline media typically show weakly-to-single scattering characteristics.

Figure 2 .
Figure 2. Summary Geological Map of Mount St. Helens Volcano, USA, adapted from United States Geological Survey (USGS) [25].Light brown shows Pleistocene-Recent age volcanic deposits and the darker brown, Oligocene-Miocene deposits.Pink areas and the dashed lines shows outcroppings of shallow subsurface plutons, whilst the red line maps out the extent of the 1980's debris flows.

Figure 2 .
Figure 2. Summary Geological Map of Mount St. Helens Volcano, USA, adapted from United States Geological Survey (USGS) [25].Light brown shows Pleistocene-Recent age volcanic deposits and the darker brown, Oligocene-Miocene deposits.Pink areas and the dashed lines shows outcroppings of shallow subsurface plutons, whilst the red line maps out the extent of the 1980's debris flows.

Figure 4 .
Figure 4. Topographical map of Mount St. Helens volcano.Receiver stations used in the study are marked by red triangles.The LP drumbeat source is identified as a red star approximately 200 m beneath the growing dome.

Figure 4 .
Figure 4. Topographical map of Mount St. Helens volcano.Receiver stations used in the study are marked by red triangles.The LP drumbeat source is identified as a red star approximately 200 m beneath the growing dome.

Figure 5 .
Figure 5. Z-component of an example seismic arrival at stations STD and CDF filtered at 3 Hz, 6 Hz, 12 Hz and 18 Hz.Calculated envelopes, from which energy density is measured, are plotted below with the analysis windows between tmin and tmax indicated as vertical bars.Note the reduction in tmax as frequency increases.

Figure 5 .
Figure 5. Z-component of an example seismic arrival at stations STD and CDF filtered at 3 Hz, 6 Hz, 12 Hz and 18 Hz.Calculated envelopes, from which energy density is measured, are plotted below with the analysis windows between t min and t max indicated as vertical bars.Note the reduction in t max as frequency increases.

Figure 6 .
Figure 6.Example histogram plots for station HSR at 15 Hz.Coefficient averages were taken from the most populous bins to minimise the effects of poor data distribution and outliers distorting results.When two or more bins have the same number of counts an average is calculated between the minimum and maximum values of those bins.

Table 1 .
Mean inversion results for the diffusion and attenuation coefficient values for each recording station with associated error.Values in brackets next stations are the total number of events used in the histogram count after NaN and infinite values have been removed.Values in brackets next to coefficient values are the total number of events used to calculate average values.

Figure 6 .
Figure 6.Example histogram plots for station HSR at 15 Hz.Coefficient averages were taken from the most populous bins to minimise the effects of poor data distribution and outliers distorting results.When two or more bins have the same number of counts an average is calculated between the minimum and maximum values of those bins.

Figure 7 .
Figure 7. Coefficient values plotted as a function of increasing frequency.Minimum and maximum of error bars are defined as the limits of the most populous histogram bin.

Figure 7 .
Figure 7. Coefficient values plotted as a function of increasing frequency.Minimum and maximum of error bars are defined as the limits of the most populous histogram bin.

Figure 8 .
Figure 8. Theoretical (dashed) vs. observed (solid) measurements for the logarithm of the averaged seismic energy density as corrected for geometrical spreading of body waves.By varying frequency the data fit can change quite drastically.

Figure 8 .
Figure 8. Theoretical (dashed) vs. observed (solid) measurements for the logarithm of the averaged seismic energy density as corrected for geometrical spreading of body waves.By varying frequency the data fit can change quite drastically.

Figure 9 .
Figure 9. Plot of diffusion coefficient values against source-receiver offset.Figure 9. Plot of diffusion coefficient values against source-receiver offset.

Figure 9 .
Figure 9. Plot of diffusion coefficient values against source-receiver offset.Figure 9. Plot of diffusion coefficient values against source-receiver offset.

Figure 10 .
Figure 10.Near-surface comparison of the diffusion coefficient (km 2 /s) at 3 Hz plotted using an inverse distance gridding method.The highest values for the diffusion coefficient are found closest to known intrusive deposits.The lowest values are found close to the volcanic edifice.Stations FL2 and CDF are found on older material and show average values.The map on the right is obtained interpolating the values obtained at different stations.

Figure 10 .
Figure 10.Near-surface comparison of the diffusion coefficient (km 2 /s) at 3 Hz plotted using an inverse distance gridding method.The highest values for the diffusion coefficient are found closest to known intrusive deposits.The lowest values are found close to the volcanic edifice.Stations FL2 and CDF are found on older material and show average values.The map on the right is obtained interpolating the values obtained at different stations.

Figure 11 .
Figure 11.Diffusion (km 2 /s) mapping gridded with a local polynomial method.3 Hz is associated with a combination of near-receiver coda wave and deeper sampling from the surface wave.Increasing depth is associated with increasing frequency, and so the strongly scattering anomaly observed between 12 Hz and 15 Hz is attributed to a rapidly broadening magma chamber beneath and slightly to the NE of the volcanic edifice.

Figure 11 .
Figure 11.Diffusion (km2 /s) mapping gridded with a local polynomial method.3 Hz is associated with a combination of near-receiver coda wave and deeper sampling from the surface wave.Increasing depth is associated with increasing frequency, and so the strongly scattering anomaly observed between 12 Hz and 15 Hz is attributed to a rapidly broadening magma chamber beneath and slightly to the NE of the volcanic edifice.

Figure 13 .
Figure 13.Attenuation (1/s) mapping using a local polynomial method.As with Figure12, results correlate with peak delay times but show a strong dependence on large scale regional trends.

Figure 13 .
Figure 13.Attenuation (1/s) mapping using a local polynomial method.As with Figure12, results correlate with peak delay times but show a strong dependence on large scale regional trends.

Table 1 .
Mean inversion results for the diffusion and attenuation coefficient values for each recording station with associated error.Values in brackets next stations are the total number of events used in the histogram count after NaN and infinite values have been removed.Values in brackets next to coefficient values are the total number of events used to calculate average values.

Table 2 .
Comparison of fit of observed data with simplified theoretical data at 9 Hz.

Table 3 .
Comparison of fit of observed data with estimated theoretical data at 9 Hz.

Table 2 .
Comparison of fit of observed data with simplified theoretical data at 9 Hz.

Table 3 .
Comparison of fit of observed data with estimated theoretical data at 9 Hz.