Amplitude‐Dependent Energy Dissipation Inferred From Spectral Analyses of Intraslab Earthquakes

While the Earth (rocky planet) can be approximated as an elastic body, it also possesses anelasticity, Q−1, which causes energy dissipation (intrinsic attenuation) during deformation. Rock deformation experiments have demonstrated that Q−1 shows a marked amplitude dependency for strain of >∼10−6, but seismic‐waves‐induced strain is generally much smaller than the threshold value. Therefore, all seismological analyses have assumed amplitude‐independent attenuation. However, theoretical models of dislocation‐induced attenuation have predicted amplitude‐dependent attenuation under high temperatures and pressures equivalent to the crustal and uppermost mantle conditions. This study investigates whether attenuation is actually independent of seismic‐wave amplitudes by a systematic analysis of a large number of spectral amplitudes of co‐located earthquake pairs with different observed amplitudes. We assume the ω2 source and circular crack models to correct for the source parameters of earthquakes. The obtained results suggest that amplitude‐dependent attenuation is required to explain the observations, which contradicts a long‐standing recognition that amplitude‐dependent attenuation is negligible for the propagation of seismic waves. Our model calculation further suggests that Q−1 is proportional to the seismic amplitude, A, whereby Q−1 ∝ A0.1–0.2, which introduces attenuation variations along a given raypath, with significantly enhanced attenuation near the hypocenter. Our result highlights the need to revisit and potentially rethink current models of the physical mechanisms influencing intrinsic attenuation.


Introduction
The Earth is generally approximated as an elastic body, such that the theory of elasticity is used to analyze seismic-wave propagation, fault slip, surface deformation, and free oscillations.However, the Earth possesses anelasticity, which causes energy dissipation (intrinsic attenuation) during seismic-wave propagation and internal deformation.Intrinsic attenuation, Q − 1 , which is defined as the ratio of the energy loss during one cycle to the total energy, has been estimated from geophysical observations, including the amplitude decay of body and surface waves (Mitchell, 1995), the broadening of free oscillation peaks (Benioff et al., 1961), and the damping of the Chandler wobble (Smith & Dahlen, 1981).
The experimental studies have revealed frequency-dependent attenuation with Q − 1 ∝ f − α at temperature conditions comparable to the upper mantle, where f is frequency and α is estimated to be 0.2-0.4(Farla et al., 2012;Gribb & Cooper, 1998;I. Jackson et al., 1992I. Jackson et al., , 2002;;Takei et al., 2014;Yamauchi & Takei, 2016).The frequencydependent attenuation can be supported by seismological observations (e.g., Lekić et al., 2009;Nakajima et al., 2013;Stachnik et al., 2004).Furthermore, the measurements of attenuation in a wide range of strain amplitudes have suggested that attenuation is constant for small strain amplitudes and is potentially enhanced with increasing amplitude (Gordon & Davis, 1968;Mavko, 1979;Stewart et al., 1983;Takei et al., 2014;Tisato & Quintal, 2014;Winkler et al., 1979;Yamauchi & Takei, 2016).However, amplitude-dependent attenuation has not yet been confirmed from seismological observations.When amplitude-dependent attenuation is considered, we can approximate attenuation as follows (e.g., Mavko, 1979): where Q − 1 L represents amplitude-independent attenuation, ε is strain amplitude, and a 1 is a positive constant.Equation 1 means that Q − 1 essentially increases proportional to the strain amplitude.With typical values of Q − 1 L and a 1 obtained in the experimental studies (Mavko, 1979;Stewart et al., 1983), strain amplitudes above which attenuation possesses amplitude dependent are estimated to be ∼10 − 6 .However, this threshold value is too large compared to seismic-wave amplitudes, because strain amplitudes of seismic waves that we observe at hypocentral distances of ∼100 km are expected to be as small as 10 − 10 -10 − 8 for M3-5 earthquakes and 10 − 7 -10 − 6 even for larger earthquakes (Figure 1) (Gomberg et al., 2001;Hill et al., 1993).In addition, the nonlinearity of attenuation with respect to strain amplitudes may disappear at pressures of >5 MPa (Winkler et al., 1979).Therefore, it is undoubtedly assumed that seismic attenuation does not depend on amplitude for seismic waves propagating in the Earth, and such amplitude-independent attenuation have been involved in major terrestrial observations, interpretations, theoretical models, and hypotheses.
Theoretical anelastic models have shown that attenuation due to dislocations would depend on amplitudes at high temperature and pressure conditions equivalent to the upper mantle (Granato & Lücke, 1956;D. D. Jackson & Anderson, 1970;I. Jackson, 2015;Karato, 1998;Karato & Spetzler, 1990).However, experimental studies that investigated dislocation-induced attenuation are limited and have not quantified amplitude-dependent attenuation  1. Dashed line represents an approximate threshold value, above which attenuation becomes nonlinear.We used a relationship A = εV p (Van Der Elst & Brodsky, 2010) with V p = 7.0 km/s for the conversion from strain amplitude, ε, in rock experiments, to seismic-wave amplitude, A. Gray area represents typical seismic-wave amplitudes observed for small to moderate-sized earthquakes.(Farla et al., 2012;Gueguen et al., 1989;Sasaki et al., 2019).Therefore, whether seismic attenuation depends on amplitude at the high pressure and temperature conditions and how it behaves during the propagation of seismic waves are still enigma.
This study investigates amplitude dependence of seismic attenuation by the systematic analysis of seismic waveforms and challenges a long-standing recognition that seismic attenuation is constant during seismic-wave propagation in the Earth.We apply the spectral ratio method to spectral amplitudes observed from co-located but different-sized earthquakes and quantify amplitude dependence of attenuation from the changes in the slope of the spectral ratio.The obtained results suggest that constant attenuation cannot explain the observations and seismic attenuation is required to depend on amplitude with the power of 0.1-0.2, which is different from the conventional assumption of amplitude-independent attenuation.

