Anomalous Attenuation of High‐Frequency Seismic Waves in Taiwan: Observation, Model and Interpretation

High resolution maps of seismic attenuation parameters in Taiwan have been obtained from a modified “Multiple Lapse Time Window Analysis” (MLTWA). At most of the stations in porous sedimentary and highly faulted areas in Taiwan, the conventional modeling of MLTWA based on the scalar theory of radiative transfer in a half‐space with isotropic scattering fails to explain the spatio‐temporal distribution of the whole S‐wave train. Using Monte Carlo simulations of wave transport, we demonstrate that this anomalous energy distribution in the coda may be modeled by multiple anisotropic scattering of seismic waves. In addition to the scattering quality factor Qsc ${Q}_{sc}$ , we introduce a parameter g $g$ (independent of Qsc ${Q}_{sc}$ ) which determines the angular redistribution of energy upon scattering (scattering anisotropy). We inverted the attenuation parameters Qsc−1 ${Q}_{sc}^{-1}$ , Qi−1 ${Q}_{i}^{-1}$ and g in three frequency bands (1–2, 2–4, and 4–8 Hz). Overall, Taiwan is more attenuating than most orogens with a mean effective scattering loss Qsc∗−1=Qsc−1(1−g) ${\left({Q}_{sc}^{\ast }\right)}^{-1}={Q}_{sc}^{-1}(1-g)$ and a mean intrinsic absorption Qi−1 ${Q}_{i}^{-1}$ of 2.5 × 10−3 and 9 × 10−3 at 1.5 Hz, respectively. Scattering loss Qsc∗−1 ${\left({Q}_{sc}^{\ast }\right)}^{-1}$ varies over more than one order of magnitude across Taiwan while absorption variations reach approximately 30%. The more attenuating zones are the Coastal Range and the Coastal Plain where scattering dominates over absorption at low frequency, and inversely at high frequency. These regions are also characterized by strong backscattering ( g $g$ < −0.85) at 1.5 Hz and rather high VP/VS ${V}_{P}/{V}_{S}$ ratio. We propose that the observed strong back‐scattering at low frequency, related to large impedance fluctuations in the crust, is induced by the presence of fluids.

. Surface traces of active faults from Shyu et al. (2016) are shown in red lines. The geology of Taiwan is characterized by five major lithotectonic units (delineated with thick dark lines): from west to east, the Coastal Plain (CP), the Western Foothills (WF), the Hsuehshan Range (HR), the Central Range (CR), and the Coastal Range (CoR). The deformation front (DF), also the boundary between the CP and WF denoted by a dashed white line, extending northward from the northern Manilia trench offshore SW Taiwan to the northern coast of Taiwan marks the frontal limit of collision-related deformation. The longitudinal Valley (LV) corresponds to the collision suture zone between the EP and PSP. The Okinawa Trough (OT) behind the Ryukyu arc-trench system is a backarc rifting basin which migrates southwestward to the Ilan Plain (IP) in northeast Taiwan. The map also shows the seismological stations used in the present study (left panel) and the distribution of earthquakes (right panel) which occurred between 1994 and 2016 with a focal depth shallower than 40 km and a local magnitude greater than 3.5 reported from the Central Weather Bureau of Taiwan.
CALVET ET AL.

