Chantal Van Dinther
email: chantal.van-dinther@univ-grenoble-alpes.fr

Ludovic Margerin

Michel Campillo

Laterally varying scattering properties in the North Anatolian Fault Zone from ambient noise cross-correlations

Keywords: Wave scattering and diffraction, Seismic noise, Coda waves, Seismic attenuation, Fault zone rheology, Numerical modelling

Intrinsic absorption and scattering properties provide us with information about the physical state and heterogeneity of the Earth's crust. These properties are usually obtained by observing the energy decay of naturally occurring earthquakes, leading to sparse spatial sampling and therefore average scattering values over a large region. The present study uses ambient noise cross-correlations to analyse the energy decay and scattering properties over a part of the North Anatolian Fault (NAF; Turkey) from the continuous records of the 73 stations of the DANA temporary array in the frequency band 0.1 -0.5 Hz. The region covered by the stations has rapidly varying geological characteristics and is highly faulted around the northern strand of the NAF. We measured in the noise correlations the space-time evolution of the energy of the coda waves. We first perform measurements in separate sub-regions. The local scattering and attenuation properties are obtained by global optimization of a 2-D solution of the radiative transfer equation for surface waves. We found that the mean free path and attenuation coefficient are considerably varying Page 2 of 37 Geophysical Journal International 2 C. van Dinther et al.

laterally with strong scattering observed in the region lying along the northern strand of NAF. The optimization provides well-constrained values for the scattering mean free path on the order of 10 km in the fault region. The mean free path is much larger (>100 km) in the neighbouring regions. We compare our global observations with a phonon based Monte Carlo simulation of scattered energy in a laterally variable scattering model. These simulations confirm the large contrast of heterogeneity between NAF and the surrounding crust and provide further constraints on the lateral extent of NAF. When sources are located inside the fault zone, we find a signature of the actual non-uniform scattering properties, observed as a concentration of energy in the fault zone for a limited amount of time. This in turn suggests that lateral variations of scattering properties should be taken into account in future monitoring studies.

properties of heterogeneities, and thereby extract detailed information on the Earth's structure and composition, one needs to quantify and distinguish between absorption and scattering, both processes contributing to coda attenuation through a relation that cannot be established in general.

One phenomenological method to investigate the relationship between the observed seismogram envelopes and the spectral structure of the random heterogeneity of the Earth is based on the scalar radiative transfer equation (RTE) (e.g. [START_REF] Wu | Multiple scattering and energy transfer of seismic waves-separation of scattering effect from intrinsic attenuation-I[END_REF]; [START_REF] Wu | Introduction: Seismic wave scattering in three-dimensionally heterogeneous earth[END_REF][START_REF] Hoshiba | Simulation of multiple-scattered coda wave excitation based on the energy conservation law[END_REF][START_REF] Hoshiba | Separation of scattering attenuation and intrinsic absorption in Japan using the multiple lapse time window analysis of full seismogram envelope[END_REF][START_REF] Hoshiba | Simulation of coda wave envelope in depth dependent scattering and absorption structure[END_REF][START_REF] Sato | Synthesis of scattered energy density for nonspherical radiation from a point shear-dislocation source based on the radiative transfer theory[END_REF][START_REF] Margerin | Radiative transfer and diffusion of waves in a layered medium: new insight into coda Q[END_REF]). In the case of a multiple-scattering medium there are two important parameters that describe the heterogeneity: the scattering mean free path and the transport mean free path * . , the reciprocal of the scattering coefficient, represents the average distance between two scattering events. * is the propagation distance required for a wave to lose memory of its initial direction. The first study taking multiple scattering into account to estimate the relative contribution between scattering and intrinsic absorption was by [START_REF] Wu | Multiple scattering and energy transfer of seismic waves-separation of scattering effect from intrinsic attenuation-I[END_REF]. The work was based on a stationary transport equation. Later, time-dependent radiative transfer theory has been introduced to describe energy propagation in randomly inhomogeneous media [START_REF] Shang | Transportation theory of multiple scattering and its application to seismic coda waves of impulse source[END_REF][START_REF] Sato | Energy transportation in one-and two-dimensional scattering media: analytic solutions of the multiple isotropic scattering model[END_REF]). The first Monte-Carlo simulations using radiative transfer for envelope synthesis were developed in parallel (e.g. [START_REF] Gusev | Monte-Carlo simulation of record envelope of a near earthquake[END_REF][START_REF] Abubakirov | Estimation of scattering properties of lithosphere of Kamchatka based on Monte-Carlo simulation of record envelope of a near earthquake[END_REF]). [START_REF] Hoshiba | Numerical basis of the separation of scattering and intrinsic absorption from full seismogram envelope: A Monte-Carlo simulation of multiple isotropic scattering[END_REF] and [START_REF] Fehler | Separation of scattering and intrinsic attenuation for the Kanto-Tokai region, Japan, using measurements of S-wave energy versus hypocentral distance[END_REF] analysed the entire S-seismogram envelopes (including the ballistic wave) to measure the ratio of scattering attenuation and intrinsic attenuation quantitatively. Their method, known as the multiple lapse-time window analysis (MLTWA), has lead to a large number of studies reporting regional values of the scattering properties (e.g. [START_REF] Mayeda | A comparative study of scattering, intrinsic, and coda Q-1 for hawaii, long valley, and central california between 1.5 and 15.0 hz[END_REF] in Hawaii, Long Valley and Central California;[START_REF] Hoshiba | Separation of scattering attenuation and intrinsic absorption in Japan using the multiple lapse time window analysis of full seismogram envelope[END_REF] in Japan; [START_REF] Jin | Separation of intrinsic and scattering attenuation in southern california using terrascope data[END_REF] in Southern California; [START_REF] Carcolé | Spatial distribution of scattering loss and intrinsic absorption of short-period S waves in the lithosphere of Japan on the basis of the Multiple Lapse Time Window Analysis of Hi-net data[END_REF] using the Hi-net stations in Japan; [START_REF] Eulenfeld | Crustal intrinsic and scattering attenuation of high-frequency shear waves in the contiguous United States[END_REF] in the contiguous United States). [START_REF] Sato | Isotropic scattering coefficient of the solid earth[END_REF] provides a comprehensive review of transport mean free paths measurements from rock sample to lithospheric scales. The majority of the these studies use records of local earthquakes with magnitudes large enough for evaluating coda properties at large lapse times. This limits the spatial coverage and therefore provides average or regional values of scattering properties instead of more specific local values.

Ambient noise cross-correlations (e.g. [START_REF] Shapiro | Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise[END_REF]) offer an attractive alternative to study the attenuation parameters on a more local scale, particularly in the lower frequency bands (f < 1 Hz). From scattering theory we know that we can reconstruct the Greens function between two stations if we correlate the coda of an impulse source, such as an earthquake [START_REF] Campillo | Long-range correlations in the diffuse seismic coda[END_REF].

Furthermore, it has been shown that we can use the coda of correlations to also reconstruct the Greens functions [START_REF] Stehly | Reconstructing Green's function by correlation of the coda of the correlation (c3) of ambient seismic noise[END_REF]. This is a good indication that the coda of correlations is analogous to the earthquake coda and contains valuable information on the propagation properties.

Wegler & Sens-Schönfelder (2007) were the first to estimate the coda attenuation Q -1 c from ambient noise auto-correlations in Japan. More recently, [START_REF] Soergel | Coda-Q in the 2.5-20 s period band from seismic noise: application to the greater alpine area[END_REF] mapped Q c in the greater Alpine region using ambient noise records in the 2.5-20 s period band and found large spatial variations of coda attenuation related to the geology. [START_REF] Hirose | A passive estimation method of scattering and intrinsic absorption parameters from envelopes of seismic ambient noise cross-correlation functions[END_REF] demonstrated that cross-correlations functions (CCF) can be used to derive scattering properties by comparison with active shot data. In their study they derived average scattering properties at the Sakurajima volcano in Japan by minimising the misfit between the energy densities measured in CCFs and the predictions of scalar RTE.

Motivated by these recent results, we use the coda of ambient noise CCFs to detect and quantify possible lateral variations of scattering properties across the North Anatolian Fault Zone. Although the mapping of these spatial variations of attenuation across the fault is interesting in its own right, it may also have some important implications for the accurate location of velocity changes based on coda wave interferometry (CWI). Indeed, geological and geophysical studies (see e.g. Ben-Zion & Sammis 2003, for a review) suggest that the assumption of uniform scattering properties underlying most CWI tomographic studies may well break down near fault zones.