Methodology
This study employs the spectral ratio method (e.g., Abercrombie et al., 2017;Mayeda et al., 2007) to investigate the amplitude dependence of seismic attenuation, whereby we analyze a pair of co-located larger (EQ-L) and smaller (EQ-S) earthquakes that have almost identical ray paths to a given seismic station (Figure 2a).If seismic attenuation is amplitude independent, then the spectral ratio of EQ-L to EQ-S would flatten above the corner frequencies ( f c ) of the two earthquakes (black dashed line in Figure 2b).Conversely, if attenuation is amplitude dependent, with higher attenuation attributed to larger amplitudes and vice versa, then the slope of the high-frequency decay in the spectra becomes steeper for EQ-L (red line in Figure 2b), and the resultant spectral ratio yields a negative slope at higher frequencies (black line in Figure 2b).The spectral ratio of larger and smaller co-located earthquakes can therefore be used to resolve the amplitude dependence of seismic attenuation if it exists, because the slope of the spectral ratio becomes steeper as the amplitude dependence of attenuation strengthens.

Theoretical Equations
The velocity amplitude spectrum (A f ) of earthquake i at station m can be expressed as follows (Scherbaum, 1990): Concept for the estimation of amplitude-dependent attenuation.(a) Schematic ray paths to a given seismic station (inverted triangle) for two co-located earthquakes.Red and blue stars represent the locations of the larger (EQ-L) and smaller (EQ-S) earthquakes, respectively.The distance between the two earthquakes is small enough to assume that the ray paths from the two earthquakes are almost identical.(b) Theoretical spectra for EQ-L (red line) and EQ-S (blue line) based on the ω 2 source model (Brune, 1970).The dashed lines represent the source spectra with constant attenuation (Q = 300 for EQ-L and EQ-S), while the solid red line represents the source spectrum with Q = 240 for EQ-L (amplitude-dependent attenuation).The EQ-L/EQ-S spectral ratios with and without amplitude dependent attenuation are shown by the solid and dashed black lines, respectively.f cL and f cS represent the theoretical corner frequencies of the EQ-L and EQ-S spectra, respectively.
where f is the frequency, S is the source spectrum, I is the instrumental response, P is the site response, B is the attenuation term, and F is a frequency-independent constant term that includes the radiation pattern.The ratio of the velocity amplitude spectra of two co-located earthquakes (EQ-L and EQ-S) at station m can be expressed as follows: Here, I and P can be canceled out since the spectral ratios are taken at the same station.The attenuation term B is expressed as follows: where α is an exponent that represents the frequency dependence of Q − 1 , Q − 1 1 im is seismic attenuation at 1 Hz along the raypath from earthquake i to station m, and t im is the travel time from earthquake i to station m.
It is usually considered that the attenuation term B can be canceled out for a pair of co-located earthquakes because the raypath from the hypocenter to the station is identical for the two earthquakes.However, if seismic attenuation exhibits amplitude dependent, then B cannot be canceled out.For amplitude-dependent attenuation, we obtain the following equation by substituting Equation 4 into Equation 3: where f c is the corner frequency and t is the travel time.Here we assume the ω 2 source model (Brune, 1970).
Examples of the spectral ratios for any pairs of earthquakes obtained via Equation 5 are shown in Figure 3a.
We calculated the theoretical f c using a circular crack model (Eshelby, 1957) with a given stress drop.We then moved the source spectrum to the left-hand side (Figure 3b) and took the logarithm of both sides (Figure 3c).We assumed t = (t L + t S )/2 because we chose co-located earthquake pairs, and divided both sides by the travel time to obtain the following equation (Figure 3d): where: and a 2 is a frequency-independent constant.The above transformation linearizes Equation 5 with respect to the differential attenuation, ∆Q − 1 m .Equation 6indicates that the slope of the spectral ratio becomes zero for amplitude independent attenuation (∆Q − 1 m = 0).
While Equation 6 can be constructed for every earthquake pair at a single station, the spectral ratio would often be unstable due to noises in waveforms.Therefore, we classified the spectral ratios that have the same range of observed time-domain amplitude ratios, A t Lm / A t Sm , to estimate ∆Q − 1 for that amplitude ratio.The amplituderatio bands were set at 0.2 intervals over the log(A t Lm / A t Sm ) range, with a moving window of 0.1.After classifying each observed spectral ratio into its corresponding amplitude-ratio bands, we corrected for the offset of each spectral ratio (Figure 3e), calculated the average spectral ratio (Figure 3f), and estimated ∆Q − 1 for a given amplitude-ratio band via a least squares fit of Equation 6 to the average spectral ratio (Figure 3g).
Since it is well-known that Q − 1 in the upper mantle depends on frequency with α = 0.2-0.4(e.g., I. Jackson et al., 1992Jackson et al., , 2002;;Lekić et al., 2009;Nakajima et al., 2013), we adopted α = 0.3.The effect of α on the results is discussed in Sections 4-5.We used a frequency range of 20-40 Hz for the spectral ratio fitting, which can avoid the effect of f c on the estimate of ∆Q − 1 (see Section 2.2).

