Petro-mineralogical controls on coda attenuation in volcanic rock samples

S U M M A R Y Seismic attenuation measurements, especially those obtained from coda decay analysis, are becoming a key data source for the characterization of the heterogeneous Earth due to their sensitivity to small-scale heterogeneities. However, the relation between the scattering attenuation measured from coda waves and physical rock properties is still unclear. The goal of this study is to identify the main petrophysical and mineralogical factors controlling coda attenuation in volcanic rocks at the laboratory scale, as a necessary step before modelling seismic waves in real volcanic media. Coda wave attenuation was estimated from ultrasonic S-wave waveforms. To quantify the heterogeneity of the rocks and link them with this attenuation parameter, we performed several categorizations of the pore and grain systems of volcanic samples. Considering that seismic attenuation in rock samples can be modelled using the framework of wave propagation in random media, a statistical analysis of shear wave velocity fluctuations was performed: this analysis gives correlation lengths ranging from 0.09 to 1.20 mm, which represents the length scale of heterogeneity in the samples. The individual evaluation of the pore space and mineral content revealed that the pores of the samples (characterized by large vesicles) have a bigger effect than the grains on the heterogeneity level. We have developed a framework where intrinsic properties of the host rocks drive seismic attenuation by correlating the petro-mineralogical characteristics obtained from image data processing and analysis, with the coda attenuation measured at ultrasonic frequencies. There is conclusive evidence that porosity alone is not the primary controller of coda attenuation: it is also changed by the alteration level (i.e. oxidation, coating of the vesicles, secondary minerals) and the size of grains and pores. Among all the parameters analysed, it appears that the pore space topology is the main contributor to scattering attenuation in the volcanic samples.

Petro-mineralogical controls on coda attenuation 1859 experimental settings to measure a quantity like Q c −1 in heterogeneous samples, where both sample size (Yoshimitsu et al. 2016) and scattering (Mavko et al. 2009) dominate coda waves. This quantity is increasingly used as a marker of seismic absorption in field-scale imaging (Sketsiou et al. 2020).
Whereas seismic experiments at the ultrasonic scale routinely focus on the information contained in the wave-arrival package, in this study we focus on the information contained in the coda waves. In the field, the attenuation of coda waves Q c −1 is given by Q c −1 = Q i −1 + Q cs −1 (Calvet et al. 2013;Shapiro et al. 2000; following the single-scattering approximation of Aki & Chouet 1975), where Q i −1 and Q cs −1 refer to the coda attenuation produced by intrinsic absorption and scattering, respectively. An intrinsic energy loss can be triggered by various mechanisms like fluid/squirt flow, internal friction, viscosity and thermal relaxation (Barton 2006); these mechanisms are mainly controlled by the presence of fluids. In the absence of intrinsic attenuation, the amplitude decay is caused by the perturbation of the seismic wavefield (change in wave-direction and/or phase) due to changes in the medium properties (i.e. inhomogeneities). Measuring the decay of seismic energy with time in small samples is challenging: reflections and conversions of coherent waves on the boundaries of the sample are likely to contribute to this decay (Zhang et al. 2014;Yoshimitsu et al. 2016); the number of wavelengths that will propagate into the sample is constrained by the sample size (Guo & Fu 2007); and phase and amplitudes of the coda waveform depend on the experimental settings (Fujisawa & Takei 2009), especially on the coupling between the source/receiver and the rock sample. It is established that changes in coupling generate uncertainties in experiments, thereby hindering the repeatability of the analyses (McKenzie et al. 1982;Subramaniyan et al. 2014).
Coherent-wave scattering in heterogeneous media depends on the ratio of the seismic wavelength λ (related to frequency and velocity) to the size of the scattering heterogeneity d. Based on this relation, scattering can be described as Rayleigh scattering (λ > d), Mie scattering (λ ≈ d), and diffusion scattering (λ < d, Mavko et al. 2009). At the laboratory scale, the diameter of the scattering heterogeneity depends either on the pore space size or on the grain size (Matsunami 1991;Nishizawa & Fukushima 2008;Liu et al. 2017). The diameter of scattering heterogeneities could represent clusters of grains (or pores) that are comparable with the wavelength (Blair 1990;Lucet & Zinszner 1992). Therefore, we could expect that at sample scale scattering can be best described by the Rayleigh and Mie regimes. How large, if any, is the effect of the grain texture on seismic attenuation?  described how attenuation varies with the grain shape/size in polycrystalline materials. Their findings demonstrate that, in the Rayleigh regime, attenuation depends on the effective volume of the grain, while in the Mie regime it depends on the grain dimension in the direction of propagation.
There is a growing interest of the rock physics community for targeting volcanic rocks to provide more reliable quantifications of: elastic properties to be used in modelling volcanic processes (Heap et al. 2020;Pistone et al. 2021); attenuation mechanisms to enhance the interpretation of volcano-seismicity (Vedanti et al. 2018;Clarke et al. 2020); and microstructural changes at ultrasonic frequencies (Durán et al. 2019;Vanorio et al. 2002). Because recent advances in microscopy techniques and digital rock physics (Andrä et al. 2013;Cnudde & Boone 2013;Markussen et al. 2019) allow quantification of the spatial heterogeneity of the samples, we can use this measured heterogeneity to forward model both direct and scattered seismic wavefield.
Using this spatial heterogeneity as a marker of velocity fluctuations can provide a realistic representation of the coda wavefield (Holliger & Levander 1992;Spetzler et al. 2002;Fukushima et al. 2003).
There are three important challenges to overcome when characterizing Q c −1 at the laboratory scale: (1) the lack of a link between petro-mineralogical and scattering parameters affecting coda decay with time, analogous to those existing for phase-dependent imaging at subduction zones and mantle scales (Holtzman et al. 2003;Karato & Weidner 2008;Wagner et al. 2008); (2) the absence of laboratory calibrations (Tisato & Madonna 2012;Subramaniyan et al. 2014;Yoshimitsu et al. 2016) specifically targeting codadecay measurements and (3) the need for a joint description of coda waves in terms of radiative transfer theory and the wave equation, similar to those performed at field scale (Przybilla & Korn 2008;Obermann et al. 2013). Overcoming these obstacles will eventually enable the quantification of the petro-mineralogical state of volcanic and geothermal environments by using their dispersive effects on seismic waves, efficiently calibrating the interpretation of imaging studies that use scattering and absorption as seismic attributes (Takahashi et al. 2007;Calvet et al. 2013;De Siena et al. 2016;Mayor et al. 2016;Prudencio et al. 2017;Sketsiou et al. 2020).
In this paper, we tackle the first two challenges. To do so, we characterized heterogeneous rocks by an extensive petrological and mineralogical description of volcanic samples in terms of their pore and grain system, using optical and electron microscopy techniques along with image processing and analyses. An experimental link between the rock properties and seismic parameters at laboratory scale was established by: (1) determining the length scale of the heterogeneities and its relationship to porosity; (2) identifying the main contributor to scattering attenuation using coda wave analysis of ultrasonic S-wave waveforms; (3) analysing how the mineralogical composition of the samples impact coda attenuation parameters; (4) evaluating the effect of the grain size and alteration level of the sample and (5) assessing the use of 3-D image analysis versus conventional 2-D microscopy images to predict spatial heterogeneity. The results were discussed considering the trade-offs reverberations and surface waves generated by the limited size of the samples produce on coda envelopes. They provide valuable insights to interpret coda attenuation parameters as a marker of petro-mineralogical properties in volcanic rocks.