Our study area is the North Anatolian Fault Zone (NAFZ; Turkey). This is a seismically active fault zone formed where the Anatolian block and the Eurasian continent meet. The ∼1200 km long dextral strike-slip fault has an east-west orientation. The region around the Izmit rupture zone is the area we are considering in current work (see Fig. 1). Here the NAFZ splays into a northern and southern strand; separating the study region into the Istanbul Zone in the north, the Armutlu Block in the centre and the Sakarya Zone in the south. The Istanbul zone is dominated by marine clastic sediments and marls, overlying a crystalline basement. Another main geological feature is the Adapazari Basin, a pull-apart sedimentary basin (S ¸engör et al. 2005). The central block mainly comprises the Armutlu and Almacik Mountains, an accretionary complex of Cretaceous and metamorphosed sediments overlying a metamorphic basement [START_REF] Yılmaz | Geological evolution of the late Mesozoic continental margin of northwestern anatolia[END_REF]. In the southern part of this block, constrained by the southern strand of the NAF, one finds a small sedimentary basin that is shallower than the Adapazari basin to the north. The southernmost part of our study region is the Sakarya Terrance. It mainly consists of limestones, red beds and other marine sediments with some metamorphic rocks towards the southern strand of the NAF, such as serpentinite, all overlying a metamorphic basement [START_REF] Yılmaz | Geological evolution of the late Mesozoic continental margin of northwestern anatolia[END_REF][START_REF] Okay | Tethyan sutures of northern Turkey[END_REF]. The northern strand of the NAF will be referred to as 'fault zone'

(FZ) in this study, because it has a wider faulted zone than the southern strand, due to the additional small normal faults. The complexity of the NAF has been documented by several recent seismological studies, thanks to the deployment of the Dense Array for North Anatolia (DANA 2012). [START_REF] Taylor | Crustal imaging across the North Anatolian Fault Zone from the autocorrelation of ambient seismic noise[END_REF] and [START_REF] Taylor | Near-surface structure of the North Anatolian Fault zone from Rayleigh and Love wave tomography using ambient seismic noise, Solid Earth[END_REF] used DANA data to provide a detailed tomographic images of seismic velocities across the fault zone. [START_REF] Gaebler | Imaging of shear wave attenuation along the central part of the North Anatolian Fault Zone, Turkey[END_REF] and [START_REF] Izgi | Crustal seismic attenuation parameters in the western region of the North Anatolian Fault Zone[END_REF] used earthquake data to infer intrinsic, scattering and total attenuation in the region. The latter study is particularly relevant to ours, because it focuses specifically on DANA data. We will therefore compare our results to the findings of [START_REF] Izgi | Crustal seismic attenuation parameters in the western region of the North Anatolian Fault Zone[END_REF] throughout the manuscript.

The manuscript consists of two main parts. In the first part, we estimate the spatial variation of intrinsic attenuation Q -1 i and scattering mean free path from ambient noise cross-correlations (Section 3) using a data regionalisation approach guided by geological considerations. In the second part we validate the inferred scattering properties by direct comparison between the data and results

of Monte Carlo simulations in media with laterally varying mean free path and intrinsic attenuation (Section 4). The next section presents the dataset and the signal processing workflow. The coda decay properties of ambient noise CCF is also briefly described.

DATA PROCESSING AND DECAY PROPERTIES OF CODA WAVES

In this section, we describe the basic processing applied to the data. We show that the coda can be effectively reconstructed and perform an empirical analysis of its decay properties.

Data pre-processing

In this study we exploit ambient noise records from the Dense Array for North Anatolia (DANA (2012); see Figure 1). This array consists of 73 three-component broadband stations, 63 of which are forming a rectangle covering an area of approximately 35 by 70 km with an average inter-station distance of 7 km. The body wave reflection study by [START_REF] Taylor | Crustal imaging across the North Anatolian Fault Zone from the autocorrelation of ambient seismic noise[END_REF] and the surface wave tomography by [START_REF] Taylor | Near-surface structure of the North Anatolian Fault zone from Rayleigh and Love wave tomography using ambient seismic noise, Solid Earth[END_REF] provide convincing evidence that microseismic noise can be used successfully at frequencies lower than 1 Hz to reconstruct empirical Green's functions from autocorrelations and CCFs with the DANA array. To extract the direct Rayleigh wave and its coda with sufficient signal-tonoise ratio, we computed CCFs using the complete 18 months of continuous data which were recorded during the period between May 2012 and October 2013. The continuous data were first divided into one-hour segments. All components were down-sampled to 25 Hz and corrected for the instrument response before removing segments containing earthquakes of magnitude ≥ 2 as inferred from a local catalog [START_REF] Poyraz | New constraints on micro-seismicity and stress state in the western part of the North Anatolian Fault Zone: Observations from a dense seismic array[END_REF]. A spectral whitening was applied to the data between 0.01 and 1 Hz, followed by 4th-order zero-phase Butterworth filtering in the 0.1-0.5 Hz frequency band. One-bit normalisation is the last step of the pre-processing and was applied to remove any remaining transient signal. We computed the full cross-correlation tensor between all pairs of stations using 1h windows of pre-processed data. For each station pair, the results from all windows were subsequently stacked to obtain the mean CCFs over the full acquisition period.

Fig. 2 shows an example of the resulting CCFs for the averaged horizontal component pairs (for all 9 separate component pairs, see Fig. S1 in Supplement). From the first arrivals we derive a velocity, which is approximately equal to 2.1 km/s, indicating that the main energy pulse is most likely composed of surface waves. Henceforth, we will assume the coda to be mostly composed of scattered surface waves. For the short inter-station distances, it is hard to differentiate between Rayleigh and Love waves. This observation is confirmed by the surface wave tomography [START_REF] Taylor | Near-surface structure of the North Anatolian Fault zone from Rayleigh and Love wave tomography using ambient seismic noise, Solid Earth[END_REF], who show that the average velocities of Love and Rayleigh waves are very similar over the study area in the 0.1-0.5 Hz frequency band. To verify convergence of the calculated CCFs we compare in Fig. 3 the envelopes of CCFs derived from 12 months of stacked data versus 18 months of stacked data. In the coda window the differences are negligible and we are thus confident that our CCFs have converged.

The later part of the envelopes however, are different, showing that the noise levels vary spatially across the study area.

Measurement of Coda attenuation

Q -1 c
A first estimate of seismic wave attenuation in the Earth's crust underneath the DANA array, for the frequency band 0.1-0.5 Hz, can be empirically obtained by measuring the energy decay of coda waves. [START_REF] Akinci | Attenuation of coda waves in Western Anatolia[END_REF] used earthquake records in higher frequency bands to derive coda-Q in Western Anatolia. Following [START_REF] Aki | Origin of coda waves: source, attenuation, and scattering effects[END_REF], we assume that the energy envelopes of noise CCFs obey the same algebro-exponential decay as earthquake data:

E(f, t) = S(f ) exp[-2π t f Q -1 c (f )] t -α (1)
where E(f, t) is the mean-squared energy envelope at lapse time t around frequency f , S(f ) is a frequency-dependent factor combining the virtual source magnitude and the site effect at the station and α is an exponent to be discussed below. E(f, t) is obtained by applying a moving smoothing window, with a length of 16 times the central period, to the squared CCF. This smoothing is sufficient to remove rapid fluctuations of the envelope while preserving the average decay properties. To further improve the stability of the measurements, we subsequently average all four horizontal components of the energy envelopes of the cross-correlation tensor after a normalization of each term has been performed at a fixed lapse-time of 100 s. In Eq. ( 1) the exponent α is a fixed parameter that depends on both the regime of scattering (from single scattering to diffusion) and the dominant wave type in the coda (body waves or surface waves). It typically varies between 1 and 2. Assuming that the ballistic and coda waves are mostly composed of surface waves, we choose α = 1. It is worth noting that in 2-D isotropically scattering media, the theoretical algebraic decay of energy is of the form t -1 in both the single-scattering and multiple-scattering regimes [START_REF] Paasschens | Solution of the time-dependent Boltzmann equation[END_REF]. Like in previous studies (e.g. [START_REF] Soergel | Coda-Q in the 2.5-20 s period band from seismic noise: application to the greater alpine area[END_REF], we estimate Q -1 c directly from the slope of the Log-Energy decay of the envelope against lapse time t with a linear least-squares method. We emphasize that the choice of α changes the absolute values of Q -1 c , but does not affect the conclusion on lateral variation on Q -1 c . We note, nevertheless, that the quality of the coda reconstruction is not uniform over the network. In particular, stations located in the south of the array show generally lower signal to noise ratios (SNR). To avoid underestimation of the decay rate due to noise contamination, we decided to adapt the duration of the measurement window to the quality of the CCF. The coda typically starts 25 s after the ballistic arrival and ends when the SNR drops below 5. We reject all CCFs for which (1) the coda duration is less than 75 s or (2) the correlation coefficient of the linear regression R 2 is less than 0.75. This provides us with good quality estimates of Q -1 c for the relatively early part of the coda. About ∼55% of the constructed CCFs pass the selection criteria, leading to 1328 Q -1 c measurements. We note that this number of Q -1 c measurements is for a maximum interstation distance of 35 km (see next section for details).

Mapping of lateral variations of

Q -1 c
To facilitate the discussion of the results, the inter-station measurements of Q c in the 0.1-0.5 Hz frequency band have been converted to 2-D maps. The study area is first discretized onto a grid of (∼8.5 km-by-12 km) pixels. At each pixel, we identify all the inter-station paths of total length smaller than 35 km that propagate through and record the corresponding values of Q -1 c . We then estimate the we show that maps derived from both causal and anti-causal parts of different component combinations are very similar. Therefore, we conclude that there is no obvious bias entailed by the possible noise directivity. Fig. 4 illustrates that Q c is relatively uniform over the entire study region. In particular, the fault zone does not distinguish itself from the surrounding crust. From this negative result, we may be tempted to conclude that the attenuation properties are uniform across DANA. This is not the case, however, as indicated by the direct observations of energy propagation presented below.

Laterally varying propagation properties

To give a direct image of the propagation properties of the region, and to detect potential local differences in properties, we have represented the energy distribution at different times deduced from the envelopes of the correlation functions. We considered virtual sources located in two very different geological contexts. The first source is located at the extreme south of the network (station DA01, Fig. 5A) while the second source is located in the immediate vicinity of the northern branch of the fault (station DA07, Fig. 5B). We measured the energy on non-normalised CCFs between these stations and all others stations of the array. The energy is averaged for all horizontal components in 8 different time windows of 30 seconds, between 5 and 75 seconds. The energy snapshots are presented in Fig. 5 after temporal correction by the average Q c . For a source south of the fault system, Fig. 5A does not show any spatial pattern after the initial flow of energy from the source. The energy distribution exhibits a speckle-like behaviour, without indication of any energy being concentrated at a particular location. On the contrary, when the source is close to the main fault, Fig. 5B shows that the energy is not spread uniformly but that it is higher inside the fault zone, especially in the eastern side, than outside the fault zone for lapse times between 25 and 65 s. The absence of such a pattern for the source in the south (Fig. 5A) indicates that the energy concentration observed in Fig. 5B along the northern strand cannot be explained by a local effect such as an amplification due to a shallow structure. One may notice that the region where we find the energy concentration is close to a kink of the main fault where intense fracturing is expected [START_REF] King | Speculations on the geometry of the initiation and termination processes of earthquake rupture and its relation to morphology and geological structure[END_REF]) and observed (Fig. 1). A dense fracture network would result in a strong scattering strength for seismic waves. No clear energy concentration is found for the scenario with a virtual source in the vicinity of the southern branch of the fault zone. These direct observations suggest specific propagation characteristics along the northern branch that are not revealed by Q c mapping. In the following, we rely on the MLTWA [START_REF] Fehler | Separation of scattering and intrinsic attenuation for the Kanto-Tokai region, Japan, using measurements of S-wave energy versus hypocentral distance[END_REF][START_REF] Hirose | A passive estimation method of scattering and intrinsic absorption parameters from envelopes of seismic ambient noise cross-correlation functions[END_REF] to evaluate quantitatively the scattering and absorption properties in different sub-regions.

MAPPING OF ATTENUATION PROPERTIES

Transport model and inversion strategy

As stated in introduction, the attenuation properties of the medium as quantified by the intrinsic attenuation Q -1 i and the scattering mean free path , are related to Q c via an unknown non-linear function.

Although the previous section did not reveal significant variations of coda-Q over the region, it does not imply that the scattering properties may not be variable. Since Q -1 c especially relies on the assumption of a simple linear decay of the Log-Energy in the time domain, it may be insensitive to finer details of the spatio-temporal distribution of energy, particularly at small spatial scales such as shown in Fig. 5. To assess the possible contrasts in scattering properties between the fault and the surrounding crust, we choose to split the study area into four sub-regions where we measure the scattering properties separately. We will refer to the sub-regions as North, 'Fault Zone' (FZ), Centre and South, following a self-explanatory naming convention (see also Fig. 1). A similar data regionalisation was employed by [START_REF] Izgi | Crustal seismic attenuation parameters in the western region of the North Anatolian Fault Zone[END_REF], but these authors did not distinguish between fault zone and centre. The definition of each sub-region results from a compromise between the following criteria. A sub-region should be: (1) small enough to ensure some uniformity in the geology and therefore in the scattering properties;

(2) large enough so that the aperture of the sub-array enables reasonable estimates of the scattering properties as further discussed below. We note that the selected stations in each sub-region are colour-coded in Fig. 1. Only station pairs for which both stations belong to the same sub-region are taken into account.

To estimate the local scattering parameters we perform a MLTWA as originally proposed for earthquake data by [START_REF] Fehler | Separation of scattering and intrinsic attenuation for the Kanto-Tokai region, Japan, using measurements of S-wave energy versus hypocentral distance[END_REF] and recently applied to ambient noise CCFs by [START_REF] Hirose | A passive estimation method of scattering and intrinsic absorption parameters from envelopes of seismic ambient noise cross-correlation functions[END_REF]. In the case of the DANA data, we measure the total energy of the four horizontal components of the CC tensor in four 15-s long time windows starting at 5 s, 25 s, 50 s and 75 s after the ballistic arrival. The choice of the first time window starting at 5 s is due to high noise levels around t =0, which is most plausibly caused by incoming teleseimic P-waves. The energy in a late time window of 15 s duration starting at lapse-time t =100 s is used to normalise the measurements for the magnitude of the virtual sources and the site effects. The time window adopted for normalization purposes generally depends on the dataset. In our case, a lapse-time of 100 s has been adopted to guarantee that the SNR is high enough and the wavefield sufficiently diffuse. The observed normalised energy densities (NEDs), E obs , are then averaged in 2-km-wide distance bins, in order to avoid bias towards one specific distance. The spatio-temporal distribution of NEDs in the four sub-regions is shown in Fig. 6. It is apparent that the spatial energy decay in the first time window is much faster inside the FZ than outside, which suggests a contrast of attenuation properties between the FZ and its environment. To confirm this interpretation, we infer the values of the scattering mean free path and the intrinsic quality factor Q i in each sub-region by comparing the observed NEDs to the predictions of a 2-D radiative transfer equation (RTE) applied to the transport of surface waves [START_REF] Sato | Energy transportation in one-and two-dimensional scattering media: analytic solutions of the multiple isotropic scattering model[END_REF]; Paasschens 1997):

P (r, t) = e -ct/l-ωt/Q i 2πr δ(ct -r) + 1 2πlct 1 - r 2 t 2 c 2 -1/2 exp l -1 ( c 2 t 2 -r 2 -ct) -ωt/Q i H(ct -r) (2)
In Eq. ( 2), P (r, t) represents the energy density at hypocentral distance r and lapse time t, for energy traveling at velocity c. The symbols δ(x) and H(x) represent, respectively, the Dirac delta function and the Heaviside step function. The first term on the RHS of Eq. ( 2) represents the direct wave contribution and is non-zero only at t = r/c. The second term models the diffuse energy forming the coda for t > c/r. Note that the dependence of the intrinsic quality factor Q i and the scattering mean free path (and in turn of the energy density P ) on the central frequency of the signal (f = ω/2π) is implicit. We assume dominance of surface waves in the direct wavefield and therefore that the coda is dominated by surface wave energy, even though it is known that there is conversion of surface wave to body wave energy and vice versa in the coda (for more details we refer the interested reader to [START_REF] Margerin | A scalar radiative transfer model including the coupling between surface and body waves[END_REF]). Even though the scattering mean free path and intrinsic absorption values would be different when taking into account this mode conversion, the focus of this study is on the first order effect of the non-uniform scattering properties, thus the assumption of surface wave dominance is sufficient.

To quantify the agreement between the observed and simulated NEDs, we introduce the following misfit function, with SR as the sum of squared residuals:

SR = 4 i=1 M j=1 log E RT E (t i , j) E obs (t i , j) 2 (3)
where E RT E denotes the energy density predicted by the 2-D transport model and the indices i, j refer respectively to the time window and the hypocentral distance bins. Note that the same normalization procedure is applied to the observed and modeled energy densities. In Eq. ( 3) the logarithm ensures that all epicentral distances and lapse-time contribute equally to the misfit. Physically, the effect of scattering is most pronounced in the early lag times. For this reason we chose not to down-weight the SR according to the variance of the measurements in Eq. ( 3), which would be conventional in a least-square fitting approach. To find the optimal value of mean free path and intrinsic attenuation, we perform a similar grid search as in [START_REF] Hirose | A passive estimation method of scattering and intrinsic absorption parameters from envelopes of seismic ambient noise cross-correlation functions[END_REF] and [START_REF] Fehler | Separation of scattering and intrinsic attenuation for the Kanto-Tokai region, Japan, using measurements of S-wave energy versus hypocentral distance[END_REF]. The search range for and Q -1 i is the same for all sub-regions with varying between 5 and 300 km with increments of 1 km and Q i varying between 60 and 200 with increments of 2. Before discussing the inversion results

in the next section we briefly recall why separation of scattering and absorption properties is made possible by MLTWA.

The basic ideas were presented in [START_REF] Fehler | Separation of scattering and intrinsic attenuation for the Kanto-Tokai region, Japan, using measurements of S-wave energy versus hypocentral distance[END_REF] perturbation located off the direct ray), we refer to [START_REF] Mayor | Sensitivity of coda waves to spatial variations of absorption and scattering: radiative transfer theory and 2-D examples[END_REF]. The fact that different time windows of the signal have quantitatively distinct sensitivities to elastic and anelastic perturbations is the key to MLTWA.

Inversion results

: absorption (Q -1 i ) and scattering (l). of the NAF we also performed the inversion, but no significant difference in neither the spatial decay rate nor the resulting and Q i were observed compared to the neighbouring central and southern

In the neighbouring sub-regions, Q i is well-constrained and the best-fitting values are slightly higher than inside the fault zone with Q i ≈ 116, 90, 106 for the north, centre and south respectively. [START_REF] Izgi | Crustal seismic attenuation parameters in the western region of the North Anatolian Fault Zone[END_REF] suggest slightly larger spatial variations of Q i at 0.75 Hz than found from ambient noise data. Our values of Q i are generally not consistent with the estimates of Q c , as we would have expected the reversed trend with higher Q i values in the south than in the north. This discrepancy may again be purely a consequence of the generally more complicated envelope shape of the data than the simple parametric form of the decay underlying Q c measurements. The misfit contours of Figure7

show that the mean free path is less well constrained outside the fault zone than inside. This is not surprising, as the mean free path ( ≥ 40km) appears to be of the same size or larger than the largest linear dimension of the sub-network. The 20% error range indicates may range from 40 km to over 300 km (300 km is the maximum value tested) in the normal crust. These values are typical of what is found worldwide [START_REF] Carcolé | Spatial distribution of scattering loss and intrinsic absorption of short-period S waves in the lithosphere of Japan on the basis of the Multiple Lapse Time Window Analysis of Hi-net data[END_REF][START_REF] Eulenfeld | Crustal intrinsic and scattering attenuation of high-frequency shear waves in the contiguous United States[END_REF][START_REF] Sato | Isotropic scattering coefficient of the solid earth[END_REF].

Effect of Velocity Model and Noise

To asses the robustness of the optimization method and to better understand the poorer constraint on we explore the effect of the surface wave velocity and of the 'noise' in the data. In each sub-region a different velocity is used, as derived from the estimated propagation time of the first arrival. Similar to the Rayleigh wave tomography results of [START_REF] Taylor | Near-surface structure of the North Anatolian Fault zone from Rayleigh and Love wave tomography using ambient seismic noise, Solid Earth[END_REF], we found the highest velocities in the centre and the lowest in the north. The actual values of the velocities, however, may be slightly different from the ones we derived via the first arrival estimate, which affects the optimization. Although marginally better fits (smaller SRs) are obtained when using lower velocities in all sub-regions, there is no significant effect on the inverted scattering properties (for more details see Supplement 3.1).

Additionally, fluctuations in the NED measurements, or 'noise', affect the optimization and the resulting scattering properties as well. This is especially the case for zones with slow energy decay with distance in combination with a short array aperture. This aspect is also explored in greater details in Supplement 3.2. The main findings from this section are that (1) Q i is always well resolved because Q i controls the typical energy ratio between the different time windows and (2) may only be well resolved when the variation of the NED in the first time window is sufficiently large over the aperture array. Otherwise, the inversion can only provide a lower bound for the mean free path.

COMPLETING THE CYCLE: COMPARING OBSERVATIONS WITH ENERGY TRANSPORT SIMULATIONS

To better understand the physical processes that play a role in the spatial and temporal energy propagation in a medium containing non-uniformly distributed scatterers, we perform Monte-Carlo simulations of 2-D energy transport in a medium with possibly variable and Q i . Technical details of the numerical implementation are presented in Appendix A. Comparing the observations with the results from the simulations provides us with insights about the scattering process in the NAFZ.