Corner Frequency Estimation
We used theoretical values for f c since this study focuses on the estimation of ∆Q − 1 .The fault radius, r, in a circular crack model can be expressed as follows (Eshelby, 1957): where M o is the seismic moment and ∆σ is the stress drop.M o is calculated from the moment magnitude, M w , using the following transformation equation (Kanamori, 1977): We used the Japan Meteorological Agency (JMA) magnitude (M JMA ) instead of M w in Equation 8 since M w is comparable to M JMA over the M3-5 range used in this study (Edwards & Rietbrock, 2009;Uchide & Imanishi, 2018).f c was calculated using the formula of Madariaga (1976): where k = 0.315 (constant) and β 0 is the shear-wave velocity near the hypocenter.We adopted the appropriate β 0 value for the focal depth of each hypocenter from the JMA 2001 velocity model (Ueno et al., 2002).
f c was finally obtained by substituting Equation 7into Equation 9 and rearranging the formulation: Since the average ∆σ of earthquakes that occurred in our study area is estimated to be ∼10 MPa (Tsuda et al., 2006), we assumed ∆σ = 10 MPa for all earthquakes to calculate f c .This assumption yields f c = 11.2Hz for an M3 earthquake, the smallest earthquake used in this study.We will discuss the validity of the assumed constant stress drop model in Section 4.3.

Data
We analyzed 457 aftershocks (M ≥ 3) of the 2003 M JMA 7.1 Miyagi-oki intraslab earthquake that occurred on 26 May in the subducting Pacific slab beneath Tohoku, Japan.Aftershocks used in this study are distributed at 60-80 km depth and occurred during the period between 27 May 2003 and 31 December 2007 (red dots in Figures 4a  and 4b).Detailed information of the origin time, latitude, longitude, depth, and magnitude of the earthquakes used in the analysis is listed in a supplementary file.
We calculated the spectral amplitudes using the P-wave signals that were recorded at 100 Hz sampling, and took the spectral ratio of each earthquake pair with an inter-event distance of ≤5 km that was recorded at a station located within epicentral distances of 200 km.We only kept the spectra that possessed signal-to-noise ratios of ≥5 over at least 90% of the 20-40 Hz band.We limited the observed EQ-S amplitudes (A t S ) to 10 − 6 -10 − 5 m/s to resolve the effect of the absolute amplitude on seismic attenuation (see Figure S1 in Supporting Information S1 for the observed amplitude distribution).
We conducted the above analyses for each of the four overlapping 1-s time widows (50% overlap), with time zero of the first window set at the P-wave onset (Figure 4c), such that we could discriminate amplitude dependence of attenuation with respect to the increasing proportion of P-coda waves (Sato, 2007;Sato & Korn, 2007).We took earthquake pairs without considering their focal mechanism similarity because our analysis exactly involved the absolute amplitudes observed at stations to classify the spectral ratios into an appropriate amplitude-ratio band.