DATA
This study focuses on volcanic samples. The dataset comes from two cored boreholes (PTA2 and KMA1) on the Big Island of Hawai'i (Jerram et al. 2019). These boreholes penetrate lava flow sequences covering a diverse range of volcaniclastic facies. We chose 25 core samples at variable depth intervals, 18 cores from the PTA2 borehole, and 7 cores from the KMA1 borehole. The specimens are cylindrical cores of around 2.54 cm in diameter and 5.08-6.35 cm in length.
We provide a summary of the petrophysical properties of the samples in Table 1 (see Supporting Information S1-S3 for a detailed description and measurement procedures). All the experiments in this study were done in dry conditions, at ambient room temperature (∼18-20 • C) and atmospheric pressure (1 bar, ∼0.1 MPa). The samples have porosities ranging from 0.1 to 50 per cent (values highly influenced by the presence of vesicles), and ultrasonic P-wave velocity from 2.5 to 5 km s -1 for the entire data set. Cross Note. The porosity, permeability and density values were provided by Volcanic Basin Petroleum Research VBPR, while the ultrasonic P-and S-wave velocities were measured as part of this study (Section 3.1).
plots between the rock physics properties (velocity, porosity, permeability and density) do not follow trends close to theoretical models ( Fig. S2), an indication of the high heterogeneity of the samples.

M E T H O D O L O G Y
To relate seismic attenuation parameters with the heterogeneity level of the volcanic rocks, the samples were analysed in terms of (1) the coda attenuation measured at ultrasonic frequencies; (2) their petro-mineralogical characteristics, via measurements of pore space geometry, 2-D textural heterogeneity, and mineral phase distribution and (3) a statistical characterization of the random media.