10.1029/2022JB025211
3 of 23 mélange complexes (Hsu, 1976). On the other hand, the northern Taiwan orogenic belt has transitioned from compressional collision into extensional collapse since the Plio-Pleistocene times (Teng, 1996). The lithospheric stretching has induced the post-collisional magmatism to produce the Northern Taiwan Volcanic Zone (NTVZ) including onshore Tatun Volcano Group (TVG) and Keelung Volcano Group (KVG) at the northern tip of Taiwan and several offshore volcanoes composed dominantly of andesitic lavas and pyroclastic rocks (K.-L. Wang et al., 1999) (Figure 1). Meanwhile, it has reactivated the ceased back-arc rifting of the Okinawa Trough (OT) and facilitated its southwestward extension onshore to the Ilan Plain (IP) in northeast Taiwan.
Over the past few decades, crustal elastic wave velocities of Taiwan have been intensively explored largely through tomographic inversion or forward modeling of massive P and S arrival times from local earthquakes and active sources (Huang et al., 2014;Y.-P. Lin et al., 2011). Though the resulting distributions of P-wave velocity ( ) , S-wave velocity ( ) , and their ratio ( ∕ ) have illuminated our understanding of the subsurface structure and tectonic evolution of the Taiwan orogen as well as the regional characterization of seismic ground motion and risk, they still remain insufficient to explain all seismic waveform signals and to infer physical properties and bulk compositions of crustal substances. As evidenced by ubiquitous coda trailing the direct ballistic waves and its exponential decay of amplitude, the crust of Taiwan being created, deformed, and destroyed through diverse geological processes exhibits the highly heterogeneous and anelastic nature. It is well known that they both contribute to seismic wave attenuation but its character and physical origin have not yet been fully investigated.
Previous attenuation models of the crust and upper mantle in the Taiwan region are mostly constrained from fitting the observed displacement spectra of direct P-and S-waves of local earthquakes at a single station or amplitude ratios at paired stations as a function of frequency to those predicted by an assumed Brune-type ( −2 ) source model (Brune, 1970). The corresponding corner frequency ( ) indicative of the finite dimension of an earthquake at which its source spectrum starts to drop abruptly, and the attenuation operator, * , defined as the integral of the reciprocal of quality factor along the path of travel, are then simultaneously determined for inferring the earthquake stress drop and dimension  and spatial distributions of and (Roth et al., 1999), respectively. Assuming frequency-independent quality factors, Y.-J. Wang et al. (2010) presented 3-D and models of the crust beneath the Taiwan orogenic belt from the tomographic inversion of * values measured at high frequencies of 2-30 Hz for P and 1-20 Hz for S. However, due to the lack of crossing ray paths and thus poor tomographic resolution in the near surface, there is neither obvious correlation found between the lateral variations of and properties of rock outcrops nor sharp change of attenuation across major geological boundaries or structural lineaments. The earlier attenuation models from the inversion of fewer * data estimated directly from the exponential decay of P-and S-wave amplitude spectra with frequency declared that the regions of lower and correlate with higher seismicity in the upper crust (<15 km depth), whereas it is opposite in the lower crust and upper mantle (K. J. Chen, 1998;K.-J. Chen et al., 1996). Cheng (2013) inverted the maximum acceleration of S-waves converted from seismic intensity reported by the Central Weather Bureau (CWB) of Taiwan for the structure, indicating the anomalously low regions extending from the upper crust to the lithospheric mantle at 75 km depth beneath the offshore volcanic island, Kueishantao, IP and NTVZ.
Despite the laboratory and field evidence for frequency-dependence of attenuation (Mavko et al., 2020), all the aforementioned studies assumed frequency-independent , overlooked the tradeoff problem between the estimated corner frequency ( ) and * , neglected the possible scattering loss causing the amplitude decay of high-frequency body waves, and implicitly attributed all the resulting variations to the anelastic properties of the crust and upper mantle varied by temperature, fluid/melt content, population of fractures or cracks, and etc. Though these models reveal some common features, such as relatively lower and in the volcanic areas of northern Taiwan, the alluvial deposit plain of southwestern Taiwan, and the high heat flow area under the SE CR, there still exist discrepancies as to details among them. For example, both Y.-J. Wang et al. (2010) and Cheng (2013) observed the significant lateral variation of attenuation across the CLF in their Q S models. However, compared to imaged in the foot-wall block of the CLF at depths ≤10 km, the former yielded relatively higher anomalies in the hanging-wall block between the CLF and Chuchih fault (CF) separating the WF and HR, while the latter resulted in comparatively lower , stronger attenuation therein.  introduced a cluster event method that substituted a conventionally single event with a cluster of nearby events pairing with stations for robust determination of magnitude-dependent and path-dependent * . Ko, Kuo, Wang, et al. (2012) applied it along with two-station spectral ratio measurements to imaging the variation across northeast Taiwan and the offshore southwestern Ryukyu subduction zone. Unlike the broadly CALVET ET AL.  (2013), a sharp, narrow zone of high anomaly which immediately abuts a low region beneath Kueishantao was found at tens to 90 km depths.
On the other hand, Chung et al. (2009) mapped the lateral variations of the coda quality factor at several narrow frequency bands between 1.5 and 18 Hz in Taiwan from local S-wave coda starting at twice the travel time of direct S. The results indicate higher in the western plain and eastern coastal areas and lower in the central mountain range at the central frequencies of 1.5 and 3 Hz but an almost reverse trend as the frequency increases. However, as demonstrated in Calvet and Margerin (2013), the measured strongly depending on the chosen lapse-time window which could lead to the biased variation of with earthquake hypocentral distance has not been cautiously investigated in this study. Zhang and Papageorgiou (2010) made a first attempt to separate the contribution of intrinsic absorption ( ) and scattering loss ( ) to total apparent S-wave attenuation ( ) by adopting a multiple lapse time window analysis (MLTWA) (Fehler et al., 1992) and found the obtained in accord with S-wave attenuation estimated by the coda normalization method. However, they only determined intrinsic and scattering at few sparsely distributed stations under the assumption of isotropic scattering. The results remain too preliminary and uncertain to shed light on the geological structure and bulk crustal properties in Taiwan.
In this work, we shall pursue the efforts put into the mapping of seismic attenuation across Taiwan with a special emphasis on the crust. While our approach is based on the classical MLTWA, our observations highlight some unusual spatial distributions of S-wave energy measured in a time window comprising the ballistic S-wave and its early coda. To interpret this original behavior, we introduce a non-isotropic scattering model which offers a simple physical explanation for the observations. Hence, in addition to the classical quality factors we are led to consider an anisotropy parameter that quantifies the tendency of seismic waves to scatter backward ( 0) or forward ( 0) . We examine both the frequency dependence and the spatial variations of the three attenuation parameters across Taiwan. In addition, we investigate the relation between the distribution of anomalous stations where non-isotropic is required to explain the data and salient geological and seismological features highlighted by previous investigations.
Our work has been made possible by the efforts of the CWB of Taiwan which, since the 1990s, has established the CWB Seismic Network (CWBSN) for monitoring regional seismic activities and routinely locating earthquakes hypocenters with manually-picked P and S arrival times. Beginning in 2012, the short-period, broadband, and strong-motion seismograph networks on more than 200 sites distributed in Taiwan and offshore islands are all equipped with high resolution 24-bit recorders and integrated to a unified data platform that has significantly increased the detection capability of local seismicity (Shin et al., 2013) and lateral resolution of crustal velocity heterogeneity across Taiwan (Huang et al., 2014). Owing to recent advances in theoretical formulation of seismic energy transport over a broad range of propagation regime and coda wave sensitivities to scattering and absorption quality factors in the anisotropically scattering crust Mayor et al., 2014), the vast amount of earthquake data set would open up new opportunities to comprehensively reassess the scattering and absorption properties of the crust in Taiwan.