Amplitude Dependence of Seismic Attenuation
We observed a marked positive correlation between the observed amplitude ratio and ∆Q − 1 for the four time windows over the amplitude ratios of 0.0 ≤ log(A t L / A t S ) ≤ 1.6 (Figures 4d and 5a).This correlation can be verified by the increase in the slope of the spectral ratio as the observed amplitude ratio increases (Figure 4d and Figure S2 in Supporting Information S1).Of note, the observed amplitude dependency of ∆Q − 1 tends to intersect the origin of the coordinate, which indicates that ∆Q − 1 becomes zero for two earthquakes with the same observed amplitudes.The positive correlation between the observed amplitude ratios and ∆Q − 1 is directly related to the increase in the EQ-L amplitude (A t L ), as A t S used in the analysis is approximately constant (10 − 6 -10 − 5 m/s) over the range of analyzed amplitude ratios (Figure 5b and Figure S1 in Supporting Information S1).The observed ∆Q − 1 values are in the range of 0-0.002, and these values can be produced by realistic attenuation values (i.e., Q − 1 L = 0.004 and Q − 1 S = 0.003 for ∆Q − 1 of 0.001).
Interestingly, ∆Q − 1 generally becomes insensitive to the amplitude ratio for the later time windows (Figure 4a).Since the absolute amplitudes in each amplitude band are almost the same even for different time windows (Figure 5b), it is likely that a larger proportion of P-coda waves weakens the amplitude dependence of ∆Q − 1 .This is supported by previous studies that coda waves generally become dominant 1− 2 s after the P-wave onset at a hypocentral distance of ∼100 km (Sato, 2007;Sato & Korn, 2007).Our results suggest that the amplitude dependence is dominant for intrinsic attenuation while it is less dominant for scattering attenuation.
Seismic stations used in this study are mostly located in the fore-arc side (blue squares in Figure 4a), where temperatures in the crust and uppermost mantle are not high enough to produce partial melting (Nakajima & Hasegawa, 2003).Therefore, we infer that amplitude-dependent attenuation arises from the seismic-waveinduced deformation of solid-state materials.It is noted that amplitude-dependent attenuation discussed in this study represents an integrated value along ray paths, which propagate in the uppermost mantle by 40-70 km and continental crust by ∼30 km (Figure 4b).Therefore, this study cannot resolve which part along a raypath exhibits amplitude-dependent attenuation.We need to carry out the same analysis using earthquakes in the continental crust to discriminate amplitudedependent attenuation in the crust from that in the upper mantle.

Effect of Source Spectra
This study used the ω 2 source model (Brune, 1970) to generate the source spectra for both EQ-L and EQ-S.However, a deviation from the ω 2 source model would produce an apparent inclination of the spectral ratio at frequencies above f c .Here we assumed that attenuation is amplitude independent, but the high-frequency decay of the source spectrum is different for EQ-L and EQ-S.This assumption yields the following spectral ratio of the two earthquakes: where n S and n L are the powers of the source spectra for EQ-S and EQ-L, respectively, and a 3 is a frequency-independent constant.We assumed the standard ω 2 source model for EQ-S with n S = 2 and evaluated the effect of n L on the high-frequency inclination of the spectral ratio.As described before, we used theoretical f c obtained from Equation 10.The right hand of Equation 11 was then fitted with the observed spectral ratios of each amplituderatio band, and the n L and a 3 values were determined via a grid search to minimize the residuals between the observed and theoretical spectral ratios.
The results show that the observed slope of the spectral ratio can be explained by a monotonic increase in the power of the source spectra of EQ-L with increasing to n L ∼ 2.6 for an M ∼ 4.5 earthquake (Figure 6 for the n L variations and Figure S3 in Supporting Information S1 for the details of the data fitting).Although it has been reported that the power of the source spectra can be variable (Abercrombie, 1995;Allmann & Shearer, 2009), no studies have previously indicated systematic and significant increases in the power of the source spectra with magnitude.Furthermore, if the differences in the power of the source spectra of the earthquake pairs are the major cause of the observed slope of the spectral ratio, we would obtain identical results for the four different time windows because the effect of the source spectra is equally included in the spectral ratio of each time window.However, we obtained the different trends of ∆Q − 1 for the four time windows (Figure 5a), whereby we conclude that possible uncertainties in the power of the source spectra have minor effect on the estimation of ∆Q − 1 .

Effect of Stress Drop Variation
This study used f c values that are theoretically calculated using a circular crack model with a constant stress drop (∆σ = 10 MPa) to correct for the observed spectral ratios.However, the stress drops of earthquakes in the subducting Pacific slab can vary by an order of magnitude or more (e.g., Kita & Katsumata, 2015;Uchide et al., 2014).We here assumed four models of normally distributed ∆σ with different average values (log(∆σ ave ) (Pa) = 6.0, 7.0, 7.5, and 8.0) and a standard deviation of 0.5 (Figure 7 for the model and results, and Figure S4 in Supporting Information S1 for the spectral-ratio fitting) and performed the same analysis to evaluate the effect of the uncertainties in f c on the spectral ratio.
The obtained results show that amplitude dependence of ∆Q − 1 is observed for each of the four stress drop values, with the amplitude dependence becoming weaker with increasing ∆σ ave .The reduced amplitude dependence for large ∆σ is due to the existence of f c for EQ-S in the 20-40 Hz frequency band  and win3), where the colors in each square denote the average absolute amplitudes for EQ-L (upper half) and EQ-S (lower half) at that amplitude band.
Figure 6.The best-fit power (n L ) of the source spectra for EQ-L to explain the observed spectral ratio without considering amplitude-dependent attenuation, where n S is fixed to be 2. Color coding is based on the average magnitude of EQ-L and EQ-S at each observed amplitude ratio.
( f c = 24 Hz for an M3 earthquake with ∆σ = 100 MPa), which would introduce an apparent negative slope of the spectral ratio.Although we often observe variable ∆σ, the average ∆σ for the aftershocks of the 2003 Miyagioki earthquake is estimated to ∼10 MPa (Tsuda et al., 2006), such that ∆σ ave = 100 MPa would be unlikely high.Therefore, we infer that the amplitude dependence of ∆Q − 1 observed in this study is robust with respect to moderate variations in the observed ∆σ.