Ultrasonic velocity experiments
In ultrasonic experiments, the seismic survey consists of generating elastic waves using an ultrasonic sensor (transducer) as the pulsing source. This pulse vibration travels inside the rock sample until it is acquired by a second transducer working as a receiver. The acquired waveform is a record of the change in time of the voltage pulse (source signal) due to its propagation through the specimen. The first step was to acquire waveforms at ultrasonic frequencies. We used a shear wave source as coda waves in the seismograms are primarily comprised of shear waves (Shapiro et al. 2000;Fukushima et al. 2003). The transducers belong to the Olympus Videoscan series. The shear-type piezoelectric crystal has a disk shape, 13 mm in diameter and a characteristic frequency of 1 MHz. The source and receiver were positioned at the ends of the central axis of the core sample. The transducers were placed in a box-holder coupled to an internal spring to guarantee constant pressure conditions, as it allows the coupling force to be independent of the assembly system ( Fig. 1). To ensure repeatability, these holders were attached to a bracket system to give stability and to hold the sample in the middle of both transducers while keeping them aligned. The pulse repetition frequency was set at 100 Hz and the pulser voltage at 400 V using an ultrasonic square wave pulser/receiver unit (Olympus 5077PR). The waves are received perpendicular to the transducer/sample interface (i.e. the longitudinal vibration travels along the vertical axis) and the output signal is displayed on a digital oscilloscope (Tektronix TDS 3012C). Each survey was conducted with identical sensors and parameters to keep consistency in the acquisition settings. Figure 2. Illustration of coda wave analysis on sample 1H. Recorded waveform after applying a Butterworth bandpass filter between frequencies 20-800 kHz (a) and its envelope (b); coda wave envelope (c); logarithm of the coda envelope and the linear fitting to obtain Q c −1 (d).

Coda wave attenuation measurements
Aki & Chouet (1975) described coda waves and modelled their energy loss as a function of frequency (f) and time (t) as: where E is the power spectra of coda waves, S(f) is the frequencydependent source factor, ∝ is a constant factor related to the geometrical spreading and Q c is the coda quality factor (inverse of coda attenuation), which may or may not be dependent on frequency. In our experiment, the frequency-dependent source factor [S(f)] and geometrical spreading are assumed to be constant for all the acquired waveforms. To obtain Q c −1 values we rearrange equation 1 by taking the natural logarithm: To compute the coda attenuation, we used the coda wave decay method: the envelope of the coda window is calculated through a smoothed Hilbert transform. Then, by cross-plotting the natural logarithm of the envelope for a central frequency of 150 kHz (left parameter in eq. 2; Q c computed at other frequencies in Supporting Information S8) versus time, we determine Q c −1 from the slope of a straight line fitting the data in a least-squares sense (Fig. 2d).
All the waveforms acquired in this study have a length of 4e-04 seconds from the origin time. To compute coda envelopes, we use a window starting 1.75e-04 s after the origin time with a length of 1.75e-04 s (Figs 2a and b). We discuss the implications of considering a different coda window as done by other studies (e.g. Guo & Fu 2007;Guo et al. 2009;Hu et al. 2018) in the Supporting Information S4, where we observed that the amplitude and frequency of secondary reflections can be disregarded, as their impact is weakened by the loss of high-frequency content in the coda (Fu et al. 2016;Figs S4.1 and S4.2).

Pore space analysis
3-D digital images of the internal microstructure (e.g. Fig. 3a), a proxy of the heterogeneities in rock samples (Andrä et al. 2013), were acquired for samples 1H, 9H, 11H, 14H and 18H, with a resolution of ∼22 μm using X-ray microcomputed tomography 'μxCT', a non-destructive acquisition technique (Cnudde & Boone 2013). Only five samples were analysed due to limited resources, the chosen five samples cover the full range of porosities (Table 1). Here, we describe the samples in terms of pore size distribution, aspect ratio and relative measurement of heterogeneity level at the pore scale (correlation length). We consider the resolution of the tomogram sufficient for these quantifications of the pore space, as the structures (pores) smaller than 0.022 mm are too small to be comparable to the size of the ultrasonic wavelength of our seismic experiments. The processing was done using the software Avizo (AVIZO 2019). The workflow can be summarized as: (1) Import the 2-D slice images generated from the micro-CT scanner (Versa XRM 500 from Zeiss-XRadia).
(3) Processing on the raw data: first crop the side-ends of the scanning to take out images with a shadow, which can compromise the intensity distribution, then apply a median filter for noise reduction.
(4) Perform manual segmentation of the data, separating voxels with CT values that refer to pores from those associated with the grains, then refine the segmentation by removing disconnected regions and smoothing to reduce the roughness in the grain-pore boundaries.
(6) Measure of pore space parameters such as equivalent diameter, sphericity, volume and area. The sphericity is expressed as: where V is the volume of a particle and A is its surface area. The equivalent diameter computes the diameter of the spherical particle of the same volume; it is given by