To facilitate the discussion, we compare two types of models: i) uniform models, where the full space has the same scattering properties everywhere (further details may be found in Supplement 4.1), and ii) non-uniform models, where space is divided into four different sub-regions including a fault zone, each with different scattering properties. The motivations for this non-uniform model are as follows (1) as discussed in section 3.2, the optimization is performed region by region, consequently there is no guarantee of global agreement between the data and the model, and (2) as we have observed in Fig. 6, the fit between the data and the model is not perfect, suggesting that there is room to improve the model of lateral variations. To maximise the ability to compare the simulations with the observations, we use the same dimensions in the simulations as in the study region. the corresponding sub-region. An example of this configuration is shown in Fig. 9. For these simulations we perform again the MLTWA and compare the results with the observations. Fig. 8 shows that for all sub-zones and all time windows the NEDs are greater in the simulations than in the observations, especially for the fault zone and the centre at early times. Taking simply the parameters derived in four uniform models and combining them to one non-uniform model seems insufficient to explain all observations in our study area. We note that the simulated NED for the t = 5 s window in the north and south appear to be constant or even increase with distance in certain epicentral distance ranges. This is related to the more complex wave propagation in a non-uniform model compared to a uniform model. One of the phenomena at work is back-scattering as shown in Fig. 12 and Fig. 13, and will be discussed in more detail in Section 4.3. We have seen in the previous sections that (i) we constrain inside the FZ but it proved to be more difficult outside, and (ii) there is a clear difference in the decay rate with distance between inside and outside the fault zone (Fig. 6). Consequently, the main focus of this section will be on explaining the observations in the fault zone. To simplify the analysis we assume uniform scattering properties outside of the fault zone from here onward (see Section S4.2 in Supplement for test on uniform Q i ). km, and outside the FZ uniform scattering properties for all sub-zones (Q i =100 and = 150 km). The difference in the models between panels C and D is the width of the FZ, 14 km for (C) and 5.5 km for (D). Although the results for the wide FZ have the correct energy decay rate with distance, using a narrower zone with strong scattering yields an almost perfect match to the observations. This is in line with the findings of the tomography studies in the region (e.g. [START_REF] Kahraman | Crustal-scale shear zones and heterogeneous structure beneath the North Anatolian Fault Zone, Turkey, revealed by a high-density seismometer array[END_REF][START_REF] Taylor | Crustal imaging across the North Anatolian Fault Zone from the autocorrelation of ambient seismic noise[END_REF], suggesting a fault zone that is perhaps not wider than 7 km down to the mantle.

C D B A C D Qi =80 ℓ=11km Qi =116 ℓ=299km Qi =106 ℓ=299km Qi =90 ℓ=225km C D

The signature of a finite width scattering zone

In the case of a structure with a fault zone of material associated with high scattering taken between regions of low scattering, the energy decay with distance depends not only on the values of , as we saw in the numerical tests of the previous section, but also on the width of the FZ. Fig. 11 shows the simulation results for fixed values of (10 and 150 km) and varying the FZ width W . We compare the behaviour of non-uniform models with that of uniform models associated with the two values of , 10 and 150 km. Note that for the non-uniform models and the uniform model with =150 km, the energies were normalised with a 15 s long time window starting at a lapse-time of 100 s, as was done in the data analysis. A first observation is that in our time and distance configuration, the non-uniform models have very similar behaviours at 50 s and almost identical at 75 s. For the latter time the results are very similar to those of the uniform model with =150 km. A first simple argument to explain this observation is that the area of the strong scattering zone becomes small in comparison with the total area covered by the scattered waves (Fig. 9). The results presented here have been normalised at 100 s, to be consistent with the MLTWA analysis, but we have verified that the absolute long time energies for the non-uniform and uniform models with = 150 km are close (see Fig. S16 in Supplement). Their relative difference follows approximately a scaling in W /t, deriving from the same simple geometrical argument. We thus observe a convergence towards the solution of a uniform model with =150 km for long time in all cases. These results thus indicate that it is difficult to detect a heterogeneity from the observation of temporal decreases at long times, which explains why our characterization with Q c yielded grossly similar results inside and outside the fault zone. The most important differences among the models are observed for the short times. Let us consider the case of the window starting 5 s after the arrival of the direct waves, which is critical in our comparison with the data. For this time window we have plotted in Fig. 11 in dotted line the uniform model decay with =10 km, the value in the band representing the fault zone. Contrary to the other models and the data, the energies were not normalised with the 100 s values for this uniform model. Based on our previous conclusions, the normalization was done with the values of the external uniform model ( =150 km) towards which the solutions of the non-uniform models tend for long times. This normalization is carried out only for interpretation purposes, while approaching the data analysis conditions.

Indeed, we know that the solutions converge towards the uniform case with =150 km, but that they show differences at finite times (see Fig. S16 in Supplement for non-normalised results). Since the FZ widths considered are small in relation to the distances travelled by the waves, and the sources are located near the FZ boundary, a significant portion of the energy escapes rapidly from the central FZ and will not be subject to intense scattering. Thus, a significant shift between the curve of the normalised =10 km uniform model and the non-uniform models is observed even at short distances.

For all widths, this shift increases with distance. Remembering that the absolute time of the energy measurement increases with distance (as the lag time is added to the ballistic wave arrival time), this can be explained by the decreasing role played by the central FZ as the total surface occupied by the diffracted waves grows, leading to the convergence of the non-uniform models towards the uniform one with =150 km for long times. We note, furthermore, that all non-uniform models exhibit a decay with distance that is larger than the one of the uniform =10 km model.

Fig. 11 illustrates that the amplitude and decay of the normalised energies vary significantly with W . In the case where W = 5.5 km, that is smaller than in the FZ, a significant shift is observed between the uniform model = 10 km and the other models at the shortest distance. The curve exhibits first a rapid decrease with distance then approaches the curve of the uniform model = 150 km, a convergence which is observed with the curves of the longer times. The amplitude is explained by the small size of the strongly diffracting FZ (5.5 km) with respect to the radius of the region sampled by the waves. Even at the shortest distance, with the time window considered (here between 6.5 and 21.5 s), the propagation of the waves is predominantly in the low-scattering environment. Hence, the position of the source in the FZ must also be taken into account when comparing amplitudes at short distances for the different geometries.

For the other values of W , the same evolution can be observed globally. For t = 5 s, the curves gradually separate to approach the asymptote at distances increasing with W . The same behaviour is observed for the other lapse-times. The evolution of the curves illustrates the transition between two limit models: strong scattering in the FZ for short times and weak for long times. Note that our configuration does not allow the first regime to be fully observed because the widths of the fault zone are small. For the windows with lapse-times larger than 5 s, the differences between the non-uniform models diminish rapidly with time or distance and the constraints on the model must be searched in a rather short time range but over a range of distances that covers the characteristic dimension/width of the scattering heterogeneity. If we consider an observation at a fixed distance, e.g. 3 km, the amplitudes for W from 10 to 20 km are very similar for t = 5 s (the waves sample the same part of the model), much more separated for 25 s (at different convergence stages), then again very close for larger lapsetimes (close to the asymptote). This numerical example shows the need for the observation of multiple lag windows at multiple distances to decipher the global view of the energy behaviour, and also the need for a complete modelling to relate it to the model parameters and W . In particular, the effects of finite size of the diffracting FZ cannot be simply expressed by a single attenuation operator for each pair (W , ) that would apply at all times.

Back-scattered Energy from the Fault Zone

In this section we pursue our investigation of the effects of a strong, narrow band of scattering on the energy evolution in space and time. The model we use is the one that seems most optimal from the previous section (with a FZ width, W , of 5.5 km and = 10 km inside, and = 150 km outside the FZ). The configuration for these simulations is different than previously. The model space is again nonuniform but the receiver line is north-south oriented. Snapshots of the resulting energy distribution, normalised energy envelopes as a function of lapse-time and the simulation configuration are shown in Fig. 12A-C, respectively. The most striking feature we observe on both the time-domain envelopes and the snapshots is a 'bump', which can be interpreted as back-scattered energy from the fault zone.

When the ballistic energy reaches the fault zone, a front of diffuse energy propagates back towards the direction of the source (Fig. 12A). Fig. 12B shows the same feature, where for stations farther away from the FZ it is easier to distinguish between ballistic energy and back-scattered energy. As we can see when the station is closest to the FZ, both back-scattered energy and ballistic energy arrive at similar times. At a lapse time of the order of 50 s, we observe an inhomogeneous distribution of energy associated with the lateral variations of scattering properties that progressively disappears at long lapse-times (Fig. 12A). If this phenomenon of strong back-scattering energy is also present in the actual data, mixture of these energies can have profound implications for the coda decay and scattering properties analysis. Unfortunately, identification of these bumps in the actual CCFs proved to be more difficult than expected.

In the map view of Fig. 13 we show four sections (A-D) perpendicular to the fault zone, with a 'virtual source' roughly in line but as far as possible from the stations to accommodate a move-out of this potential back-scattered energy. Sections A-C have permanent station KO07 (grey star) and section D has permanent station KO06 (blue star) as their virtual source. The simulation (Fig. 12B) discussed above is most similar to section B of Fig. 13. Additionally, the colour codes for the envelopes are kept similar. It is expected for the back-scattered energy to arrive later in the stations farther away from the fault zone. Unlike the predictions from the simulations, we fail to observe back-scattered energy in the data (Fig. 13). The simple model with a strong narrow band of scattering (Fig. 12) seems to make predictions that are impossible to discern in the data. This suggest that the fault zone is far more complex than the simple conceptual model proposed in this paper. Other medium properties such as the density and seismic velocity may be significantly different inside the fault zone versus outside. This can cause a high acoustic impedance contrast and therefore trap energy within the fault zone [START_REF] Ben-Zion | Characterization of fault zones[END_REF]. In our simulations we have not taken other medium properties than intrinsic attenuation and scattering mean free path into account. It also seems reasonable to think that there are variations of scattering strength along the strike of the fault.

Fault zone coda energy concentration

As described in Section 2.4, when a virtual source is inside the northern branch of the fault the scattered energy appears to localise inside the fault at a distance larger than 20 km and for lapse times between virtual source is in the southern most part of the study area (Fig. 5A). We observe that the relative energy spot is slowly propagating away from the source. This suggests that the energy concentration could be related to the diffuse energy in the fault zone where the small mean free path implies a slow expansion with respect to the surrounding medium. To verify this further we performed an additional simulation, where the source is put in the middle of the fault zone, and the receivers are placed in a half-circle over the fault zone with the virtual source at its centre. The source-receiver distance is therefore constant, 24 km, and the total of 49 receivers are placed with 1 km increment of distance to the median line of the fault zone. This allows for relative energy analysis without the source receiver distance dependence in the region around the fault zone. The scattering properties (uniform Q i =100, =10 km inside and =150 km outside the fault zone) and width of the fault zone (∼5 km) are similar to the optimal model described in the previous section. Fig. 14A shows clearly the diffuse energy concentrated in the fault zone. At 35 s, the energy is actually minimum in the region of the fault, although a small secondary maximum is visible in the centre of the fault zone. For this early time, the dominant effect is the strong scattering attenuation in the fault zone while the diffuse energy in the fault is relatively weaker than in the surroundings.

At 50 s the diffuse waves have a higher energy in the fault zone than in the surrounding region. This dominance lessens at 100 s when the scattered wave field is increasingly homogenised. The negative halo just outside the fault zone of Fig. 14A cannot be observed in the data (Fig. 5B). This may be because the data is not sampled densely enough, with the array having a 7 km inter-station spacing.

Another contributing factor could be the oversimplification of our model. The actual fault zone may have a more gradual change in scattering properties than our model that has an abrupt large scattering contrast. A comparison of the energy as a function of time for one receiver outside the fault zone 1 ) and one inside the fault zone (r 2 ) is shown in Fig. 14B, where we can observe in the inset that the energy density of r 2 exceeds that of r 1 after 40 s. The location of the receiver outside the fault zone is chosen close to the fault zone but outside the negative halo. The simulation is reproducing the observation of a concentration of energy in the fault zone during a limited period of time starting much later than the passing of the ballistic waves and their possible local resonances. This effect, and particularly the fact that it is late and transient, cannot be explained by simple wave amplification effects and is therefore the signature of the actual non-uniform properties of scattering. The stronger scattering attenuation towards the eastern part of the fault zone reported by Izgi et al. ( 2020) is yet another observation supporting the idea of a non-uniform strong zone of scattering. Although we see similar behaviour between the data and the model, the match between the two is far from perfect quantitatively, likely due to the complexity of the fault zone and surrounding geology. While our model certainly oversimplifies the reality, it offers an explanation for the late concentration of energy, which is consistent with the geological observation of a narrow heterogeneous fault zone. In the southern branch of the NAF we did not find any evidence of a similar transient local increase of relative energy. This is in with the findings of Section 3.2 i.e. the lack of significantly different scattering properties for this branch of the fault zone to the surrounding southern region. The absence of both concentrated energy and distinct scattering properties in the southern branch may be explained by the smaller damage zone, e.g. due to the lower activity of this fault (e.g. Aslan et al. 2019, and references therein).