Effect of Site Response
It is known that there are nonlinear site responses when the input seismicwave amplitudes are large (Assimaki et al., 2008;Bonilla et al., 2011;Sawazaki et al., 2006).If site responses show nonlinearity, we may observe an apparent slope in the spectral ratio.However, we could rule out the effect of these nonlinear site responses because we obtained different amplitude dependences of ∆Q − 1 for the four analyzed time windows, where almost identical seismic-wave amplitudes were used in the four time windows (Figure 5b and Figure S2 in Supporting Information S1) and the effect of nonlinear site responses would appear equally in each time window.Therefore, we conclude that the observed amplitude dependence of ∆Q − 1 does not arise from differences in the source parameters of the analyzed earthquake pairs or the nonlinearity of the site responses, but rather represents the actual amplitude dependence of intrinsic attenuation.

Effect of α
In this study, we utilized α = 0.3.However, previous seismic observations reported a wide range of α from 0.1 to 0.7 in the crust and the upper mantle (e.g., Nakajima et al., 2013;Stachnik et al., 2004;Takahashi, 2012;Thouvenot, 1983).We therefore evaluated whether the use of the other α values still result in the amplitude dependence of Q − 1 .The calculation results demonstrated that the amplitude dependence of Q − 1 is robust, even though the absolute ∆Q − 1 values vary with α (Figure S5 in Supporting Information S1).However, ∆Q − 1 becomes very large (∼0.01), for example, for α = 0.6 with 1.4 ≤ log(A t L / A t S ) ≤ 1.6, and the ∆Q − 1 of 0.01 requires attenuation pairs, such as (Q − 1 L , Q − 1 S ) of (0.01, 0) and (0.02, 0.01).These attenuation values are too large and cannot be supported by the seismic observations (Liu et al., 2014;Nakajima et al., 2013;Wang et al., 2017).In addition, such large amplitude dependence of ∆Q − 1 could probably be detected when we analyze spectra of individual seismic waveforms with large amplitude variations.Therefore, we infer that value of α is likely less than 0.5, even though the optimal value of α cannot be constrained only from this study.

Model Description
When intrinsic attenuation exhibits amplitude dependent, attenuation gradually changes along the raypath, depending on the amplitude of seismic waves at each point.We here formulate the amplitude dependence of attenuation as a function of seismic-wave amplitudes and quantify how seismic attenuation varies with amplitude.The frequency of seismic wave considered here is 1 Hz.
We assume that the P-wave signal radiating from hypocenter to seismic station only decays via geometrical spreading and intrinsic attenuation.The amplitude (A d ) of seismic waves at hypocentral distance d (km) can be reproduced using the initial amplitude (A 0 ) that radiated from the hypocenter (d = 0) as follows: where V p is the P-wave velocity (km/s) and Q is the quality factor at 1 Hz.
Microscopic rock deformation is often observed as a power-law creep, whereby the strain rate is proportional to the nth power of the stress (Bai et al., 1991;Mei & Kohlstedt, 2000;Sasaki et al., 2019).Here we assume that intrinsic attenuation is proportional to the nth power of the seismic amplitude, where intrinsic attenuation at distance d can be expressed as a function of A d as follows: where C is a phenomenological constant that controls the absolute value of attenuation.Once A 0 is given, we can calculate A d and Q − 1 d iteratively via Equations 12 and 13 at each distance d (0 ≤ d ≤ D) along a given raypath (Figure 8), where D is the hypocentral distance.It is noted that we assume that the crust and uppermost mantle possess the same characteristics of amplitudedependent attenuation.

Parameter Search
We first considered the ranges of n = 0.05-0.50(0.05 interval) and C = 0.001-0.300(0.001 interval).We then searched for the optimal initial amplitudes of EQ-L and EQ-S at a hypocenter, A 0 L and A 0 S , respectively, for each n-C pair to reproduce the average amplitudes of EQ-L and EQ-S observed at a station (colors in Figure 5b).We set D of 87 km, which is calculated from the average focal depth of earthquakes (71 km) and average epicentral distance to stations (50 km) analyzed in this study (Figure 8).V p at each point was calculated by the linear interpolation of the one-dimensional velocity structure of the JMA2001 model (Ueno et al., 2002).This procedure repeated for each amplitude-ratio band of the spectral ratio, which is set at 0.2 intervals of log(A t L / A t S ) range with a moving window of 0.1.
After we obtain the optimal A 0 L and A 0 S values for individual amplituderatio bands, we calculated Q − 1 d for each of EQ-L and EQ-S via Equations 12 and 13 at a 1-km interval along a raypath.Then we averaged Q − 1 L and Q − 1 S values for the entire raypath and calculated ∆Q − 1 from the path-averaged We finally calculated the residuals between the observed and theoretical ∆Q − 1 values over the range of 0.0 ≤ log(A t L / A t S ) ≤ 1.6 and determined the optimal pair of n and C values that best account for observed trend of amplitude dependence of attenuation.We applied this procedure only to the first time window (win1) (see Figure 4c) because later time widows contain a large proportion of coda waves but we did not involve the scattering effect in the calculation of seismic-wave amplitude in Equations 12 and 13.