Data and Observations
In this study, we make a first selection of earthquakes with a focal depth less than 40 km that occurred between 1994 and 2016. A total of 17,149 events with a local magnitude greater than 3.5 are extracted from a catalog provided by the CWB (Figure 1). The maximal magnitude in our data set is 7.07. Waveforms are collected on broadband, short-period and accelerometric networks operated by the CWBSN and BATS (Broadband Array in Taiwan for Seismology, Institute of Earth Science, Academia Sinica). A total of 236 stations are used (see Figure 1). We select data verifying the following criteria in the frequency band 2-4 Hz: (a) the hypocentral distance is smaller than 100 km (b) S-wave coda duration is greater than 75 s with a signal-to-noise ratio (SN) greater than 4; (c) the average intensity decreases monotonically in the coda window from the maximum of the envelop. This last criterion eliminates undesired signals such as aftershocks in the coda. Our final data set contains about 330,000 high-quality waveforms.
Next, the spatio-temporal distribution of the energy of the whole S-wave train is analyzed using the Multiple Lapse Time Window method developed by Fehler et al. (1992) and Hoshiba (1993). At every station, 5 of 23 we collect events with a hypocentral distance less than 60 km and a focal depth smaller than 15 km. The recorded waveforms are first filtered in three frequency bands ([1-2], [2][3][4], and [4-8] Hz) using bandpass Butterworth filters of order 3. Contrary to previous studies (Carcolé & Sato, 2010;Hoshiba, 1995), we were unable to exploit higher frequencies. The primary reason is the severe degradation of the signal-to-noise ratio in the S-wave coda, which entails the rejection of the vast majority of waveforms according to the aforementioned selection criterion (2).
Following the MLTWA procedure introduced by Fehler et al. (1992), we compute the instantaneous energy of the signals by summing the squares of the ground velocities on three components. Next, energy integrals denoted by ( = 1, 2, 3) are calculated in three consecutive time windows of 15 s duration. The first window starts 1 s before the S-wave arrival time in order to properly take into account the whole energy of the ballistic wave. To avoid contamination by the P-wave coda, the signal is set to zero for times smaller than ( − 1) s. To correct for the source magnitude and local site amplification, we compute the fourth energy integral, U 4 , in a 15 s duration time window starting 60 s after the origin time of the earthquake. We subsequently normalize the energy in the windows ( = 1, 2, 3) by the energy in window 4 , and correct for the geometrical spreading with a factor 4 2 . This last operation compensates for the divergence of 1 as the hypocentral distance approaches 0. In spite of all reasonable care taken to select the data and remove source and site effects, our measurements are sometimes contaminated by outliers. To circumvent this difficulty, we employ a simple denoising technique to the data: (a) energy ratios are sorted as a function of hypocentral distance; (b) the trimmed mean of the energy ratios is computed in moving a window of length seven samples and a stride of one sample. In effect, the minimum and the maximal values are removed before the mean is computed. Because the window is relatively short, the procedure does not introduce too much spatial smearing. Furthermore, stations with less than 20 measurements are systematically discarded. Figure 2 shows, in a semi-logarithmic scale, the normalized energy ratios, 4 2 ∕ 4 , as a function of hypocentral distance at two stations located in distinct tectonic units in Taiwan: CHK in the CoR and STY in the CR. At low frequency (1-2 Hz), the energy distributions at CHK and STY are clearly distinct. STY exhibits a "classical" behavior as previously reported in a number of studies based on MLTWA from the literature (see, as an archetypal example of MLTWA, the measurements in Japan by Carcolé and Sato (2010)): (a) at small hypocentral distance , the energy ratio between the first and the second window is greater than the energy ratio between the second and the third window, so that we have, in general, log 10 ( 1∕ 2) > log 10 ( 2∕ 3) ; (b) The normalized energy ratio in the first window U 1 decreases with hypocentral distance, whereas the ratios 2,3 show the opposite behavior. By contrast, CHK exhibits an "anomalous" behavior: (a) the energy is much more "uniformly" distributed in the coda: the energy ratio between the first and the second window is close to the energy ratio between the second and the third window, namely, log 10 ( 1∕ 2) ∼ log 10 ( 2∕ 3) ; (b) the three energy ratios log 10 ( 4 2 ∕ 4 ) show a clear curvature with a positive slope at small hypocentral distance. Most of the stations located in the WF and the CoR exhibit an "anomalous" energy distribution similar to CHK. At higher frequency (>4 Hz), we note that most stations exhibit the more "classical" energy distribution in the coda, indicating a change in the propagation regime.

