Moment magnitude estimates for central Anatolian earthquakes using coda waves

A proper estimate of moment magnitude, which is a physical measure of the energy released at an earthquake source, is essential for better seismic hazard assessments in tectonically active regions. Here a coda wave modeling approach that enables the source displacement spectrum modeling of the examined event was used to estimate moment magnitudes of central Anatolia earthquakes. To achieve this aim, three-component waveforms of local earthquakes with magnitudes $\mathrm{2.0}\le M_{\mathrm{L}}\le \mathrm{5.2}$ recorded at 69 seismic stations, which were operated between 2013 and 2015 within the framework of the Continental Dynamics–Central Anatolian Tectonics (CD–CAT) passive seismic experiment, were utilized. An inversion on the coda wave traces of each selected single event in the database was performed in five different frequency bands between 0.75 and 12 Hz. The resultant moment magnitudes (M_w coda) exhibit a good agreement with routinely reported local magnitude (M_L) estimates for the study area. Apparent move-out that is particularly significant around the scattered variation of M_L–M_w coda data points for small earthquakes (M_L < 3.5) can be explained by possible biases of wrong assumptions to account for anelastic attenuation and seismic recordings with a finite sampling interval. Finally, I present an empirical relation between M_w coda and M_L for central Anatolian earthquakes.

1 Introduction

Robust and stable knowledge of source properties (e.g., moment magnitude estimates) is crucial in seismically active countries such as Turkey for a better evaluation of seismic hazard potential as this highly depends on the establishment of reliable seismicity catalogs. Moreover, accurate information on source parameters could be important when developing regional attenuation properties.

Conventional types of magnitude scales (M_L, m_b, M_S), such as the empirically derived results of using direct wave analyses, can be biased due to various effects such as source radiation pattern, directivity, and heterogeneities along the path since they may cause drastic changes in direct wave amplitude measurements (e.g., Favreau and Archuleta, 2003). Several early studies depending on the analysis of local and/or regional coda envelopes have indicated that coda wave amplitudes are significantly less variable by a factor of 3 to 5 compared to direct wave amplitudes (e.g., Mayeda and Walter, 1996; Mayeda et al., 2003; Eken et al., 2004; Malagnini et al., 2004; Gök et al., 2016). In fact, local or regional coda waves are usually considered to be generally composed of scattered waves. These wave trains can be simply explained by the single-scattering model of Aki (1969), which has been proven to be virtually insensitive to any source radiation pattern effect in contrast to direct waves due to the volume-averaging property of the coda waves sampling the entire focal sphere (e.g., Aki and Chouet, 1975; Rautian and Khalturin, 1978). In Sato and Fehler (1998) and Sato et al. (2012) an extensive review of the theoretical background of coda generation and advances in empirical observations and modeling efforts can be found.

There have been several approaches used for extracting information on earthquake source size via coda wave analyses. These approaches can be mainly divided into two groups. The first group of studies can be considered the parametric approach and essentially employs a coda normalization strategy in which measurements require a correction for empirically derived quality factors representing seismic attenuation parameters (e.g., intrinsic and scattering). In this case, the adjustment of final source properties is achieved with the help of some reference events whose seismic moments are previously estimated based on waveform inversion methods. For the forward generation of synthetic coda envelopes, either single-backscattering or more advanced multiple-backscattering approximations are used. An example for this group is an empirical method originally developed by Mayeda et al. (2003) to investigate seismic source parameters such as energy, moment, and apparent stress drop in the western United States and in the Middle East. They corrected observed coda envelopes for various influences, for instance path effect, S-to-coda transfer function, site effect, and any distance-dependent changes in coda envelope shape. Empirical coda envelope methods have been successfully applied to different regions with complicated tectonics such as northern Italy (e.g., Morasca et al., 2008), Turkey and the Middle East (e.g., Eken et al., 2004; Gök et al., 2016), and the Korean Peninsula (e.g., Yoo et al., 2011).