Textural heterogeneity
The analysis of the grain system was conducted on images derived from polished thin sections, based on the hypothesis that textural heterogeneity controls seismic scattering. These polished thin sections were made from a slice of the rock samples cut from one end  of all the core plugs. To describe their mineralogy, the samples were imaged using a Zeiss Axio Imager M2m optical imaging system under plane-polarized light conditions (Fig. S5). The 2-D heterogeneity was quantified using the textural arrangement, following the procedure describe by Mukerji & Prasad (2007). The spatial texture heterogeneity is quantified by three parameters estimated from the intensity images (relative to X-ray absorption): coefficient of variation (CV), correlation length (CL) and anisotropy ratio (AR). The procedure is described as follow and illustrated in Fig. 4 for sample 1H: (1) The CV (analogue to the 2-D textural heterogeneity) was assessed as the ratio between the standard deviation and the mean of the intensity images (Figs 4a and b).
(2) To obtain the CL we first computed the autocorrelation (ACF) of the entire image, next we extracted 180 radial profiles of the ACF (Fig. 4c). The correlation length was estimated for each profile as the lag value where the autocorrelation drops to 1/e of its maximum (Fig. 4d). Their average gave the CL.
(3) The AR was measured as the ratio of the maximum to the minimum correlation length over all the profiles.

Phase distribution
The composition and morphology of the samples were measured on the polished thin sections using an automated mineralogy mapping (TIMA: Tescan Integrated Mineral Analyzer acquisition system) to test if changes in chemical compositions affect seismic attenuation parameters. The backscattered electron (BSE) images contain information on the atomic number and X-ray emission of the element present, which allocates their chemical composition. The initial mineral segmentation output was reprocessed for quality control and to assign minerals to unidentified phases in the measured data. The processing consisted of analysing the unknown spectra in terms of the elements that must be present in the mineral and their characteristic atomic concentrations (e.g. Fig. S6). This allowed us to compute phase abundance and analyse various mineral texture properties and grain size distributions from the mineral mapping image (Fig. S5, Tables S6.1 and S6.3). The grain size is defined as the diameter of an equivalent circle with the same area as the detected grain. To report the grain size (Gs) in μm the TIMA software (TIMA 2019) use where A is the number of pixels inside the particle and P the pixel spacing.

Statistical characterization of a random medium
Our assumption is that scattering attenuation in the data set depends on the contrast between the elastic properties of the mineral grains and pores. At the sample scale, these contrasts are randomly distributed and can be described as a random medium. In such a medium the level of heterogeneity can be characterized by statistical functions (e.g. Gaussian, Exponential and von Karman), computing the spatial autocorrelation function ACF of their statistical velocity functions (Holliger & Lavander, 1992Sato et al. 2012). We thus use differences in grain size and mineralogical composition as a proxy for velocity fluctuations in the samples. We analysed the shear wave velocity fluctuations in the rock samples following a modified version of the methodology applied by Sivaji et al. (2002) for compressional waves, and Fukushima et al. (2003) for shearwaves to characterize the rock heterogeneities. Our workflow can be summarized as follow: (1) We use as the initial image the one corresponding to the mineral mapping (Fig. 5a). A square section of 256 mm 2 of the input image was transformed into a greyscale intensity image using ImageJ software (Fig. 5b).
(2) We compute the average shear wave velocity using the percentage of minerals in the samples and their velocities (Tables S6.1 and S6.2). The velocity values were selected from literature assuming isotropic minerals. The zeolites (observed on the thin sections during the microscopy analysis) show different levels of alteration in different samples, and the velocity for this phase group is then chosen within a range, where the chosen value depends on the level of mineral alteration; (3) We substitute the intensity data of each mineral with its shear-wave velocities (Table S6.2) to create a velocity image (Fig. 5c). For the vesicles we used air elastic properties (i.e. Vs = 0 km s -1 ).
(4) We compute velocity fluctuations by subtracting the average velocity.
(5) We apply a median filter to smooth the mapping from abrupt changes between small grains and pores, which are going to be averaged as part of the sample matrix. Then the main contribution in the fluctuations comes from particles (grains or pores) with size bigger than the median grain size in the sample.
(6) We generate 30 1-D profiles evenly distributed within the image (Fig. 5d).  Fig. 2d. Samples have been grouped in colour-zones based on the alteration level and the presence of secondary minerals. (b) Light microscope images for samples 1H, 5H, 9H, and 23H. The distribution of pore sizes (c) and pore sphericity (d) are shown for the five samples analysed using the 3-D digital images. A sphere has sphericity equal to 1, thus here 'irregular' corresponds to values lower than 0.25, 'low' to values between 0.25 and 0.5, 'medium' to values between 0.5 and 0.75 and 'high' to values higher than 0.75.
(7) We calculate the probability spectral density function (PSDF) by Fourier transforming the profiles and estimated the averaged PSDF (Fig. 5e).
(8) The ACF is obtained by the inverse Fourier transform of the averaged PSDF. The resultant ACFs are best modelled by the 1-D von Karman ACF: where x is the spatial lag, a is the correlation length (scale-length of the heterogeneity), ε is the rms of the fractional velocity fluctuations, (h) and Bh are the gamma function and modified Bessel function of the second kind of order h, h is the Hurst number rounding between 0 and 1. Our data fit for h = 0.10. The value of a and ε are estimated by fitting this function (eq. 6) to the computed ACF (Fig. 5f).