CONCLUSIONS

This paper demonstrates that a high scattering limited fault zone is a reasonable first order model for the North Anatolian FZ in the area of the Izmit rupture zone. A narrow fault zone of ∼5.5 km width (± 50 %) with a scattering mean free path of the order of 10 km, surrounded by a medium with 150 km scattering mean free path, has been obtained by analysis of ambient noise cross-correlations in the 0.1-0.5 Hz frequency band. We verified that our model reproduces a concentration of energy in the fault zone when the source is inside, similar to what we observe in the data. Currently, ambient noise monitoring methods rely on the assumption that scattering properties are the same across fault zones.

Here, we have demonstrated that this is not the case. Therefore, laterally varying scattering properties need to be taken into account in future monitoring efforts.

From a methodological perspective we can conclude that one can use the coda of correlations for attenuation studies. Although we know that scattering properties vary between different types of crust, e.g. between volcanic and normal crust, variations on such short scale have not been reported yet. Our results have direct implications for applications such as monitoring velocity changes in the Earth's crust. In particular, our study suggests that a good estimate of the spatial variations of scattering properties is key to correctly recover the temporal variations of seismic properties. Future improvements of our scattering model of the NAFZ should be focused on incorporating lateral variations along the fault zone, because the data revealed that there is even more complexity than described by our model.

A larger array, with an extension towards the east, is necessary to capture this complexity. on the right-hand side is over all the directions in the plane. We remark that if e(t, r, k) solves Eq. (A.1), it also solves the following equation where the effective scattering mean free time τ e is constant: where δ(k, k ) is the delta function on the unit circle, τ nc (r) is the scattering mean free time for the null or delta scattering events and τ -1 e = τ (r) -1 + τ nc (r) -1 . In effect, all that we have done is to add the same term τ nc (r) -1 e(t, r, k) on both sides of Eq. (A.1). But the benefit is in fact immense because in the fictitious medium with null collisions, the scattering mean free path is constant, so that the most basic implementations apply straightforwardly. Note that there is a certain degree of arbitrariness in the definition of τ e . Assuming that τ is bounded in the domain of interest, we may simply take τ -1 e = sup τ (r) -1 .

(∂ t + ck • ∇ + τ -
In practice, these ideas may be very easily implemented as follows. The transport process is simulated in a medium where the mean free path e = cτ e is constant, so that the simple exponential probability law described above still applies. At each collision point r c , we must decide whether a true or a fictitious scattering event occurs. To carry out this task, we select a uniformly distributed random number in the interval ]0, 1[ and compare it to the local scattering conversion rate. If: < τ (r c ) -1 /τ -1 e (A.3) then a genuine scattering occurs and we proceed as usual with the selection of a new propagation direction. Otherwise, we have a null collision and the propagation direction is unchanged up to the the next collision. The main computational overload comes from the fact that we are in fact simulating more scattering events than actually occur in the real world. But this is largely counter-balanced by the considerable simplification of the free path length selection as well as the absence of complicated geometrical tracking of the particle.

  A) Map of the study region showing the DANA network (triangles) and the main faults (red; from Emre et al. (2018)). The selection of stations for the four sub-regions used in the MLTW analysis are indicated by the coloured triangles; triangles with a coloured outline belong to two sub-regions as indicated in the legend. The array is composed of 6 columns of stations, labeled A to F from west to east. Each column comprises 11 stations, numbered 1 to 11 from south to north. B) Geological map of the the wider study region (adjusted from Taylor et al. (2019) and Akbayram et al. (2016)). The Adapazari and Pamukova basin, are indicated by AB and PB, respectively.

Figure 2 .

 2 Figure 2. Bin-stacked cross-correlations between all pairs of North-components recorded at DANA in the 0.1-0.5 Hz frequency band. The dominant arrivals show energy propagating at velocities that are typical of surface waves.

Figure 3 .

 3 Figure3. Envelopes of ambient noise cross-correlations obtained by stacking 12 months of data (red dotted line) versus 18 months of data (solid black line). As a reference, the grey solid (dashed) lines show the energy averaged over all four station pairs for 18 (12) months stacked data, respectively. The two dashed vertical grey lines show the coda window, preceded by ballistic waves and followed by noise. The four panels correspond to a representative sample of station pairs across the array; A and B show the envelopes for a station pair in the north and centre, respectively. C and D both show the envelopes for station pairs in the south.

Figure 4 .

 4 Figure 4. A) Q c map derived from the causal part of the energy envelope of CCFs averaged over horizontal components. It shows the arithmetic mean of Q c for all rays that cross the cell with a minimum of 5 rays. B) Ray coverage C) Relative standard deviation in percents.

Figure 5 .

 5 Figure 5. Snapshot of seismic energy distribution derived from ambient noise cross-correlations. The energy values have been normalised for display purposes only. Assuming an average coda-Q of 150 in the study area and an average frequency of 0.3 Hz, the energy per time window is multiplied by e 2πf t/Qc . A) The virtual source is located in the south at station DA01. B) The virtual source station DA07 (depicted by a yellow star) is located inside the fault zone. The colour map indicates the normalised energy. The moving black box highlights the concentrated energy in the fault zone, slowly propagating towards the eastern side of the fault zone.

Figure 6 .

 6 Figure 6. Normalised energy density (NED) in four sub-regions located inside the FZ (B) and outside the FZ (A, C, D). Individual markers show the observed NED. The lines represent the NEDs modelled by the 2-D RTE using the best fitting parameters indicated at the top of each plot. The errorbars represent ± 1 STD of the individual measurements per bin. The colours correspond to four time windows, with lighter colours for increasing lapse times (5 s, 25 s, 50 s, and 75 s after the ballistic waves).

Figure 7 .

 7 Figure 7. Residuals and optimal values for the grid search in the 4 sub-regions of NAF. A), C) and D) for the regions outside the FZ, north, centre and south respectively; B) for inside the FZ. The green colour indicates the values of normalised sum of the residuals (SR; as in Eq. 3), with the darkest colour for the smallest misfit. The values are normalised w.r.t. the minimum SR. The best fitting values are indicated by the red crosses ('X').

Figures 6

 6 Figures 6 and 7 show the results of the optimization procedure for the scattering mean free path and the intrinsic absorption in each sub-region. The comparison of the predictions of the best-fitting model with the data indicates that a simple 2-D RTE is sufficient to capture the general spatio-temporal energy distribution across DANA, provided that different attenuation parameters and velocities are

Figure 8 .

 8 Figure 8. Comparison of spatio-temporal energy evolution between observations (markers) and simulations (lines) for all sub-regions: A) north, B) FZ, C) centre and D) south. The colours indicate the different 15 s-long time-windows, starting at 5 s, 25 s, 50 s and 75 s respectively. The results correspond to the initial non-uniform model, based on the scattering properties derived from the regionalisation approach and schematically shown in the lower right corner.