Best-Fit Parameters
The amplitude dependence of ∆Q − 1 can be reproduced in which the optimal values exist in 0 ≤ C ≤ 0.1 and 0.05 ≤ n ≤ 0.4 (Figure 9a), with smaller n values (<0.2) yielding ∆Q − 1 trends that are more comparable to the observations (Figure 9b and Figure S6 in Supporting Information S1 for results with other α values).However, the theoretically calculated absolute Q − 1 values for n = 0.05 become larger than 0.01 over the amplitude range comparable to seismic amplitudes (10 − 6 ≤ A(m/s) ≤ 10 − 2 ) (purple line in Figure 9c), and these Q − 1 values are much larger than the values (0.001-0.007) observed beneath the fore-arc side of northeastern Japan (Liu et al., 2014;Nakajima et al., 2013;Wang et al., 2017).Therefore, we infer that Q − 1 exhibits an amplitude dependence with n = 0.1-0.2, which yields Q − 1 values comparable to the observations.We performed the same calculations for other values of stress drop and found that the optimal n-C values both decrease when we assume large stress drop values (Figure S7 in Supporting Information S1).Red and blue stars represent the hypocenter of EQ-L and EQ-S, respectively.We adopted the hypocenter depth to be 71 km and the epicentral distance to station to be 50 km, both of which correspond to the average values for earthquakestation pairs used in this study.The hypocentral distance is 87 km.A raypath is calculated with the 1D velocity structure of JMA2001.Colors of the line denote P-wave velocity at each depth.Note that the hypocenters of EQ-L and EQ-S are separated by 3 km, which is the average distance of pairwise earthquakes used in this study, to show how co-located earthquakes are actually distant, but this separation is not included in the model calculation.
Journal of Geophysical Research: Solid Earth 10.1029/2023JB028492 TERO AND NAKAJIMA Figure 9d shows an example of the inferred variations in seismic-wave amplitudes and Q − 1 along a raypath for two typical amplitudes (10 − 4.5 and 10 − 5.5 m/s) observed at station when n = 0.1 and C = 0.017.We assumed the hypocentral distance of 87 km and involved the effect of geometrical spreading and amplitude-dependent attenuation.Our calculation suggests that Q − 1 varies by ∼60% along the raypath, depending on variations in seismic-wave amplitudes.Figure 10a shows attenuation variation along a 100-km long raypath, which is normalized by attenuation at station.Attenuation at hypocenter is ∼1.6 and ∼2.6 times larger than that at station for n = 0.1 and n = 0.2, respectively.It is noted that normalized attenuation is always unity along a raypath if attenuation does not depend on seismic-wave amplitudes.
Figure 10b summarizes the maximum and average normalized attenuation along ray paths for four hypocentral distances of 50, 100, 150, and 200 km.The maximum attenuation exhibits a marked positive correlation with the hypocentral distance and becomes ∼1.8 and ∼3 times larger than that at station for n = 0.1 and n = 0.2, respectively.This indicates that the assumption of amplitude-independent attenuation significantly underestimates intrinsic attenuation as it approaches to hypocenter.In contrast, the path-average normalized attenuation does not largely change for the four hypocentral distances but is enhanced by 10%-30% compared to attenuation with no amplitude dependency.The enhanced path-average attenuation could be involved in precise estimations of source parameters such as corner frequency of earthquakes and in the calculation of strong ground motion for disaster mitigation.c) Amplitude dependence of Q − 1 at 1 Hz, which is calculated via Q − 1 (f),(A) = CA n for the five highlighted n-C pairs in panel (a).Dark gray bands denote the amplitude ranges of seismic waves that were not included in this study, and light gray band denotes inferred amplitude ranges near hypocenter after correcting for the effect of geometrical spreading and intrinsic attenuation from the observed amplitudes (Figure S1 in Supporting Information S1 for the observations and Figure S8 in Supporting Information S1 for the model calculations).Pink band denotes the acceptable range of Q − 1 (0.001-0.007) estimated from seismological observations (Liu et al., 2014;Nakajima et al., 2013;Wang et al., 2017).(d) Examples of changes in seismic-wave amplitudes along the raypath for the hypocentral distance of 87 km with two typical amplitudes (10 − 5.5 m/s and 10 − 4.5 m/s) observed at a station.Color coding is based on the Q − 1 values at each point, which is calculated using n = 0.1, C = 0.017, and the seismic-wave amplitude at each point.