Coda attenuation at sample scale
The coda attenuation values in the samples range from 0.016 to 0.031, corresponding to coda quality factors (Q c ) between 32 and 61. Given the uncertainties associated with each measurement, coda attenuation (Fig. 6a) shows no simple dependency on porosity values. While coda attenuation seems to decrease with porosity, the lack of a simple polynomial attenuation-porosity trend suggests that, apart from porosity, other rock physical properties control attenuation in the volcanic data set being analysed. In addition to the porosity level, the pore space was characterized by its pore size distribution and sphericity level (Figs 6c and d). Pore size refers to the equivalent diameter of the pores of a given volume, and sphericity is a measure of how spherical the pores are. Characterizing the pore space is important as, for instance, samples 1H and 5H have similar porosity levels but different coda attenuation values and textural features (Fig. 6b). The discrepancies in their coda attenuation could be related to: (1) the differences in pore shape, the pore sphericity level is higher in sample 1H than in 5H, where pores follow a preferential direction; and (2) the fact that sample 1H has minerals of a size comparable to the pores, while sample 5H has a more homogeneous matrix of minerals characterized by smaller grain sizes. On the other hand, the sizes and shapes of pores and minerals in samples 9H and 23H are similar. Their different coda attenuation values could thus be explained by (1) the presence of secondary minerals, coating the boundaries of the vesicles in sample 9H, and (2) the moderate alteration (note the darker colour produced by the oxidation) of sample 23H. The dependence of coda attenuation on intrapore minerals could be related to changes in shear compliance (Barton 2006) produced by the intrinsic properties of these secondary minerals layering the pores. A robust link is expected for the cases in which the alteration level is controlled by secondary minerals with high clay content, as the relationship between clay content and (total) attenuation confirmed by Klimentos & McCann (1990) for compressional waves in sandstones. As the attenuation of clay minerals is highly influenced by water saturation and temperature (Biryukov et al. 2016), further experiments with fluid interaction could help to establish causality between secondary minerals (alteration) and coda attenuation.

Mineralogical composition of the samples
Five main groups of mineral phases were classified in the samples: feldspars, amphiboles, pyroxenes, olivine and zeolites. The remaining minerals represented less than 5 per cent of the sample composition (others, Fig. 7a). All the samples are matrix-supported, with a matrix mainly composed of minerals of the amphibole and feldspar groups, with some samples having a high concentration of either olivine or pyroxenes (or both). Zeolites with different alteration levels are encountered in some samples, mainly filling (or coating) the vesicles, but also as a detached mineral. The mineral concentration and grain size vary sample by sample (Figs 7b and 8). There is no clear correlation between the distribution of the phases and the coda attenuation levels; we argue though that the samples with median grain size larger than pore size (coincidental with the lowest porosity in the data set) and strong variations in mineral In the secondary y-axis, we observe the ratio between the grain size (for the largest 10 per cent) and the samples wavelength λ (computed at 150 kHz in red and at the dominant frequency recorded in the ultrasonic seismograms in green) reported in percentages. (b) grain size distribution for sample 1H, the horizontal axis shows midpoint ranges; bars indicate the mass per cent (area multiplied by density) of minerals. (c) Dependency of coda attenuation on grain size ratio, coloured by porosity level in the sample. (d) display a zoom of (c) for grain ratios between 2 and 6 to appreciate the distribution.
concentration (e.g. 18H, 2H and 17H, Fig. 7b) show the highest attenuation values (Q c −1 > 0.028). Once computed the mineral mapping of the sample, it was possible to calculate the equivalent diameter of the grains, described as grain size. The biggest grain size in the data set is less than 1.65 mm. The P80 (percentile 80) of the grain sizes is lower than 0.50 mm in most of the samples, this means that 80 per cent of the grains in the samples are smaller than 0.50 mm ( Fig. 8a; Table S6.3). This represents around 2.5-5 per cent the size of the wavelength at a central frequency of 150 kHz at which the attenuation values were computed (9-20 mm). Large grains correspond to the second peak of the bimodal distribution shown in Fig. 8(b) and about 20 per cent of the total grains, not included in the P80 statistic. What effect could these large grains have on the elastic parameters? We compare the ratio between the size of the large grains in the sample and their P80 (grain size ratio, Fig. 8c) to assess the influence of grain size differences in the samples on coda attenuation. Only a third of the samples have a ratio higher than 4. The samples with the highest ratio (1H and 19H) show low-to-average coda attenuation. When comparing the size of the large grains with the wavelength (dots in Fig. 8a), we consider that these grains influence the wave propagation only in samples in which the grain size is at least 5 per cent of the wavelength (dash line); for the rest, they will act as a matrix with average elastic properties.