The second type of approach depends on estimating source and structural properties through a joint inversion technique. This technique employs a simultaneous optimization of source-, path-, and site-specific terms via a fitting procedure between a physically derived synthetic coda envelope and an observed coda envelope within a selected time window that includes both the observed coda and direct S-wave parts. Although the conventional coda normalization method essentially relies on a correction for undesired effects of the source and site amplifications, it may fail for small events with a shorter coda. This mainly stems from random seismic noise that dominates the coda, which does not satisfy the requirement of the homogeneous distribution of energy in space. In the present study, I avoid this shortcoming by involving source excitation and site amplification terms in the inversion process. To achieve this, radiative transfer theory (RTT) is employed for analytic expression of synthetic coda wave envelopes. The method was originally developed by Sens-Schönfelder and Wegler (2006) and successfully tested on local and regional earthquakes ($\mathrm{4}\le M_{\mathrm{L}}\le \mathrm{6}$) detected by the German Regional Seismic Network. Further, it has been applied to investigate source- and frequency-dependent attenuation properties in different geological settings such as the upper Rhine Graben and Molasse basin regions in Germany, western Bohemia–Vogtland in Czechia (Eulenfeld and Wegler, 2016), the entire United States (Eulenfeld and Wegler, 2017), and the central and western North Anatolian Fault Zone (Gaebler et al., 2018; Izgi et al., 2018). A more realistic earth model with anisotropic scattering conditions was considered earlier by Gusev and Abubakirov (1996) and yielded peak broadening effects of direct seismic wave arrivals. This approach that examines the propagation of P-wave elastic energy and the effect of conversion between P- and S-wave energies was later used in Zeng et al. (1991), Przybilla and Korn (2008), and Gaebler et al. (2015).

In the current work I present source spectra as the output of a joint inversion of S-wave and coda wave parts extracted from 487 local earthquakes with magnitudes $\mathrm{2.0}<M_{\mathrm{L}}<\mathrm{4.5}$ detected in central Anatolia. The approach used here employs an isotropic acoustic RTT approach for the forward calculation of synthetic coda envelopes. Gaebler et al. (2015) have observed that modeling results from isotropic scattering were almost comparable with those inferred from relatively more complex elastic RTT simulations with anisotropic scattering conditions. The use of a joint inversion technique is advantageous since it is insensitive to any potential bias, which could be introduced by external information, i.e., source properties of a reference that is obtained separately from other methods for calibration. This is mainly because of the fact that I utilize an analytical expression of a physical model involving source- and path-related parameters to describe the scattering process. Moreover, the type of optimization during joint inversion enables estimates for source parameters of relatively small-sized events compared to the one used in coda normalization methods.

2 Regional setting

The present tectonic setting of Anatolia and the surrounding regions is mainly the outcome of northward-converging movements among the African, Arabian, and Eurasian plates. To the west, the subducting African Plate with slab rollback dynamics beneath Anatolia along the Hellenic Trench has led to back-arc extension in the Aegean and western Anatolia, while compressional deformation to the east around the Bitlis–Zagros suture was explained by collisional tectonics (e.g., Taymaz et al., 1990; Bozkurt, 2001) (Fig. 1). Central Anatolia is located between an extensional regime to the west due to subduction and a compressional regime to the east due to collisional tectonics. There are several fault systems responsible for ongoing seismic activity in the region. The major fault zone, the Central Anatolian Fault Zone (CAFZ) (Fig. 2), which primarily represents a transtensional fault structure with a small amount of left-lateral offset during the Miocene (e.g., Koçyiğit and Beyhan, 1998), can be considered a boundary between the carbonate nappes of the Anatolide–Tauride block and the highly deformed and metamorphosed rocks in the Kırşehir block. To the northwest of the CAFZ, the Tuz Gölü Fault Zone (TGFZ) (Fig. 2), which is characterized by a right-lateral strike-slip motion with a significant oblique-slip normal component, appears to be collocated with the Tuz Gölü Basin sedimentary deposits as well as the crystalline rocks within the Kırşehir block (e.g., Çemen et al., 1999; Bozkurt, 2001; Taymaz et al., 2004; Çubuk et al., 2014). At the southwest tip of the study region, the EAFZ generates large seismic activity that can be identified by a rather complicated seismotectonic setting: predominantly left-lateral strike-slip motion that is well correlated with the regional deformation pattern and with existing local clusters of thrust and normal faulting events on NS- and EW-trending subsidiary faults, respectively (Bulut et al., 2012). Such complicated behavior explains kinematic models (e.g., Riedel shear, anti-Riedel shear models) of the shear deformation zone evolution (Tchalenko, 1970). It connects to the NAFZ at the Karlıova Triple Junction (Bozkurt, 2001) and to the south splits into various segments near the Adana Basin (Kaymakci et al., 2006) (Fig. 2). Toward the south, the EAFZ reaches the Dead Sea Fault Zone (DSFZ) that has a key role in accommodating northward relative motions of the Arabian and African plates with respect to Eurasia.