4. 1 Figure 9 .

 19 Figure 9. Example of the configuration for the simulations. Depending on the specific case, uniform or nonuniform, the model is divided in sub-regions where a source and an east-west receiver line is placed with receiver spacing 3.5 km. The specific configuration shown above is used for the non-uniform case with different fault zone widths. The different widths are realized by shifting the southern boundary (dashed red line). The fault zone boundary towards the north (right side) is fixed, as is the position of the receiver line w.r.t. the northern boundary and the source.

Figure 10 .

 10 Figure 10. Comparison of spatio-temporal energy evolution between observations (markers) and simulations (lines) for the FZ only. The upper two panels correspond to simulations in uniform models with A) = 10 km and Q i = 100, B) = 150 km and Q i = 100. The lower two panels are for a simple model with =10 km and Q i = 80 inside the FZ, and = 150 km and Q i = 100 outside the FZ. C) for a wide fault zone, of 14 km and D) for a narrow FZ of 5.5 km. The colours indicate the different 15 s wide time windows, starting at 5 s, 25 s, 50 s and 75 s respectively.

Fig. 10

 10 Fig.10shows four panels with observation and simulation results inside the FZ. We first discuss the simple uniform models: panel A for a model space where Q i =100 and = 10 km, and panel B where Q i =100 and = 150 km. Both models do not match the observations to a desirable degree, especially in the first time window. For the model with strong scattering (A), the spatial energy decay in the first time window of the simulations is similar to the observations but the values themselves

Figure 11 .

 11 Figure 11. Comparison of spatio-temporal energy evolution of simulations using uniform and non-uniform models with different FZ width. The dashed and dotted lines show NEDs for uniform models with =150 and 10 km, respectively. The NEDs of the non-uniform models are indicated by the continuous lines and markers, for different widths of the fault zone, W , as indicated in the legend. The colours indicate the different 15 s wide time windows, starting at 5 s, 25 s, 50 s and 75 s respectively.

Figure 12 .

 12 Figure 12. A) Snapshots of energy propagation, showing the effect of a strong scattering zone. The graphs show the energy density as function of distance (from right, north, to left, south) at different times. The red dashed lines indicate the fault zone. B) Normalised energy density versus time for four different receivers (7 km apart) in the same simulation. C) Simple non-uniform model used for the simulation exhibiting a narrow fault zone, with W = 10 km, = 10 km inside the FZ, = 150 km outside the FZ and Q i = 100 everywhere. The receiver line is oriented north-south, in line with the source in the north, 15 km from the first receiver (yellow star). In panels B and C, colours are used to encode the distance between the station and the fault zone.

Figure 13 .

 13 Figure 13. Observations of energy envelopes for station pairs oriented perpendicular to the fault zone. The same colour code as in Fig. 12B-C is used to indicate the distance between the stations and the fault zone. The map view shows sections A-D. For section A-C, the virtual source is the grey star (permanent station KO07). For section D, the virtual source is the blue star (permanent station KO06). The distance to the virtual sources are indicated on the map. Note: envelopes have been re-scaled with respect to the distance, d, between virtual source and station (20 × d) to enhance the visibility of energy move-out, i.e. the envelopes are shifted on the log-scale to separate them.

Figure 14 .

 14 Figure14. Simulation of energy entrapment in the fault zone for a non-uniform model with a source inside the fault zone. A) On a linear scale the energy versus distance from the centre of the fault zone is shown, before (t =35 s ), during (t =50 s), and after (t =100 s) energy entrapment. The black and red circles denote station r 1 and r 2 respectively, as indicated on panel C. B) On a linear scale the energy versus time at two receivers, r 1 (outside the FZ; black solid line) and r 2 (at the centre of the FZ; red dash-dotted line). The inset shows the energy density ratio between the two receivers (r 2 /r 1 ) in the time window indicated by the grey box. C) Configuration of the simulation. Black triangles represent the 49 receivers, the distance to the middle of the FZ between the stations increments by 1 km. All stations are at equal distance from the source (yellow star).

  1 e )e(t, r, k) = τ (r) -1 2π p(k, k )e(t, r, k )dk + τ nc (r) -1 2π δ(k, k )e(t, r, k )dk , (A.2)

Geophysical Journal International

ACKNOWLEDGMENTS

We thank Frederik Tilmann and an anonymous reviewer for their useful comments. We thank Cécile Lasserre and Marco Bohnhoff for useful discussions. We acknowledge the support from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No 742335, F-IMAGE). This study uses the continuous recording data of the Dense Array for Northern Anatolia (DANA), which is part of the Faultlab project, a collaborative effort by the University of Leeds, Bogazii University Kandilli Observatory and Earthquake Research Institute (BU-KOERI) and Sakarya University. For the earth quake catalog we would like to thank Niyazi Turkeli.

We remind the reader that Monte-Carlo simulations designate a class of stochastic methods of solution of radiative transport equations. In this approach, energy is represented by discrete seismic phonons who undergo a sequence of collisions upon which their direction of propagation changes according to the laws of scattering. In a medium with uniform scattering properties, the distance d between two collisions (also called free path length) is governed by an exponential probability law of the form p(x < d < x + dx) = -1 exp(-x/ )dx with the scattering mean free path. For greater details, we refer the reader to the literature as summarized in [START_REF] Shearer | The global short-period wavefield modelled with a Monte Carlo seismic phonon method[END_REF] and Sato et al. (2012). To simulate the transport of seismic energy in media exhibiting spatial variations of scattering and absorption, we have employed the so-called method of null or delta-collisions [START_REF] Lux | Monte Carlo particle transport methods: neutron and photon calculations[END_REF]. This is an exact method of simulation that maps the transport process from a medium where attenuation properties vary spatially onto a medium where they are constant. Let us remark first that, in the framework of Monte-Carlo simulations, absorption may be treated as a scattering process that reduces the energy of a particle by a factor B, with B the local value of the albedo. Hence, to get a grasp on the method it is sufficient to treat the case of non-uniform scattering properties, which is also most relevant to our applications.

Consider a 2-D transport equation with a spatially varying mean free path (r). The equation governing the transport of the energy is (e.g. Paasschens 1997):


M Olteanu

¢¡ ¤£

J Rynkiewicz
email: rynkiewi@univ-paris1.fr

B.Maillet ¨ § ¥¡ ¦£

Non-linear Analysis of Shocks when Financial Markets are Subject to Changes in Regime

Violent turbulences are often striking the financial markets and an Index of Market Shocks (IMS) was recently introduced in the attempt of quantifying these turbulences. Regime switching linear models have already been used in modelling the conditional volatility of returns. In this paper we propose a description of the IMS with hybrid models integrating multi-layer perceptrons and hidden Markov chains. After sudying the prediction performance of these models, we focus on the series separation and the index behaviour subject to the hidden states.

Introduction

The financial "crashes" or "crisis" we hear speaking of in the media, are translated on the markets by large price movements. Whether it is the Great Depression started in 1929, the one day 22.6% drop of US stocks in October 1987 or the World Trade Center attacks in September 2001, the markets have experienced powerful turbulences since their creation. Until recently, there was no way of quantifying these crisis or measuring how large they were. Traditional measures of risk, such as the implied or the historical volatility, have some drawbacks in being considered as suitable tools.

Following the initial approach of Zumbach and al (2000), [START_REF] Maillet | An Index of Market Shocks Based on Multiscale Analysis[END_REF] introduce the Index of Market Shocks (IMS), an easily computable measure that quantifies the importance of market movements. The index was designed in the attempt of having a scale for measuring the extent of financial crisis, and also as a tool for comparing the intensity of strongly agitated periods. Briefly, its construction is a financial translation of the Richter scale for measuring earthquakes intensity. The first step is to compute volatilities, from the highest to the lowest frequency available. Then, a Principal Components Analysis is performed and an increasing function of the "dissipated energy" (the volatility, in this case) is applied to the resulting factors.

The aim of this paper is to seek for a good method of modelling and forecasting the Index of Market Shocks values. The basic idea is to find a hidden Markov chain -neural net hybrid model that would provide, besides good data fitting and forecasts, a state separation that describes the market behaviour in normal times and in crisis periods.

Autoregressive Hidden Markov Chain Models

Before giving the results, we make a brief description of the autoregressive models with hidden Markov chains.

Financial series historics, taken over long periods, show important breaks in the series behaviour. These breaks are due to several reasons, such as bankruptcies and panics on the markets, changes in governments policies or wars. One way to model a time series, while taking into account these "regime" changes, is to use hidden Markov chains.

