Moment magnitude estimates for Central Anatolian 1 earthquakes using coda waves 2

8 Abstract. Proper estimate of moment magnitude that is a physical measure of the energy released at earthquake 9 source is essential for better seismic hazard assessments in tectonically active regions. Here a coda wave 10 modeling approach that enables the source displacement spectrum modeling of examined event was used to 11 estimate moment magnitude of central Anatolia earthquakes. To achieve this aim, three component waveforms 12 of local earthquakes with magnitudes 2.0 ≤ M L ≤ 5.2 recorded at 69 seismic stations which have been operated 13 between 2013 and 2015 within the framework of the CD-CAT passive seismic experiment were utilized. An 14 inversion on the coda wave traces of each selected single event in our database was performed in five different 15 frequency bands between 0.75 and 12 Hz. Our resultant moment magnitudes (M W -coda) exhibit a good 16 agreement with routinely reported local magnitude (M L ) estimates for the study area. Apparent move-out that is, 17 particularly, significant around the scattered variation of M L -M W -coda data points for small earthquakes 18 (M L <3.5) can be explained by possible biases of wrong assumptions to account for anelastic attenuation and of 19 seismic recordings with finite sampling interval. Finally, we present an empirical relation between M W -coda and 20 M L for central Anatolian earthquakes.


Introduction
The 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 establishment of reliable seismicity catalogs.Moreover, accurate information on source parameters could be important when developing regional attenuation properties.
Conventional type of magnitude scales (M L , m b , M S ) as the result of empirically derived 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).Instead 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 that 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 have 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 study on the theoretical background of coda generation and advances of empirical observations and modelling efforts can be found in details.
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 as the parametric approach and essentially employs 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, adjustment of final source properties are achieved with the help of some reference events whose seismic moments are previously estimated based on waveform inversion methods.For forward generation of synthetic coda envelopes, either single-backscattering or more advanced multiple-backscattering approximation are used.An example to 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 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 method have been successfully applied to different regions with complicated tectonics such as northern Italy (e.g.Morasca et al., 2008), Turkey and Middle East (e.g.Eken et al., 2004;Gök et al. 2016); or Korean Peninsula (e.g.Yoo et al., 2013).
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 physically derived synthetic coda envelope and 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 the 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 homogeneous distribution of energy in space.
In the present study, we avoid this shortcoming by involving source excitation and site amplification terms in the inversion process.To achieve this, the 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 (4 ≤ Ml ≤ 6) detected by the German Regional Seismic Network.Further it has been applied to investigate source and frequency dependent attenuation properties of different geological settings, i.e., Upper Rhine Graben and Molasse Basin regions in Germany and western Bohemia/Vogtland in Czechia (Eulenfeld and Wegler, 2016); entire United States (2017); central and western North Anatolian Fault Zone (Gaebler et al., 2018;Izgi et al., 2018).A more realistic earth model in which anisotropic scattering conditions were earlier considered by Gusev and Abubakirov (1987) yielded peak broadening effects of the 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 and Aki (1991), Przybilla and Korn (2008), Gaebler et al. (2015).
In the current work I present source spectra as the output of a joint inversion of S-and coda waves parts extracted from 487 local earthquakes with magnitudes 2.0 < ML < 4.5 detected in central Anatolia.The approach used here employs isotropic acoustic RTT approach for 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 we utilize an analytical expression of physical model involving source, and path related parameters to describe the scattering process.Moreover the type of optimization during joint inversion enables the estimates for source parameters of relatively small sized events compared to the one used in coda-normalization methods.

Regional Setting
Present tectonic setting of Anatolia and surrounding regions have been mainly the outcome of the northward converging movements among Africa, Arab, and Eurasian plates.To the west, the subducting African plate with a slab roll-back dynamics beneath Anatolia along 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 the subduction, and a compressional regime to the east due to the 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 as 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, 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 et al., 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 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 nearby 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 Arabian and African Plates with respect to Eurasia.