Characterization of the heterogeneous medium
We quantified the heterogeneity level in the samples by three different approaches: the measure of 2-D textural heterogeneity, 3-D pore space analysis, and statistical characterization of the random media. The 2-D textural heterogeneity parameters (CV, CL and AR) cannot be directly correlated to attenuation values (Fig. 9a). The samples have a coefficient of variation CV between 40 and 50 per cent (which we consider as the heterogeneity level), except for sample 18H, in which CV is 29.8 per cent (this data tip was excluded from the plot to highlight variations for the rest of the data set). The anisotropy ratio varies from 1.38 to 1.85, except for sample 4H for which the anisotropy reaches 6.88 (again not included in the plot; this high value can be related to the presence of vesicles with a diameter as large as ∼8 mm - Fig. S5). The correlation Petro-mineralogical controls on coda attenuation   length goes from 0.04 to 0.40 mm. These values are an average of the textural mineralogy present in a 256 mm 2 segment of the thin sections. To understand how variable this quantification is at different scales, we divided the segment into 256 subsegments of 1 mm 2 (Fig. 9b). The range of values changes sample by sample, especially for CV, a marker of the heterogeneity of the data set. AR and CL have lower uncertainty than CV except for samples 1H, 11H, 5H and 7H.
The spatial autocorrelation of the pore space was best fitted by an exponential function where x is the spatial lag, a the correlation length and ε the standard deviation of the spatial fluctuations of the pores. Then we estimated the correlation length from the exponential fitting to the observed autocorrelation function (ACF - Fig. 10a). This correlation length is a relative measurement of the heterogeneity level associated with the presence of pores (Fig. 10b). Only five samples were analysed, the results will have to be tested on a wider data set to determine a dependency of coda attenuation on correlation lengths.
The statistical analysis of the random media consisted of determining the correlation length a of the samples, this time considering both pores and minerals, and building up the estimations from the shear-wave velocity fluctuations. The values of ε range from 8 to 13 per cent while a range from 0.09 to 1.20 mm (Table S7). This relatively small variation in ε indicates that the heterogeneities in the dataset have similar intensity (as was noted from the values of CV) but with distinctive characteristic scale lengths. The statistical method gave us by far the best correlation with the rock's physical properties (porosity). Samples with porosity lower than 30 per cent show that porosity increases linearly with increasing correlation length (R = 0.72), which translates as well into an increase of the scattering strength (Fig. 11).
For samples with porosity higher than 30 per cent, parameters other than porosity have a greater effect and the linear relationship is lost. For instance, if we compare samples 4H, 6H and 11H (which have similar porosity levels): (1) the anomalously high a value for sample 4H is due to the size of the vesicles; (2), the anomalouslylow a value for sample 6H is related to the presence of grains with sizes comparable to the pore sizes and (3) the lowest a value of sample 11H is associated to the irregular pore shapes and the presence of secondary minerals filling the vesicles. On the other hand, why does the sample with the highest porosity has such a small correlation length? The pores in sample 21H have high sphericity, the composition of the sample has one predominant phase mineral (80 per cent amphibolite), and the median grain size is 0.064 mm. Therefore, despite its high porosity, sample 21H is one of the most homogenous samples in the data set: porosity level on its own is not enough to estimate the level of heterogeneity in the sample.

D I S C U S S I O N S
Numerous studies have intended to differentiate between scattering and intrinsic (anelastic) attenuation. We consider that the estimated coda attenuation at ultrasonic frequencies in the studied samples is dominated by scattering attenuation based on the following assumptions: (1) the coda wave field is a product of the multiple scattering inside the samples (i.e. driven by the physical properties of the rock); (2) intrinsic attenuation is mainly related to the geometrical spreading which we consider relative equal between the samples (given the laboratory setting used for the experimental acquisition and the similar geometry of the samples); (3) as the system was not disturbed by external forces, the changes observed in the waveforms are create by changes in the wave propagation and (4) there is consistency of internal surfaces forces between grains from sample to sample, then other mechanisms that could contribute to anelastic attenuation (converting seismic energy into heat) at dry conditions (e.g. Frictional sliding between grain boundaries, grain contact adhesion) are neglected: they have an insignificant weight over the coda attenuation calculated in this study.