Comparison With Experimentally Derived Amplitude Dependence
Experimental studies under low-pressure conditions have shown a marked amplitude-dependent attenuation for ε > ∼10 − 6 (e.g., Mavko, 1979;Stewart et al., 1983;Winkler et al., 1979), and ε of 10 − 6 corresponds to seismicwave amplitudes of 10 − 3 -10 − 2 m/s with the relation of A = εV p (Van Der Elst & Brodsky, 2010).However, these seismic-wave amplitudes are two or three orders of magnitude larger than typical amplitudes that we observe for M3-5 earthquakes (Figure 1 and Figure S1 in Supporting Information S1).Therefore, the amplitude-dependent attenuation derived from experimental studies suggests that attenuation is constant for seismic waves radiated from small to moderate earthquakes.
We here evaluate whether the experimentally derived amplitude dependence of attenuation can explain the observations.In the calculation, we assume two sets of Q − 1 L and a 1 values in Equation 1, (Q − 1 L , a 1 ) = (10 − 2 , 10 − 3 ) and (10 − 3 , 0.02), which are taken from typical values obtained in the experimental studies (Mavko, 1979;Stewart et al., 1983).We also consider one optimal set, (Q − 1 L , a 1 ) = (10 − 2.76 , 10 − 4.62 ), that best explains the observation with Equation 1.The calculated results demonstrate that the experimentally derived amplitude-dependent attenuation cannot reproduce the observed amplitude dependence of ∆Q − 1 (see red, orange, and blue lines in Figures 11b and 11c).Therefore, we conclude that seismic attenuation possesses a weak power-law amplitude dependent over an amplitude range of at least 10 − 6 ∼10 − 2 (m/s), suggesting that attenuation depends on amplitude over the entire raypath even for seismic waves radiated from small to moderate earthquakes (see Figure 10 and Figure S8 in Supporting Information S1).

Conclusion
Conventional data analyses and theoretical calculations, including seismic wave propagation and free oscillations of the Earth, have assumed that the energy dissipation is amplitude independent.This study carried out a careful and systematic spectral analysis of co-located M3-5 earthquakes with different observed amplitudes and revealed that seismic attenuation in the crust and the uppermost mantle is amplitude dependent of seismic waves with 0.1-0.2power of the amplitude for the amplitude and frequency ranges of 10 − 6 -10 − 2 m/s (strain of 10 − 10 -10 − 6 ) and 20-40 Hz, respectively.Although it may be debatable whether amplitude-dependent attenuation can be observed for deformation with smaller amplitudes and lower frequencies such as the propagation of tele-seismic waves and

Figure 1 .
Figure 1.Attenuation (Q − 1 ) versus strain amplitude (ε) plot for Q − 1 L = 0.0058 and a 1 = 0.0013 (see Figure 2 in Tisato & Quintal, 2014) in Equation1.Dashed line represents an approximate threshold value, above which attenuation becomes nonlinear.We used a relationship A = εV p (Van Der Elst & Brodsky, 2010) with V p = 7.0 km/s for the conversion from strain amplitude, ε, in rock experiments, to seismic-wave amplitude, A. Gray area represents typical seismic-wave amplitudes observed for small to moderate-sized earthquakes.

Figure 3 .
Figure3.Examples of step-by-step processing of the spectral ratios at the amplitude-ratio band of 0.9 ≤ log(A t L / A t S ) ≤ 1.1.(a) The spectral amplitudes of EQ-L and EQ-S are calculated, and their ratios are taken, whereA f L / A f S = R m .(b)The source spectra on the right side of Equation 5 are moved to the left side.(c) The common logarithm of the spectral ratios in panel (b) is taken.(d) The spectral ratios in panel (c) are divided by the average travel time of EQ-L and EQ-S.(e) The spectral ratios in panel (d) are shifted along the y-axis so that the mean value of each spectral ratio becomes zero at 20− 40 Hz.(f) All the spectral ratios in panel (e) are stacked, and the average spectral ratio is calculated (black line).(g) The right-hand side of Equation 6 is fitted to the black line at 20− 40 Hz using the least squares method (green line).Note that the vertical scale is different for each panel.

Figure 4 .
Figure 4. (a) Distributions of seismic stations (blue squares) and earthquakes (red dots) used in this study.Red triangles indicate active volcanoes.Stations are colored according to the number of analyzed earthquake pairs.(b) Vertical cross-section of hypocenters between 38.5°N and 39.5°N.Gray dots represent M JMA ≥ 2 earthquakes that occurred during the period 2003-2021.(c) Example of waveforms recorded at station N.YMDH for two co-located earthquakes with M JMA 4.1 (EQ-L) and M JMA 3.1 (EQ-S).Brackets below the P-wave signal denote the 1-s time windows (one noise and four signal windows) that were used to calculate the amplitude spectra.(d) Examples of the spectral ratio fits for four different amplitude-ratio bands of the win1 time window.Gray thin lines represent the observed spectral ratios, and black line represents the average spectral ratio.Green line represents the best-fit theoretical spectral ratio over the 20-40 Hz frequency band.The best-fit ∆Q − 1 value, its standard deviation, and the number of spectral ratios is provided at the lower right in each panel.