Modeling the Anomalous Distribution of Seismic Energy
Conventional modeling of MLTWA is based on the scalar theory of radiative transfer (RT) in a half-space containing isotropic scatterers. In this model, the spatial distribution of seismic energy around frequency f depends on the scattering and absorption quality factors ( ) and ( ) only. To illustrate the unusual distribution of seismic energy at CHK, we show in Figure 3 the a posteriori fits of MLTWA measurements obtained with the conventional isotropic multiple-scattering model (gray lines). This model provides an excellent match to the data at station STY but fails at capturing the curvature of the spatial energy distribution near the source at CHK.
Several authors have previously argued that non-isotropic scattering is important to model the envelope shape of seismograms near the ballistic arrivals (Calvet & Margerin, 2013;Gaebler et al., 2015;Gusev & Abubakirov, 1996;Heller et al., 2022;Hoshiba, 1995;Zeng, 2017). Hence, it is natural to hypothesize that the assumption of isotropic scattering may be at the origin of the discrepancy between the observations and predictions of the RT model at station CHK. To confirm the plausibility of this argument, we have incorporated a non-isotropic scattering pattern into Monte-Carlo simulations of the multiple-scattering process, following previous works of Calvet and Margerin (2013). We note that for elastic waves, both forward and backward scattering are possible depend-  Figure S1 in Supporting Information S1. Blue, green and red dots correspond to the logarithm of the energy ratio in the first, second and third window, respectively, corrected for the geometrical spreading 4 2 , where is the hypocentral distance. In the map, surface traces of active faults from Shyu et al. (2016) are shown in red lines and the major lithotectonic units are delineated with thick dark lines. 7 of 23 ing on the nature of the heterogeneities and the frequency regime (Wu & Aki, 1985). A relatively simple and realistic analytical model of scattering non-isotropy which handles the two situations is provided by the so-called Henyey-Greenstein function (Henyey & Greenstein, 1941): In Equation 1, the unit vectors ′ and represent the directions of propagation before and after scattering, respectively. The parameter is the mean value of the cosine of the scattering angle and bounded between −1 and 1. Negative (resp. positive) values of correspond to preferentially backward (resp. forward) scattering, while isotropic scattering is obtained by setting = 0 . As shown by Margerin and Nolet (2003), the Henyey-Greenstein model may be understood as a limiting case of scattering in an inhomogeneous random medium characterized by a Von-Kármán power spectrum, when the Hurst exponent tends to 0. Physically, controls the roll-off of the power spectrum at large wavenumbers so that → 0 corresponds to random media that are rich in small-scale fluctuations. We refer the reader to Chapter 2 in Sato et al. (2012) for an in-depth discussion of heterogeneity power spectra.
The impact of non-isotropic scattering on the spatial distribution of seismic energy contained in the time window of ballistic waves + the first 15 s of the coda (referred to as the first time window hereafter) is illustrated in Figure 4. These simulations are carried out for scattering anisotropy that varies from predominantly forward ( = 0.85) to strongly backward ( = −0.95) . The isotropic case ( = 0) is also shown for reference. The values of intrinsic and scattering attenuation, as quantified by inverse quality factors, are fixed at the values that are typical of station CHK as summarized in Table 2 (isotropic model). It is customary to quantify the scattering attenuation with a length scale l termed scattering mean free path and defined as: with the S-wave velocity. It represents the attenuation length of coherent waves or the typical distance traveled freely by the waves between two scattering events.
In Figure 4, we observe that the main effect of forward scattering is to increase uniformly the normalized energy in the first time window, as compared to the isotropic case. To understand at least qualitatively this effect, it is necessary to introduce yet another quantity known as the transport mean free path and defined as: The transport mean free path is a crucial length scale as it governs the energy level in the late coda. Indeed, from first principles of multiple-scattering theory, it can be shown that the transport of coda wave energy is governed by a diffusion equation at long lapse-time, with a diffusion constant given by = * ∕3 (in 3-D space) (see Sato (2019)   8 of 23 excitation is proportional to −3∕2 . Furthermore it has been proposed that * is a better length scale than to quantify the attenuation of ballistic waves, as it takes into account the fact that only the energy scattered at large angles contributes to the attenuation of direct waves observed on seismograms (Sato, 1982;Wu, 1982). Hence, for a given value of the scattering mean free path, forward scattering simultaneously increases the energy contained in the ballistic waves and reduces the coda level at long lapse-time. When combined, these two effects explain the strong increase of the normalized energy contained in the first time window in the case of forward scattering, as compared to the isotropic case.
Following the same line of reasoning as in the forward-scattering case, we expect backward scattering to simultaneously increase the attenuation of ballistic waves and the coda level at long lapse-time. As anticipated, we observe in Figure 4 that backscattering tends to reduce the normalized energy in the first time window. Furthermore, backscattering changes the curvature of the curve representing the energy in the first time window as a function of hypocentral distance from straight to concave. This effect is also clearly observed in the data at station CHK. Hence, strong backscattering of the waves appears as a good candidate to explain the anomalous seismic energy distribution at station CHK. This hypothesis is confirmed by the examination of Figure 3, which shows the bet-fitting non-isotropic models obtained by a least-squares inversion of the data where we let −1 , −1 and vary (see the next section for implementation details). On the one hand, at station STY, anisotropic scattering is not necessary to explain the data and the best-fitting value of is close to 0. On the other hand, at station CHK, the introduction of strong back-scattering ( = −0.99) improves dramatically the fit between the multiple-scattering model and the observations. This result supports the use of the Henyey-Greenstein model to represent the scattering anisotropy of seismic waves. In the next section we present an inversion method to simultaneously retrieve absorption, scattering attenuation and anisotropy in a MLTWA approach.

Inversion Method
Following Carcolé and Sato (2010), we implement a Levenberg-Marquardt algorithm to minimize a misfit function defined as the sum of the squared differences between the base 10 logarithms of observed and modeled normalized energy ratios: where the index runs over all earthquakes and the superscripts and refer to model-predicted and observed data, respectively. To guide the inversion, we compute the partial logarithmic derivatives using the approach of Takeuchi (2016) for the last two parameters. To compute the sensitivity with respect to , we use the basic ideas of the differential Monte-Carlo method (Lux & Koblinger, 2018), which allows to compute the seismic energy in two media with slightly different anisotropy parameters (unperturbed) and + (perturbed) in a single simulation. It suffices to note that for a given sequence of scattering events, the energy detected in the coda in the perturbed and unperturbed media are related by: In Equation 5, each factor of the product is equal to the probability to be scattered from direction to direction +1 in the perturbed medium, divided by the probability of the exact same event in the unperturbed medium. Each factor should therefore be interpreted as a weight which restores a posteriori the correct frequency of occurrence of the scattering from direction to direction +1 in the perturbed medium. Applying a Taylor expansion to Equation 5, we arrive at the following formula for the partial derivative of the energy with respect to the anisotropy parameter : where the logarithmic derivative of the Henyey-Greenstein function is given by: It should be noted that in practice Equation 6 has to be averaged over a large number of realizations to obtain a sufficiently accurate estimate of both the energy and its partial derivative. Empirically, we found that the simulation of 4 × 10 6 independent trajectories was sufficient for our purposes.
Since the Levenberg-Marquardt algorithm is a local optimization method, it is important that the starting model be sufficiently close to the best-fitting model in the parameter space, in order to avoid a local minimum trap. We therefore adopt a two-step inversion scheme to estimate the attenuation parameters −1 , −1 and g at every station. First, −1 and −1 are estimated for 19 fixed values of ranging from −0.95 to 0.85 with an increment Δ = 0.1 . Second, starting from the best solution obtained at step 1, the three parameters ( , , ) are optimized with the aid of the Levenberg-Marquardt algorithm, thanks to the implementation provided by the GNU Scientific Library (Galassi et al., 2002). The first inversion step also allows us to study the possible trade-off between the inverse quality factors and the anisotropy parameter. Figure 5 illustrates step 1 of our inversion approach and the results obtained at stations STY and CHK, which display distinct energy distributions in the coda. A posteriori fits of the MLTWA measurements (dots) for the best anisotropic (black solid lines) and isotropic (gray dashed lines) models are shown on the left. For the simulations, the S-wave velocity is deduced from the experimental traveltime curves constructed at each station from the seismic bulletin provided by the CWT. Table 1 provides the best-fitting values of −1 , −1 and for these two stations. The misfit function, the effective scattering attenuation (1∕ * ) and the absorption ( −1 ) (obtained after step 1) are also plotted as a function of together with the typical uncertainties. Following our discussion on the transport mean free path, the effective scattering loss is defined as 1∕ * = −1 (1 − ) . This quantity is thought to be a better estimate of the scattering attenuation of direct waves than −1 since it takes into account the fact that the energy scattered around the forward direction is not entirely lost. From the uncertainties shown in Figure 5, it is clear that absorption −1 is better constrained than the scattering parameters −1 and . We also remark that −1 is strongly correlated with , in contrast with −1 .  Figures S2 to S14 in Supporting Information S1. Energy ratio measurements (dots), the best-fitting isotropic (dashed lines) and anisotropic (solid lines) models are shown on the left panels, where blue, green and red dots correspond to energy ratios in the first, second and third window, respectively. The station locations are indicated on the map furthest to the right: red lines are surface traces of active faults from Shyu et al. (2016) and the five major lithotectonic units are delineated with thick dark lines. The middle panels show the misfit as a function of the scattering anisotropy ( ), while the right panels indicate 1∕ (black solid lines) and 1∕ * (gray solid lines) for the best models as a function of the scattering anisotropy (g) with standard errors delineated by dashed lines. Station STY is characterized by a "classical" distribution of the seismic energy. The misfit function decreases with g in the range [−0.95, 0.0] and increases rather slowly for positive . We do not observe a significant difference between the a posteriori fits for anisotropic models with → 1 and a posteriori fit for the isotropic multiple-scattering models. In other words, all models with −0.4 can explain equally well the MLTWA measurements at STY. For g > −0.4, we observe on the one hand that −1 is almost constant, regardless of the level of anisotropy. On the other hand −1 increases with g indicating that strong forward-scattering is compatible with the data but a higher scattering level is required.
Unlike STY, station CHK shows an "anomalous" distribution of the seismic energy, characterized by the downward curvature of the normalized energy as a function of hypocentral distance near = 0 km. The corresponding misfit function increases with over the whole range [−0.95, 0.85]. A posteriori fits for models with strong back-scattering ( ≃ −1) are significantly improved, in particular at short hypocentral distances for the first time window. To summarize, an anomalous distribution of the energy in the coda is better explained by a multiple-scattering model with strong back-scattering ( −0.8) , while a "classical" energy distribution in the coda can be equally well explained by isotropic or forward-scattering models. In the latter case, we note an increase of scattering attenuation −1 with the anisotropy parameter . The equivalence of forward-scattering and isotropic scattering models in the case of a "normal" energy distribution was previously pointed out in the context of sedimentary basins by Gaebler et al. (2015).
Results at all available stations in the three frequency bands [1, 2], [2,4], and [4,8] Hz are summarized in Table 4. In Figure 6, we compare intrinsic and scattering attenuation ( −1 and * −1 ) for isotropic ( = 0) and anisotropic ( ≠ 0) models. We recall that = * for isotropic scattering ( = 0) . We also compute the misfit ratio, 2 iso ∕ 2 aniso , to compare quantitatively the goodness of the fit to the data for anisotropic and isotropic models.

of 23
At low frequency (1-4 Hz), the data fits of anisotropic models are systematically better ( iso∕ aniso > 1.5) when the optimal value of is strongly negative (typically −1 < < −0.8 ). In that case, the isotropic models systematically underestimate the absorption (by a factor ∼1.5) and the effective scattering attenuation (by a factor 2-3) compared to the best anisotropic model with negative . We also observe that at most of the stations where the best-fitting model has positive , the improvement of the fit to the data is modest compared to isotropic models. The scattering attenuation (1∕ ) found in isotropic models is generally close to the effective scattering attenuation (1∕ * ) deduced from anisotropic scattering models. This result agrees with the findings of Gaebler et al. (2015). We also find that the absorption retrieved from the data does not vary significantly between isotropic and anisotropic models with 0 .
At high frequency (f = 6 Hz), there is a significant reduction of the number of stations at which the inversion can be performed, due to the overall strong level of attenuation. As a consequence of the degraded data coverage, the inversion of the scattering anisotropy is more delicate. Nevertheless, it is worth noting in Figure 6 that models with positive g exhibit stronger effective attenuation than their isotropic counterparts (1∕ * (anisotropic) > 1∕ (isotropic)) . This behavior is also visible to a lesser extent at lower frequencies for models with 0 and weak attenuation. In these models, the transport mean free path can be of the same order or even larger than the typical hypocentral distance. In that case, energy transport is likely sensitive to details of the scattering that are not encapsulated in the parameter . Indeed, complete information on the scattering pattern requires the knowledge of all the moments of the cosine of the scattering angle. In summary, in the weak scattering regime, the effective scattering attenuation 1∕ * is not sufficient to describe the transport of seismic energy. Sufficiently detailed knowledge of the scattering pattern appears to be required. The inequivalence of isotropic and anisotropic models also applies to the strong backscattering regime. In the next section, the inversion results are further explored by analyzing the spatial variations of attenuation and scattering anisotropy across Taiwan. , that is, the relative contribution of scattering versus absorption, and scattering anisotropy ( ) Table 4) are assigned to the location of each station as previously done by Carcolé and Sato (2010). This simple mapping approach can be justified on the basis of the sensitivity analysis of Mayor   . These authors have demonstrated that the sensitivity of the intensity in the coda to attenuation parameters is greater at the location of the station and the source. As MLTWA is a single station approach, there is a redundancy of paths sampling the vicinity of the seismic station. Our results are therefore expected to be representative of the attenuation properties in the vicinity of each station. We will provide in this section a description of the most striking features observed on the various maps. In the case of absorption −1 and scattering * −1 , the color scale has been specifically adapted to better highlight the lateral variations: white color corresponds to the average attenuation for Taiwan (Table 2), while dark red and dark blue correspond to the minimum and maximum values obtained in each frequency band, respectively.