Scaling law of scattering attenuation
Following the scattering classification by Aki & Richards (1980), we projected the relation between ka and kL at a central frequency of 150 kHz (Fig. 12a), in which k is the wavenumber (k = 2π f/V, where f is the frequency and V-the S-wave velocity), a is the correlation length, and L is the length of the sample. Note that L is relatively constant (around 5 cm) in the data set, as we were looking to recreate similar conditions for the acquisition of the ultrasonic waves. Larger samples could allow a larger mapping of the scattering inside the sample. However, different interactions with the borders of the samples would be registered within coda waves, making our analysis more difficult.
In this study, we have quantified three different length scales of heterogeneity: when a corresponds to the correlation length from the pore space analysis (blue in Fig. 12a), the medium perturbation falls into the wide scattering regime: 0.1 < ka < 1 and, at such frequencies, it can be described by effective medium theory (Carcione 2015;Igel 2016). When a relates to the mineral grains then ka < 0.1: thus, if textural heterogeneity controls scattering, the sample could be defined as a homogeneous body (orange in Fig. 12a), with weak scattering for wavelengths corresponding to a frequency of 150 kHz. If the correlation length a is the one computed from the statistical fluctuations of the velocity in the media (yellow in Fig. 12a), the scattering in the samples varies between weak and moderate and can thus be modelled by using either effective medium theory or an equivalent homogeneous medium. As this is the method showing the lowest uncertainties, this result shows that at the sample scale we cannot work with a single scattering regime.
The need to work between scattering regimes is confirmed by using the ratio between the wavelength λ and the size of the heterogeneity d (Fig. 12b), a common approach in the rock physic community (Blair 1990;Mavko et al. 2009). Most of the coda attenuation measurements in the samples were taken in media that are described by the Mie scattering regime (less than λ/d ≈ 2π ). In sample 17H λ>10d, meaning that the heterogeneous rock behaves as a homogenous medium, in which the waves propagate through the sample without scattering on pores and grains; in that case, the contribution of scattering on the coda attenuation is low. In Fig. 12(b), λ is fixed at a frequency of 150 kHz: increasing this frequency results in the scatters moving towards Mie scattering (lower λ/d ratio).

Coda attenuation versus heterogeneity level
The primary hypothesis in this study was that 'as porosity is a direct expression of heterogeneity level in rock samples, and scattering attenuation is a measure of the heterogeneity in the media, then scattering attenuation is controlled by the porosity in the samples'. However, as we have demonstrated in this paper, many other factors influence the attenuation parameters at the laboratory scale, such as the shape and size of the pore space (Figs 6 and 10), the grain dimensions ( Fig. 8; also reported by , the mineral composition (Figs 7 and 11), and the presence of secondary minerals (Figs 6 and 11;consistent with Best et al. 1994). It is still impossible to establish the exact level of contribution of each factor. A similar level of complexity exists when correlating the velocity of elastic waves with petrophysical properties (Vedanti et al. 2018;Garia et al. 2019). Nonetheless, we have evidence that pore space is the main contributor to scattering (Figs 6 and 11). The correlation length from the statistical analysis provided the most reliable connection between elastic parameters and rock physics (Fig. 11). This result is consistent with previous studies that measured velocity fluctuations at the field scale (Holliger & Levander 1992).
The scattering in the data set is best defined by the Mie scattering regime ( Fig. 12; Liu et al. 2017), a result already obtained when analysing shales (Hu et al. 2018). This outcome is contrary to the notion that the Rayleigh regime is sufficient to describe scattering in rock samples at the laboratory scale, allowing to describe seismic waveforms in samples by analogy to the small-scattering assumption used in the far-field . In seismic exploration, the typical wavelength for shear waves is around 150 m at a dominant frequency of 20 Hz, then the scattering regime for a given scale of inhomogeneities in the media fluctuates based on the resolvable frequency content. If wavelength and scatterers have a similar dimension (ka and/or kL becomes larger), ray theory cannot be applied, and techniques specific to computational wave propagation and random media theory play a major role in modelling seismic energy ( Fig. 12a; discussed also in Igel 2016).

Laboratory calibrations for the acquisition of coda waves
Researchers have recognized and tackled the need for appropriate laboratory settings to measure seismic attenuation (e.g. Sivaji et al. 2002;Fukushima et al. 2003;Best et al. 2007;Guo & Fu 2007;Fujisawa & Takei 2009;Tisato & Madonna 2012;Tisato & Quintal 2013;Subramaniyan et al. 2014). Most of these works have targeted the attenuation of direct wave packages, but there is still a significant lack of laboratory guidance for measuring coda wave attenuation. One major drawback of ultrasonic experiments is the high uncertainty associated with the quantification of amplitudedependent parameters measured from the tail portion of the seismogram Yoshimitsu et al. 2016). As there is a need for calibrating stochastic attenuation parameters on real rock samples, we have assessed the best experimental setting to acquire suitable waveforms for performing reliable analyses on coda decay (Supporting Information S3 and S4). Despite these efforts, the uncertainty of the coda attenuation estimation is relatively high (Fig. 6a).