Let us consider ¢¡ ¤£ ¦¥ §£ © the observed time series and let £ ¥ §£ © be a homogeneous Markov chain defined by its state space "! $# ¡ &% $' (' )' (% # &0 21 , 3 54 76 98 and the 3 A@ B3 transition matrix C with D FE HG I QP R S£ ¢T ¡ VU XW S£ Y aF¥ , U % b dc % e' )' (' )% 3 . If we suppose, with no loss of generality, that the chain state space is the canonical basis of f 0 and we note g £ ¢T ¡ 7 S£ ¢T ¡ ih V S£ ¢T ¡ W & S£ ¥ , then an autoregressive hidden Markov chain model has the following form:p

S£ ¢T ¡ 7C q S£ sr tg £ ¢T ¡ ¡ £ ¢T ¡ vu xw $y ¢¡ £ £ © S T ¡ ¥ sr Sw y £ ¢T ¡ where ¡ £ £ © eT
¡ defines the vector ¢¡ ¤£ © eT ¡ % e' )' (' )% ¡ ¤£ §¥ , u xw $y 4 ! $u % $' (' )' )% u ¦ i1 is an autoregressive function of order , w $y I4 ! % e' )' (' )% ¦ d1 is a real strictly positive number, e £ ¦¥ §£ © are independent, identically distributed f g h % c ¥ .

Hamilton's study [START_REF] Maillet | An Index of Market Shocks Based on Multiscale Analysis[END_REF] on the GNP (Gross National Product) growth is one of the first applications of the HMC models, where the series is supposed to behave as a mixture of linear autoregressive models. But one can also consider nonlinear autoregressive functions, such as the multilayer perceptrons and thus, we get the so called hybrid MLP-HMC models (multilayer perceptron -hidden Markov chain models).

Research of a MLP-HMC Model

The data used in the analysis are daily values of the IMS, computed for the MSCI (Morgan Stanley Capital International) World Equities Index, from June 25th, 1998 until July 25th, 2002 (1066 observations).

We are investigating whether the series is "regime changing" and whether we can distinguish between different market behaviours by using a hybrid MLP-HMC model. The linear model and the one hidden-layer perceptron that were onsidered in a preliminary study showed that the significant lags were 1, 2, 4, 5 and 6. We fixed this input vector and we let the number of experts and hidden units vary up to three.

In the end, two configurations were selected and further investigated, on the basis of two criteria : the architecture with the smallest mean squared error and the one with has the "best" transition matrix. The "best" transition matrix is chosen such that its trace divided by its dimension is the closest to one and thus, for the corresponding configuration, we have the most powerful segmentation of the series. These "empirical" criteria were chosen because of the impossibility of using the usual likelihood ratio test (one of the regularity conditions is not fulfilled, that is the nonsingularity of the information matrix).

Prediction Performance

The first step of this study was to compare the estimated MLP-HMC model with the results of a linear ARMA and of a multi-layer perceptron.

After investigating all ARMA(p,q) models, for p and q in a specified range, an AR(6) minimizing the BIC criterion was selected. The research of the multi-layer perceptron is done using the REGRESS software [START_REF] Rynkiewicz | Modèles hybrides intégrant des réseaux de neurones artificiels à des modèles de chaînes de Markov cachées : applications à la prédiction des séries temporelles[END_REF], which performs an automatic research through all possible configurations with a hidden layer. The "best" model is chosen to minimize the BIC criterion and the "Statistical Stepwise" [START_REF] Cottrell | Neural Modeling for Time Series : a Statistical Stepwise Method for Weight Elimination[END_REF] algorithm is used to eliminate the non significant connexions, once the "dominating" perceptron was found. The final selected model has two hidden units, twelve parameters and the input vector goes up to the sixth lag of time.

The summary in the next table shows that, for this series, there is no interest in considering MLP-HMC in forecasting purposes as there is no improvement when compared to the linear model or to the perceptron. The global mean squared error is computed on the whole series (1060 observations), while the training and test errors are computed when splitting the series (800 values and, respectively, 260).

Model

Global Next, we will focus on studying the state separation of the series. Since the transition matrix does not suggest an interesting separation and since the table above implies overfitting, we will not further investigate the model which minimizes the mean squared error and we will only consider the models which maximize the trace of the transition matrix. We chose to study a two-state model and a three-state model.

A Two State Model

For the two-state model, we have the following estimated transition matrix :

C ¡ h ' ¢ h ¤£ ¤h ¥£ h ' c §¦ h ¥¤© h ' h ¢ c ¢ ¤ h ' £ ¤© ¢ ¥ ©
The first expert has three hidden units with an associated noise variance of 0.43045, while the second is a linear model (0 unit) with an associated variance of 0.05637.

In Figure 1, we plotted the conditional probabilities of the first expert for Russia, Brazil and United States. We note that they are close to one when the market passes through a crisis period.

The first example concerns the russian crisis in August 1998 for which the main causes were the strong ruble devaluation and the country debt default, combined with the echos of the Eastern Asia crisis in 1997. On the graphic, one can see that the peaks correspond to the key moments of the crisis : August 11th, the stock market is at its lowest level, while the short-term bond market crashes, August 17th, the Central Bank widens the ruble/dollar exchange rate corridor and a 90-day moratorium on foreign debt payments is announced, August 27th, Russia's main foreign exchange market, the MICEX, shuts down every currency trade.

The next example concerns the Brazilian crisis in January 1999, mainly due to the important fiscal deficit and the high inflation rate. Here too, we see the decisive moments : January 7th, the governor of Brazil's third biggest state announces a 90-day moratorium on the debt payement to the federal government, January 13th, the finance minister resignes and the real/dollar currency band is lowered, January 15th, the Central Bank abandons its defense of the real, despite the holding of important foreign reserves. Finally, the last example we propose here is the September 2001 crisis. Even before the terrorist attacks, the american economy showed signes of recession. On the graphic, we see the two moments : September 6th, the stock prices drop due to the bad news on the economy health, September 13th, the delayed opening of the markets after September 11th.

A Three State Model

The research of a three state model having the most powerful segmentation of the series selected an architecture with the following estimated transition matrix :

C v ¡ ¢ h ' £ ¤£ ¤© ¨h ' h h ¦ © ¤£ h ' h ¥ £ ¤£ h ' h h ¤£ © ¤© h ' £ h ¥£ Sc h ' c c §¦ ¥ h ' c h ¤£ h Sc h ' c ¥£ ¦ ¤Sc h ' £ © c £ c ¦ §
The first expert is a three hidden units perceptron, while the two others have one hidden unit. The associated estimated variances are 0.42687, 0.04422 and 0.275537.

Studying the results for this model is less obvious than for the two-state one. We chose to describe the behaviour of the shock index subject to the three experts by a Kohonen map classification. We considered as input variables the IMS and the conditional probabilities of the three perceptrons. In Figure 2, the first two graphics contain the variables means (normalized values) in the clusters, which are homogeneous. We see that the high values of the index (strongly agitated periods on the markets) are associated with the first expert, while the small ones correspond to the second one, and this opposition is obvious on the map. As for the intermediary values of the index, we may expect them to go with the third experts. This is only partially true : we have, indeed, mean values of the index for which the conditional probability of the third expert is the largest, but also cases of mixture between the first and the third or between the second and the third expert.

In Figure 3 and4 we have the results of an hierarchical classification performed on the map. If we keep three clusters, we note that the index values evoluate through clusters, from small to high values, the first and the last cluster are homogeneous and correspond, respectively, to calm periods and to crisis moments. The middle cluster is less homogeneous because it integrates the intermediate values that are associated to the third expert, but to the mixtures of experts mentioned above as well. This is even more obvious if we cut the classification tree at five clusters : we do not change the first and the last cluster, but the middle one is splitted according to different expert combination.

An immediate application of the Kohonen classification would be constructing a portfolio selection strategy, by identifying the calm periods and the turbulences on the values of the computed IMS.

Conclusion

Using a hidden Markov chain model for this kind of financial series does not improve the prediction results obtained with classical linear or non-linear models. The interest in using them is related to the state separation they provide. Thus, one may, on one hand, identify key dates in market evolution and, on the other hand, derive strategies for portfolio selection using the experts conditional probabilities and the index values.

Figure 1 :

 1 Figure 1: Conditional probabilities of the first expert (Russia, Brazil, U.S.)

Figure 2 :

 2 Figure 2: Kohonen map for the three state model (IMS,e1,e2,e3)

Figure 3 :

 3 Figure 3: Clustering of the Kohonen map in three clusters

Figure 4 :

 4 Figure 4: Clustering of the Kohonen map in five cluster