Figure 5 .
Figure 5. (a) Estimated ∆Q − 1 with respect to the observed amplitude ratios for four time windows.Vertical bars represent one standard deviation of ∆Q − 1 .(b) ∆Q − 1 for two time windows (win1and win3), where the colors in each square denote the average absolute amplitudes for EQ-L (upper half) and EQ-S (lower half) at that amplitude band.

Figure 7 .
Figure 7. (a) The distribution of assumed stress drops, where we adopted n of 6, 7, 7.5, and 8.The standard deviations of log∆σ ave are assumed to be 0.5.(b) The results from different values of log∆σ ave .Colored circles represent the estimates with different log∆σ ave and white circles denote our results for win1 shown in Figure 5a.

Figure 8 .
Figure 8. Schematic diagram of the model calculation.Red and blue stars represent the hypocenter of EQ-L and EQ-S, respectively.We adopted the hypocenter depth to be 71 km and the epicentral distance to station to be 50 km, both of which correspond to the average values for earthquakestation pairs used in this study.The hypocentral distance is 87 km.A raypath is calculated with the 1D velocity structure of JMA2001.Colors of the line denote P-wave velocity at each depth.Note that the hypocenters of EQ-L and EQ-S are separated by 3 km, which is the average distance of pairwise earthquakes used in this study, to show how co-located earthquakes are actually distant, but this separation is not included in the model calculation.

Figure 9 .
Figure 9. Summary of model calculation.(a) Distribution of the residuals (grayscale) between the observed and theoretical ∆Q − 1 values for the win1 window.The optimal solutions for five n-C pairs are shown by colored stars.(b) Comparison of the observed (white circles) and theoretical (colored lines) ∆Q − 1 values, which were calculated using the five highlighted n-C pairs in panel (a).Line colors correspond to those of the stars in panel (a).(c) Amplitude dependence of Q − 1 at 1 Hz, which is calculated via Q − 1 (f),(A) = CA n for the five highlighted n-C pairs in panel (a).Dark gray bands denote the amplitude ranges of seismic waves that were not included in this study, and light gray band denotes inferred amplitude ranges near hypocenter after correcting for the effect of geometrical spreading and intrinsic attenuation from the observed amplitudes (FigureS1in Supporting Information S1 for the observations and FigureS8in Supporting Information S1 for the model calculations).Pink band denotes the acceptable range of Q − 1 (0.001-0.007) estimated from seismological observations(Liu et al., 2014;Nakajima et al., 2013;Wang et al., 2017). (d) Examples of changes in seismic-wave amplitudes along the raypath for the hypocentral distance of 87 km with two typical amplitudes (10 − 5.5 m/s and 10 − 4.5 m/s) observed at a station.Color coding is based on the Q − 1 values at each point, which is calculated using n = 0.1, C = 0.017, and the seismic-wave amplitude at each point.

Figure 10 .
Figure 10.(a) Examples of attenuation variation along a raypath, which is normalized by attenuation at a seismic station.Circles and triangles denote the normalized attenuation calculated with the optimal values of n = 0.1 and 0.2, respectively, and are colored by amplitudes at each point.Amplitudes at hypocenter (D = 0 km) are adjusted so that amplitudes observed at station (D = 100 km) become ∼10 − 5 m/s.(b) Normalized attenuation as a function of hypocentral distances (D = 50, 100, 150, and 200 km) for the optimal values of n = 0.1 and n = 0.2.Colored symbols denote the maximum attenuation along the raypath, which corresponds to attenuation at hypocenter, and the colors represent amplitudes at hypocenter.White symbols with error bars represent the path-average normalized attenuation and its standard deviation.The results for n = 0.1 and n = 0.2 are slightly shifted to the horizontal direction to avoid the overlap of symbols and error bars.The observed amplitudes are adjusted to ∼10 − 5 m/s at each hypocentral distance.

Figure 11 .
Figure 11.(a) Amplitude-dependent attenuation for three sets of parameters (Q − 1 L , a 1 ) with Equation 1 (red, orange, and blue lines) and the best-fit model to our observations with Equation 13 (black line).The black line represents the observationbased optimal function of Q − 1 for n = 0.1.See text for details.Seismic-wave amplitudes, A, are calculated with A = εV p (Van Der Elst & Brodsky, 2010) for V p = 7.0 km/s.(b) Comparisons between experimentally predicted and observed ∆Q − 1 .Black circles represent the observed ∆Q − 1 values.Colors of lines are the same as in panel (a).(c) The same as in panel (b) but with the y-axis in a logarithmic scale.