Scattering
Compared to values reported worldwide (Sato, 2019), the mean scattering level is globally high in Taiwan, particularly at low frequency (see Table 2). Furthermore, we note that the mean value of scattering attenuation * −1 decreases significantly with frequency, over about one order of magnitude between 1.5 and 6 Hz ( Table 2).
This strong frequency dependence is also observed on single-station estimates (Table 4). Figure 7 exhibits huge lateral variations of scattering, about 2 orders (resp. 1 order) of magnitude at low (resp. high) frequency. Despite the relatively small geographical area, scattering fluctuations in Taiwan are one order of magnitude greater than those observed in Japan or in Korea (Carcolé & Sato, 2010;Rachman et al., 2015). In Table 3, similar tectonic units were grouped together to compare their average level of seismic attenuation. In the 1-4 Hz frequency band, scattering is consistently strong around active faults in the CP and northern WF, as well as in alluvial deposits of the IP and Pingtung Plain (south of the WF). The LV and the CoR, composed of exotic volcanic arc rocks and associated mélange materials, are also characterized by strong scattering. In comparison, the HR, CR, and central WF exhibit relatively low scattering strength. In these regions, we find that the scattering is, respectively, 4 times and 3 times less intense in the 1-2 and 2-4 Hz frequency bands, as compared to neighboring tectonically active or volcanic regions (Table 3). In the highest frequency range (4-8 Hz), scattering is on the whole much weaker. The HR (except in the vicinity of the geological boundary fault on the eastern edge) and south-central WF are still marked by the lower values of scattering, while the eastern flank of the HR, CR and CoR are globally the most attenuating areas. Even if the coverage is sparse, active faults are consistently marked by rather low scattering at high frequency.