Upscaling rock physics properties of volcanic rocks.
Upscaling results between nanometres and micrometres, millimetres or kilometres has its trade-off. In this study, we assume that quantifications performed at the micrometre scale (image analysis) hold at the core scale (mm). For example, correlation length values show low variations when measured in subsamples of 1mm 2 in comparison to measures done over wider areas (256 mm 2 , Fig. 9). Similar results were obtained when measuring porosity using different methods (Fig. S1.1). We infer that measures done on 2-D planes are a reliable representation of the 3-D core. There are more systematic methods to assess if the size of the sample over which the properties are computed is large enough (referred to as REV Representative Elementary Volume, e.g. Singh et al. 2020) but this is beyond the scope of this paper; furthermore, no single method has been validated for all rock physics parameters.
Seismic measurements at the laboratory scale are taken at ultrasonic frequencies because the sensitivity of ultrasonic waves allows us to detect pores, grains, and cracks at a scale of less than half a millimetre. While the transducer used in the experiments emitted a signal with a characteristic frequency of 1 MHz, the output signal had a central frequency of around 150 kHz. This is half of the received dominant frequency for shear waves in sandstones for the same type of source transducers (Fu et al. 2020). Given the velocity ranges of the samples and their average length, we could only record one to three wavelengths per sample: therefore, we could not map individual heterogeneities of the order of the grain size (max 1-2 mm, but on average ∼0.1-1 mm). Grains and (micro)pores act as a group, turning into effective clusters with dimensions comparable to the wavelength (Lucet & Zinszner 1992). We acquired waveforms using compressional-wave transducers of higherfrequency (5 MHz) targeting smaller wavelengths, but the dominant frequency observed in the volcanic data set was again lower than 200 kHz.
A similar loss of high-frequency information is a central problem in volcanic imaging exploration of sub-basaltic reservoirs using direct wave phases and amplitudes (Gallagher & Dromgoole 2007;Eide et al. 2018). Coda attenuation modelling at the field scale shows that the coda wavefield is sensitive to lateral changes in scattering regime (from Rayleigh to Mie) that depend on the ratio between wavelengths and scale of heterogeneity (e.g. at Campi Flegrei caldera: De Siena et al. 2013 andDe Siena et al. 2017).
Pore space could thus be the main controller of coda wave techniques, widely used to image and monitor volcanoes. In order to upscale our results, it will be necessary to model coda waves at both field and sample scales. While this can be done analytically, for example in the Rayleigh regime and for more homogeneous samples , wave equation modelling and radiative transfer techniques including boundary conditions (De Siena et al. 2013;Obermann et al. 2013) are likely necessary in the volcanic case. Then, similar experiments and modelling will have to be performed in fluid-saturated volcanic samples to properly characterize the effect of fluids on coda waves.

C O N C L U S I O N S
We characterized volcanic samples in terms of their rock physics properties, pore space topology, texture and mineralogy to deduce their relationship with coda attenuation measurements. The parameters measured from the rocks (velocity, porosity and density) as a function of those measured from the scattered wavefield (coda attenuation) do not follow simple polynomial trends. Our findings suggest that: (1) The correlation length estimated from statistical functions gives a good approximation of the scale length of the heterogeneities. The intensity of the velocity fluctuations represented by the rms of the velocity fluctuations is quite similar within the data set. The differences in these heterogeneous rocks are defined by their scale length. For samples with porosity lower than 30 per cent, there is a linear trend of increasing correlation length with increasing porosity.
(2) There is no conclusive evidence to quantify the individual contribution of each of the analysed parameters in coda attenuation. However, the results point to pore space as the primary driver of seismic attenuation in volcanic samples, given that the heterogeneities triggering the perturbation of the ultrasonic seismic wavefield mainly correspond to the size and shape of the vesicles.
(3) The relationship between porosity and coda attenuation is likely influenced by mineral alteration in the samples (i.e. secondary minerals, coating of vesicles and oxidation), the pore size and shape (sphericity level) and the presence of grains with sizes similar to or larger than the pore size.
(4) The mineral phases and size of the grains determine if we can assume the sample as homogenous. When the grain size is not comparable with the vesicles size, or when the largest grains are smaller than 5 per cent of the wavelength size, the medium becomes an agglomeration of grains (independently of their mineralogy), which act as a bulk matrix with averaged elastic properties and where the wave travels without interference. Otherwise, the mineral grains and their mineralogy would render the matrix inhomogeneous.
(5) We consider that parameters computed from 3-D X-ray images and those computed from 2-D petrographic images are appropriate to characterize the heterogeneity level of the data set. However, the textural analysis was ineffective in establishing a relationship with coda attenuation.
The analysis of volcanic samples undertaken here strengthens the idea that more laboratory-scale studies should focus on complex samples (e.g. volcanic rocks), which better represent field conditions when trying to model the full seismic wavefield. Further modelling work will help define the individual contribution of the rock physics properties analysed in this study to the attenuation parameters. This new understanding will guide the identification of controlling parameters that can be upscaled to model geophysical observations in the field.

A C K N O W L E D G E M E N T S
We thank Dr John Millett (Volcanic Basin Petroleum Research VBPR) for providing the rock samples, Dr Maxim Lebedev (Curtin University) for helping us with the acquisition of the μCT images, and Prof Andrew Putnis for his insights in the mineralogical description of the samples. We also thank the Aberdeen-Curtin Alliance (http://aberdeencurtinalliance.org/) for providing funding for this research. Part of this study was undertaken using the TIMA instrumentation (ARC LE140100150) at the John de Laeter Center, Curtin University. This work was supported by resources provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia.

DATA AVA I L A B I L I T Y
The data underlying this paper are available in Mendeley Data, at http://dx.doi.org/10.17632/ysmt4kcjtn.2. Results-related data are available in the paper and its online supplementary material.