Frequency-dependent body-Q and coda-Q in Karlıova Triple Junction and its vicinity, eastern Turkey

* Correspondence: alperdemirci@comu.edu.tr


Introduction
One of the key considerations in earthquake hazard-based studies of a seismically active region is knowing the correct seismic quality factor which controls the mechanism of attenuation and amplitude decay of the corresponding seismic energy that propagates through the medium with distance. The commonly-used attenuation unit is a dimensionless quality factor (Q), which expresses the decay of wave amplitude during wave propagation in the medium (Knopoff, 1964), explained by the ratio of seismic energy lost (ΔE) to harmonic total energy (E) per cycle of vibration: (Q= -2πE / ΔE). Since being introduced by Aki and Chouet (1975), many descriptive and evaluation studies of seismic wave attenuation, serving several purposes, have been performed on different tectonic environments worldwide. Besides providing crustal information about the corresponding region (Hoshiba, 1993;Bianco et al., 2002), the general use of a quality factor in seismic attenuation satisfies the need for assessing earthquake hazard (Pulli, 1984) and the calculation of local magnitudes, spectral source parameters, and synthetic seismograms (Abercrombie, 1997).
Lower Q values indicate strong dissipation of seismic energy in a regional sense. Theoretically, an infinite (or excessive) value of the attenuation factor corresponds to almost fully elastic behavior in the medium (Lay, 2015). Hence, it can be clearly understood that tectonically stable regions are represented by high Q values, whereas regions having high seismicity show relatively lower Q values. As well as large-scale tectonic movements and crustal deformations such as warping, folding, and faulting (Kumar et al., 2014), the regional Q may be also influenced by the chemical compounds of units, crustal impedance, and temperature changes as well as fluid saturation in the medium (Hauksson and Shearer, 2006).
Determination of Q has been widely performed via different wave phases (P, S, Lg, and Coda) observed on seismograms and titled Qp, Qs, QLg, and Qc (Yoshimoto et al., 1993;Kim et al., 2004;Sharma et al., 2008;Ma'hood et al., 2009). Direct phases of P and S waves are propagated through a path between the source and seismic station, but coda waves are commonly accepted as back-scattered seismic waves that form as a result of the interaction between several heterogeneities and body waves (Herraiz and Espinosa, 1987) and include information about deeper parts of the lithospheric medium (Aki and Chouet, 1975). In other words, while Qp and Qs, determined using direct body waves, address shallow crustal medium attenuation characteristics, Qc contains information from the deeper lithospheric medium.
The coda normalization method (CNM) developed by Aki (1980) can be performed using one station for multievents, or vice versa, in order to estimate the Qs factor by normalization of direct S waves to coda waves. Yoshimoto et al. (1993) reported that direct P waves were also suitable for the technique and named it the extended coda normalization method (ECNM). ECNM is based on the fact that spectral amplitudes of the direct phases of P, S, and coda waves display a distinct decay rate by the increasing source-station distance. Moreover, recent developments related to the propagation and amplitude characteristics of coda and body waves have led to an increase in attenuation studies of related seismic phases through the Earth's lithosphere, resulting in a better knowledge of attenuation parameters in various regions in the world having different tectonic and geological features. Studies carried out in different tectonic environments using the ECNM have both spread the use of the method worldwide and also allowed different tectonic features to be compared globally. Thereupon, ECNM was used by many researchers such as Yoshimoto et al. (1993) for Kanto, Japan;Ma'hood et al. (2009) for East Central Iran; Padhy (2009) for Bhuj, India; Kumar et al. (2014) for Kinnaur, Himalayas; Mukhopadhyay et al. (2016) in the Aswan reservoir, Egypt; and Castro et al. (2008Castro et al. ( , 2009 in volcanic regions of the Apennines (Italy) and Sonora (Mexico). These studies were conducted on volcanic or highseismicity environments in crustal-scale and addressed earthquakes at a maximum hypocentral distance of 160 km in high frequency ranges (up to 63 Hz).
Coda attenuation is also commonly represented through the frequency-dependent Qc factor. The single backscattering model (SBM) of coda waves (Aki and Chouet, 1975) is used to estimate Qc values. The main aspect of the model is that coda wave trains are backscattered S waves. Coda waves are assumed to consist of single back-scattered waves; in other words, it is commonly accepted that S waves are scattered only once due to the heterogeneity of the crust and upper mantle. The assumption is that the dimensions of the randomly-distributed scattering elements are larger than the wavelength. This was explained in detail by Wu and Aki (1988), who classified the perturbation of wave velocity as a function of the dimension of the scatterers in terms of various propagation mediums identified by the wavelength and intensity of heterogeneity in the region. They also concluded that coda wave studies examining the scattering produced by heterogeneities vary from 0.1 km to 10 km in the crust and upper mantle.
Studies related to coda attenuation have been attempted by numerous researchers throughout Turkey and worldwide. Aki and Chouet (1975) Akıncı and Eyidoğan (1996) for the Erzincan region in eastern Turkey, Sertçelik (2012) for the East Anatolian Fault Zone (EAFZ) and its five subregions, Sertçelik and Güleroğlu (2017) on the easternmost segments of the North Anatolian Fault Zone (NAFZ), and Rahimi et al. (2009Rahimi et al. ( , 2010 on Mt. Sabalan and the Alborz region in Iran. These studies noted the relationship between the Qc values they obtained and the degree of heterogeneity in the corresponding area. In this study, Qp, Qs, and Qc variations and frequencydependent functions for each factor were evaluated using a dense dataset of local earthquakes recorded in the vicinity of Karlıova Triple Junction (KTJ) (Figure 1). The major geological and tectonic features of this region were briefly discussed. KTJ is one of the most active regions in terms of seismicity in Turkey and has been subject to numerous hazardous earthquakes in the historical and instrumental periods (Ambraseys and Finkel, 1995). Thus, the obtained frequency-dependent attenuation functions, which constitute the main output of the present study, may play a key role as one of the required input parameters in seismic risk assessment, early warning (Lior et al., 2015) ground motion modelling (Akıncı et al., 2010) and site attenuation (kappa) modelling (Van Houtte et al., 2018) studies. These functions have also been frequently used for the exploration of geothermal resources and oil and natural gas reservoirs (Klimentos, 1995;Jyothi et al., 2017) in recent years.
The study area was divided into three subregions in consideration of the tectonic branches of the KTJ, namely, the EAFZ, NAFZ, and Varto Fault Zone (VFZ). ECNM and SBM were performed to derive the body wave attenuation factors (Qp and Qs) and coda wave attenuation factor (Qc), respectively. Waveform data were extracted for 111 regional earthquakes that occurred between January 2010 and December 2017 and were recorded by six threecomponent broadband seismic stations operated by the Kandilli Observatory and Earthquake Research Institute (KOERI) network 1 . The local magnitude, epicenter coordinates, and event depths were used according to calculated parameters by the routine processes performed by the KOERI network 2 . The focal depths were found to be shallower than 10 km and the local magnitudes were between 2.8 and 6.1. After a detailed analysis of the digital seismograms, taking into consideration sufficient coda length, a total of 1164 high-quality waveforms were processed.

Tectonic and geological outline of the study area
Eastern Anatolia, having a complex set of dynamic and tectonic regime characteristics, stands out among the most significant instances of a continental collision on Earth. The westward propagation of the Anatolian block, also resulting in intense seismicity, probably started 12 Ma ago (Dewey et al., 1986;McQuarrie et al., 2003) due to the collision of the Arabian and Eurasian plates along the Bitlis-Zağros thrust zone. The tectonic escape of Anatolia, starting from KTJ, is accommodated by the NAFZ and EAFZ (Şengör, 1980). Reilinger et al. (2006) reported that the eastern part of the Anatolian plate, or just the western part of the KTJ, had a differential motion with a velocity of 20 mm/year with respect to the Eurasian plate.
The NAFZ, which forms the northern border of the westward-propagating Anatolian block as a first-order tectonic structure, continues as a convex line represented by a right lateral strike-slip movement. It is approximately 2000 km long, starting from the KTJ and continuing to approximately 100 km south of the Black Sea coastline (Armijo et al., 1999) (Figure 1a). Earthquake focalmechanisms also show right lateral transtensional motions around the Yedisu Segment, situated in the western part of the KTJ or eastern part of the NAFZ (Figure 1b), which is taken into account in this study.
Although it has also strike-slip deformation characteristics and a similar seismicity rate, the EAFZ is not as well-known as the NAFZ (Aksoy et al., 2007) Figure 1. a) Location of the area investigated (modified from Bozkurt, 2001), b) shaded relief map, active faults (Şaroğlu et al., 1992) and focal mechanism solutions of some events (Kalafat et al., 2009) in the region, c) geological map of KTJ and its surroundings (modified from Dilek and Sandvol (2009) Figure 1a).
To the east-southeast and the third branch of the KTJ, the 50-55 km continuation of the NAFZ is called the Varto Fault Zone (Gürboğa, 2016). This continuation is approximately parallel to the Yedisu segment of the NAFZ, but the fault mechanisms represent a dominantly reverse-thrust character (Figure 1b). It contains four or more subparallel segments within the fault zone. Şaroğlu (1985) reported that the Varto Fault intercepts the Varto Caldera to the east of Karlıova (Figure 1c). The structural data obtained in six different fault segments and splays are significant proof that the kinematic background of the fault was subjected to a series of shortening and expansion regimes (Karaoğlu et al., 2017). Additionally, the end of the VFZ is also regarded as the termination point of the NAFZ (Gürboğa, 2016) (Figure 1). Two earthquakes (Mw = 6.8 and Mw = 6.2) caused heavy damage to the VFZ on 19 and 20 August 1966 (Ambraseys and Zatopek, 1968). Büyüksarac and Bektaş (2018) reported that during the 1966 Varto events, earthquake intensity decreased rapidly, especially in areas north of the epicenter. They pointed to the heavy damage in a village located 13 km northeast of Varto, whereas in another village 23 km northeast with similar construction characteristics, hardly any damage was observed. These findings may relate to the high attenuation capacity of the Varto region.
More recent geological studies showed some of the primary characteristics of structures located between the EAFZ and NAFZ and the effect of these faults on the main fault zones. Additionally, recent modelling studies have pointed out complex faulting (secondary faults between the NAFZ and EAFZ, and EAFZ and VFZ) due to strain partitioning in the vicinity of the KTJ (Şengör, 2014;Sançar et al., 2015;Zabcı et al., 2015;Seyitoğlu et al., 2019). Due to the spatially based approach of the method used, a number of different station-source distances are required for investigation on a regional scale. Hence, the effects of these small-scale secondary faults have been ignored so far, on the assumption that secondary faults in the proximity of main fault zones have similar slip characteristics.
An overview of the geological context of the KTJ and its circumjacent region ( Figure 1c) shows that all tectonic structures located between the EAFZ and NAFZ were active since Plio-Quaternary in the region. The East Turkish High Plateau (ETHP) (Dewey et al., 1986), surrounded by the Bitlis massif to the south and eastern Pontides to the north, occupies a large part of eastern Anatolia ( Figure 1c). The Bitlis massif is tectonically underlain by a late Cretaceous-early Tertiary ophiolitic melange in the south and this melange also tectonically overlies the sedimentary deposits of Arabian foreland (Dilek and Sandvol, 2009). The ETHP has the basement units of ophiolites and ophiolitic melanges (Dilek and Sandvol, 2009). A large part of the NAFZ branch of the KTJ is located on the ETHP. The KTJ and its surroundings, including many Plio-Quaternary volcanic activities that occurred after the initial Arabia-Eurasia collision between 8 and 6 Ma (Dilek and Sandvol, 2009), is mostly covered by volcanic rocks (Figure 1c). These volcanic formations consist mainly of intercalated lavas with subordinated ignimbrites and sedimentary layers, aged from 6.9 + 0.9 to 1.3 + 0.3 Ma (Innocenti et al., 1982;Keskin et al., 1998;Karaoğlu et al., 2016). Karaoğlu et al. (2018) also revealed the presence of magma chambers located at shallow depths between 2 and 5 km below the Varto volcano following field observations and analytical results. A common opinion is that the strato-volcanoes are highly heterogeneous structures. The unconsolidated volcanic deposits, accumulated over the years and led to high seismic impedance contrast, such as lahars, ash falls, lava, and pyroclastic flows, may be attributed to the increasing effect of scattering (Wegler and Lühr, 2001) and intrinsic attenuation (Prudencio et al., 2017). Considering the fact that the dense spatial distribution of volcanic units with high attenuation capacity and presence of the Varto volcano in the study area also increase heterogeneity, especially for the Varto Fault which intercepts the Varto Caldera; it can be emphasized that the outputs of the analyses are largely dependent on geological factors.

Data
A dataset based on six broadband stations recording 111 small-to moderate-sized local seismic events throughout the area between coordinates 39.50 °E to 42.00 °E and 38.00 °N to 40.00 °N was analyzed. The study area was divided into three subregions as the best way to represent each fault zone in terms of earthquake distribution. These subregions are being represented by different epicentral distribution ellipses and are shown in Figure 2. Meirova and Pinsky (2014) suggested the use of earthquake data having an epicentral distance smaller than 120 km in order to bypass aliasing effects due to the additional energy caused by surface waves which appear immediately after the arrival of the S phase. Considering that the waveforms may be affected in whole or in part by surface waves, the maximum station-event distance was therefore chosen so as not to exceed 120 km. The local magnitude range is between 2.8 and 6 (90% of events satisfy 3 < M L < 5 and 6% have M L ≥ 5) with shallow (<10 km) focal depths (Table 1).
The average Vs and Vp velocities for the crustal medium were taken as ~3.3 km/s and 5.8 km/s respectively in the analyses (Vanacore et al., 2013).
The stations in the study were selected to represent each of the tectonic elements which form the KTJ area. The dataset of the KARO station at the intersection of the fault zones was also used for the analysis of all three regions. Moreover, the VRTB station in the VFZ, the YEDI and ERZN stations in the NAF region, and the data of the PTK and BNGB stations in the EAF region were used in addition to the KARO station. Therefore, while the NAF and EAF were represented by three stations, the VFZ region was analyzed using two stations ( Figure 2). The data of broadband seismometers (Güralp CMG-40T) have a flat velocity response between 0.05 and 20 Hz and were recorded with a 0.01 s sampling rate and suggested Nyquist frequency of 50 Hz. Considering the Nyquist frequency and instrumental response of the data used, and the frequency content of regional seismic waves, analyses were performed up to a maximum frequency of 18 Hz.

SBM for estimation of Qc (coda-Q)
The SBM of coda waves (Aki and Chouet, 1975) was used to estimate Qc values. The observed coda wave amplitude on a seismic record consists of varied components such as the source, site, instrument, and medium effects. The decay rate of coda amplitude in a specific frequency range is due to medium effects including geometrical spreading and attenuation but is independent of source, site, and radiation effects (Aki, 1969).
The coda amplitude of the filtered seismogram for a lapse time (t) in a given frequency range is associated with the Qc parameter and is expressed as follows: ( 1) where S(f) indicates the source factor of the coda wave at a corresponding frequency independent of the time and radiation pattern. The geometrical spreading factor (α) was taken as 1 (Sato and Fehler, 1998) because it concerns the backscattered body waves. In this case, Equation (1) can be reshaped after logarithmization as: ( 2) and (3) By calculating the slope (slp) of the linear least-square fit of , the Qc value can be determined easily by Equation (4) below: (4) Using an eight-pole Butterworth filter with two-thirds bandwidth, the root-mean-square (RMS) amplitudes of the band-passed data were calculated by applying a moving time window of 1 s for all central frequencies.
Additionally, in order to avoid loss of energy in the RMS amplitudes due to a decrease in the number of data points, relatively narrower moving time window lengths were  chosen at higher frequencies. Pulli (1984) emphasizes that local events allow the determination of coda decay at lapse times from 25 to 60 s. Hence, the logarithmic amplitudes vs. t for four different lapse time lengths (20, 30, 40, and 50 s) were plotted and a least-squares approximation was performed to calculate the slope of the fitted straight line for corresponding central frequency. Data representing a signal-to-noise ratio of greater than 3, and a best-fit line having a correlation coefficient of more than 0.8 for each central frequency, were both taken into account. Signalto-noise ratios were calculated by dividing the RMS amplitudes of the 10-s signal data portion in the middle of the coda window into the portion before the first arrival of P waves as noise. The lapse time windows starting at twice the travel time of the S-wave onset were used here (Aki and Chouet, 1975;Rautian and Khalturin, 1978;Aki, 1980) ( Figure 3). It is an important challenge of extracting P wave coda part based upon the fact that P coda waves could be overlapped by S waves at low source-receiver distances; thus, in such studies, S coda waves are generally preferred rather than P coda waves (Sertçelik and Güleroğlu, 2017). Therefore, the frequency dependence of the quality factor (Q c ) that satisfies the power-law approximation with a function of , where is the quality factor at 1 Hz and n corresponds to frequency dependency, was calculated separately for all subregions using S coda waves.

Qp and Qs (body-Q) estimation via ECNM
The ECNM (Yoshimoto et al., 1993) was used to estimate Qs and Qp in the crustal medium. This method is based on the fact that the proportionality between the coda, S, and P wave spectral amplitudes shows a decay rate independent of the epicentral distance to the receiver. Due to the normalization process to the coda amplitudes, the instrument, site, and source effects are also removed.
Qp and Qs can be calculated by using earthquake data at a varied range of epicentral distances. Due to the spatially based approach of the method, the main criterion is the number of events for single station analyses in seismically active regions; whereas, in weak seismic zones, the required data for analysis can only be acquired by a large number of seismic stations located at different distances from the event epicenter. The following equation was used to calculate the quality factors (Q p,s ) for the corresponding seismic wave phase: where Ap,s (f,r) is the half values of peak-to-peak amplitudes of the direct P or S wave transformed for 1.28-s window length and r indicates the source-receiver distance. Studies using a single horizontal component pointing out the amplitude similarity (e.g., Tusa and Gresta, 2008;Bora and Biswas, 2017) were taken into account and additional trial calculations were performed separately for two horizontal components in order to be on the safe side. It was observed that the obtained results were very close to each other. Therefore, EW components were selected for the Qs analyses and vertical component records were used for the 2010 Geometrical spreading (r α ) was taken as r -1 for epicentral distances smaller than twice the Moho depth (h m~4 5 km for the investigation area (Çınar and Alkan, 2017)). For greater distances, the geometrical spreading term (r α ) of Equation (5) took the form of (r.hm) α and α was taken as -0.5 (Herrman and Kijko, 1983;Ma'hood et al., 2009;Kumar et al., 2014;Singh et al., 2012). Qp and Qs values were calculated with the slope of the best correlated linear regression line, which represents a decrement of the coda normalized RMS amplitudes of P and S wave envelopes along increasing epicentral distances. The regression process was performed on seven different frequency bands with central frequencies ranging from 1.5 Hz to 18 Hz. For the corresponding central frequency, Qp,s values were calculated using Equation (6). The regression line for each central frequency is shown in Figures 4 and 5. (6) Similar to coda-Q, using the slope of regression lines calculated at corresponding central frequencies, the frequency dependence of the quality factor (Qp,s) that satisfies the power-law approximation with a function of, where Q 0 is the quality factor at 1 Hz and n corresponds to frequency dependency, was calculated separately for all subregions using S coda waves.
Within this theoretical framework, the present study contains frequency-dependent Q relations that were !"# ! !,! inferred using Q estimates for each central frequency of both P and S waves. Naturally, the mean value of the Qs/ Qp ratios for each tectonic branch of the KTJ was also calculated. This ratio should be calculated specifically in such studies because it is a comparable parameter on a global scale and contains valuable information about the heterogeneity and degree of fluid saturation of the investigated tectonically-active region.

Results and discussion
In the scope of the present study, the Qp and Qs variations representing the attenuation characteristics for the KTJ region were determined using the ECNM. The study area was separated into three subregions considering the tectonic branches of the junction area. Qp and Qs for all the studied subregions were estimated using the fitted curve slope of logarithmic values of the RMS amplitude ratio of the associated seismic phase to the coda phase versus hypocentral distance for seven different frequency bands. The data were subjected to least-squares linear regression and best fit lines for different frequency bands, which were estimated in order to obtain Qp and Qs values and their frequency dependencies, as shown in Figures 4  and 5, respectively. The linear regression lines denote an increase in attenuation factors with increasing central frequency.
The minima and maxima of Qp were 7 at 1.5 Hz for the VFZ and 739 at 18 Hz for the NAFZ, respectively. Qs values were also observed between a minimum of 26 for VFZ and a maximum of 1259 for NAFZ (Table 2). Comparatively, the VFZ seems to have a significantly higher attenuation than the other branches of KTJ. The graphs of the frequency-dependent power-law fit, Q(f) = Q 0 f n , (where Q 0 is the attenuation quality factor at 1 Hz and n is the frequency dependency parameter) of Q P and  Table 2. Observations suggest that there is a strong relationship between quality factor and frequency. The Qp and Qs functions have a similar frequency dependency and it can be clearly seen that P-phases are attenuated more rapidly than S-phases (Qs/Qp >1) for the entire frequency range in the whole studied area ( Figures  4 and 5 and Table 2). The Qs/Qp ratio was also calculated in order to explain the observed values of Q and to reveal possible attenuation mechanisms. As a result, these ratios were found to be greater than unity (>1) at all central frequencies in the range of 1.5-18 Hz, indicating that scattering is a key attenuation phenomenon that controls the distance-dependent amplitude decrement of the body phases radiated in the KTJ region. The observation of such high Qs/Qp ratios in many regions can be associated with high scattering effects in relation to the high level of lateral heterogeneity (Hough and Anderson, 1988).
Additionally, at all central frequencies from 1.5 to 18 Hz except 3 Hz, Qs/Qp ratios observed on the VFZ are greater than the theoretical ratio of 2.41 calculated by Sato (1984); however, these ratios are less than 2.41 for the NAFZ and EAFZ over the entire frequency range. Qs/Qp ratios yield results greater than unity, as already calculated in volcanic areas (e.g., Mt. Etna -Patane` et al., 1994) and in high seismicity regions for the upper crustal medium (e.g., Bianco et al., 1999;Hough et al., 1999). Moreover, Vassiliou et al. (1982) investigated the attenuation characteristics of sedimentary rocks in a laboratory environment and concluded that partially saturated rocks have a Qs/Qp ratio greater than 1. In this respect, it may be suggested that a partially fluid-saturated upper crustal structure exists for the entire study area.
Several studies of seismic attenuation have been carried out in several volcanic areas and contributed significantly to geological and geophysical investigations. For example, Del Pezzo et al. (1986) calculated frequency-dependent Qs variation and remarked on the differences with other tectonic regions; these were attributed to the presence of magma effects. Del Pezzo et al. (2006) also revealed the attenuation effect on volcanic Mount Vesuvius using S-coda envelopes and calculated the Q value as ~10 at a central frequency of 2 Hz. Vila et al. (1997) reported low Q values, smaller than 40, obtained using P-wave dispersion analysis at Campi Flegrei, suggesting that the region has a very low Q. Additionally, the seismic attenuation structure of the volcano Rabaul in Papua New Guinea was studied by Gudmundsson et al. (2004) using regional earthquake data and the spectral ratio method to estimate the Qp and Qs factors. It was revealed that near-surface materials with low Q values (15-20) inside the caldera have comparatively more attenuation capacity than those outside the caldera. In another study, Castro et al. (2008) evaluated S-wave attenuation in a volcanic region of the Apennines (South Italy) via the distance-based spectral decay method and described high attenuation characteristics in the upper crust which they defined in the functional form of Qs = (18.8) f (1.7) ( Table 3).
The low QP and QS values obtained in this study correspond to those of seismically-active and volcanicallyactive (or previously active) regions in the world. In particular, the low Qp and Qs values (<4 Hz) obtained for VFZ are the lowest among all previous results. The fact that other studies presenting very low Q values similar to VFZ had been performed in areas with high volcanic activity (Del Pezzo et al., (1995) (Mt. Etna); Castro et al., (2009) (Sonora Volcanic Region); and Castro et al. (2008) (Apennine Volcanic Region)) suggest that the high  (Figures 6a and 6b) conducted in seismically active regions on body-Q. The results also show similar trends and similar values, in the range of 1.5-18 Hz ( Figure 6). Coda attenuation quality factor (Qc) is another significant parameter that explains the degree of regional tectonic activity (Jin and Aki, 1988). Many studies have mentioned the relevance of Q 0 , n, and seismicity in various regions. Van Eck (1988) correlated the frequency dependence of coda attenuation factor and seismicity level while Jin and Aki (2005) observed high n and low Q values in all regions identified as having high tectonic activity in Japan. Additionally, Yun et al. (2007) and Dasovic et al. (2012) presented results supporting these findings for South Korea and on the Pannonian-Dinaridic contact zone.
Although the Q values were determined by different methods, they usually produce comparable body and coda Q values (e.g., Aki, 1980). In the present work, the body-Q values obtained with ECNM are lower than the Qc values for all corresponding central frequencies. Del Pezzo et al. (1987) described these numerical differences between body-Q and coda-Q as being caused by the propagation of wave paths. Additionally, they stated that in the assumption of homogeneous half-space, the body and coda-Q values should coincide exactly, but in reality, the increase of Q values with depth and attenuation in the crust is stronger than in the mantle. Thus, comparatively high coda-Q values are thought to be due to the effect of depth, especially for longer coda durations.
In addition to the characterization of body wave attenuation, the attenuation values of the coda waves were also determined at the same central frequencies for four coda lapse time (20, 30, 40, and 50 s) windows.
According to results for the corresponding coda window lengths, the calculated average values of Q 0 for VFZ with standard deviations varied from 14.6 ± 4.3 for a 20 s lapse time window to 45 ± 11 for a 50 s lapse time window. These values are 96 ± 22 to 133 ± 38 and 53 ± 12 to 94 ± 17 for NAFZ and EAFZ, respectively. The decreasing trend of n can be expressed as 1.2 ± 0.1 to 1.0 ± 0.1 in VFZ, 1.0 ± 0.1 to 0.8 ± 0.09 in NAFZ, and 0.95 ± 0.11 to 0.86 ± 0.10 in EAFZ for a lapse time window ranging from 20 to 50 s. The results given above show the strong lapse time dependency between Q 0 and n by an increase in Q 0 and a decrease in n with ascending lapse time length (Table 4). Sertçelik (2012) and Sertçelik and Güleroğlu (2017) also reported a decreasing n value with increasing lapse time in the EAFZ and eastern part of the NAFZ. The potential impact of source effects on the coda waves can be minimized by using longer lapse times, because the Qc may be affected by the focal mechanism (Rautian and Khalturin, 1978).
For a longer lapse time, the average volume from which coda waves are recorded is also larger and the Qc determined includes the impact of a greater area. Woodgold (1994) noted that if the lapse time is higher than or equal to 30 s, significant seismic energy is transmitted to the mantle. This may indicate that Qc values are time-dependent and increase over a longer lapse time. Giampiccola et al. (2002) noted that the dependence of frequency-Qc decreases when heterogeneity decreases    ( along with depth. Sato (1984) also showed that the overall decrease in the frequency dependency (n) value with increasing lapse time points to the short-wavelength component dominance of inhomogeneity for larger lapse times; thus, large-scale inhomogeneity becomes less dominant in deeper parts. Taking into account the results of Woodgold (1994), 20, 30, 40, and 50 s lapse time analyses were performed. In addition, a 30 s lapse time analysis was selected for comparison with the Qc analyses performed in and around the study area. The results obtained are shown in detail in Figures 7 and 8.
The obtained curves of Qc vs. frequency for all three subregions present similar attenuation characteristics to previous studies of seismically-active regions in Turkey and its surroundings (Table 5 and Figure 8). It is clearly seen that attenuation characteristics of seismic waves in the KTJ region exhibit similar frequency-dependent relations to previous studies performed in the vicinity of KTJ and to seismically-active volcanic regions on the Arabian plate in Iran (Table 5). Furthermore, the frequency-Qc curves obtained by using similar frequency ranges and lapse time  Consequently, keeping in mind that volcanoes and fault zones increase the degree of heterogeneity of a crustal medium and lead to a high seismic wave scattering (Zieger et al., 2016), and considering the effect of the Varto volcano in the VFZ with its widespread volcanic deposits, it should be mentioned that a dense heterogeneity directly affects the seismic wave attenuation in such a complex environment. As a major result of the study, this can be considered as the main reason for obtaining the lowest Qp, Qs, and Qc values in VFZ out of the whole KTJ area in all frequency ranges. On the other hand, NAFZ stands out as having relatively high Q values. These high Q values for all seismic phases in NAFZ exhibit relatively low attenuation characteristics, particularly at higher frequencies (f ≥ 3 Hz) compared to other branches of the KTJ. The fact that most of the events occurred on Paleocene ophiolitic melange units and that there are less volcanic deposits than in the other regions may explain the high values on the NAFZ.

Conclusions
In the present paper, frequency-dependent functions representing attenuation characteristics based on P, S, and coda wave phases for the KTJ region were obtained based on ground motion data recorded by the regional seismic network. The obtained functions may play a key role in modelling ground motion in order to achieve better assessment of future earthquake hazards in the region. While many studies have been carried out to discover the coda wave attenuation functions in the study area and its near vicinity, there is no research which shows the attenuation characteristics of the body waves in such a detailed manner. Therefore, this paper is the first to reveal fault zone-based attenuation characteristics together with P, S, and coda wave phases. The obtained frequencydependent attenuation functions for each branch of the KTJ are as follows: Varto Fault Zone (VFZ); · Qc = (23 ± 6) f (1.14±0.11) (30 sec), Qp = (4 ± 1.4) f (1.2±0.13) and Qs = (8 ± 2) f (1.3±0.1) Eastern part of North Anatolian Fault Zone (NAFZ); · Qc = (106 ± 34) f (0.93±0.12) (30 sec), Qp = (22 ± 2) f (1.16 ± 0.1) and Qs = (43 ± 7) f (1.18 ± 0.07) East Anatolian Fault Zone (EAFZ); · Qc = (67 ± 13) f (0.92 ± 0.1) (30 sec), Qp = (25 ± 9) f (0.96 ± 0.14) and Qs = (62 ± 10) f (0.81 ± 0.06) Consequently, the obtained frequency-dependent Q functions could be used as the basic inputs for many future scientific and engineering studies, such as regional earthquake risk analyses, strong motion simulations, and the most accurate determination of earthquake spectral source parameters. Moreover, obtained Q-body attenuation models may also contribute to further researches which will discuss the site-specific attenuation (kappa) models, especially at higher frequencies. Low Q values across the study area may be associated with the abundant distribution of volcanic deposits throughout the region. NAFZ is represented by the highest Q values, while the lowest Q values were obtained in the VFZ. These low values are considered to be caused by the dense crustal heterogeneity of Varto Volcano to the north of the VFZ.

Acknowledgments
The author would like to thank the editor and anonymous reviewers for their valuable contributions. Graham Lee is also thanked for proofreading the paper. Some figures in the manuscript were created using the Generic Mapping Tool (Wessel and Smith, 1995).