Figure 1 Major tectonic features of Turkey and the adjacent areas. The plate boundary data used here are taken from Bird (2003). Subduction zones are black, continental transform faults are red, continental rift boundaries are green, and spreading ridge boundaries are yellow. NAFZ, EAFZ, and DSFZ are the North Anatolian Fault, East Anatolian Fault, and the Dead Sea Fault, respectively.

Figure 2 Epicentral distribution of all local events selected from the study area in the KOERI catalog. Gray circles represent earthquakes with poor quality that are not considered for the current study, while black indicates the location of local events with good quality. Red circles among these events are 487 events used in coda wave inversion since they are successful at passing the quality criteria for a further preprocessing procedure.

3 Data

The present work utilizes three-component waveforms of local seismic activity detected at 72 broadband seismic stations (Fig. 2) that were operated for 2 years between 2013 and 2015 within the framework of a temporary passive seismic experiment, Continental Dynamics–Central Anatolian Tectonics (CD–CAT) (Portner et al., 2018). I benefit from revisited standard earthquake catalog information that is routinely released by the Kandilli Observatory and Earthquake Research Institute (KOERI) (publicly available at http://www.koeri.boun.edu.tr, last access: 15 February 2019) to extract waveform data for a total of 2231 examined events with station–event pair distance less than 120 km and focal depths less than 10 km. Most of the detected seismic activity in the study area is associated with several fault zones in the region, i.e., the EAFZ, CAFZ, DSFZ, and TGFZ. Here I note that the use of only local earthquakes is to exclude possible biases, which may be introduced by Moho-boundary-guided Sn waves. Upper crustal earthquakes with less than 10 km focal depths are preferred in this study to exclude the effect of relatively large-scale heterogeneities on coda wave trains. Additionally, I performed a visual inspection over all waveforms to ensure high-quality waveforms. The final event number was reduced to 1193. Selected station and event distributions can be seen in Fig. 2.

Observed waveforms were prepared at five different frequency bands with central frequencies at 0.75, 1.5, 3.0, 6.0, and 12.0 Hz via a Butterworth band-pass filtering process. In the next step, I applied a Hilbert transform to filtered waveform data in order to obtain the total energy envelopes. An average crustal velocity model was used to predict P- and S-wave onsets on envelopes and then based on this information: (i) the noise level prior to the P-wave onset was eliminated; (ii) the S-wave window was determined starting at 3 s prior to and 7 s after S-wave onset, as this allowed for the inclusion of all direct S-wave energy; and (iii) starting at the end of the S-wave window, a coda window of 100 s at maximum was determined. The length of coda windows can be shorter when the signal-to-noise ratio (SNR) is less than 2.5 or when there are coda waves from two earthquakes (e.g., because of an aftershock sequence) within the same analysis window, which can cause another rise instead of a decline in the envelope. I omit earthquakes with less than 10 s of coda length from the database. Taking into account these criteria, coda waveforms extracted from 6541 source–receiver pairs were used for further data processing.

4 Method

I adopted an inversion procedure that was originally developed by Sens-Schönfelder and Wegler (2006) and later modified by Eulenfeld and Wegler (2016). The forward part, which involves the calculation of energy density for a specific frequency band under the assumption of an isotropic source, is expressed in Sens-Schönfelder and Wegler (2006) as follows:

$$\begin{array}{}\text{(1)}& E_{\mathrm{mod}}\left(t,r\right)=WR\left(r\right)G(t,r,g)e^{-bt},\end{array}$$

where W gives the source term and is frequency dependent. R(r) indicates the energy site amplification factor and b is the intrinsic attenuation parameter. $G(t,r,g)$ represents Green's function, which includes the scattered wave field as well as the direct wave and its expression, and is given by Paasschens (1997) as follows:

$$\begin{array}{ll}{\displaystyle }G\left(t,r,g\right)& {\displaystyle }=e^{(-v_{\mathrm{0}}t{g}_{\mathrm{0}})}\left[{\displaystyle \frac{\mathit{\delta }\left(r-v_{\mathrm{0}}t\right)}{\mathrm{4}\mathit{\pi }{r}^{\mathrm{2}}}}+{\left({\displaystyle \frac{\mathrm{4}\mathit{\pi }{v}_{\mathrm{0}}}{\mathrm{3}{g}_{\mathrm{0}}}}\right)}^{-\frac{\mathrm{3}}{\mathrm{2}}}{t}^{-\frac{\mathrm{3}}{\mathrm{2}}}\right.\\ {\displaystyle }& {\displaystyle }\times {\left(\mathrm{1}-{\displaystyle \frac{r^{\mathrm{2}}}{v_{\mathrm{0}}^{\mathrm{2}}{t}^{\mathrm{2}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{8}}}K\left(v_{\mathrm{0}}t{g}_{\mathrm{0}}{\left(\mathrm{1}-{\displaystyle \frac{r^{\mathrm{2}}}{v_{\mathrm{0}}^{\mathrm{2}}{t}^{\mathrm{2}}}}\right)}^{\frac{\mathrm{3}}{\mathrm{4}}}\right)\\ \text{(2)}& {\displaystyle }& {\displaystyle }\left.H\left(v_{\mathrm{0}}t-r\right)\right].\end{array}$$

Here the term within the Dirac delta function represents the direct wave and the other term indicates scattered waves. v0 describes the mean S-wave velocity, while g0 is the scattering coefficient.

Possible discrepancies between predicted (Eq. 1) and observed energy densities for each event at each station with N_ij time samples (index k) in a specific frequency band can be minimized using

$$\begin{array}{}\text{(3)}& \mathit{ϵ}\left(g\right)=\sum _{i,j,k}^{N_{\mathrm{S}},N_{\mathrm{E}},N_{ij}}{\left(\ln E_{ijk}^{\mathrm{obs}}-\ln E_{ijk}^{\mathrm{mod}}\left(g\right)\right)}^{\mathrm{2}}.\end{array}$$

Here, the numbers of stations (index i) and events (index j) are shown by N_S and N_E, respectively. Optimization of g will be achieved by fulfilling the following equality:

$$\begin{array}{}\text{(4)}& \ln E_{ijk}^{\mathrm{obs}}=\ln E_{ijk}^{\mathrm{mod}}\end{array}$$

or

$$\begin{array}{}\text{(5)}& \ln E_{ijk}^{\mathrm{obs}}=\ln G\left(t_{ijk},r_{ijk},g\right)+\ln R_i+\ln W_j-b{t}_{ijk}.\end{array}$$

Equation (5) simply defines an overdetermined inversion problem with ${\sum }_{i,j}{N}_{ij}$ equation systems and with $N_{\mathrm{S}}+N_{\mathrm{E}}+\mathrm{1}$ variables, and thus b, R_i, and W_j can be solved via a least-squares technique. ϵ(g) can be defined as a sum over the squared residuals of the solution. As can be seen from Eq. (1), there is an obvious trade-off between R_i and W_j, which I can manage by fixing the geometrical mean of R_i to 1 (ΠR_i=1). Equation (1) also implies a rather moderate trade-off between W_j and b. The trade-off between g and other inverted parameters is usually small since this parameter is fixed through the energy ratio of the direct S wave and the level of the coda waves (Gaebler et al., 2018).

Eulenfeld and Wegler (2016) present a simple procedure to perform the inversion.

• i.

  Calculate Green's functions through the analytic approximation of the solution for 3-D isotropic radiative transfer (e.g., Paasschens, 1997; Sens-Schönfelder and Wegler, 2006) by using fixed scattering parameters, and minimize Eq. (5) to solve for b, R_i, and W_j via a weighted least-squares approach.

• ii.

  Calculate ϵ(g) using Eq. (3).

• iii.

  Repeat (i) and (ii) by selecting different g to find the optimal parameters g, b, R_i, and W_j that finally minimize the error function ϵ.

In Fig. 3 an example for the minimization process that was applied at five different frequency bands is displayed for one selected event at recorded stations of the CD–CAT project.

Figure 3 An example from the inversion procedure explained in Sect. 4. Here coda envelope fitting optimization is performed on band-pass-filtered (4–8 Hz) digital recordings of an earthquake (9 April 2014, M_w coda 3.2) extracted for seven seismic stations that operated within the CD–CAT array. The large panel displays the error function ε as a function of g₀. The thick blue cross represents the optimal value of g=g₀. The small panels show the least-squares solution of the weighted linear equation system for the first six guesses and the optimal guess for g₀. The dots and gray curves indicate the ratio between energy (E^obs) and the Green's function (G) obtained for direct S waves and observed envelopes at various stations, respectively (please note that during this optimization process envelopes are corrected for the obtained site corrections R_i). The slope of the linear curve in each small panel yields −b in relation to the intrinsic attenuation. The linear curve has an intercept of W, representing source-related terms in the right-hand side of Eq. (5).

The minimization described above for different frequencies will yield an unknown spectral source energy term, W_j, as well as a site response, R_i, and attenuation parameters, b and g, that will satisfy optimal fitting between observed and predicted coda wave envelopes. An example for this fitting can be seen in Fig. 4. The present study deals with the frequency dependency of W_j since this information can later be useful to obtain the source displacement spectrum and thus seismic moment and moment magnitudes of analyzed earthquakes using the formula of the S-wave source displacement spectrum for a double-coupled source in the far field, which is given by Sato et al. (2012):

$$\begin{array}{}\text{(6)}& \mathit{\omega }M\left(f\right)=\sqrt{{\displaystyle \frac{\mathrm{5}{\mathit{\rho }}_{\mathrm{0}}{v}_{\mathrm{0}}^{\mathrm{5}}W}{\mathrm{2}\mathit{\pi }{f}^{\mathrm{2}}}}},\end{array}$$

where W indicates the radiated S-wave energy at a center frequency f, while v₀ and ρ₀ represent the mean S-wave speed and medium density, respectively.

Figure 4 (a) Results of the inversion for the 9 April 2014 M_wcoda 3.2 earthquake: sample fits between observed and calculated energy densities in the frequency band 0.5–1.0 Hz are given for six different stations (see the upper right corner for event ID, station name, and distance to hypocenter). Note that light blue curves represent the observed envelope. Smoothed observed and calculated envelopes in each panel are presented by blue and red curves, respectively. Blue and red dots show the location of the average value for observed and calculated envelopes within the S-wave window, respectively. (b) The same as in (a) obtained in the frequency band 4.0–8.0 Hz.

The relation between the obtained source displacement spectrum and seismic moment value was described earlier in Abercrombie (1995) by

$$\begin{array}{}\text{(7)}& \mathit{\omega }M\left(f\right)=M_{\mathrm{0}}{\left(\mathrm{1}+{\left({\displaystyle \frac{f}{f_c}}\right)}^{\mathit{\gamma }n}\right)}^{-\frac{\mathrm{1}}{\mathit{\gamma }}},\end{array}$$

where n is related to the high-frequency falloff and γ is known as a shape parameter that controls the sharpness of the spectrum at corner frequencies between the constant level M₀ (low-frequency part) and the falloff with f⁻ⁿ (high-frequency part). Taking the logarithm of Eq. (7) gives

$$\begin{array}{}\text{(8)}& \ln \mathit{\omega }M\left(f\right)=\ln M_{\mathrm{0}}-{\displaystyle \frac{\mathrm{1}}{\mathit{\gamma }}}\ln \left(\mathrm{1}+{\left({\displaystyle \frac{f}{f_c}}\right)}^{\mathit{\gamma }c}\right).\end{array}$$

Equation (8) describes an optimization problem in which the observed source displacement spectrum data (left-hand side) can be inverted for four unknown source parameters, M₀, γ, n, and f_c (right-hand side), in a simultaneous least-squares inversion scheme. Finally, moment magnitude, M_w, can be calculated from the modeled source parameters and seismic moment, M₀, using a formula given by Hanks and Kanamori (1979):

$$\begin{array}{}\text{(9)}& M_{\mathrm{w}}={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\log }_{\mathrm{10}}{M}_{\mathrm{0}}-\mathrm{10.7}.\end{array}$$

5 Results and discussion

5.1 Coda wave source spectra

Figure 5 displays observed values of source spectra established by inserting the inverted spectral source energy term W at each frequency in Eq. (6) for all analyzed events. Each curve in this figure represents the model spectrum estimate based on the inversion procedure described in the previous section. Modeled spectrum characteristics computed for 487 local earthquakes whose geographical distribution is presented in Fig. 2 suggest, in general, that I was able to obtain a typically expected source displacement spectrum with a flat region around the low-frequency limit and a decaying behavior above a corner frequency.

Figure 5 All individual observed (black squares) and predicted (gray curve) source displacement spectra observed at 72 stations from 487 local earthquakes in central Anatolia.

Owing to the multiple-scattering process within small-scale heterogeneities that make coda waves gain an averaging nature, the variation in coda amplitudes due to differences in the source radiation pattern and path effect is reduced (Walter et al., 1995; Mayeda et al., 2003). Eulenfeld and Wegler (2016) found that the radiation pattern has only a minor influence on the S-wave coda, while it might disturb attenuation models inferred from direct S-wave analyses unless the station distribution relative to the earthquakes indicates good azimuthal coverage.

Conventional approaches (e.g., Abercrombie, 1995; Kwiatek et al., 2011) to estimate source parameters, such as corner frequency, seismic moment, and high-frequency falloff through fitting of observed displacement spectra observed at a given station in an inversion scheme, could be misleading since these methods usually (i) assume a constant value of attenuation effect (no frequency variation) defined by a factor exponent ($-\mathit{\pi }ft{Q}^{-\mathrm{1}}$) over the spectrum (ii) and assume an omega-square model with a constant high-frequency falloff parameter, n=2. Following Sens-Schönfelder and Wegler (2006) and Eulenfeld and Wegler (2016), however, I estimate attenuation parameters (intrinsic and scattering) separately within a simultaneous inversion procedure in which the high-frequency falloff parameter varies. This is fairly consistent with early studies (e.g., Ambeh and Fairhead, 1991; Eulenfeld and Wegler, 2016) wherein significant deviations from the omega-square model (n>3) were reported, implying that the omega-square model as a source model for small earthquakes must be reconsidered in its general acceptance. It has been previously observed that the source spectra, especially for large earthquakes, could be better explained by models of two corner frequencies (e.g., Papageorgiou and Aki, 1983; Joyner, 1984; Atkinson, 1990). Recently, Denolle et al. (2016) observed that a conventional spectral model of a single-corner frequency and high-frequency falloff rate could not explain P-wave source spectra of thrust earthquakes with magnitudes M_w 5.5 and above. Instead, they suggested the double-corner frequency model for large global thrust earthquakes with a lower-corner frequency related to source duration and with an upper-corner frequency suggesting a shorter timescale unrelated to source duration, which exhibits its own scaling relation. Uchide and Imanishi (2016) reported that similar differences from the omega-square model would also be valid for smaller earthquakes by using a spectral ratio technique that involves empirical Green's function (EGF) events to avoid having complete knowledge of path and site effects for shallow target earthquakes (M_w 3.2–4.0) in Japan. The source spectra for many of the target events in their study suggested a remarkable discrepancy from the omega-square model for relatively small earthquakes. They explained such differences by incoherent rupture due to heterogeneities in fault properties and applied stress, the double-corner frequency model, and the possibility of a high-frequency falloff exponent value slightly higher than 2. In my case, the smallest event was with M_w coda larger than 2.0; thus, I had no chance to make a similar comparison to that of Eulenfeld and Wegler (2016). However, high-frequency falloff parameters varied from n=0.5 to n=4. A notable observation in the distribution of n was n=2 or n=2.5, which would be better explained for earthquakes with M_w coda > 4.0, whereas the smaller magnitudes exhibited a more scattered pattern of variation in n (Fig. 7). Eulenfeld and Wegler (2016) claimed that the use of separate estimates of attenuation or correction for the path effect via an empirically determined Green's function would be a better strategy to invert station displacement spectra for source parameters. This is mainly because for smaller earthquakes (with n>2), in particular, assuming an omega-square model can distort the estimates of corner frequency and even seismic moment, especially in regions where Q is strongly frequency dependent. Thus, independent estimates of Q during station displacement spectra inversions for source parameters must be taken into account or the influence of paths such as attenuation must be removed via empirically determined Green's functions (Eulenfeld and Wegler, 2016).

Figure 6 Scatter plot between local magnitudes (M_L) of analyzed events with coda-wave-derived magnitudes (M_w coda) of the same events. The outcome of a linear regression analysis yielded an empirical formula (e.g., Eq. 10) to identify the overall agreement represented by the gray straight line. Yellow and red dashed lines indicate the upper and lower limit of linearly fitting to that scatter.

5.2 Coda-wave-derived magnitude vs. M_L catalog magnitude

A scatter plot between catalog magnitudes based on local magnitudes (M_L) and my coda-derived magnitudes (M_w coda), which are inferred from the resultant frequency-dependent source displacement spectra and thus seismic moment (e.g., Eq. 9), is shown in Fig. 6. Such a comparison suggests an overall coherency between both types of magnitudes. This implies that a very simple model of a first-order approximation for S-wave scattering with an isotropic acoustic radiative transfer approach can be efficient to link the amplitude and decaying character of coda wave envelopes to the seismic moment of the source.

Figure 7 Same scatter plot displayed in Fig. 6. Here the color code indicates the estimated high-frequency falloff parameter for each inverted event.

In the present study, a linear regression analysis performed between M_w coda and M_L magnitudes (Fig. 5) resulted in an empirical formula that can be employed to convert local magnitudes into a coda-derived moment magnitude calculation of local earthquakes in this region:

$$\begin{array}{}\text{(10)}& M_{\mathrm{w}-\mathrm{coda}}=\mathrm{1.1655}\pm \mathrm{0.0337}\times M_{\mathrm{L}}-\mathrm{0.7085}\pm \mathrm{0.0128}.\end{array}$$

Bakun and Lindh (1977) empirically described the linear log seismic moment–local magnitude relation between seismic moments (Mo) and local magnitudes (M_L) for earthquakes near Oroville, California. Several other studies have also investigated an optimum relation between M_w and M_L by implementing linear and/or nonlinear curve-fitting approaches. Malagnini and Munafò (2018) proposed that two different linear fits separated by a crossover M_L=4.31 could represent M_L–M_w data points obtained from earthquakes in the central and northern Apennines, Italy. Several coefficient-of-regression analyses in their fits account for the combined effects of source scaling and crustal attenuation as well as regional attenuation, focal depth, and rigidity at the source. Goertz-Allmann et al. (2011), for instance, introduced a hybrid type of scaling relation that is linear below M_L 2 and above M_L 4 with a quadratic relation in between ($\mathrm{2}\le M_{\mathrm{L}}\le \mathrm{4}$) for earthquakes in Switzerland detected between 1998 and 2009. Edwards and Rietbrock (2009) employed a second-order polynomial equation to relate local values routinely reported in Japan Meteorological Agency (JMA) magnitudes and moment magnitudes. More recently, using multiple spectral ratio analyses, Uchide and Imanishi (2018) estimated relative moment magnitudes for the Fukushima Hamadori and northern Ibaraki prefecture areas of Japan and reported a quadratic form of correlation between JMA magnitudes and moment magnitudes. The resultant empirical curve in Uchide and Imanishi (2018) implied a considerable discrepancy between the moment magnitudes and the JMA magnitudes, with a slope of 1∕2 for microearthquakes suggesting possible biases introduced by anelastic attenuation and the recording by a finite sampling interval.

The apparent move-out in Fig. 5 and Eq. (10) presumably stems from the use of different magnitude scales for comparison. Conventional magnitude scales such as M_L, inferred from phase amplitude measurements, are seemingly sensitive to attenuation and 2-D variation along the path (Pasyanos et al., 2016). Unlike local magnitude scales, seismic-moment-based moment magnitude (M_w) essentially represents a direct measure of the strength of an earthquake caused by fault slip and is estimated from a relatively flat portion of source spectra at lower frequencies that can be less sensitive to near-surface attenuation effects. The consistency between coda-derived moment magnitude and local magnitude scales for the earthquakes with M_w coda > 3.0 indicates that my nonempirical approach successfully worked in this tectonically complex region. This observation is anticipated for relatively large earthquakes, since more energy will be characteristic at lower frequencies. I observed a similar type of consistency in early studies that investigate source properties of local and regional earthquakes based on empirical coda methods with simple 1-D radially symmetric path correction (e.g., Eken et al., 2004; Gök et al., 2016). Coda-wave-derived source parameters were obtained with high precision in Mayeda et al. (2005), Phillips et al. (2014), and Pasyanos et al. (2016) following the use of 2-D path-corrected station techniques to consider amplitude–distance relationships. Observable outliers in Fig. 5 for events with less than M_w 3.5, however, can be attributed to possible biases in local magnitude values taken from the catalog or small biases in my intrinsic ($Q_{\mathrm{i}}^{-\mathrm{1}}$) and scattering ($Q_{\mathrm{s}}^{-\mathrm{1}}$) attenuation terms. Another possible contribution to such a mismatch might be associated with the influences of mode conversions between body and surface waves or surface-to-surface wave scattering that are not restricted to low frequencies (< 1 Hz) (Sens-Schönfelder and Wegler, 2006).

6 Conclusions

This study provides moment magnitude estimates as a direct physical measure of seismic energy for local earthquakes with magnitudes $\mathrm{2.0}\le M_{\mathrm{L}}\le \mathrm{5.2}$ recorded at 69 seismic stations in central Anatolia. The source displacement spectra were obtained following the application of a coda wave modeling procedure that employs a simultaneous optimization of source-, path-, and site-specific terms by fitting physically derived synthetic coda envelopes and observed coda envelopes. Radiative transfer theory was used for the analytic expression of synthetic coda wave envelopes. Overall consistency between M_w coda and M_L suggests that my nonempirical approach successfully worked in this tectonically complex region. Variation of the high-frequency falloff parameter indicated that for smaller earthquakes (n>2) assuming an omega-square model can distort the estimates of corner frequency and even seismic moment, especially in regions where Q is strongly frequency dependent. Since the present study mainly focuses on source properties of local earthquakes in the study area, scattering and intrinsic attenuation properties that are other products of my coda envelope fitting procedure will be examined in detail within a future work. Finally, a linear regression analysis resulted in an empirical relation developed between M_w coda and M_L, which will be a useful tool in the future to quickly convert catalog magnitudes into moment magnitudes for local earthquakes in the study area.