Data
The present work utilizes three-component waveforms of local seismic activity detected at 72 broadband seismic stations (Fig. 2) that have been operated for 2 years between 2013 and 2015 within the framework of a temporary passive seismic experiment, the Continental Dynamics-Central Anatolian Tectonics (CD-CAT) (Portner et al., 2018).We benefit from revisited standard earthquake catalogue information that is routinely released by the Kandilli Observatory and Earthquake Research Institute (KOERI) (publicly available at http://www.koeri.boun.edu.tr) 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 to several fault zones in the region, i.e., the EAFZ, CAFZ, DSFZ, TGFZ, etc.Here we 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 effect of relatively large-scale heterogeneities on coda wave trains.Additionally, we performed a visual inspection over all waveforms to ensure high-quality waveforms.Our final event number reduced to 1193.Selected station and event distributions can be seen in Figure 2.
Observed waveforms were prepared at 5 different frequency bands with central frequencies at 0.75, 1.5, 3.0, 6.0, 12.0 Hz via a Butterworth band-pass filtering process.In the next step, we applied 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 Pwave onset was eliminated (ii) S-wave window was determined starting at 3s prior to and 7 s afterwards S-wave onset as this allowed to include all direct S-wave energy, (iii) starting at the end of the S-wave window, a coda window of 100s at maximum was determined.Length of coda windows can be shorter when 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.
We omit the earthquakes with less than 10 s of coda length from our database.Taking into account of these criteria, finally coda waveforms extracted from 6541 source-receiver pairs were used for further data process.

Method
We 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 calculation of energy density for a specific frequency band under assumption of an isotropic source, is expressed in Sens-Schönfelder and Wegler (2006) as follows: (1) where W gives source term and it is frequency dependent.R(r) indicates the energy site amplification factor and b is intrinsic attenuation parameter.(, , ) represents Green's function that includes scattered wave field as well as direct wave and its expression is given by Paasschens (1997) as follows: Here the term within Dirac delta function represents direct wave and other term indicates scattered waves.v0 describes the mean S-wave velocity while g0 is the scattering coefficient.
Possible discrepancy 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: Here, the number 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 following equality: Equation 5 simply define an overdetermined inversion problem with number equation systems and with N S + N E + 1 variables and thus b,  !, and  !can be solved via a least-squares technique.  can be defined as sum over the squared residuals of the solution.As can be seen from equation 1 that there is an obvious trade-off between  ! and  !, which we can manage by fixing the geometrical mean of  ! to 1 (Π != 1).Equation 1 also implies rather moderate trade-off between  ! and b.Trade-off between g and other inverted parameters are 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 recipe 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 equation 5 to solve for b,  !, and  !via a weighted least-squares approach.
(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.
Minimization described above for different frequencies will yield unknown spectral source energy term,  ! as well as site response, R i and attenuation parameters, b, and g that will satisfy optimal fitting between observed and predicted coda wave envelopes.Example for this fitting can be seen in Figure 4.The present study deals with frequency dependency of  !since this information can be later useful to obtain source displacement spectrum and thus seismic moment and moment magnitudes of analyzed earthquakes using the formula of the Swave source displacement spectrum for a double-couple source in the far-field, which is given by Sato et al. (2012): where W indicates the radiated S-wave energy at a center frequency f while  ! and  !represent the mean Swave speed and medium density, respectively.
The relation between the obtained source displacement spectrum and seismic moment value was earlier described in Abercrombie (1995) by: where n is related to the high-frequency fall-off and γ is known as shape parameter that controls the sharpness of spectrum at corner frequency between the constant level M 0 (low frequency part) and the fall-off with f −n (high frequency part).Taking the logarithm of equation 7 gives: Eq. 8 describes an optimization problem where the observed source displacement spectrum data (left-hand side) can be inverted for four unknown source parameters,  !, , n, and  !(right-hand side) in a simultaneous leastsquares inversion scheme.Finally moment magnitude, M W can be calculated from modeled source parameters, seismic moment, M 0 using a formula given by Hanks and Kanamori (1979): that we were able to obtain typically expected source displacement spectrum with a flat region around the low frequency limit and a decaying behaviour above a corner frequency.

Coda wave source spectra
Owing to the multiple-scattering process within small scale heterogeneities that makes coda waves gain an averaging nature, the variation in coda amplitudes due to differences in source radiation pattern and path effect are reduced (Walter et al., 1995;Mayeda et al., 2003).Eulenfeld and Wegler (2016) found that radiation pattern would have only a minor influence on the S-wave coda while it might disturb attenuation models inferred from the direct S-wave analyses unless the station distribution relative to the earthquakes indicates a good azimuthal coverage.
Conventional approaches (e.g.Abercrombie, 1995;Kwiatek et al., 2011) to estimate source parameters such as corner frequency, seismic moment, high-frequency fall-off 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 exp (−πftQ −1 ) over the spectrum, (ii) and assume omega-square model with a constant high-frequency fall-off parameter, n=2.Following Sens-Schönfelder and Wegler ( 2006) and Eulenfeld and Wegler (2016), however, we estimate attenuation parameters (intrinsic and scattering) seperately within a simultaneous inversion procedure in which high-frequency fall-off parameter varies.This is fairly consistent with early studies (e.g.Ambeh and Fairhead, 1991;Eulenfeld and Wegler, 2016) where 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.
Earlier it has been well-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 conventional spectral model of a single-corner frequency and high-frequency fall-off rate could not explain P wave source spectra of thrust earthquakes with magnitude Mw 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 time scale unrelated to source duration, which exhibits its own scaling relation.Uchide and Imanishi (2016) reported similar differences from the omega-square model would be valid also for smaller earthquakes by using spectral ratio technique that involves empirical Green's function (EGF) events to avoid having a 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 omegasquare 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 possibility of a high-frequency falloff exponent value slightly higher than 2. In our case, the smallest event was with M W -coda larger than 2.0, thus we had no chance to make a similar compared to that of Eulenfeld and Wegler (2016).
However, high-frequency fall-off parameters varied from n=0.5 to n=4.A notable observation in the distribution of n was n=2 or n=2.5 would be better explained for earthquakes with M W -coda >4.0 whereas the smaller magnitudes exhibited more scattered pattern of variation in n (Figure 7).Eulenfeld and Wegler (2016) claimed that the use of separate estimates of the attenuation or correction for path effect via emprically determined Green's function would be better strategy in order to invert station displacement spectra for source parameters.This is mainly because smaller earthquakes (with n>2), in particular, assuming 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 path such as attenuation must be removed via empirically determined Green's functions (Eulenfeld and Wegler, 2016).

Coda wave-derived magnitude vs. M L catalogue magnitude
A scatter plot between catalogue magnitudes based on local magnitudes (M L ) and our coda-derived magnitudes (M W -coda) that are inferred from resultant frequency dependent source displacement spectra and thus seismic moment (e.g.Eq. 9) is shown in Fig. 6.Such 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 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.
In the present study, a linear regression analyses 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 coda-derived moment magnitude calculation of local earthquakes in this region: More recently, using multiple spectral ratio analyses Uchide and Imanishi (2018) estimated relative moment magnitudes for the Fukushima Hamadori and the northern Ibaraki prefecture areas of Japan and reported a quadratic form of correlation between JMA magnitudes and moment magnitudes.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.
Apparent move-out in Fig. 5 and Eq. 10, presumably stems from the use of different magnitude scales for comparison.Conventional magnitudes scales such as M L , mb inferred from phase amplitude measurements are seemingly sensitive to attenuation and 2D 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 relatively flat portion of source spectra at lower frequencies that can be less sensitive to the 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 our non-empirical 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.We observed similar type of consistency in early studies that investigate source properties of local and regional earthquakes based on emprical coda methods with simple 1-D radially symmetric path correction (e.g.Eken et al., 2004;Gök et al., 2016).Coda waves-derived source parameters were obtained with high-precision in Mayeda et al. (2005), Phillips et al. (2014), Pasyanos et al. (2016) following the use of 2-D path-corrected station techniques to consider the amplitude-distance relationships.Observable outliers in Figure 5, for the events with less than Mw 3.5, however, can be attributed to the either possible biases on local magnitude values taken from the catalogue or small biases on our intrinsic ( !!! ) and scattering (  !!! ) attenuation terms.One another possible contribution to such mismatch might be associated to the influences of mode conversions between body and surface waves or surface-to-surface wave scattering that are not restricted to low frequencies (<1Hz) (Sens-Schönfelder and Wegler, 2006).

Conclusions
This study provides moment magnitude estimates as a direct physical measure of the seismic energy for local earthquakes with magnitudes 2.0 ≤ M L ≤ 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 envelope and observed coda envelopes.The Radiative Transfer Theory was used for analytic expression of synthetic coda wave envelopes.Overall consistency between MW-coda and ML suggests that our non-empirical approach successfully worked in this tectonically complex region.Variation of high-frequency fall-off parameter indicated that for smaller earthquakes (n>2) assuming 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 our coda envelope fitting procedure will be examined in details within a future work.Finally, a linear regression analysis resulted in an empirical relation developed between MW-coda and ML, which will be a useful tool in the future to quickly convert catalogue magnitudes into moment magnitudes for local earthquakes in the study area.

Data and resources
The python code used for carrying out the inverse modeling is available under the permissive MIT license and is distributed at https://github.com/trichter/qopen.We are grateful to the IRIS Data Management Center for maintaining, archiving and making the continuous broadband data used in this study open to the international scientific community.The KOERI is specially thanked for providing publicly open local seismicity catalogues.14    19

Figure 5
Figure 5 displays observed values of source spectra established by inserting 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 Figure 2 suggest, in general, 10)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.Beside this several other studies investigated to find an optimum relation between M W and M L by implementing linear and/or nonlinear curve-fitting approaches.Malagnini and Munafò (2018) proposed two different linear fits separated by a crossover M L =4.31 could represent M L -M W data points obtained from earthquakes of 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 source.Goertz- Allmann et al. (2011), for instance, introduced hybrid type of scaling relation that is linear below M L 2 and above M L 4 and a quadratic relation in between (2 ≤ M L ≤ 4) for earthquakes in Switzerland detected between 1998 and 2009.Edwards and Rietbrock (2009) employed a second-order polynomial equation to relate local magnitudes routinely reported in the Japan Meteorological Agency (JMA) magnitude and moment magnitude.

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

Figure 2 :
Figure 2: Epicentral distribution of all local events selected from the study area in the KOERI catalogue.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 quality criteria of further pre-processing procedure.

Figure 3 :
Figure 3: An example from the inversion procedure explained in chapter 3.Here coda envelope fitting optimization is performed on band-pass filtered (4-8Hz) digital recordings of an earthquake (2014 April 09, M W -coda3.2) extracted for 7 seismic stations that operated within the CD-CAT array.Large panel at the lower left-hand side displays the error function ε as a function of g 0 .Thick blue cross here represent the optimal value of g = g 0 .Other small panels at upper and right-hand side show the least-squares solution of the weighted linear equation system for the first 6 guesses and optimal guess for g 0 .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 notice that during this optimization process envelopes are corrected for the obtained site corrections R i ).The slope of linear curve at each small panel yields -b in relation to the intrinsic attenuation.The linear curve has an intercept of W representing source related terms at the right-hand side of equation 5.

Figure 4 :
Figure 4: a) Results of the inversion of the 2014-April-09, M W -coda3.2 earthquake: Sample fits between observed and calculated energy densities in the frequency band 0.5-1.0Hz are given for 6 different stations (see upper right corner for event ID, station name, and distance to hypocenter).Note that light blue curves represent observed envelope.Smoothed observed calculated envelopes in each panel are presented by blue and red curves, respectively.Blue and red dots exhibit 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.0Hz.

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

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