Absorption
Similar to scattering attenuation, the mean level of absorption is rather high in Taiwan and slightly decreases with frequency from −1 = 9.37 × 10 −3 at 1.5 Hz to −1 = 3.22 × 10 −3 at 6 Hz (Table 2). Figure 8 illustrates the rather strong lateral variations of absorption, with an amplitude of a few tens of % at a scale of a few tens of kilometers. Yet, it is worth pointing out that these variations are modest compared to those shown on the scattering attenuation maps. The spatial distribution of absorption depends on the frequency and the amplitude of the variations slightly decreases with frequency ( Table 2). As shown on the 1-2 Hz map, stations located in the LV and CoR affiliated with the accreted Luzon arc and forearc of the PSP affinity, composed mainly of andesitic volcanoclastics, flysch deposits, and ophiolite-bearing mélange, show higher absorption than those located in the CR and HR that constitute the metamorphic slate belt and schist complex (see also Table 3). Sedimentary materials located in highly tectonized area in the CP and WF, that is, in the vicinity of active faults, exhibit on average stronger absorption than sedimentary materials in the IP or near the Peikan Basement High. On the 2-4 Hz map, we observe grossly the same features except for the IP where absorption is now almost at the same level as in the CoR/LV. The absorption signature of active faults in the WF is also less clear. Even if the spatial coverage is much more sparse at high frequency (f > 4 Hz), the spatial distribution of high-frequency absorption does not seem to be controlled by the surface geology: we rather observe an east-west gradient. Figure 9 shows the relative contribution of scattering and absorption, also termed albedo

Albedo
, in the three frequency bands [1,2], [2, 4], and [4, 8] Hz. The albedo is an important quantity for seismic hazard assessment. For the same level of total attenuation, a small value of albedo implies that attenuation will reduce the amplitude of the whole signal. To the contrary, a large value of albedo is expected to promote the transfer of energy CALVET ET AL.

10.1029/2022JB025211
13 of 23 from direct to coda waves. In this case the reduction of the amplitude of the direct waves can be associated to an increase in the duration of the ground motion. On average, * 0 decreases as the frequency increases. We observe that scattering is universally dominant ( * 0 > 0.5 ) at 1-2 Hz in areas where scattering is strong around active faults distributed in the CP, WF, LV, CoR, and IP while scattering and absorption are equivalent in the HR and CR. At higher frequency (f > 4 Hz), absorption becomes globally dominant in the whole taiwanese crust (see also Tables 2 and 3).   Figure 10 shows the spatial variations of scattering anisotropy ( ) across Taiwan in the three frequency bands ([1, 2], [2, 4], and [4, 8] Hz). When −0.5 < < 0.5 , scattering is mostly isotropic and the corresponding values of g are shown in white. Red (resp. blue) colors correspond to predominantly backward (resp. forward)-scattering. The size of the colored disks depends on the ratio between the misfit function ( 2 ) for an isotropic and an anisotropic model. 2 ∕ 2 > 1.5 indicates that an anisotropic model significantly improves the fit to the data. As previously discussed, anisotropic models with strong back-scattering ( −0.8) better explain the data than isotropic ones in certain areas. At low frequency (1-2 Hz), g is globally negative in the crust of Taiwan (Table 2). The zones of strong back-scattering ( −0.8) are mainly located close to active faults in CP, WF and IP. Stations in LV and CoR are also globally characterized by negative values of g. By contrast, crystalline massifs in CR and HR exhibit higher values of scattering anisotropy ( −0.5) . On average slightly increases with frequency (see also Table 2), but the spatial distribution of is almost similar in the two first frequency bands. At higher frequency (f > 4 Hz), the spatial distribution of is not spatially coherent and does not show clear relations with tectonic structures or geological units.

