Evaluating the applicability of a screen di ﬀ raction approximation to local volcano infrasound

Atmospheric acoustic waves from volcanoes at infrasonic frequencies ( „ 0.01–20 Hz) can be used to estimate important source parameters for hazard modeling, but signals are often distorted by waveﬁeld interactions with topography, even at local recording distances ( ă 15 km). We present new developments toward a simple empirical approach to estimate attenuation by topographic di ﬀ raction at reduced computational cost. We investigate the applicability of a thin screen di ﬀ raction relationship developed by Maekawa [1968, doi: https://doi.org/10.1016/0003-682X(68)90020-0 ]. We use a 2D axisymmetric ﬁnite-di ﬀ erence method to show that this relationship accurately predicts power losses for infrasound di ﬀ raction over an idealized kilometer-scale screen; thus providing a valida-tion for the scaling to infrasonic wavelengths. However, the Maekawa relationship overestimates attenuation for realistic volcano topography (using Sakurajima Volcano as an example). The attenuating e ﬀ ect of di ﬀ raction may be counteracted by constructive interference of multiple reﬂections along concave volcano slopes. We conclude that the Maekawa relationship is insu ﬃ cient as formulated for volcano infrasound, and suggest future modiﬁcations that may improve the prediction capability.


Introduction
Erupting volcanoes produce atmospheric acoustic waves dominant in the infrasonic frequency range ("0.01-20 Hz) that can be used to detect, locate, and characterize ongoing eruptions, and to estimate important eruption source parameters [e.g. De Angelis et al. 2019;Johnson and Ripepe 2011;Matoza et al. 2019]. However, outdoor propagation effects complicate direct interpretation of source processes. For example, recordings of volcano infrasound have been used to estimate volume flux [e.g. Iezzi et al. 2019a;Johnson and Miller 2014;Kim et al. 2015], vent radius [Muramatsu et al. 2018], and crater geometry [Johnson et al. 2018a;Johnson et al. 2018b], but wavefield interactions with topography lead to variable source estimates from signals recorded at different azimuths [e.g. Iezzi et al. 2019a;Kim and Lees 2014;Lacanna and Ripepe 2013;Maher et al. 2020;McKee et al. 2014].
In particular, diffraction over the crater rim and other topographic barriers [e.g. Kim and Lees 2011;Lacanna et al. 2014] distorts waveform character and reduces acoustic power compared to predictions based on spherical spreading over flat ground. The effect of Corresponding author: smaher@ucsb.edu diffraction over topography on acoustic power can be accounted for in full wavefield numerical simulations [Fee et al. 2017b;Iezzi et al. 2019b;Kim et al. 2015;Lacanna and Ripepe 2013;Matoza et al. 2009], but a rapid analytical approach would be useful in some scenarios. For example, rapid correction of observed signals is desirable when time constraints prohibit the use of numerical modeling, such as when planning station locations for the deployment of a temporary infrasound network during a rapid eruption response. It could also be advantageous when infrasound is generated by valley-confined surficial mass movements [e.g. Allstadt et al. 2018] or eruptions from a crater with rapidly evolving morphology [e.g. Ortiz et al. 2018].
The simplest expression for acoustic power loss to diffraction assumes a single point of diffraction over a thin screen, as proposed by Maekawa [1968]. This relationship was empirically derived from experiments in which attenuation was estimated as the difference in sound pressure levels behind a rigid screen and the predicted values for free field propagation. Later studies incorporated the thin screen approximation into expressions for thick barriers [Fujiwara et al. 1977] and wedges [Maekawa and Osaki 1985]. These and various other numerical approaches [e.g. Hadden and Pierce 1981;Pierce 1974] deal exclusively with audible-frequencies (20 Hz-20 kHz) and require complicated calculations and knowledge of the barrier geometry relative to the source and receiver.
The thin screen approximation has been discussed in the context of volcano infrasound by Lacanna and Ripepe [2013] and Ishii et al. [2020], but the accuracy of predictions made according to the Maekawa [1968] relationship have not been evaluated. Lacanna and Ripepe [2013] used a 2D finite-difference simulation to show that synthetic amplitudes recorded on a topographically-obstructed side of the crater rim of Stromboli Volcano are attenuated by 8 dB compared to those recorded at the same distance on a less topographically-obstructed side of the crater. Ishii et al. [2020] computed theoretical pressure fields in the vicinity of thin screen approximations to topography at Sakurajima Volcano, where the 2D topography profiles are rotated to match the elevations of source and receiver, then replaced by flat ground with a thin screen at the height and distance of the tallest obstruction. Ishii et al. [2020] calculate idealized pressure fields as a sum of two contour integrals representing geometrical and diffraction aspects for monopole radiation in cylindrical polar coordinates [Li and Wong 2005;Macdonald 1915].
The results of Ishii et al. [2020] suggest that accounting for the attenuation by the thin screen provided a better prediction of the relative amplitudes at several stations than a 1/r spreading correction alone. However, Ishii et al. [2020] did not compare their attenuation results to the empirical relationship proposed by Maekawa [1968], which could provide simpler and more rapid estimates than their analytical approach. Furthermore, the results of Ishii et al. [2020] neglect full-wavefield effects associated with realistic topography (e.g. focusing of reflections down concave slopes) and may therefore overestimate the true losses.
We seek to directly evaluate the applicability of the empirical relationship for thin screen diffraction [Maekawa 1968] to infrasonic wavelengths in general and to volcano acoustics in particular. We use numerical modeling to generate synthetic waveforms for full wavefield propagation over flat ground and multiple barrier types, including Sakurajima topography, to obtain attenuation estimates directly comparable to the Maekawa [1968] experiments. Our approach is conceptually similar to that of Le Pichon et al. [2012], who used large numerical simulations to derive a semiempirical predictive relationship for infrasound attenuation due to the influence of stratospheric winds on the acoustic wavefield during propagation over a global scale. Analogously, we use numerical simulations to provide progress towards a semi-empirical relationship for infrasound attenuation due to the effect of topographic diffraction on wavefield propagation at a local scale. Our study differs from Le Pichon et al. [2012] in that we use synthetic results to investigate the applicability of an existing empirical relationship rather than to propose a new empirical relationship. These results represent a first step towards the development of an accurate predictive relationship and directions for future research are discussed in Section 9.4. This manuscript is organized as follows. Section 2 describes the thin screen diffraction approximation that is investigated here in the context of volcano infrasound. Section 3 describes Sakurajima Volcano and the observational dataset that is used in this study. Section 4 describes our primary methodology for numerically simulating infrasound propagation. Section 5 presents the main results of our simulations. Section 6 compares the results of our simulations to those from a different methodology that has been previously used with Sakurajima [Kim and Lees 2014]. Section 7 examines attenuation by diffraction in the frequency domain. Section 8 compares our synthetic results and thin screen predictions to observational data from Sakurajima. Section 9 discusses the implications of the study, limitations of our methodology, and directions for future work. Section 10 concludes the manuscript with the main findings and interpretations of the study.

Thin screen diffraction approximation
To first order, a barrier to acoustic propagation can be approximated by a screen with height and length much greater than its thickness. The acoustic shielding effect of a thin screen was quantified in a seminal study by Maekawa [1968], who investigated the efficiency of a thin screen to reduce sound pressure for general noise reduction applications (e.g. road noise). Maekawa [1968] recorded sound pressure behind a 2 m tall screen during a controlled indoor experiment and used the results to calculate the sound reduction compared to free field propagation. The shielding effect of the screen causes acoustic power losses during diffraction which can be related to the Fresnel number N [ Pierce 1981]. The Fresnel number is the extra diffracted path length normalized by half a wavelength: where R t is the length of the diffracted path, R d is the length of the direct path to a receiver (line-of-sight slant distance) and λ is wavelength. The orange and green lines over the volcano topography profiles in Figure 1 show example line-of-sight and diffracted paths, respectively, where the diffracted path is computed as the shortest possible source-receiver path over obstructions. Maekawa [1968] presents an empirical chart of insertion loss as a function of Fresnel number which has become a standard reference for first-order noise reduction estimates in the field of audible acoustics [e.g. Ekici and Bougdah 2003;Hohenwarter 1990;Plotkin et al. 2009;Vu 2007]. Lacanna and Ripepe [2013] present a

Figure :
Map of Sakurajima Volcano in southern Japan with m resolution around the crater and m resolution around the perimeter (interpolated to m resolution). Contours are marked at meters and the meter contour is marked in bold. The active vent is in Showa crater (upward triangle). Downward triangles represent infrasound sensors deployed for eight days in [Fee et al. ]. Profiles at right show topography between the vent and each sensor. Dash-dotted and solid lines show the direct and shortest diffracted raypaths, respectively, with lengths corresponding to R d and R t in Equation , respectively. logarithmic fit to this chart which was originally introduced by Tatge [1973] and is commonly used in engineering applications: where IL pT is the "predicted" insertion loss, in decibels, with subscript "T" indicating the study [i.e. Tatge 1973]. This expression carries the advantage of fitting the entire range of Fresnel numbers considered in the chart; however, the accuracy of the fit is reduced at N ă 1, because the abscissa is adjusted to make the fit a straight line in the range´0.3ă N ă 1 [Maekawa 1968]. Note that negative Fresnel numbers indicate regions with direct line of sight to the source, but close enough to the shadow zone to feature influence of diffraction effects.
Since the source-receiver geometries at Sakurajima mainly feature N ă 1, greater accuracy in this range is desirable. We therefore follow the approach of Yamamoto and Takagi [1992], who provide independent fits in the ranges´0.3 ă N ă 1 and N ě 1, with a maximum deviation of 0.5 dB in the range´0.3 ă N ă 1: for´0.3<N<1 10 log 10 pN q`13 for N ě 1. ( It is important to note that while the Maekawa [1968] chart is commonly used as a first-order noise reduction estimate regardless of source-receiver geometry [e.g. Ekici and Bougdah 2003;Ettouney and Fricke 1973;Plotkin et al. 2009;Scholes et al. 1971;Vu 2007;Young et al. 2015], Maekawa [1968] states that it is intended for cases where the source and receivers are in the air. This condition allows direct comparison between diffracted measurements and free-field estimates since the influence of ground reflection can be avoided. However, in volcano infrasound studies, the sources and receivers are commonly on the ground. A chart relating N to attenuation for the case when sources and receivers are on the ground is provided by Fehr [1951] and can be fit by the equation: IL pF " 10`10 log 10 pN q. (4) We find that Equation 4 generally predicts the attenuation for our synthetic volcano infrasound results better than Equations 2 and 3; however, the original paper by Fehr [1951] unfortunately does not provide explanation or justification for the relation. The chart is nonetheless reproduced in Figure 3 of Maekawa [1968] and was commonly cited before Maekawa's empirical relation became popular [e.g. Ettouney and Fricke 1973;Jonasson 1972;Purcell 1957;Scholes et al. 1971].
The predicted insertion loss given by Equations 2, 3, and 4 can be tested against experimental results described by modeled insertion loss (IL m ) [Lacanna and Ripepe 2013]: where p s is peak pressure behind an obstacle and p r is peak pressure at the same distance without the obstacle. Here, p s is peak pressure of the first wave cycle of a synthetic signal recorded behind a rigid obstacle using a finite-difference time-domain (FDTD) method for linear infrasound propagation [Section 4; de Groot-Hedlin 2017]. The synthetic peak pressure p r is recorded at the same receiver location but with a flat lower boundary. We mainly use peak waveform pressures instead of Fourier amplitudes at a specific frequency because signal frequency content may vary inconsistently with different source processes such that a single frequency component may not be dominant in every event. However, in Section 8 we compare Fourier amplitudes at 0.5 Hz between observations, predictions, and synthetics to allow direct comparison to results from Ishii et al. [2020]. We hereafter refer to Equations 2, 3, and 4 as "predicted" insertion loss and to Equation 5 as "modeled" insertion loss in order to express the difference between predictions and our results from numerical modeling. We hypothesize that modeled insertion loss will be well approximated by predicted insertion loss if the attenuating effect of the barrier in question (e.g. wedge, Sakurajima topography) is similar to that of a thin screen.

Sakurajima Volcano
Sakurajima is a highly active andesite-dacite stratovolcano in southern Kyushu, Japan, near the southern rim of the Aira Caldera ( Figure 1). The 20 km 3 edifice consists of lava flows and pumice fall deposits from at least 12 Plinian eruptions [Aramaki 1984]. From 1955 to present the eruptive activity has been characterized by ash-rich Vulcanian explosions [Ishihara 1985]. During the periods 1955-2006 and 2017-present, eruptions occurred from the summit vent of Minamidake, whereas from 2006 to 2017, eruptions occurred from a vent in Showa Crater on the southeast flank of Minamidake [Iguchi et al. 2013].
Due to the volcano's frequent activity ("1000 events per year since 2009) and proximity to the city of Kagoshima, Sakurajima is one of the world's bestinstrumented volcanoes [e.g. Iguchi et al. 2013] and a common target for infrasound studies. Acoustic waves generated by the explosive Vulcanian eruptions feature large amplitudes [449 Pa at 2.3 km; Fee et al. 2014;Matoza et al. 2014] and are sometimes visually observed as shock-like changes in luminance [Ishihara 1985;Yokoo and Ishihara 2007]. These studies and others [e.g. Aizawa et al. 2016;Cimarelli et al. 2016;Garcés et al. 1999;Iguchi et al. 2013;Kawakatsu et al. 1992;Maryanto et al. 2008;Miwa and Toramaru 2013;Morrissey et al. 2008;Smith et al. 2018;Tameguri et al. 2002;Uhira and Takeo 1994;Yokoo et al. 2009] illustrate the value of Sakurajima as an archetypal research target.

Finite-difference modeling
Our primary methodology involves numerical simulations of linear infrasound propagation with the FDTD code developed by de Groot-Hedlin [2017]. The model features a 2D cylindrical coordinate system in a sourcereceiver plane that is axisymmetric about the left boundary. This quasi-2D geometry models amplitude decreases with range identically to spherical spreading from a point source in 3D Cartesian space (1/r decay without topography). In a 2D Cartesian model, the amplitude decrease with range would instead reflect a line source (1{ ? r decay). References in this manuscript to "2D" modeling refer to the axisymmetric method of de Groot-Hedlin [2017] rather than a 2D Cartesian model. The simulations are advanced in time with finite differences over two time-steps using the modified midpoint method of Press et al. [1996], yielding second order accuracy in time of the solution. The program uses a staggered grid in which pressure is updated in cell centers while particle velocities are updated at cell boundaries. This staggered grid approach provides second order accuracy in the spatial derivates of the finite differences [Ostashev et al. 2005]. The lower boundary features rigid stair-step topography while the top and right boundaries feature absorbing perfectly matched layers [Berenger 1994].
The main results described in Section 5 are obtained from simulations performed in a 12ˆ12 km model space with a constant sound speed of 340 m s´1 and an ambient density of 1.229 kg m´3. We require a minimum of 20 grid cells per wavelength at 5 Hz, yielding grid spacings of 3.4 m and time steps of 2.4 ms. The accuracy of synthetic frequency components is generally considered acceptable for discretization up to 10 grid nodes per wavelength [Taflove and Hagness 2005]. In these simulations the discretization criterion of 10 nodes per wavelength is met up to a frequency of 10 Hz, so synthetic results can be considered accurate up to 10 Hz. However, we chose to only interpret frequencies up to 5 Hz to maintain confidence in the numerical accuracy of the results. Note that in Section 6 we present simulations with a maximum source frequency of 1 Hz and a constant sound speed of 349 m s´1 to facilitate direct comparison to 3D modeling by Kim and Lees [2014].
We run simulations with a variety of lower boundary conditions including flat ground, an idealized thin screen 800 m in height, (Figure 2, third row), a 1 km The height and distance of the thin screen is designed to scale the Maekawa [1968] experiment to infrasonic wavelengths and to create overlap in Fresnel numbers with the Sakurajima topography profiles. The geometries of the wedge and wide barrier are designed to create overlap in Fresnel numbers, facilitating direct comparison of the resulting attenuation to that of the screen and Sakurajima. We calculate Fresnel numbers for the synthetic stations assuming a wavelength of 340 m, corresponding to "1 Hz peak frequency in the synthetic signals, and the homogeneous sound speed of 340 m s´1.
For volcano topography we use Sakurajima Volcano as a case study because it is a common target for volcano infrasound research (see Section 3) and features rugged topography that complicates the acoustic wavefield [Johnson and Miller 2014;Kim and Lees 2014;Lacanna et al. 2014;McKee et al. 2014]. We use 2D topography profiles for each of five infrasound sensors deployed at Sakurajima in July 2013  ( Figure 1). We place the acoustic source inside Showa crater at the left boundary and propagate the wavefield over topographic profiles corresponding to the azimuths to the sensors in the 2013 deployment. For the lower boundary we use topography profiles from a 10m resolution digital elevation model (DEM) covering a broad area over Kagoshima Bay and a smaller 1-m resolution DEM covering the summit region and Showa crater. We fuse the two DEMs together after interpolating the 10-m resolution DEM to 1-m resolution over a rectangular mesh using a bivariate spline approximation. While our approach may be generalized to any topographic profile, the results presented here follow from previous work on the 2013 dataset Johnson and Miller 2014;McKee et al. 2014; and may inform future studies that use the same data.
Synthetic peak pressures are recorded at 100 m spacing behind each obstruction ( Figure 2) and compared to those recorded over flat ground to estimate IL m (Equation 5). We use peak pressure to calculate IL m such that the direct wavefront is measured rather than laterarriving reflections. Our approach allows us to thoroughly compare modeled insertion loss to predicted insertion loss by providing synthetic pressures at low computational cost for numerous source-receiver geometries representing a range of Fresnel numbers.

Results
The primary results of our study are the comparisons between predicted insertion loss from Maekawa's empirical relationship (Equation 3) and modeled insertion loss (Equation 5) based on the peak pressures of synthetics from 2D axisymmetric FDTD simulations. Fig-ure 2 shows wavefield snapshots for five simulations. The compressional source is centered at an altitude of 660 m, the height of Showa crater. In the flat ground case, most acoustic energy is concentrated at the wavefront and the singular ground reflection behind it. In the Sakurajima KUR case, multiple reflections down the concave slope (see also Figure 1) lead to wavefield complexity and focusing, as also noted previously by Lacanna et al. [2014]. In the cases of the thin screen, wedge, and square barrier, diffraction reduces amplitudes behind the obstruction compared to the portions of the wavefield propagating higher than the obstruction. Reflections from the mirrored obstruction across the axisymmetric left boundary appear in the last one or two snapshots, but these phases do not influence the peak pressure insertion loss results because only the first arriving wave cycles are used in the analysis. Figure 3A shows IL m results all synthetic receivers over each lower boundary type. For 0.1 ă N ă 1 the thin screen results (red circles) show good agreement (ă1 dB difference) with Equation 3 (black line). For the square barrier (cyan squares) and wedge (yellow triangles), the attenuation is underpredicted and overpredicted, respectively, by Maekawa's empirical relationship by a few dB, as expected from previous studies [e.g. Fujiwara et al. 1977;Maekawa and Osaki 1985;Pierce 1974]. Attenuation at synthetic Sakurajima receivers (diamonds in Figure 3A) is significantly overpredicted by Equation 3 and even negative at KUR at low Fresnel numbers (N ă 0.2), indicating gains in power relative to the flat ground case, rather than losses, due to topographic focusing effects. The volcano results are better predicted by Equation 4 (grey line), which is formulated for sources and receivers on the ground, but which only covers N ą 0.1.
As a guide, Figure 3B shows how the Fresnel number of each receiver varies with range (distance from the left boundary). Fresnel numbers generally decrease with distance as the difference between R t and R d decreases; however, at close ranges over the Sakurajima profiles, Fresnel numbers can increase with range as the propagation path becomes more obstructed by the crater rim.

Comparison of 2D cylindrical to 3D
Cartesian modeling To corroborate the results of our 2D numerical modeling in an axisymmetric geometry, we run additional simulations for direct comparison to previous 3D Cartesian finite-difference modeling with Sakurajima topography [Kim and Lees 2014]. As in the simulations presented by Kim and Lees [2014], we implement a homogeneous atmosphere with 349 m s´1 sound speed, a maximum source frequency of 1 Hz with a discretization of 30 nodes per wavelength. Synthetics are recorded at the five stations of the July 2013 de- We additionally run new 2D cylindrical and 3D Cartesian simulations that approximate flat-ground propagation for the same source-receiver distances so that attenuation by topography may be estimated using Equation 5. In the 2D cylindrical case, the flat ground simulation features a flat surface at sea level, with the source on the ground at the left boundary. Receivers are placed on the ground at distances corresponding to the slant distances (R d ) of the Sakurajima stations, measured from the edge of the source radius. We calculate Fresnel numbers for the synthetic stations assuming a wavelength of 873 m, corresponding to "0.4 Hz peak frequency in the synthetic signals, and the sound speed 349 m s´1. We approximate the source function of this simulation by placing a receiver 10 m outside the source radius. The source radius is 1386 m and the model dimensions are 12ˆ12 km.
Both methods use a staggered grid scheme [Ostashev et al. 2005] such that synthetic pressures are recorded at grid cell centers. This means that receivers cannot be placed directly on the ground and are instead located 5 meters above the surface in these simulations. When the height of a receiver is large compared to its range and the acoustic wavelength this may introduce undesired interference between the direct and groundreflected waves which would be in-phase on the surface. This effect is largely insignificant for the wavelengths considered here and at receivers more than a few hundred meters from the source. However, in the 3D Cartesian case, the source function receiver is placed 10 meters from the point source and thus is subject to phase-shifting in the ground-reflected wave. For this reason we approximate flat ground propagation in the 3D Cartesian case by recording pressure in the free field and multiplying the amplitudes by two. This approach accurately models the constructive interference of the in-phase direct and ground-reflected waves in the flat ground case by accounting for the image sources inherent to flat ground propagation [Kim et al. 2012;Morse and Ingard 1986]. For the 2D cylindrical case, the source function receiver is placed at a significant distance from the center of the source (1396 meters) compared to the 5.5 m distance between the ground and the receiver, making negligible the interference of ground-reflected waves. Figure 4 compares waveforms and power spectral density (PSD) estimates for the 2D cylindrical and 3D Cartesian simulations described above. In the 2D cylindrical cases the waveforms exhibit bipolar main pulses with relatively symmetrical compression and rarefaction phases ( Figure 4A), while the PSDs feature broad peaks around 0.4 Hz (Figures 4C, 4D). Waveforms from 2D cylindrical simulations with topography are generally time-delayed and reduced in amplitude compared to the simulations with flat ground ( Figure 4A). In contrast, the 3D Cartesian simulations feature strongly asymmetrical waveforms with single main compression phases ( Figure 4B) and nearly white power spectra (Figures 4E, 4F). The differences in waveform shape can be largely attributed to the different source functions. The 2D axisymmetric models uses a spatially distributed Gaussian pulse at time zero, while the 3D model uses a time-varying function (Blackman-Harris window) at a point-source. Figure 5 compares the distribution of peak amplitudes and modeled attenuation (IL m ; Equation 5) for the 2D cylindrical and 3D Cartesian synthetic cases described above. In both the 2D and 3D cases, peak am-plitudes for the flat-ground approximations decay linearly with distance from the source function receiver, as expected for spherical spreading ( Figure 5A and 5B, respectively). When topography is included, the amplitudes are reduced at each station in the 2D case (Figure 5A), whereas in the 3D case, the amplitudes are increased at ARI, KOM, and KUR ( Figure 5B). These differences are reflected in the IL m results ( Figure 5C), which feature positive attenuation at every station in the 2D cases but only at HAR and SVO in the 3D cases. However, both modeling approaches produce the same relative amplitude and distribution across the network (e.g. ARI is always more attenuated than KUR and less attenuated than HAR). The attenuation predicted by IL pY is clearly greater than the IL m results from both 2D and 3D modeling.
While both the 2D cylindrical and 3D Cartesian models present validated solutions to the wave equation [de Groot-Hedlin 2017;Fee et al. 2017b;Iezzi et al. 2019b;Kim et al. 2015;Kim and Lees 2014;Maher et al. 2020], several differences between the methods may account for the larger amplitude reductions in the 2D topography simulations compared to the 3D simulations. In the 2D model the source function is a spatially distributed Gaussian pulse, whereas in the 3D case the source is a time-varying Blackman-harris window at a single grid node. These different sources clearly pro-   duce wavefields with different frequency content as exemplified by the 0.4 Hz peak in the 2D flat synthetics ( Figure 4C) in comparison to the nearly white spectrum in the 3D free field synthetics ( Figure 4E). Since attenuation increases with frequency (see Section 7), we speculate that the less attenuated low-frequency components (<0.5 Hz) in the 3D model may contribute to the heightened amplitudes when compared to the 2D model. Our comparison between 2D and 3D simulations is inherently limited by these different source functions; however, adapting the methods to achieve equivalent source implementations is outside the scope of this study.
Additionally, the 2D model operates in an axisymmetric geometry, assuming azimuthal symmetry in topography around the source-receiver plane, whereas the 3D model uses a Cartesian geometry that allows for more realistic topography. Wavefield interactions with topography outside the 2D azimuthal plane may account for some discrepancies between the methods. For example, the line-of-sight propagation path for KUR follows a drainage channel (Figure 1) that may allow for constructive interference of reflections from the valley walls. This presumably increases recorded amplitudes at KUR compared to both spherical spreading over flat ground and to the 2D model, which does not account for the channel. Ray tracing [e.g. Blom 2020; Blom and Waxler 2012;Green et al. 2012] could potentially be used to investigate the influence of offpath reflections on the waveform; however, this is outside the scope of this study. Regardless, the preliminary comparison presented here corroborates the main finding of our 2D simulations, which is that attenuation of infrasound at Sakurajima is overpredicted by the screen diffraction approximation.

Frequency Dependence of Attenuation
While Equations 2 and 5 estimate attenuation for time domain amplitudes at fixed wavelengths (e.g. 340 m in Section 5 and 349 m in Section 6) the thin screen approximation also enables analysis of the frequencydependence of attenuation. We hypothesize that, if the thin screen approximation holds true for volcano infrasound, attenuation of 2D and 3D synthetics for Sakurajima topography will consistently increase with frequency to match the logarithmic relationship proposed by Maekawa [1968]. We investigate this aspect by first defining a new frequency-dependent variable for the Fresnel number`N˘, where the wavelength in Equation 1 is written as a function of ambient sound speed pc 0 q and a frequency vector pf q: This reformulated Fresnel number can be substitued into Equation 3 to express the empirical relation of Maekawa [1968] in the frequency domain´x IL p¯: x IL pY " # 5`8|N | 0.438 for´0.3 ăN ă 1 10 log 10`N˘`1 3 forN ě 1.
Finally, the predictions can be compared to estimates from numerical modeling´x IL m¯b y substituting the waveform amplitudes in Equation 5 for power spectra: wherep s is a PSD estimate (unit Pa 2 Hz´1) for a synthetic waveform from a simulation with topography at the lower boundary, andp r is PSD for a synthetic waveform from a simulation with a flat lower boundary. We examine the frequency dependence of predicted attenuation for five Sakurajima stations and modeled attenuation for the 2D and 3D simulations with Sakurajima topography described in Section 6. Figure 6 shows results for each station. As expected, x IL p increases as a logarithmic function of frequency, with y-intercepts controlled byN . Predicted attenuation increases by approximately 5 dB from 0.1 to 1 Hz, about half an order of magnitude in units of Pa 2 Hz´1. Similarly, x IL m for the 2D simulations generally increases with frequency; however, the values are lower than predicted by x IL p , and the increases are smaller (e.g. 2 dB increase from 0.1 to 1 Hz at SVO). Values of x IL m for 3D simulations are also lower than x IL p ; however, different patterns are observed for stations with relatively lowN (ARI, KUR and KOM) than for stations with higherN (HAR and SVO). For example, at KUR, 3D x IL m decreases with frequency and becomes negative above 0.1 Hz ( Figure 6C), whereas at SVO the 3D x IL m is relatively constant at approximately 3 dB below 1 Hz ( Figure 6E). Artificial increases of high-frequency energy in finite-difference models may be produced by spurious reflections from the corners in the stair-step topography [e.g. de Groot-Hedlin 2004]; however, this effect may be considered negligible at frequencies for which topography is finely discretized. We ran additional 2D and 3D simulations with finer discretization (5 m grid spacing) and observed no changes in the frequency range of interest ("0.1-1 Hz). We therefore interpret the decreased attenuation with frequency in the 3D model as a result of wavefield interactions with topography outside the source-receiver plane, such that they are not captured in the 2D model. This finding provides further evidence for the influence of full wavefield effects in counterbalancing power losses by diffraction.
The results of this frequency-dependent analysis provide further evidence to suggest that the thin screen approximation overestimates the diffraction losses for volcano infrasound. However, several limitations are inherent to this approach. The synthetic spectra are highly dependent on the source function, which differs significantly between the 2D and 3D methods (see Figure 4), and may be changed to model different source types [e.g. Iezzi et al. 2019b;Kim et al. 2012]. Furthermore, the approach cannot distinguish between the effect of diffraction and other wavefield-topography interactions such as reflections and scattering, which may also vary with frequency [e.g. Embleton 1996;Pierce 1981] and with topography outside the 2D azimuthal plane. Numerical artifacts may also influence the spectra. For example, spurious reflections may occur from stair-step boundaries in the discretized topography [e.g. de Groot-Hedlin 2004]. However, we consider this effect negligible in our case since we use a relatively fine grid spacing of 30 nodes per wavelength at 1 Hz, in comparison to the minimum of 10 nodes per wavelength [Taflove and Hagness 2005]. Further work to investigate these numerical limitations is outside the scope of this study; however, our results for the cases and frequencies considered suggest that attenuation by diffraction over topography is less than predicted by Maekawa's empirical relationship.

Comparison of Synthetic Results to Observed Data
In this study we primarily investigate numerical simulations because they allow direct estimation of wavefield changes for different lower boundary conditions; however, the ultimate goal of evaluating the thin screen approximation is to rapidly estimate attenuation in observed signals. Here we assess the applicability of our results to 30 eruption events recorded by five stations at Sakurajima during 18-26 July, 2013 ]. These events represent a subset of 74 short-term average/long-term average (STA/LTA) detections made by Matoza et al. [2014]. Events in this subset feature high signal-to-noise ratios (SNR) as evidenced by three or more stations with power spectra above the median of network-averaged noise spectra in the 0.1-to 10-Hz frequency band [e.g. Figure 7B; details in Maher et al. 2020]. The signals vary in waveform character from sustained low-amplitude jet noise to high-amplitude explosions (up to 449 Pa at 2.4 km) with asymmetrical shock-like waveforms Matoza et al. 2014]. To accurately reproduce the variety of observed signal types would require multiple modifications to the synthetic source functions, which is outside the scope of this study; however we observe a general agreement between synthetic and observed data as described here. Figure 7 shows an example of one event out of 30 and features an explosion followed by lower amplitude oscillations. The waveform shapes lie in between the symmetry of the 2D synthetics ( Figure 4A) and the strong asymmetry of the 3D synthetics ( Figure 4B). The shape of the power spectra agree with synthetics in that a broad peak is present below 1 Hz ( Figure 7B). The ob-served peak pressures (squares) are overpredicted by the 1/r relationship at each station, reflecting anomalously high amplitude at KUR. The amplitude decay for this event is also compared to the 2D and 3D simulation results, where the synthetics are scaled up to match the observation at KUR (circles and diamonds in Figure 7C). The 2D synthetic peak pressures follow 1/r at KOM and SVO while the 3D synthetics are reduced, more closely matching the observations. While Figure 7 exemplifies one event in the dataset, the amplitudes of all the events are summarized in Figure 8A. The boxplots show the distribution of observed peak amplitudes normalized to the amplitudes at KUR (p KUR ), while markers show the corresponding values for 2D synthetics (blue circles), 3D synthetics (green diamonds), thin screen predictions (red squares) and spherical spreading predictions (black line). The broad ranges in observed amplitudes reflect variability in source and propagation conditions across the 30 events recorded over an eight day period. These process potentially include variable source directionality [e.g. Iezzi et al. 2019b;Jolly et al. 2017;Matoza et al. 2013], eruption style [e.g. Fee et al. 2014;Matoza et al. 2014], winds [e.g. de Groot-Hedlin 2017; Kim et al. 2018;Lacanna et al. 2014;Sabatini et al. 2016], temperature inversions [e.g. Fee and Garcés 2007], and changes in vent/crater geometry [e.g. Fee et al. 2017a;Ortiz et al. 2018]. Our synthetics and predictions therefore do not provide a direct representation of each individual event; however, the general agreement in relative amplitude distribution validates our methodologies against the observations. Figure 8B shows the distributions considering the Fourier amplitude at 0.5 Hz (p w ) rather than peak pressure. Figure 8B may be compared to Figure 3c of Ishii et al. [2020], although the events considered are different, and screen predictions are based on Equations 2 and 5 rather than pressure field calculations. Use of the Fourier amplitude decreases the spread of observed values at each station and more clearly shows the overestimation of ARI, HAR, SVO values by the 1/r prediction.
The values from 3D synthetics fall within the peak pressure range of observations at all stations, while the values from other methods fall within the observation ranges at ARI, KOM, and SVO, but are overpredicted at HAR (largest N ). The predicted values according to the screen approximation generally agree with the values obtained by Ishii et al. [2020], although the value at KOM is slightly larger than the 1/r value. This finding shows that our method based on the empirical relationship of Maekawa [1968] gives similar relative amplitudes arrived at by Ishii et al. [2020]; however, we argue that it is misleading to normalize the amplitudes by KUR in this way. By making the values relative to KUR rather than to the source, the normalization obscures the overestimation of attenuation by the screen approximation in the case of volcano topography (Fig-ure 3A). Furthermore, the values obtained according to the screen approximation do not appear to provide an appreciable improvement in matching the observations when compared to 1/r or results from numerical modeling (Figure 8).

Discussion
We aim to develop a correction scheme for improved accuracy in volcano-acoustic source parameter estimates by evaluating a rapid method for predicting diffraction losses. A thin screen approximation to diffraction has been proposed in the context of volcano infrasound as an improvement over the 1/r spreading correction alone [Ishii et al. 2020]; however, the predictive ability of the original empirical relationship for diffraction losses [Maekawa 1968] has not been evaluated for volcano infrasound. In this study we use FDTD modeling [de Groot-Hedlin 2017] to show that Maekawa's original empirical relationship (Section 2) successfully predicts power losses due to infrasound diffraction over an idealized kilometer-scale thin screen. This finding validates the scaling of Maekawa's relation to large length scales and infrasonic wavelengths. Furthermore, the empirical relationship yields comparable relative amplitude distributions across the Sakurajima network to the approach taken by Ishii et al. [2020]. However, the absolute amplitudes predicted by the screen approximation overestimate the attenuation in the case of realistic volcano topography. We interpret this finding as a consequence of reflections from topography counteracting losses from diffraction.

Application of the thin screen approximation to infrasonic wavelengths
The IL m results from our scaled thin screen simulation ( Figure 3A) agree well with the empirical relationship of Maekawa [1968], indicating that Equation 3 is a reasonable approximation for infrasonic wavelengths recorded on the ground in the case of a kilometerscale thin screen. Furthermore, this agreement suggests that the 2D axisymmetric geometry of our numerical model space is adequate to describe losses for a thin screen. Moreover, the insertion loss results for a wedge and wide barrier generally agree with the empirical prediction (ă5 dB difference), suggesting that the screen approximation gives a reasonable first-order estimate of infrasound diffraction losses for a variety of simple long barrier types when receivers are on the ground. Equation 3 therefore provides a rapid method to estimate acoustic power losses for infrasound diffraction over simple barriers without resort to time-or computationally-intensive numerical simulations.

Figure :
[A] Example of observed waveforms for one event at Sakurajima at : 6: 6. 6 on / / (UTC). Amplitudes are normalized by the maximum at ARI ( 6 Pa) for ease of plotting.
[B] Multitaper PSD estimates for the waveforms in Figure A (original waveform amplitudes). The grey shaded region indicates the noise range across the network as defined by the th and th percentiles of networkaveraged PSD curves for % overlapping hourly time-windows during the deployment [Maher et al. ]. The dashed grey line is the median noise condition ( th percentile).
[C] Peak pressure vs distance for the waveforms in Figure A (squares). The peak pressures for D simulations with topography (circles), D simulations with topography (diamonds). The synthetic and predicted amplitudes are scaled to match the observed KUR amplitude. Black line represents expected amplitude decay from spherical spreading based on KUR amplitude. This example is one out of events used to construct the boxplots in Figure 8.

Implications for modeling diffraction of volcano infrasound
The IL m results from our 2D simulations with Sakurajima Volcano topography are systematically smaller than predicted by IL pY , with differences on order of 10 1 dB (Figure 3). These findings are corroborated by 3D simulations with Sakurajima topography, which feature even smaller IL m than the 2D results (Figure 5C). In the 3D simulations IL m is negative at stations ARI, KUR, and KOM, indicating gains in amplitudes relative to flat ground propagation ( Figure 5C). This suggests that this predictive empirical relationship is not adequate to describe losses to diffraction over complex volcano topography. The losses are better predicted by a different relation proposed by Fehr [1951] and represented by Equation 4, as shown by grey lines in Figures 3A and 5C. This relation is relevant in that it is meant for scenarios with sources and receivers on the ground; however, it is not very useful in that it is given no justification by Fehr [1951], and is only applicable at N ě0.1, whereas volcano infrasound recording geometries may have lower Fresnel numbers. Our conclusions contrast with those of Ishii et al. [2020], who suggest that the observed relative amplitude distribution at Sakurajima is better modeled when accounting for the effect of a thin screen than when only geometrical spreading is considered. We found comparable relative amplitude distributions to Ishii et al. [2020] using our method ( Figure 8); however we consider it misleading to normalize the screenpredicted amplitudes relative to station KUR in this way. Our synthetics show that attenuation at KUR is nearly zero or even negative ( Figure 5C), whereas the screen approximation would predict 8 dB attenuation for the corresponding Fresnel number. Direct use of the Maekawa [1968] empirical relationship results in overprediction of losses to diffraction in the case of volcano infrasound. Consequently we consider the screen approximation inappropriate, even if normalization by KUR causes the relative amplitude distribution to be comparable to observations. Furthermore, while Ishii et al. [2020] claim that the predictions from the screen approximation offer an improvement over predictions based on 1/r spreading for  (Figure 8). They used the standard deviation distance (S " |X´µ|{σ , where X is the predicted amplitude, µ is the mean of the observations and σ is the standard deviation of the observations) to show that S values are generally lower for the screen prediction than 1/r. We calculated S for the observations in Figure 8 and found larger values for the screen method than the 1/r method at stations ARI, KOM, and SVO (Table 1). While this analysis is clearly limited by small sample size (30 events) and the use of only four stations, it illustrates the lack of appreciable improvement of the screen method compared to 1/r.
We conclude that the primary cause of relatively small IL m values for volcano topography is the focusing of the acoustic wavefield by multiple reflections down concave-upward slopes (e.g. Figure 2). Focusing has been previously postulated at Sakurajima to explain amplitude variations with azimuth [e.g. Kim and Lees 2014;Lacanna et al. 2014]. Ishii et al. [2020] also observed that amplitudes from 31 Vulcanian explosions during August to September 2017 were consis- tently higher than expected from spherical spreading at one station (JIG). Notably, Ishii et al. [2020] used KUR as a reference station to estimate relative amplitude distrubtions based on 1/r, yet the topographic profile for KUR has one of the most concave slopes in the network. Out of the stations considered in this study, KUR also has the lowest IL m values, indicating the least attenuation losses to diffraction. Acoustic focusing at this station may explain why the KUR-relative amplitudes are larger than observed amplitudes at several stations including ARI, HAR, and SVO. We interpret that the concavity of the slope facilitates constructive interference of reflected waves which counteract amplitude losses from diffraction.

Limitations of our methodology
Our results suggest that the diffraction of infrasound over volcano topography produces significantly less attenuation than predicted for a thin screen with equivalent Fresnel number; however, several limitations in our methodology must be acknowledged. Firstly, our numerical simulations are limited in their treatment of atmospheric propagation effects which may increase or decrease power losses depending on the conditions. We use a homogeneous atmosphere to isolate the effect of diffraction and allow direct comparison to the empirical relationship of Equation 3, but in real outdoor propagation the sound speed generally decreases with altitude [Embleton 1996;Pierce 1981]. This vertical gradient enhances upward refraction of the wavefield and leads to reduced acoustic power at the receivers (greater IL m ) compared to the homogeneous atmosphere case. Lacanna et al. [2014] showed that wind and atmopsheric temperature gradients contributed to spatial and temporal deviations from 1/r amplitude decay at Sakurajima over a 21-month period. Additionally, atmospheric variability can change the propagation conditions over short time periods [e.g. Chunchuzov et al. 2011;de Groot-Hedlin 2017;Fee and Garcés 2007;Green et al. 2012;Iezzi et al. 2018;Johnson et al. 2012;Matoza et al. 2009;Ortiz et al. 2018;Sanderson et al. 2020]. For example, changes in atmospheric conditions in the boundary layer (lowermost "3 km of troposphere) caused up to "10 dB pressure differences in recordings of controlled chemical explosions at local distances [<10 km; Kim et al. 2018]. Consequently, our comparison between synthetic and predicted amplitudes to observations (Figure 8) is limited by the unmodeled atmospheric variability.
Secondly, the rigid topography used in our simulations may introduce minor variations in the synthetic waveforms compared to the observations. Real volcano topography can feature ground of finite-impedance (e.g. loose volcanic tephra layers) that may absorb some acoustic energy or convert it to seismo-acoustic coupled waves [e.g. Matoza et al. 2009] and consequently increase IL m values. A heterogenous distribution of impedance around the volcano could contribute to variations in the relative amplitudes across the network. This effect is commonly assumed negligible for infrasonic wavelengths and ă20 km path lengths [Matoza et al. 2019], but warrants further consideration since volcanic landscapes often feature loosely-consolidated near-surface materials (e.g. tephra layers or soil). Additionally, our simulations feature stair-step topography which may lead to spurious reflections and artificial diffraction side-lobes [de Groot-Hedlin 2017]; however, these are much smaller in amplitude than the compressional first arrivals that we use to calculate IL m . We also ran additional simulations with finer discretization (5 m grid spacing in both the 2D and 3D models) and observed no changes in the frequency range of interest ("0.1-1 Hz), indicating negligible influence of artifacts arising from topography discretization.
Finally, we acknowledge that the 2013 infrasound network at Sakurajima is not an optimal configuration for thorough diffraction analysis, since there are only five receivers and they are distributed azimuthally rather than radially. We attempted to improve this in our simulations by placing synthetic receivers at 100 m intervals on the topographic profiles ( Figure 2) such that we are able to represent a range of Fresnel numbers ( Figure 3). However, for comparison to observations we are limited to the five stations and 30 events (Figure 8). Further investigation of the thin screen diffraction approximation should target a volcano with more events and stations, preferably with a radial line configuration, such as at Yasur Volcano Matoza et al. 2017]. Our focus on Sakurajima has allowed di-rect comparison to the results of Ishii et al. [2020] in addition to previous studies that have addressed topographic effects at this volcano [e.g. Kim and Lees 2014;Lacanna et al. 2014;McKee et al. 2014;Yokoo et al. 2014].

Directions for future work
Our results suggest that the thin screen diffraction approximation is inadequeate to describe power losses for volcano infrasound, providing a first step towards the development of an analytical method for this purpose. Computationally-intensive numerical simulations remain the most accurate tool for estimating these losses, but future research could work towards a more rapid predictive relationship. For example, we only considered power losses over five topographic profiles at Sakurajima Volcano, but a larger simulation set involving multiple volanoes and azimuths could be used to derive a new semi-empirical relationship between N and IL m . Such an analysis could incorporate regional (15-250 km) to remote (ą250 km) recording distances to generalize the relationship for applications beyond local deployments. We note that in addition to N , it may also be useful to consider some measure of topographic roughness [e.g. 2D Fourier power spectral analysis of topography; Perron et al. 2008;Richardson and Karlstrom 2019], since multiple small barriers could reduce recorded acoustic power without increasing the Fresnel number.
Finally, our conclusions draw attention to the significance of acoustic focusing during propagation down concave slopes. This phenomenon has been previously postulated at Sakurajima [e.g. Lacanna et al. 2014] and seems to counteract the losses due to diffraction. The focusing effect could be further investigated with a geometrical acoustics approach such as ray tracing [e.g. Blom and Waxler 2012;Green et al. 2012]. A comprehensive correction scheme for topographic effects should account for the most significant processes; however, the relative contributions from diffraction, focusing, scattering, and reflections are currently unclear. Further analysis of these effects is required to achieve a rapid correction scheme to account for topography without costly numerical simulations.

Conclusions
We used synthetic pressure results from a 2D FDTD method for linear infrasound propagation in an axisymmetric geometry to numerically model the acoustic power losses to diffraction over kilometer-scale barriers including a thin screen, a wedge, a wide barrier, and Sakurajima Volcano topography. We found that the modeled acoustic power losses for the kilometer-scale screen are well-predicted by an empirically-derived re-lationship [Maekawa 1968] for diffraction of audible sound (20-20,000 Hz) over thin meter-scale screens, thus affirming the applicability of the relationship to infrasonic wavelengths. The power losses are also well predicted to first order for diffraction over the wedge and wide barrier, but the empirical relationship overpredicts losses for Sakurajima topography by "10 1 dB. We corroborate these results by comparison to 3D simulations with Sakurajima topography [Kim and Lees 2014], which also feature smaller attenuation values than predicted by Maekawa's empirical relationship. The relative amplitude distributions for 2D synthetics, 3D synthetics, and thin screen predictions generally agree with data from 30 eruption events at Sakurajima, but the screen predictions do not provide an appreciable improvement over the commonly used correction for geometrical spreading.
We conclude that the thin screen approximation proposed by Maekawa [1968] is an inappropriate model for diffraction by volcano topography because attenuation is overestimated in the topography case. This conclusion contrasts with that of Ishii et al. [2020], who conclude that their implementation of a thin screen approximation improves predictions of relative amplitude distributions at Sakurajima compared to geometrical spreading alone. We hypothesize that the true attenuation is smaller than predicted by Maekawa's empirical relationship because constructive interference of reflections down the concave volcano slopes (focusing) causes amplitude increases that counteract losses by diffraction.
Full wavefield numerical simulations such as finitedifference methods remain the most appropriate approach to account for topography [Iezzi et al. 2019a;Kim et al. 2015;Kim and Lees 2014;Lacanna et al. 2014;Lacanna and Ripepe 2013;Maher et al. 2020]. However, an analytical or empirically-derived method that depends only on readily-available topographic information is desirable for rapid interpretations of infrasound data in a rapid deployment or network planning scenario. Future work towards this goal could benefit from the development of a new semi-empirical relationship based on a large suite of numerical simulations with numerous volcanoes and a variety of recording distances.