Scattering Anisotropy
To summarize, we observe that volcanic materials in LV are characterized by strong scattering and absorption. This property of volcanos is broadly compatible with the literature (Carcolé & Sato, 2010;Del Pezzo et al., 2001; 10. Scattering anisotropy ( ) maps in three frequency bands: results of multiple lapse time window analysis are plotted at the location of the station. Models with strong backscattering (resp. strong forward scattering) are shown in red (resp. blue). Isotropic models are in white. The size of the circles correspond to the misfit ratio between isotropic and anisotropic models (χ iso /χ aniso ). Pink lines are surface traces of active faults from Shyu et al. (2016) and the five major lithotectonic units are delineated with thick dark lines (see Figure 1 for more details).

Comparison With Others Studies
Y.-J. Wang et al. (2010) have previously determined the attenuation structure of Taiwan, from t* measurements on P-and S-wave spectra. In their study, the authors hypothesize that seismic attenuation is (a) frequency independent, (b) mainly controlled by anelastic processes (effects of temperature and/or fluid saturation). In the first 10 km, they observe that central and northern CR are globally less attenuating than southern CR, CoR/ LV, and CP. Northern WF shows lower attenuation than HR while the southern WF (mainly the Pingtung Plain) and Southern CR are more attenuating than the surrounding regions. Y.-J. Wang et al. (2010) proposed that the strong attenuation anomaly in southern CR may be ascribed to the lower crust exhumation associated with an orogen-parallel extrusion. They also systematically observed a sharp contrast of the total attenuation across the Chelungpu, Kaoping and Chaochou faults. Our results are rather different. First, Figures 7-10 and Tables 2 and 3 clearly show the frequency dependence of attenuation: (a) the mean total attenuation varies over about one order of magnitude from 3.467 × 10 −2 at 1.5 Hz to 4.95 × 10 −3 at 6 Hz, (b) the frequency dependence is stronger for scattering than for absorption. Second, scattering is generally dominant in CP, WF, LV, and CoR (see Table 3), so that * -measurements cannot be directly interpreted as absorption only. Third, even if our spatial coverage is not perfect, most actives faults in the WF (from north to south) are characterized by both strong scattering and absorption at low frequency. In particular, the Chelungpu Fault, as well as the southern part of the Kaoping fault, are mainly marked by strong absorption and scattering (at low frequency) both in the hangingwall and the footwall of the fault. By contrast, the Chaochou Fault is rather characterized by low scattering and low absorption. Lastly, northern WF is systematically more attenuating than HR, and CR has overall both low scattering and low absorption compared to surrounding regions. Consequently, the absence of an attenuation anomaly in southern CR may not support the collisional tectonic model with a restricted subduction of the Eurasian crust discussed by Y.-J. Wang et al. (2010).
The seismic attenuation structure of Taiwan is clearly original compared to other regions of the world. In Figure 11, we have reported 1∕ * and 1∕ in the frequency band 1-2 Hz for two kinds of geological structures which correspond to both extremes in terms of seismic attenuation: volcanoes and old orogens. Volcanoes are characterized by strong scattering To make the comparison meaningful, we have distinguished between the main geological units of Taiwan as follows: crystalline massifs (HR and CR) are shown in blue, sedimentary materials (CP, IP, and WF) are in green, and areas affected by volcanism (mainly LV/CoR) are in yellow. Globally, the values of scattering and absorption in Taiwan are very close to the one of volcanoes. Sediments in CP and IP, highly tectonized areas in WF, and volcanic materials in LV have a similar level of scattering as volcanoes while absorption is close to the higher bound of −1 for volcanoes. By contrast, taiwanese crystalline massifs (CR and HR) are more attenuating by one order of magnitude (both in scattering and absorption) than old orogens. In Figures 11  and 12, we also compare taiwanese attenuation data to other regions (Japan, North America and Korea).  As Taiwan and Japan are both in a subduction environment, one might expect similar levels of attenuation. Figures 11 and 12 shows that: (a) in the frequency band [1-8] Hz, absorption is systematically stronger in Taiwan than in Japan (at least by a factor 2 at low frequency f < 2 Hz), North America (in the whole frequency range) and Korea (except in the first frequency band); (b) scattering attenuation in Taiwan is one order of magnitude higher than in Japan, North America and Korea. Carcolé and Sato (2010) have also applied a MLTWA approach to derive absorption and scattering maps, but they have used an isotropic scattering model (as well as Eulenfeld and Wegler (2017) for North America and Rachman et al. (2015) for Korea). As discussed in the previous section, an isotropic model may systematically underestimate attenuation (both absorption and scattering) compared to an anisotropic model with strong back-scattering. Consequently, distinct hypothesis on scattering anisotropy may partially explain some of the differences between Taiwan and Japan (or North America and Korea). Of course, this argument holds only in the case where the spatio-temporal distribution of energy is anomalous, an hypothesis for which we have no direct confirmation so far in the case of Japan. It is also noteworthy in Figure 11 that the spatial variations of attenuation in Taiwan are much larger than in Japan which covers a larger area -compared to Taiwan-but with a comparable diversity of geological structures.

On the Distribution of Fluids in the Taiwanese Crust
Huang et al. (2014) estimated the velocity structure ( and ) of the taiwanese crust. They obtained a reliable ∕ model constrained at the surface (from 30-60-500 m depth) using both seismic data and logging data from the Engineering Geological Database for TSMIP (EGDT). From the surface to 10 km depth, the main plains and basins (Taipei Bassin, Western Plain and Pingtung Plain), the LV, and the highly tectonically active regions in WF and CP show low and and high ∕ ratios ( Figure 13). These regions also correlate with areas of strong absorption. High ∕ ratio, low velocities and strong absorption are usually interpreted in terms of fluid content as many mechanisms of intrinsic attenuation are related to the presence of fluids within cracks and/or pores in crustal rocks (see Mavko et al. (2020) for a review). The Pingtung Plain, characterized by low velocities, high ∕ and strong attenuation, is associated with a rapid deposition environment and accumulates unconsolidated coastal and estuarine sediments over more than 5 km thickness. It is reasonable to assume that such  (Lacombe et al., 2003;Sens-Schönfelder & Wegler, 2006;Sens-Schönfelder et al., 2009) and volcanoes by red diamonds (Del Pezzo et al., 2001;Mayeda et al., 1992;Ugalde et al., 2010;Wegler, 2003). Minimum and maximum values for absorption and scattering in Japan (Carcolé & Sato, 2010), Korea (Rachman et al., 2015) and North America (Eulenfeld & Wegler, 2017) are also shown.   (Table 3) are indicates by thick lines. Results for Japan (Carcolé & Sato, 2010), Korea (Rachman et al., 2015) and North America (Eulenfeld & Wegler, 2017) are also shown.
CALVET ET AL.
10.1029/2022JB025211 20 of 23 materials contain a significant amount of fluids. Similarly, high ∕ , low velocities and strong attenuation in LV may also reflect the presence of fluids or melts related to the volcanic activity. These areas are also characterized by strong scattering (Figure 10) but more importantly by strong back-scattering ( −0.5) at low frequency additional information on the nature of the crustal heterogeneities. Takemura and Yoshimoto (2014) previously demonstrate a strongly heterogeneous structure in the Kanto area (Japan) interpreted as the result of the dehydration of the subducting oceanic crust of the PSP. They assume that these strong small-scale velocity heterogeneities may be related to a random distribution of fluid which may induce low propagation velocities, high -ratio, strong attenuation and generate wide-angle scattering. Wu and Aki (1985) demonstrate that an elastic perturbation may be decomposed into a sum of "velocity-type" and "impedance-type" perturbations. A "velocity-type" heterogeneity corresponds to an object marked by a contrast of velocity without a contrast of impedance with the environments and vice-versa for an "impedance-type" heterogeneity. In the low frequency regime, "impedance-type" perturbations scatter seismic energy in the backward direction producing a negative . By contrast, "velocity-type" perturbations scatter seismic energy in the forward direction resulting in a positive . Fluid inclusions correspond to strong "impedance-type" perturbations. As a consequence, negative-values of at low frequency, combined with other observations such as a low velocities, high ∕ ratio and high absorption, could also mark the presence of fluids in LV, CP, WF and in the main sedimentary basins. It is also worth noting that g globally increases with frequency in these areas. This observation is compatible with the hypothesis of the presence of fluids since scattering theory predicts g that becomes positive at higher frequency even for "impedance-type" heterogeneities.

Conclusions
Using a MLTWA approach and a multiple anisotropic scattering model, we have produced the first crustal attenuation maps in Taiwan that distinguishes between scattering and anelastic attenuation. Despite a very limited geographical extent, the Taiwanese crust shows a strong spatial variability of seismic attenuation. Across the island, scattering varies over more than one order of magnitude and absorption fluctuations are of the order of a few tens of percent. The spatial distribution and the amplitude of the two attenuation mechanisms vary very strongly with frequency, but scattering is overall dominant at low frequencies (f < 2 Hz) while anelasticity dominates at higher frequencies. Our results confirm that attenuation maps previously obtained from * measurements cannot be directly interpreted as absorption only. The frequency dependence in scattering may be interpreted in terms of the scale and nature of the heterogeneities (volcanic/magmatic or clastic materials, melts, etc…), and the frequency dependence of absorption probably depends on the composition and distribution of melts or fluids. In this study, we demonstrate that scattering anisotropy -quantified by the mean cosine of the scattering angle -is needed to properly describe the energy distribution in the S-wave coda of local earthquakes at stations located in sedimentary, tectonized and volcanic areas. also strongly varies with frequency. At low frequency, we observe that strong back-scattering ( −0.5) is generally correlated with strong scattering, strong absorption, low velocities and high ∕ ratio. As back-scattering is observed only for impedance fluctuations, it suggests that the parameter could be used as a marker of the presence of fluids. Hence, our maps possibly highlight the spatial distribution of fluids in the Taiwanese crust in particular in the vicinity of active faults.
Our study also demonstrates that the widely-used multiple isotropic scattering model may in some cases be too simplistic. In particular, neglecting the anisotropy of scattering leads to an underestimation of the attenuation (both in absorption and in scattering). Even if scattering anisotropy significantly improves the agreement between the multiple-scattering model and observations in Taiwan, this anisotropic model is probably still too simplistic given the complexity of the Taiwanese crust. First, the vertical stratification of the velocity and the depth variation Figure 13. Comparison between the scattering anisotropy at 1.5 Hz (colored circles) and the velocity ratio ∕ at the surface (background color on the right) (Huang et al., 2014). Models with strong backscattering (resp. strong forward scattering) are shown in red (resp. blue). Isotropic models are in white. The size of the circles correspond to the misfit ratio between isotropic and anisotropic models (χ iso /χ aniso ). Pink lines are surface traces of active faults from Shyu et al. (2016) and the five major lithotectonic units are delineated with thick dark lines (see Figure 1 for more details).

of 23
of the Moho as a function of the station location have not been considered in this work. Margerin et al. (1998); Margerin (2017) have demonstrated that the partial confinement of seismic energy in a crustal waveguide gives rise to energy leakage which mimics the effect of absorption and potentially biases the estimation of anelasticity. Second, our MLTWA approach only gives access to the lateral variations of attenuation without constraint on the depth structure. In the spirit of Ogiso (2019), future works should explore in more details the sensitivity of the energy ratio in the coda to the vertical stratification and lateral variations of both scattering loss and absorption.