# Imaging the crustal and upper mantle structure of the North Anatolian Fault: A Transmission Matrix Framework for Local Adaptive Focusing

Rita Touma, Arthur Le Ber, Michel Campillo, Alexandre Aubry

► To cite this version:

Rita Touma, Arthur Le Ber, Michel Campillo, Alexandre Aubry. Imaging the crustal and upper mantle structure of the North Anatolian Fault: A Transmission Matrix Framework for Local Adaptive Focusing. Journal of Geophysical Research: Solid Earth, 2023, 128 (11), pp.e2023JB026704. 10.1029/2023JB026704. hal-04027628v2

HAL Id: hal-04027628

https://hal.science/hal-04027628v2

Submitted on 17 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Imaging the crustal and upper mantle structure of the North Anatolian
Fault: A Transmission Matrix Framework for Local Adaptive Focusing

Rita Touma,¹,² Arthur Le Ber,¹ Michel Campillo,² and Alexandre Aubry¹

¹Institut Langevin, ESPCI Paris, PSL University, CNRS, Paris, France

²ISTerre, Université Grenoble Alpes, Maison des Géosciences, BP 53, F-38041 Grenoble, France

(Dated: October 17, 2023)
Abstract

Imaging the structure of major fault zones is essential for our understanding of crustal deformations and their implications on seismic hazards. Investigating such complex regions presents several issues, including the variation of seismic velocity due to the diversity of geological units and the cumulative damage caused by earthquakes. Conventional migration techniques are in general strongly sensitive to the available velocity model. Here we apply a passive matrix imaging approach which is robust to the mismatch between this model and the real seismic velocity distribution. This method relies on the cross-correlation of ambient noise recorded by a geophone array. The resulting set of impulse responses form a reflection matrix that contains all the information about the subsurface. In particular, the reflected body waves can be leveraged to: (i) determine the transmission matrix between the Earth's surface and any point in the subsurface; (ii) build a confocal image of the subsurface reflectivity with a transverse resolution only limited by diffraction. As a study case, we consider seismic noise (0.1-0.5 Hz) recorded by the Dense Array for Northern Anatolia (DANA) that consists of 73 stations deployed for 18 months in the region of the 1999 Izmit earthquake. Passive matrix imaging reveals the scattering structure of the crust and upper mantle around the NAFZ over a depth range of 60 km. The results show that most of the scattering is associated with the Northern branch that passes throughout the crust and penetrates into the upper mantle.

# I. Introduction

The North Anatolian Fault zone (NAFZ) is one of the major continental right-lateral strike slip faults, and forms a border between the Eurasian continent and the Anatolian block. With an extremely well developed surface expression, it is one of the most active faults in the Eastern Mediterranean region [1, 2]. It is over 1600 km long and extends from eastern Turkey in the east to Greece in the west and, historically, has been subject to many destructive earthquakes [3, 4]. The seismic activity of such large faults constitutes a continuous hazard/threat to the surrounding regions and big cities, especially Istanbul city located to the West of the fault.

Faults are well defined at the surface by the localized deformation and displacement delineating the fault traces, but their deep structure remains poorly understood [5]. The understanding of such major fault systems and seismic hazard requires a characterization of the geometrical and seismic properties of the crust and upper mantle. A large number of geological and geophysical studies have discussed the complexity of fault zones and their relation with their deep roots [4]. They
are not only confined in the mid crust; indeed, models suggest that they penetrate deep into the crust and extend to the upper Mantle [6]. If so, faults develop into shear zones, corresponding to a volume of localized deformation accounting for the relative displacement of the tectonic blocks.

Seismic imaging techniques, especially reflection, refraction, and tomographic methods, constitute a very powerful tool to characterize fault zones and report the variation of the properties of the crust and the upper mantle. They rely on the study of wave propagation inside the Earth that is governed by the density and elastic properties of the rocks. To properly probe the medium, waves should be generated by a dense distribution of seismic sources. Conventional seismic exploration techniques use either earthquakes as seismic sources, or explosions and vibrators to generate seismic waves in regions with weak seismicity. Because of the limitations in the earthquake distribution and high cost of active methods, there is a need for alternative imaging approaches that would not rely on any coherent source. In the 2000's, the extraction of deterministic information about the Earth structure from ambient seismic noise revolutionized the field of seismology [see e.g. 7]. It was shown that the cross-correlation of diffuse waves or ambient seismic noise recorded at two stations provides an estimate of the Green's function between those two stations [see e.g. 8]. The reflection response of the medium is then retrieved and can be applied to build tomographic or structural images of the Earth. Because ambient noise is dominated by surface waves, their Green's function component can be easily extracted [9]. It has been demonstrated that, under energy equipartition, body wave reflections can also be retrieved from ambient seismic noise cross-correlations [10–12]. Body waves contain valuable information on the structure of the medium in depth and can be investigated to obtain high-resolution images of the crust and the mantle [13].

Faults are usually imaged indirectly through strong velocity contrasts in tomographic profiles [14], or through the offset of geological layers observed in reflectivity images [15]. However, tomographic images exhibit a relatively poor resolution, while reflection imaging methods are strongly sensitive to the available velocity model. Interestingly, a reflection matrix approach has been recently proposed to cope with these issues. Originally developed in acoustics [16, 17] and optics [18, 19], this approach has been recently applied to passive seismology [20, 21]. By considering high frequency seismic noise (10-20 Hz), high resolution images of complex areas, such as volcanoes [22] and fault zones [21], have been obtained over a few km depth. In this paper, we aim to characterize the crustal structure of the NAFZ at a much larger scale (until 60 km depth). To that aim, a lower frequency bandwidth (0.1-0.5 Hz) has been considered. At the corresponding wavelengths, the subsurface reflectivity can be considered as continuous rather than
being seen as a discrete distribution of scatterers as in our previous works [21, 22]. As we will see,
this continuous reflectivity can be exploited to enable a local and adapted auto-focus on each part
of the subsurface image, thereby showing an important robustness to the inaccuracy of the initial
wave velocity model.

Seismic matrix imaging is based on the passive measurement of the reflection matrix $\mathbf{R}$ asso-
ciated with a network of geophones. It contains the set of impulse responses between each pair
of geophones extracted from cross-correlations of seismic noise. Based on the available velocity
model, a focused reflection (FR) matrix is built by applying a redatuming process to $\mathbf{R}$ [16, 20].
This matrix contains the impulse response between virtual sources and receivers synthesized inside
the medium. In the following, matrix “input” and “output” will refer to virtual sources (downgoing
waves) and receivers (upgoing waves), respectively. This FR matrix is powerful as it first provides
an image of the subsurface reflectivity by considering its diagonal elements (i.e when virtual source
and receiver coincide); this is the *so-called* confocal image. Moreover, its off-diagonal elements
allow a local quantification of aberrations in the vicinity of each virtual source. Those aberrations
correspond to the imperfections of the image induced by the mismatch between the wave velocity
model and the actual seismic velocity distribution in the subsurface. In contrast with previous
works [21, 22], a multi-layered wave velocity model is here considered rather than just an homo-
geneous model. This more sophisticated description of seismic wave propagation enables a better
time-to-depth conversion in the confocal image and a better focusing process. Nevertheless, the
FR matrix still highlights residual aberrations that result from the mismatch between the velocity
model and the actual velocity distribution. The fluctuations of wave velocities actually induce
phase distortions on the focused wave-fronts that result in a blurry image of the NAFZ subsurface.

To overcome these detrimental effects, the FR matrix can be first projected in a plane wave
basis. By exploiting the angular input-output correlations of the reflection matrix, phase distortions
of the incident and reflected wave-fronts can be identified and compensated. This is the principle
of the CLASS algorithm (acronym for closed-loop accumulation of single scattering), originally
developed in optical microscopy [18, 23, 24]. Applied for the first time to seismology in the present
study, CLASS successfully compensates for spatially-invariant aberrations and will be shown to
clearly improve the confocal image of the NAFZ subsurface.

Nevertheless, high-order aberrations subsist and are addressed through the distortion matrix
concept in a second step. Originally introduced in ultrasound imaging [17] and optical mi-
croscopy [19], this operator contains the phase distortions of the incident and reflected wave-fronts
with respect to the propagation model. It was recently exploited in passive seismology in order
to image the San Jacinto Fault zone scattering structure that exhibits a sparse distribution of scat-
terers [21]. Here, we apply it in a new scattering regime since the NAFZ subsurface exhibits, in
the frequency range under study, a continuous reflectivity distribution made of specular reflectors
and randomly distributed heterogeneities. In this regime, a local time reversal analysis of the dis-
tortion matrix can be performed in order to retrieve the transmission matrix between the Earth's
surface and any point of the subsurface [25]. This transmission matrix is a key tool since its phase
conjugate provides the optimized focusing laws that need to be applied to the reflection matrix in
order to retrieve a diffraction-limited image of the subsurface. While most conventional reflection
imaging techniques are strongly sensitive to the available velocity model, the reflection matrix
approach is robust with respect to its limitations. An approximate velocity distribution is actually
sufficient since a time-reversal analysis of seismic data enables a local and adapted auto-focus on
each part of the subsurface image.

To image the crustal structure of NAFZ, we use data from the Dense Array for Northern Anatolia [26] that was deployed over the western segment of the fault, in the latest rupture region during the 1999 Izmit (M = 7.6) and Düzce (M = 7.2) earthquakes [27, 28]. The dense array was installed temporarily between May 2012 and October 2013. It consists of 73 3-component broadband seismometers, 66 stations arranged along 11 east-west lines and 6 North-South lines forming a rectangular grid and covering an area of 35 km by 70 km with a nominal inter-station spacing of ~7km (Fig. 1a). Seven additional stations were deployed east of the rectangular array in a semi-circle shape. In this region, the fault splits into two major strands: the northern (NNAFZ) and southern (SNAFZ) strands (Fig. 1b). The northern strand, where most of the continuous deformation occurs according to geodetic studies [2, 29], has been subject to a series of major earthquakes in the last century, among them the 1999 Izmit Earthquake. On the contrary, the latest rupture of the southern branch dates back to the fifteenth century [30]. The fault delineates three tectonic blocks (Fig. 1b): (i) the Istanbul Zone (IZ) situated North of the northern branch, (ii) the Sakarya zone (SZ) situated to the South of the southern branch and (iii) the Armutlu-Almacik crustal block (AA) located in the center, between the two fault strands [31–33]. Differences in crustal composition and properties between these blocks have been reported. Strong velocity contrasts were found across the fault strands by several tomographic studies [34–37] and full waveform inversion studies [38, 39]. Low velocity zones are found below the surface traces of the SNAFZ and NNAFZ [36, 37]. The crust of Istanbul and Armutlu Blocks is characterized by high velocities
while SZ shows relatively low velocities [35–37, 40].

The present study reveals the 3D scattering structure of the medium below this major fault. It does not only image planar interfaces, but provides a direct insight on the heterogeneities that mainly sit in the vicinity of the strands. The observed results complement previous studies conducted in the region. A step in the Moho is detected below the Northern branch, and several sub-Moho structures are observed in the North confirming that the northern branch penetrates in the upper mantle. The southern strand does not have a strong signature in the scattering profiles.

## II. Passive Seismic Matrix Imaging

### A. Reflection matrix in the geophones basis

To apply matrix imaging, we used the ambient seismic noise recorded at DANA (see Fig. 1) to compute the cross-correlation functions of horizontal EE component over the 18 months of recording period. The choice of the EE component is made because it displays a better signal-to-noise than the NN component. With this choice of body wave component, the waves being dealt with are mostly shear waves that have been reflected in depth. First, the data were down-sampled at 25 Hz and corrected from instrument response. Then, the data were split into one-hour windows. Each window is band-pass filtered between 0.1 and 0.5 Hz after applying a spectral whitening between 0.01 and 1 Hz [41]. The cross-correlation between each pair of stations is computed over one-hour windows and finally stacked to obtain the mean cross-correlation function with time lags ranging from -35 to +35 s. The causal and the anti-causal parts of the cross-correlations are then summed in order to improve the convergence towards the Green's functions between seismic stations. Although, considering seismic noise in a higher frequency range would allow, in principle, to improve the resolution of the images, matrix imaging requires the Nyquist criterion to be fulfilled: The inter-station distance (7 km) shall be of the order of a half-wavelength. Considering a S-wave velocity $c_0 = 1700$ m/s near the surface, this criterion led us to choose the 0.1 – 0.5 Hz frequency range ($\lambda = 5.7$ km at the central frequency, with $\lambda$ the wavelength at the Earth surface). The ambient noise energy in the frequency band considered in this study comes from the secondary microseisms (5 – 10 s period band) produced in the ocean [42–44] and constitutes one of the most energetic parts of the seismic noise.

The symmetric cross-correlations can be stacked in a time-dependent response matrix $R_{ss}(t)$.
<figure><img src="image_2.png" /></figure>

FIG. 1. Study region and location of DANA array [26]. (a) Map of the study region and location of DANA array geophones (black triangle. The surface traces of the NAFZ are represented by the red lines [45]. The blue line indicate the location of the cross sections represented in Figs. 4a and 8a. (b) Geological map of the region (adjusted from [40] and [46]). The major geological blocks are represented: Istanbul zone (IZ) in the North, Armutlu-Almacik (AA) block in the center and Sakarya zone (SZ) in the South. The Adapazari and Pamukova basin location are indicated by AB and PB, respectively. (c) Ambient noise EE cross-correlation filtered between 0.1 and 0.5 Hz. The correlations between pairs of stations located South of the SNAF and having at least an angle of 45° with the East-West direction are plotted. These cross-correlograms are stacked over different distances. The predominant surface wave contribution and shear wave echoes induced by planar reflectors in the subsurface are highlighted by red and green lines, respectively.

One element $R(s_i, s_j, t)$ of this matrix corresponds to the impulse response between geophones located at positions $s_i$ and $s_j$. In other words, $R(s_i, s_j, t)$ contains the seismic wave-field recorded at receiver $s_i$ if a pulse was emitted by the virtual source $s_j$ at time $t = 0$. In the following, we will thus refer to columns and lines of the response matrix as input and output wave-fields, respectively.

Figure 1c shows the impulse responses between pairs of stations located South of the SNAF and having at least an angle of 45° with the East-West direction because of the EE polarization
considered in this study. These responses are stacked for different inter-station distances. Figure 1c is dominated by surface-wave energy travelling at about 3000 m/s. Given the horizontal polarization of seismic waves considered here, these surface waves correspond to Love waves. Interestingly, some nearly vertical weaker events are also observed at larger times. Due to their horizontal polarization, they correspond to shear waves reflected by the subsurface heterogeneities in depth.

Another way to highlight the surface and bulk wave components is to investigate the response matrix $\mathbf{R}_{ss}$ in the frequency domain. To this aim, a temporal Fourier transform is applied to $\mathbf{R}_{ss}(t)$. For each frequency $f$ in the bandwidth of interest (0.1-0.5 Hz), a monochromatic matrix $\overline{\mathbf{R}}_{ss}(f)$ is obtained. The different wave components can then be discriminated by a plane wave decomposition of the output wave-fields, such that

$$ \overline{\mathbf{R}}_{sk}(f) = \overline{\mathbf{R}}_{ss}(f) \times \mathbf{P}_0, \qquad (1) $$

where the symbol $\times$ stands for the standard matrix product. $\mathbf{P}_0 = [P_0(\mathbf{s}, \mathbf{k})]$ is the Fourier transform operator that connects each geophone's position $\mathbf{s}$ to the transverse wave vector $\mathbf{k} = (k_x, k_y)$ of each angular component of the wave-field:

$$ P_0(\mathbf{s}, \mathbf{k}) = \exp(i\mathbf{k} \cdot \mathbf{s}), \qquad (2) $$

where the symbol $\cdot$ denotes the scalar product. Figure 2a shows the result of this plane wave decomposition by displaying the mean angular distribution of the output wave-field at $f_0 = 0.2$ Hz. More precisely, this distribution is displayed as a function of the ratio between spatial frequencies $k_x/(2\pi)$ and $k_y/(2\pi)$ and frequency $f$. In that representation, surface waves emerge along a circle whose radius should correspond to the slowness $c_L^{-1}$ of Love waves, with $c_L \sim 3000 \text{ m.s}^{-1}$ [40]. On the other hand, reflected bulk waves are distributed over a disk of radius $c_0^{-1}$.

Figure 2a clearly reveals: (i) The contribution of Love waves with a dominant intensity lying along the y-direction; (ii) A bulk wave component arising at small spatial frequencies and corresponding to the nearly vertical echoes already highlighted by Fig. 1c. The latter contribution is a priori induced by extended reflectors in depth. On the contrary, diffuse scattering can generate reflected waves over a larger distribution of angles. It can therefore account for the incoherent background observed in Fig. 2a. This background could result from the averaging of the random
<figure><img src="image_3.png" /></figure>

FIG. 2. Apparent slowness of seismic echoes contained in the response matrix. (a,b) Plane wave decomposition of the Earth's response matrix $R_{ss}$ at output ($f_0 = 0.2$ Hz, Eq. 1) averaged over the set of green geophones displayed in panels (e,f), respectively. (c,d) Plane wave decomposition of the filtered response matrix $R'_{ss}$ at output ($f_0 = 0.2$ Hz, z=25 km, Eq. 9) averaged over the set of green geophones in panels (e,f), respectively. The dashed circles correspond to apparent velocities of 5000 m/s (green), 3000 m/s (red) and 2000 m/s (blue). Main echoes associated with Love waves (red ellipse), vertical shear waves (green ellipse) and off-axis shear waves (blue circle) are also highlighted in panels (b,d).

speckle pattern exhibited by the angular distribution of the reflected wave-field when only a few sources are considered (Fig. 2b). Nevertheless, it is difficult at this stage to be more affirmative since an imperfect convergence of noise correlations could also lead to such a random wave-field.

A redatuming process is thus required to enhance the weight of scattered shear waves from seismic noise correlations and image the reflectivity of the deep structures around the NAFZ.

## B. Redatuming process

An image of the medium reflectivity can be obtained by applying a double focusing operation to $R_{ss}$ [20, 21]. It consists in back-propagating the response measured at the surface into wave-fields below the surface as if there were sources and receivers inside the medium. This is similar to the
"wave-field extrapolation" concept that forms the basis of the migration process [47]. It requires performing beamforming operations both at emission and reception. On the one hand, focusing in emission consists in applying appropriate time delays to the emitted sources so that waves constructively interfere and focus on one point inside the medium. Physically, this operation amounts to synthesizing a virtual source inside the medium. On the other hand, focusing in reception is carried out by applying proper time delays to the received signals so that they can constructively interfere. As in emission, this focusing operation can be seen as the synthesis of a virtual receiver inside the medium. This operation is known as "redatuming" in seismology [48] and consists of virtually moving sources and receivers from the surface to the medium below (Fig. 3a). Generally, an image of the sub-surface is built by considering the response of virtual source and receiver placed at the same location (Fig. 3f). On the other hand, the principle of matrix imaging consists in decoupling both locations [16].

In the following, the reflection matrix will be expressed in three different bases: (i) the geophones basis where the matrix $R_{ss}$ represents the cross-correlations between all pairs of stations located at $s(x_s, y_s, 0)$, (ii) the focused basis corresponding to the location $r(x, y, z)$ of virtual sources and receivers synthesized by the focusing operations and in which the image of the medium is built, and (iii) the spatial Fourier basis $k = (k_x, k_y)$ that will be first used for wave-field extrapolation and then for aberration correction.

### C. Propagator from the geophones to the focused basis

The focusing operations described in section II B provide the FR matrix that plays a pivotal role in matrix imaging. We now show how this matrix can be obtained through simple matrix operations.

Mathematically, the response between virtual sources and receivers is obtained from the response matrix at the surface through the Green's functions, that describe the propagation between each geophone and each point inside the medium using a wave velocity model. Switching between bases can be easily achieved by simple matrix products in the frequency domain. We first define the plane-wave propagator $T_0(z, f)$, that enables a direct projection of the response matrix from the geophones' basis to the focused basis $\rho = (x, y)$ at each depth $z$. Each monochromatic response matrix $\bar{R}_{ss}(f)$ can be projected in the focused basis both at input and output by applying appropriate phase shifts associated with downgoing waves at input and upgoing waves at output
to provide the FR matrix $\overline{\mathbf{R}}_{\rho\rho}$ (Fig. 3a). Under a matrix formalism, this operation can be written as follows:

$$ \overline{\mathbf{R}}_{\rho\rho}(z, f) = \mathbf{T}_{0}^{\dagger}(z, f) \times \overline{\mathbf{R}}_{ss}(f) \times \mathbf{T}_{0}^{*}(z, f). \quad (3) $$

where the symbols * and † stands for phase conjugate and transpose conjugate, respectively.

A model for the wave velocity distribution inside the medium is required. In this study, and since the horizontal EE cross-correlation functions are considered, only an estimate of the S-wave velocity is required. Unlike [20] and [21] that considered a homogeneous P-wave velocity model, here a layered 1-D S-wave velocity model is used for the focusing process. A combination of two models derived by [49], for the first 5 km, and by [50], for deeper layers is displayed in Table I. Compared to an homogeneous model, such a layered model will allow better time-depth conversion and will limit the aberration level in the subsurface image. Note, however, that our wave propagation model will not account for multiple reflections that could, in principle, take place between the interfaces of the different layers.

<table><caption>TABLE I. 1-D S-wave velocity model. $c_i$ versus depth following [50] and [49].</caption><thead><tr><th>Layer #i</th><th>Depth (km)</th><th>c<sub>i</sub> (m s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>0</td><td>0 - 1</td><td>1700</td></tr><tr><td>1</td><td>1 - 3.5</td><td>2500</td></tr><tr><td>2</td><td>3.5 - 14</td><td>3200</td></tr><tr><td>3</td><td>14 - 26</td><td>3500</td></tr><tr><td>4</td><td>26 - 40</td><td>3600</td></tr><tr><td>5</td><td>40 - 60</td><td>4300</td></tr></tbody></table>

In a layered medium, the forward and backward extrapolation of the reflection matrix can be performed through the decomposition of the wave-field into plane waves [47]. Indeed, plane waves can be easily extrapolated by applying a simple phase shift. To that aim, we define the spatial transfer function, $F_i(\mathbf{k}, f)$, that models the ballistic propagation of shear waves through the $i^{th}$ layer:

$$ F_i(\mathbf{k}, f) = \begin{cases} \exp(-iq_i \Delta z_i) & \text{for } \sqrt{k_x^2 + k_y^2} < 2\pi f / c_i \\ 0 & \text{otherwise} \end{cases}, \quad (4) $$

with

$$ q_i = \sqrt{\left(\frac{2\pi f}{c_i}\right)^2 - k_x^2 - k_y^2}, $$

the vertical component of the wave vector $\mathbf{p}_i = (k_x, k_y, q_i)$ in the $i^{th}$ layer, $c_i$ the wave velocity in
the $i^{th}$ layer of our model and $\Delta z_i$ its thickness. To propagate the plane waves from the surface to depth $z$, we define the wave-field extrapolator, $F(z, f) = [F(k, z, f)]$, as the product of the spatial transfer functions of the $N$ layers above the considered depth as follows:

$$F(k, z, f) = \begin{cases} \exp(-iq_N(z - z_N)) \prod_{i=1}^{N-1} F_i(k, f) & \text{for } \sqrt{k_x^2 + k_y^2} < 2\pi f / c_N \\ 0 & \text{otherwise.} \end{cases}, \quad (5)$$

where $z_i$ is the depth at which starts the $i^{th}$ layer. The phase propagator, $T(z, f) = [T(s, \rho, z, f)]$, can finally be expressed as follows:

$$T_0(z, f) = [P_0 \circ F(z, f)] \times P'_0. \quad (6)$$

where the symbol $\circ$ refers to the Hadamard product (i.e. element wise matrix multiplication). In the equation 6, the term-by-term product arises because wave propagation in the plane wave basis is modeled by the scalar product of each plane wave component $P_0(k, s)$ with the overall spatial transfer function $F(k, z)$. The matrix multiplication stands for the inverse Fourier transform that enables us to project the propagated wave-field from the plane wave basis to the focused basis. $P'_0 = [P'_0(\rho, k)]$ is actually the Fourier transform operator linking the focused and plane wave bases. It connects the transverse wave vector $k = (k_x, k_y)$ of each plane wave to the transverse coordinates $\rho = (x, y)$ of each focusing point:

$$P'_0(\rho, k) = \exp(i\mathbf{k} \cdot \rho). \quad (7)$$

To avoid aliasing during the change of basis between the plane wave and focused bases, a Shannon criterion should be respected. The transverse wave components, $k_x$ and $k_y$, are required to fulfill the following condition:

$$\sqrt{k_x^2 + k_y^2} < 2\pi f / c_0. \quad (8)$$

The resolution $\delta k$ of the Fourier plane is conditioned by the size of the array $\mathcal{D} = 50$ km such that $\delta k = 2\pi/\mathcal{D}$. By properties of the Fourier transform, the transverse resolution $\delta \rho_0$ in the focal plane, that corresponds to the distance between the focusing points $\mathbf{r}$, is chosen to be $\sim \lambda/2 \sim 2$ km to circumvent spatial aliasing in the focused basis.
As shown by Eq. 5, wave components of spatial frequencies larger than the wavenumber $2\pi f/c_i$ cannot be transmitted through the $i^{th}$ layer. As a consequence, redatuming acts as a low-pass filter in the spatial frequency domain. To illustrate this phenomenon, one can back-project the focused reflection matrix in the geophone basis,

$$ \overline{\mathbf{R}}'_{ss}(f) = \mathbf{T}_0(z, f) \times \overline{\mathbf{R}}_{\rho\rho}(z, f) \times \mathbf{T}_0^\top(z, f), \quad (9) $$

where the superscript $\top$ stands for matrix transpose. From $\overline{\mathbf{R}}'_{ss}(f)$, one can investigate the angular decomposition of the output wave-fields as previously done for the original response matrix $\overline{\mathbf{R}}_{ss}$ in Figs. 2(a) and (b) (Eq. 1). The result is displayed in Figs. 2(c) and (d). The comparison with their original counterparts highlights the low pass-filter operated by redatuming: The surface wave component is discarded and only the shear waves associated with spatial frequencies $k < 2\pi f/c_N$ are kept. Figure 2(d) displays: (i) a low-spatial frequency component associated with specular reflectors; (ii) several off-axis bright spots associated with peculiar single scattering events at depth $z$ and; (iii) a diffuse background that is difficult to interpret at this stage since it can be due to random single scattering, multiple scattering or noise resulting from the imperfect convergence of cross-correlations towards the Green's function. To enhance the single scattering contribution with respect to the other undesirable components for imaging, the idea is now to perform a time gating operation to enhance the single scattering contribution.

### D. Broadband focused reflection matrix

Equation 3 simulates focused beamforming for both downgoing (input) and upgoing (output) shear waves at each frequency. Each spectral component of the wave-field can then be recombined to provide a broadband focused reflection matrix $\overline{\mathbf{R}}_{\rho\rho}(z)$:

$$ \overline{\mathbf{R}}_{\rho\rho}(z) = \int_{f_1}^{f_2} df \overline{\mathbf{R}}_{\rho\rho}(z, f), \quad (10) $$

with $f_1 = 0.1$ Hz and $f_2 = 0.5$ Hz. This operation amounts to performing an inverse Fourier transform at time $t = 0$ over the frequency band [0.1 0.5] Hz. Note that this inverse Fourier transform is only performed over positive frequencies in order to have access to both the amplitude and phase of the wave-field. The time origin here corresponds to the ballistic time in the focused
<figure><img src="image_4.png" /><figcaption>FIG. 3. Focused reflection matrix. (a) The response matrix $K$ is projected onto a focused basis at each depth $z$ (Eq. 10), thereby synthesizing a set of virtual sources ($\rho_{in}$) and receivers ($\rho_{out}$) scanning laterally the field-of-view. In presence of fluctuations in the seismic velocity, focused waves are distorted while travelling from the surface to the plane, thereby enlarging and distorting the virtual geophones. (b) This effect gives rise to a off-diagonal spreading of backscattered energy in the focused reflection matrix $R_{\rho\rho}(z)$ shown here at depth $z = 25$ km. (c) The corresponding intensity profile, averaged over whole the field-of-view, provides the so-called RPSF $I$ (Eq. S2). The white circle represents the diffraction-limited transverse resolution ($\delta\rho_0 \sim 6$ km) at the considered depth. (d) The ideal RPSF that would be obtained in absence of aberrations is shown for comparison (e) Confocal image $\mathcal{I}$ (Eq. 13) built from the diagonal of $R_{\rho\rho}$. The white box represents the dimensions of the rectangular array of geophones and the red lines represent the NAFZ fault traces at the surface. (f) The confocal image corresponds to a simultaneous focusing process at input and output ($r_{in} = r_{out}$). In panels (b)-(e), the color scale refers to the scattering intensity. It is normalized by the maximum value of the scattering energy at the considered depth.</figcaption></figure>

basis. Equation 10 thus corresponds to a time gating operation that tends to select singly-scattered waves associated with a scattering event in the focal plane. Each coefficient $R(\rho_{out}, \rho_{in}, z)$ of $R_{\rho\rho}(z)$ contains the complex wave-field that would be recorded by a virtual geophone located at $r_{out} = (\rho_{out}, z)$ if a virtual source at $r_{in} = (\rho_{in}, z)$ emits a pulse of length $\delta t = \Delta f^{-1}$ at the central frequency $f_0$, with $\Delta f = f_2 - f_1 = 0.4$ Hz.

The FR matrix can be expressed theoretically as follows [16, 21, 51],

$$R_{\rho\rho}(z) = H^{\top}(z) \times \Gamma(z) \times H(z), \qquad (11)$$
which yields, in terms of matrix coefficients,

$$R(\rho_{\text{out}}, \rho_{\text{in}}, z) = \int d\rho H(\rho, \rho_{\text{out}}, z)\gamma(\rho, z)H(\rho, \rho_{\text{in}}, z). \qquad (12)$$

The matrix Γ describes the scattering process in the focused basis. In the single scattering regime, this matrix is diagonal and its coefficients correspond to the subsurface reflectivity γ(ρ, z) at depth z. H(z) is the focusing matrix whose coefficients H(ρ, ρ<sub>in/out</sub>, z) correspond to the point spread functions (PSFs) of the redatuming process. These PSFs represent the spatial amplitude distribution of the focal spots for each focusing point r<sub>in/out</sub> = (ρ<sub>in/out</sub>, z). They thus account for the lateral extent of each vitual source/detector at r<sub>in/out</sub>.

An example of the broadband FR matrix R<sub>ρρ</sub> is shown at depth z = 25 km in Fig. 3b. R<sub>ρρ</sub> is a four-dimension matrix concatenated in 2D as a set of blocks [20]. If the wave velocity model was correct, the PSFs of the redatuming process would be close to be point-like [H(ρ, ρ<sub>in/out</sub>, z) ≃ δ(ρ - ρ<sub>in/out</sub>), with δ the Dirac distribution] and the focused reflection matrix R<sub>ρρ</sub> almost diagonal [R(ρ<sub>out</sub>, ρ<sub>in</sub>, z) ≃ γ(ρ<sub>in</sub>)δ(ρ<sub>in</sub> - ρ<sub>out</sub>), see Eq. S1]. Here, the backscattered energy is far from being concentrated along the diagonal of R<sub>ρρ</sub>, which is a manifestation of the gap between our layered wave velocity model (Table I) and the real shear wave velocity distribution in the subsurface.

### E. Confocal image

Nevertheless, one can try to build an image of the medium reflectivity at each effective depth z by considering the diagonal elements of the FR matrix, i.e where the virtual sources and receivers coincide (ρ<sub>in</sub> = ρ<sub>out</sub> = ρ<sub>c</sub>, see Fig. 3c). It yields the so-called *confocal image*:

$$\mathcal{I}(\rho_c, z) = R(\rho_c, \rho_c, z). \qquad (13)$$

Fig. 3e shows the resulting 2D image I at z = 25 km retrieved from the diagonal of the FR matrix in Fig. 3b. By injecting Eq. S1 into the last equation, I can be expressed as the transverse convolution between the medium reflectivity and the confocal PSF H<sup>2</sup>, such that:

$$\mathcal{I}(\rho_c, z) = \int d\rho \gamma(\rho, z) H^2(\rho, \rho_c, z). \qquad (14)$$

Such an image is thus a reliable estimator of the reflectivity at depth z only if the wave veloc-
ity model is close to reality. In this ideal case, the spatial extent of the PSF is only limited by diffraction and the transverse resolution is given by [52]:

$$ \delta\rho_0 = \lambda/(2 \sin \theta) \qquad (15) $$

where $\theta = \arctan(D/2z)$ is determined by the size of the array $D = 50$ km, and corresponds to the maximum angle under which a focusing point sees the geophones' array.

By stacking the confocal image computed at each depth $z$, a 3D image of the reflectivity can be obtained. The cross-section at Lon 30.37° is displayed in Fig. 4a. 2D confocal images at $z = 15, 30$ and $40$ km are also shown in Fig. 4b. Unlike the transverse resolution, the axial resolution $\delta z$ is limited by the frequency bandwidth: $\delta z \sim c/\Delta f \sim 8.7$ km, with $c$ the shear wave velocity at the considered depth. The sections of the 3D image displayed in Figs. 3e, 4a and 4b show a greater reflectivity in the central part of the field-of-view, i.e right below the geophones' array, but no direct correlation can be found between the image and the location of the fault strands. In fact, as we will see now, lateral wave speed heterogeneities strongly degrade the transverse resolution of the redatuming process and induce strong aberrations in the confocal image. This image is thus not a reliable estimator of the medium reflectivity at this stage and cannot be interpreted.

## F. Quantification of aberrations

The FR matrix can provide more than a confocal image since its off-diagonal elements can lead to a quantification of aberrations. To that aim, a relevant observable is the distribution of the backscattered intensity around a common midpoint point $(\rho_m, z)$ as a function of the relative position, $\Delta\rho = \rho_{out} - \rho_{in}$, between the input and output focusing points [16, 21]:

$$ I(\Delta\rho, \rho_m, z) = |R(\rho_m - \Delta\rho/2, \rho_m + \Delta\rho/2, z)|^2. \qquad (16) $$

In the following, we will refer to this quantity as the reflection point spread function (RPSF).

To express this quantity theoretically, we first make a local isoplanatic approximation in the vicinity of each point $(\rho_m, z)$ [51]. Isoplanicity means here that waves which focus in this region are assumed to have travelled through approximately the same areas of the medium, thereby undergoing identical phase distortions. The PSF can then be considered to be spatially invariant within this local region. Mathematically, this means that, in the vicinity of each point $(\rho_m, z)$, the
<figure><img src="image_5.png" /><figcaption>FIG. 4. Original confocal image of NAFZ. (a) Vertical North-South cross-section at 30.37°E. The North-South profile is oriented perpendicular to the fault traces. The location of the profile is shown in Fig. 1a. The locations of the southern (SNAF) and northern (NNAF), and the major crustal blocks (SZ: Sakarya zone, AA: Armutlu-Almacik and IZ: Istanbul zone) are labeled. The color scale refers to the scattering intensity. It is normalized by the maximum value of the scattering energy inside the volume. (b) Depth slices retrieved from the 3D scattering volume at z = 15, 30 and 40 km with (c) their corresponding RPSFs.</figcaption></figure>

spatial distribution of the PSF, $H(\rho, \rho_{in/out}, z)$, only depends on the relative distance between the point $\rho$ and the focusing point $\rho_{in/out}$. This leads us to define a local spatially-invariant PSF $H_L$ around each common mid-point $(\rho_m, z)$ such that:

$$H(\rho, \rho_{in/out}, z) = H_L(\rho - \rho_{in/out}, \rho_m, z). \qquad (17)$$

Under this local isoplanatic assumption, the RPSF can be derived analytically in different scattering regimes. On the one hand, for large reflectors such as horizontal interfaces between geological units, the medium reflectivity can be assumed as locally constant and the RPSF is given by (see Supplementary Section S1):

$$I(\Delta\rho, \rho_m, z) = |\gamma(\rho_m, z)|^2 \left| H_L \circledast H_L \right|^2 (\Delta\rho, \rho_m, z). \qquad (18)$$
where the symbol ⊗ denotes convolution. On the other hand, for diffuse scattering, the medium reflectivity can be considered, in first approximation, as randomly distributed. Under that assumption, the mean RPSF is then proportional to the convolution between the incoherent output and input local PSFs, independently from the medium's reflectivity [16] (see Supplementary Section S1):

$$ \langle I(\Delta\rho, \rho_m, z) \rangle \propto \left[ |H_L|^2 \stackrel{\Delta\rho}{\otimes} |H_L|^2 \right] (\Delta\rho, \rho_m, z) \quad (19) $$

where the symbol $\langle \cdots \rangle$ stands for an ensemble average. Whatever the scattering regime, the spatial extension of the RPSF is thus roughly equal to the lateral dimension of the PSF. If we assume a Gaussian PSF, this equality is strict. The RPSF is thus a direct indicator of the focusing quality and its spatial extent directly provides an estimation of the local transverse resolution of the confocal image.

Fig. 3c displays the RPSF averaged over the whole field-of-view at depth $z = 25$ km. For sake of comparison, Fig. 3d shows the ideal (i.e diffraction-limited) RPSF that would be obtained in absence of aberrations. The comparison between Figs. 3c and d highlights the impact of aberrations resulting from the mismatch between the wave velocity model of Table I and the real wave speed distribution. Indeed, the full width at half maximum $w$ of the intensity profile is increased by a factor $\sim 6$ compared to its diffraction-limited value (white circle in Fig. 3d, Eq. 15) at depth $z = 25$ km. This effect explains the blurred aspect of the confocal image displayed in Fig. 3e at the same depth. The impact of aberrations is also illustrated in Fig. 4c that displays the depth evolution of the RPSF inside the Earth. As with the diffraction-limited resolution (Eq. 15), the transverse extension of the RPSF also increases with $z$ but it shows a much larger extension.

In the following we will show how matrix imaging can restore an optimal resolution for this image.

### III. Exploiting the input-output angular correlations of the wave-field: The CLASS algorithm

In order to compensate for aberrations, the reflection matrix can be first projected in the plane wave basis:

$$ \mathbf{R}_{kk}(z) = \mathbf{P}'_0 \times \mathbf{R}_{\rho\rho}(z) \times \mathbf{P}'_0^\top, \quad (20) $$

Using Eq. 7, the last equation can be rewritten, in terms of matrix coefficients, as a double spatial
Fourier transform:

$$R(\mathbf{k}_{\text{out}}, \mathbf{k}_{\text{in}}, z) = \sum_{\rho_{\text{out}}} \sum_{\rho_{\text{in}}} e^{-i\mathbf{k}_{\text{out}} \cdot \rho_{\text{out}}} R(\rho_{\text{out}}, \rho_{\text{in}}, z) e^{-i\mathbf{k}_{\text{in}} \cdot \rho_{\text{in}}} \quad (21)$$

Each coefficient of the matrix $R_{kk}(z) = [R(\mathbf{k}_{\text{out}}, \mathbf{k}_{\text{in}}, z)]$ thus contains the medium response between input and output transverse wave vectors $\mathbf{k}_{\text{in}}$ and $\mathbf{k}_{\text{out}}$. By injecting Eq. 11 into Eq. 20, the matrix $R_{kk}$ can be expressed as follows:

$$R_{kk}(z) = T(z) \times \Gamma(z) \times T^{\top}(z), \quad (22)$$

where $T(z) = P'_{0} \times H(z)$ is the transmission matrix describes plane wave propagation between the focused and the plane wave bases. Its coefficients $T(\mathbf{k}, \rho, z)$ correspond to the angular decomposition of the wave-field produced at the Earth surface for a point-like virtual source located at $(\rho, z)$. This matrix is critical for imaging since its inversion can provide a direct access to the subsurface reflectivity, without relying on a precise wave velocity model.

As a first step towards the estimation of $T$, we can go one step further in the isoplanatic approximation (Eq. 17) by assuming a full transverse-invariance of the PSF across the field-of-view. This leads us to define a laterally-invariant PSF $H_I$, such that:

$$H_I(\rho - \rho_{\text{in/out}}, \rho_m, z) \approx H_I(\rho - \rho_{\text{in/out}}, z). \quad (23)$$

This strong assumption means that wave speed heterogeneities are modelled by a phase screen of transmittance $\tilde{H}_I = [\tilde{H}_I(\mathbf{k}, z)]$ in the plane wave basis, such that $T(z) = \tilde{H}_I(z) \circ P'_{0}$, where $\tilde{H}_I(\mathbf{k}, z) = \int d\rho H_I(\rho, z) e^{-i\mathbf{k} \cdot \rho}$ is the Fourier transform of the spatially-invariant PSF $H_I(\rho, z)$. The aberration transmittance $\tilde{H}_I$ grasps the phase distortions undergone by downgoing and upgoing wave-fields during their travel between the Earth surface and the focal plane at effective depth $z$.

Under this full isoplanatic approximation (Eq. 23), a theoretical expression of $R_{kk}$ can be derived in the single scattering regime [17]:

$$R(\mathbf{k}_{\text{out}}, \mathbf{k}_{\text{in}}, z) = \tilde{H}_I(\mathbf{k}_{\text{in}}, z) \tilde{\gamma}(\mathbf{k}_{\text{in}} + \mathbf{k}_{\text{out}}, z) \tilde{H}_I(\mathbf{k}_{\text{out}}, z), \quad (24)$$

where $\tilde{\gamma}(\mathbf{k}, z) = \int d\rho \gamma(\rho, z) \exp(-i\mathbf{k} \cdot \rho)$ is the 2D Fourier transform of the medium's re-
flectivity $\gamma(\rho, z)$ at each depth $z$. In absence of aberrations ($\tilde{H}(\mathbf{k}, z) \equiv 1$), the reflection matrix expressed in the plane wave basis exhibits a deterministic coherence along its antidiagonals ($\mathbf{k}_{\text{in}} + \mathbf{k}_{\text{out}} = \text{constant}$, see Fig. 5c): $R(\mathbf{k}_{\text{in}}, \mathbf{k}_{\text{out}}, z) = \tilde{\gamma}(\mathbf{k}_{\text{in}} + \mathbf{k}_{\text{out}}, z)$ [53, 54]. This peculiar property is a manifestation of a phenomenon called the memory effect in wave physics [55, 56]. When an incident plane wave ($\mathbf{k}_{\text{in}}$) illuminates a scattering medium, it gives rise to a reflected wave-field ($\mathbf{k}_{\text{out}}$) that exhibits a speckle feature (Fig. 5a) due to the random interference between partial waves induced by each scatterer lying at depth $z$. When this incident plane wave is rotated by a certain angle ($\mathbf{k}_{\text{in}} + \Delta\mathbf{k}$), the reflected wave-field is tilted in the opposite direction ($\mathbf{k}_{\text{in}} - \Delta\mathbf{k}$, see Fig. 5b). This correlation between the downgoing and upgoing wave-fields accounts for the deterministic coherence along the antidiagonals of $\mathbf{R}_{\mathbf{k}\mathbf{k}}$ (Fig. 5c). However, in the present case, this property is not checked because of the phase distortions undergone by downgoing and upgoing wave-fields induced by wave speed heterogeneities. Mathematically, this is accounted by the phase screen $\tilde{H}_I(\mathbf{k}, z)$ in Eq. 24 that breaks the correlation between coefficients lying along the same antidiagonal (Fig. 5d).

The principle of the CLASS algorithm [18, 23, 24] consists in restoring this coherence by applying a phase correction, $\exp[-i\phi_C(\mathbf{k}, z)]$, to the input and output of $\mathbf{R}_{\mathbf{k}\mathbf{k}}$:

$$R^{(C)}(\mathbf{k}_{\text{out}}, \mathbf{k}_{\text{in}}, z) = e^{-i\phi_C(\mathbf{k}_{\text{out}}, z)} R(\mathbf{k}_{\text{out}}, \mathbf{k}_{\text{in}}, z) e^{-i\phi_C(\mathbf{k}_{\text{in}}, z)}, \quad (25)$$

where $\phi_C$ is the estimator of the aberration phase law whose phase conjugate maximizes the coherence along the antidiagonals of $\mathbf{R}_{\mathbf{k}\mathbf{k}}$. To compute $\phi_C$, the first step is to perform a coherent sum of $\mathbf{R}_{\mathbf{k}\mathbf{k}}$'s coefficients along each of its antidiagonals [see Fig. 5e]:

$$C(\mathbf{k}_{+}, z) = \sum_{\mathbf{k}_{\text{out}}} R(\mathbf{k}_{\text{out}}, \mathbf{k}_{+} - \mathbf{k}_{\text{out}}, z) \quad (26)$$

with $\mathbf{k}_{+} = \mathbf{k}_{\text{in}} + \mathbf{k}_{\text{out}}$. As shown in Supplementary Section S2, $C(\mathbf{k}_{+}, z)$ is a rough estimator for the spatial frequency spectrum $\tilde{\gamma}(\mathbf{k}_{+}, z)$ of the medium reflectivity (see Supplementary Section S2). The second step consists in performing the Hadamard product (element-wise product) between the phase conjugate of the vector $\mathbf{C}$ and the matrix $\mathbf{R}_{\mathbf{k}\mathbf{k}}$ (Fig. 5f):

$$R'(\mathbf{k}_{\text{out}}, \mathbf{k}_{\text{in}}, z) = R(\mathbf{k}_{\text{out}}, \mathbf{k}_{\text{in}}, z) C^*(\mathbf{k}_{\text{out}} + \mathbf{k}_{\text{in}}, z). \quad (27)$$
The last operation amounts to compensate for the phase fluctuations of the reflectivity spatial frequency spectrum $\tilde{\gamma}$ in $R_{kk}$. The isoplanatic phase distortion $\tilde{H}_I$ is finally estimated by summing the columns of the compensated matrix $R'_{kk}$ (see Fig. 5f)

$$ \phi_C(\mathbf{k}_{\text{out}}, z) = \arg \left[ \sum_{\mathbf{k}_{\text{in}}} R'(\mathbf{k}_{\text{out}}, \mathbf{k}_{\text{in}}, z) \right], \qquad (28) $$

The phase conjugate of the resulting wave-front, $\exp[-i\phi_C(\mathbf{k}, z)]$, tends to realign in phase the coefficients lying on the same antidiagonal of $R_{kk}$ (Eq. 25). As shown in Supplementary Section S2, this phase realignment in the plane wave basis is equivalent to a maximization of the confocal intensity in the focused basis.

The corresponding FR matrix $R_{\rho\rho}^{(C)}$ can be deduced as follows:

$$ R_{\rho\rho}^{(C)}(z) = P_0^{\dagger}(z) \times R_{kk}^{(C)}(z) \times P_0^{*}(z) \qquad (29) $$

A corrected confocal image is extracted from the diagonal of $R_{\rho\rho}^{(C)}(z)$ and displayed in Fig. 6a at depth $z = 25$ km. It should be compared with the original image shown in Fig. 3e. While the latter one displays a random-like feature, the corrected image reveals a greater reflectivity in the North that can be correlated with the expected damage around the Northern branch of the fault. The comparison between these two images illustrates the benefit of the correction process. The gain in resolution can be assessed by averaging the RPSF (Eq. S2) over the whole field-of-view (see Fig. 6c). Compared to the original RPSF displayed in Fig. 3c, we can notice that a large component of the off-diagonal energy has been brought back to the confocal lobe (white circle). The resolution $w$ is reduced from 40 km to 8 km but it is still larger than the diffraction-limited resolution ($\delta\rho_0 \sim 6$ km at the considered depth). A diffuse component subsists and can be explained by the spatially-varying residual aberrations, $\delta\tilde{H}_L(\mathbf{k}, \rho_m, z)$, that have not been compensated by the CLASS algorithm, such that $\delta\tilde{H}_L(\mathbf{k}, \rho_c, z) = \tilde{H}_L(\mathbf{k}, \rho_c, z)e^{-i\phi_C(\mathbf{k},z)}$.

A local compensation of higher order aberrations is thus required. This issue is handled in the following section by investigating the reflection matrix and its distorted component between the focused and plane wave bases.
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 5. Principles of the memory effect and CLASS algorithm. (a) When an incident plane wave of wave vector $\mathbf{p}_{\text{in}}$ insonifies a scattering medium, the reflected wave-field exhibits a random speckle pattern. (b) When this incident wave-field is rotated by an angle $\theta$, the reflected wave-field is shifted by the opposite angle $-\theta$: This is the so-called memory effect. (c) This phenomenon results in a deterministic coherence along the antidiagonals ($\mathbf{k}_{\text{in}} + \mathbf{k}_{\text{out}} = \text{constant}$) of the reflection matrix $\mathbf{R}_{\text{kk}}$ expressed in the plane wave basis. The phase of a matrix $\mathbf{R}_{\text{kk}}$ is displayed for sake of illustration. This matrix has been obtained from an ultrasound experiment performed on a medium of random reflectivity (acoustic phantom) in the conditions described by [16]. (d) Same matrix as in (c) but in presence of aberrations. Each complex coefficient of $\mathbf{R}_{\text{kk}}$ is here represented by blue arrows using a Fresnel diagram. (e) The first step of CLASS (Eq. 26) consists in a coherent sum of coefficients lying along the same antidiagonal to estimate the spatial frequency spectrum $\tilde{\gamma}$ of the medium reflectivity. Each coefficient of the resulting vector $\mathbf{C}$ is represented by a red arrow in the complex plane. (f) The second step of CLASS consists in a compensation of $\mathbf{R}_{\text{kk}}$ by $\mathbf{C}^*$ to compensate for the phase of reflectivity spectrum $\tilde{\gamma}$ (Eq. 27). The coefficients of the resulting matrix $\mathbf{R}'_{\text{kk}}$ are depicted with orange arrows. A sum over lines or columns of $\mathbf{R}'_{\text{kk}}$ provides the isoplanatic aberrated wave-front $\exp(i\phi_c)$ represented by purple arrows (Eq. 28).</figcaption></figure>
<figure><img src="image_8.png" /></figure>

FIG. 6. Aberration correction process at $z = 25$ km. (a) Confocal image obtained after applying the conjugate of (b) the CLASS phase law $\phi_C$ computed at this depth. (c) RPSF obtained after CLASS correction. (d) Confocal image obtained after performing four iteration steps of the distortion matrix process. The red lines represent the NAFZ fault traces at the surface. The yellow dashed lines delineate the regions over which a local aberration phase law $\phi(k, r)$ has been estimated. (e) Corresponding input aberration phase laws $\phi_{in}(k, r)$ obtained at the end of the process. The correlation coefficients between the corresponding aberration transmittances and the central one are displayed below each phase mask. (f) RPSF at the end of the matrix imaging process.

#### IV. Matrix approach for adaptive focusing: The local distortion matrix

The distortion matrix $D$ was already introduced in ultrasound [17, 25], optics [19, 57] and seismology [21]. Several applications proved the efficiency of this matrix in overcoming aberrations and improving the image quality. Recent works in seismology [21] and optics [19] have shown that for certain scattering regimes (specular reflectors or sparse scattering), there was a one-to-one association between the eigenstates of $D$ and the isoplanatic patches present in the field-of-view. Here, this property does not hold because the NAFZ subsurface exhibits a continuous but fluctuating reflectivity (see Supplementary Section S3). In this scattering regime, local distortion matrices should be considered over restricted areas in which the isoplanatic hypothesis is ideally
fulfilled [25, 57].

In this section, the distortion matrix concept is applied to the CLASS FR matrix obtained in the previous section for compensation of spatially-distributed aberrations. The process is outlined by five steps: (i) projection of the CLASS FR matrix at output into the plane wave basis (Fig. 7a), (ii) the realignment of the reflected wave-fronts to form a distortion matrix $D = [D(k_{out}, \rho_{in}, z)]$ (see Fig. 7b), (iii) the truncation of $D$ into local distortion matrices $D'(r_p)$, (iv) the singular value decomposition of $D'(r_p)$ to extract a residual aberration phase law for each point $r_p$ and build an estimator $\hat{T}$ of the transmission matrix $T$ (Fig. 7c); (v) the phase conjugation of $\hat{T}$ to correct for output residual aberrations (Fig. 7d). All of these steps are then repeated by exchanging output and input bases.

## A. The distortion matrix

The output of the CLASS algorithm is a focused reflection matrix $R_{\rho\rho}^{(C)}$ that still exhibits laterally-varying aberrations. To assess these residual aberrations, the first step is to choose a basis in which the distortion of the CLASS wave-front is the most spatially-invariant. In a horizontally multi-layered medium such as NAFZ, the plane-wave basis is the most adequate since plane waves are the propagation invariants in this geometry. A plane-wave projection is performed at the output of the CLASS FR matrix $R_{\rho\rho}^{(C)}$:

$$R_{k\rho}^{(C)}(z) = P'_{0} \times R_{\rho\rho}^{(C)}(z), \qquad (30)$$

$R_{k\rho}^{(C)}(z) = [R^{(C)}(k_{out}, \rho_{in}, z)]$ connects each input focusing point $r_{in} = (\rho_{in}, z)$ to the CLASS wave-field in the plane wave basis (Fig. 7a).

The CLASS wave-field can be understood as a sum of two components: (i) a geometrical component described by the reference matrix $P'_{0}$, containing the ideal wave-front generated by a source at $r_{in}$ according to the propagation model described in Table I (dashed black curves in Fig. 7a); (ii) a distorted component due to spatially distributed aberrations that subsists after the CLASS procedure. The latter component refers to the residual phase distortions that should be isolated from the CLASS wave-field in order to be properly compensated. This can be done by subtracting the ideal wave-front that would be obtained in absence of aberrations (i.e the geometrical component) from each CLASS wave-front induced by each input focusing wave at $r_{in}$. Such operation can be
<figure><img src="image_9.png" /><figcaption>FIG. 7. Local aberration correction. (a) One-side plane wave decomposition of each CLASS matrix yields the reflected wave-front associated with each focusing point $r_{in}$. (b) By removing the geometrical curvature of each reflected wavefront (dashed line in a), one can study the phase distortions over each isoplanatic patch identified by their midpoint $r_p$. This operation amounts to realigning the wavefronts in each isoplanatic area as if they were generated by input focal spots virtually shifted to the origin. (c) SVD of each distortion matrix yields an aberration phase law $\phi_D$ for each spatial window by combining coherently each focal spot to synthesize a virtual coherent reflector. The set of aberration phase laws forms the estimator $\hat{T}$ of the transmission matrix. (d) The phase conjugate of $\hat{T}$ provides the focusing laws to compensate for phase distortions for each patch $r_p$.</figcaption></figure>

expressed mathematically via a Hadamard product between $R_{k\rho}^{(C)}(z)$ and $P_0'$. It yields the residual distortion matrix $D(z)$:

$$D(z) = R_{k\rho}^{(C)}(z) \circ P_0', \quad (31)$$

The matrix $D(z)$ connects any input virtual source $r_{in}$ to the residual distortion exhibited by the CLASS wave-field expressed in the plane wave basis (Fig. 7b). By removing the geometrical component of the CLASS wave-field, spatial correlations are highlighted between distorted wave-fields
induced by neighbour virtual sources $\mathbf{r}_{\text{in}}$ [19]. Such correlations are a manifestation of a spatial invariance of residual aberrations over areas generally referred to as isoplanatic patches [25].

## B. Local distortion matrices

Our strategy is to divide the field-of-view into a set of overlapping regions (Fig. 6d). Each region is defined by a central midpoint $\mathbf{r}_p = (\rho_p, z_p)$ and a spatial extension $L$. For each region, the local residual D-matrix is defined as:

$$D'(\mathbf{k}_{\text{out}}, \rho_{\text{in}}, \mathbf{r}_p) = D(\mathbf{k}_{\text{out}}, \rho_{\text{in}}, z_p)W(\rho_{\text{in}} - \rho_p), \quad (32)$$

where $W(\rho)$ is a spatial window function such that $W(\rho) = 1$ for $|x| < L$ and $|y| < L$, and zero elsewhere. Ideally, wave-front distortions should be invariant over each region, meaning that the virtual sources $\mathbf{r}_{\text{in}} = (\rho_{\text{in}}, z)$ associated with each region belong to the same isoplanatic patch. However, in practice, this hypothesis is not fulfilled. The isoplanatic length actually scales as the typical transverse dimension over which the wave velocity fluctuates. On the one hand, the dimension $L$ of the window function should therefore be reduced to cover the smallest isoplanatic region as possible in order to provide a local and sharp measurement of aberrations. On the other hand, it should also be large enough to include a sufficient number of realizations of disorder in order to unscramble the effect of aberrations from the medium's reflectivity [25]. To reach a good estimate of the aberration phase law, the number of input focusing points in each region should be one order of magnitude larger than the number of resolution cells mapping the CLASS focal spot (Fig. 6c) [17]. This is why the initial CLASS step was important to initiate the aberration correction process and reduce the extension of the focal spots before a local and finer compensation of residual aberration by means of the D-matrix concept. The area covered by the CLASS focal spot being $20 \times 14 \text{ km}^2$ (Fig. 6c), the extent of the window is chosen to be $55 \times 55 \text{ km}^2$.
### C. Singular value decomposition

Assuming local isoplanicity in each spatial window $W_L$ (Eq. 17), the coefficients of each distortion matrix $D'(\mathbf{r}_p)$ matrix can be expressed as follows [17]:

$$D'(\mathbf{k}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, \mathbf{r}_p) = \delta \tilde{H}_L(\mathbf{k}_{\text{out}}, \mathbf{r}_p) \underbrace{\int d\boldsymbol{\rho} \gamma(\boldsymbol{\rho} + \boldsymbol{\rho}_{\text{in}}, z) \delta H_L(\boldsymbol{\rho}, \mathbf{r}_p) e^{i\mathbf{k}_{\text{out}} \cdot \boldsymbol{\rho}}}_{\text{virtual source}} \quad (33)$$

This equation can be seen as a product between two terms: the output aberration transmittance and a virtual source term modulated by the medium's fluctuating reflectivity $\gamma(\boldsymbol{\rho}, z)$. The goal is now to unscramble these two terms in order to get a proper estimation of the aberration transmittance $\delta \tilde{H}_L(\mathbf{k}_{\text{out}}, \mathbf{r}_p)$ at each point $\mathbf{r}_p$.

In practice, this can be done through a singular value decomposition (SVD) of each local distortion matrix $D'(\mathbf{r}_p)$:

$$D'(\mathbf{r}_p) = U(\mathbf{r}_p) \times \Sigma(\mathbf{r}_p) \times V(\mathbf{r}_p)^{\dagger} \quad (34)$$

where $\Sigma$ is a diagonal matrix containing the real positive singular values $\sigma_i$ in a decreasing order $\sigma_1 > \sigma_2 > \dots > \sigma_N$. $U(\mathbf{r}_p)$ and $V(\mathbf{r}_p)$ are unitary matrices whose columns, $U_i(\mathbf{r}_p) = [U_i(\mathbf{k}_{\text{out}}, \mathbf{r}_p)]$ and $V_i(\mathbf{r}_p) = [V_i(\mathbf{r}_{\text{in}}, \mathbf{r}_p)]$, correspond to the output and input singular vectors, respectively. The physical meaning of this SVD can be intuitively understood by considering the asymptotic case of a point-like input focusing beam: $\delta H_L(\boldsymbol{\rho}, \mathbf{r}_p) = \delta(\boldsymbol{\rho})$. In this ideal case, Eq. 33 becomes: $D'(\mathbf{k}_{\text{out}}, \boldsymbol{\rho}_{\text{in}}, \mathbf{r}_p) = \delta \tilde{H}_L(\mathbf{k}_{\text{out}}, \mathbf{r}_p)\gamma(\boldsymbol{\rho}_{\text{in}}, z_p)$. Comparison with Eq. 34 shows that, a first approximation, $D'(\mathbf{r}_p)$ is of rank 1. The first output singular vector $U_1(\mathbf{r}_p)$ yields the residual aberration transmittance $\delta \tilde{H}_L(\mathbf{r}_p)$ while the first input singular vector $V_1(\mathbf{r}_p)$ directly provides the medium reflectivity over the spatial window $W_L$.

However, despite the CLASS correction, the input PSF $\delta H_L$ remains far from being point-like (Fig. 6c). The spectrum of $D'(\mathbf{r}_p)$ then displays a continuum of singular values but the first eigenstate of $D'(\mathbf{r}_p)$ is still of interest. $V_1(\mathbf{r}_p)$ corresponds to a rough estimate of the medium reflectivity that allows realignment in phase for each input focal spot. Therefore, the SVD process allows the synthesis of a coherent virtual reflector that can be leveraged for the estimation of the residual aberration transmittance $\delta \tilde{H}_L(\mathbf{r}_p)$ (Fig. 7c). More precisely, this is the normalized output singular vector $U_1(\mathbf{r}_p) = [U_1(\mathbf{k}_{\text{out}}, \mathbf{r}_p)/|U_1(\mathbf{k}_{\text{out}}, \mathbf{r}_p)|]$, that constitutes a relevant estimator for
$\delta \tilde{\mathbf{H}}_L(\mathbf{r}_p)$ [25]. The estimator of the transmission matrix is then given by the Hadamard product:

$$ \hat{\mathbf{T}}_{\text{out}}(\mathbf{r}_p) = \mathbf{P}'_0 \circ e^{i[\phi_D^{(\text{out})}(\mathbf{r}_p) + \phi_C^{(\text{out})}(z_p)]} \quad (35) $$

with $\phi_D^{(\text{out})}(\mathbf{r}_p)$, the phase of $U_1(\mathbf{r}_p)$. The phase conjugate of $\hat{\mathbf{T}}_{\text{out}}$ provides the focusing laws to compensate for the output phase distortions over each patch (Fig. 7d). The same method can be repeated by exchanging the focused and Fourier bases between input and output in order to estimate the transmission matrix $T_{\text{in}}$ [25]. The whole process is iterated once to refine the estimation of $T_{\text{out}}$ and $T_{\text{in}}$.

## D. Transmission matrix estimator

The input phase laws obtained at the end of the aberration correction process are displayed in Fig. 6e for the central regions of the field-of-view highlighted in Fig. 6d. Although they show some similar features (in particular the low spatial frequency components), they also display some differences that are quantified by the correlation coefficient between the different phase masks with the central one. The value of this coefficient is reported below each phase mask. This correlation coefficient goes from 0.86 for the closest spatial windows to 0.39 for the furthest ones. One can also notice that clear differences in the phase behaviour can be observed between the north, center and the south of the field-of-view. The presence of these lateral differences is consistent with the three geological blocks in the region (Fig. 1b). The latter observation together with the correlation coefficient value show the importance of estimating a different phase law for each area and justifies the implementation of a local aberration correction process.

## E. Local compensation of spatially-distributed aberrations

Using $\hat{\mathbf{T}}_{\text{out}}$ and $\hat{\mathbf{T}}_{\text{in}}$, a corrected FR matrix can be finally obtained:

$$ \mathbf{R}_{rr}^{(D)}(z) = \hat{\mathbf{T}}_{\text{out}}^{\dagger} \times \mathbf{R}_{kk}(z) \times \hat{\mathbf{T}}_{\text{in}}^{*} \quad (36) $$

The corresponding confocal image and RPSF are displayed at $z = 25$ km in Figs. 6(d) and (f). The comparison with their CLASS counterparts [Figs 6(a) and (c)] shows the importance of the local D-matrix analysis. The diffuse background is clearly reduced and the RPSF is nearly similar
to its ideal value (Fig. 3d), with almost all the backscattered energy contained in the white circle accounting for the diffraction limit. The residual background in Fig. 6f is probably associated with high-order aberrations whose coherence length (isoplanatic area) is smaller than the size *L* of the window function *W*(*ρ*).

<table><caption>TABLE II. Confocal gain and resolution at each step of the aberration correction process.</caption><tr><td>Correction steps</td><td>0</td><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td></tr><tr><td>Correction type</td><td>0</td><td>CLASS</td><td>CLASS</td><td>D</td><td>D</td><td>D</td><td>D</td></tr><tr><td>Correction side</td><td></td><td>Output</td><td>Input</td><td>Output</td><td>Input</td><td>Output</td><td>Input</td></tr><tr><td>Confocal gain (dB)</td><td></td><td>2.41</td><td>5.17</td><td>7.1</td><td>9</td><td>9.16</td><td>9.26</td></tr><tr><td>Resolution *w* (km)</td><td>40</td><td>20</td><td>8</td><td>7</td><td>6</td><td>6</td><td>6</td></tr></table>

To be more quantitative, a confocal gain can be computed from the intensity ratio between the corrected (Fig. 6a and d) and initial (Fig. 3e) images. The transverse resolution can also be estimated from the full width at half maximum *w* of the RPSF. The confocal gain and the resolution are reported in Table II at each step of the aberration correction process for depth *z* = 25 km. Strikingly, the transverse resolution is enhanced by a factor ~ 7 compared with its initial value and the confocal intensity is increased by more than 9 dB. These values highlight the benefit of matrix imaging for in-depth probing of NAFZ at a large scale.

In the next section, the 3D image of the medium around the NAFZ is now revealed by combining the images derived at each depth. A structural interpretation is then provided in light of previous studies on the NAFZ.

## V. 3D structure of the NAFZ

The previous sections have shown the process for a local compensation of phase distortions. Performing this correction process at each depth allows to uncover a well-resolved 3D image of the subsurface.

Figure 8a shows a North-South cross-section from the final 3D image. This cross-section is chosen at the same location as the one in Fig. 4a and crosses the two fault strands. It also spans the three geological units: Istanbul zone (IZ), Armutlu-Almacik (AA) and Sakarya zone (SZ). The scattering generated by the heterogeneities of the medium induce a decrease of the backscattered energy with depth. Consequently, a drop of amplitude is observed in the 3D images. In order to compensate for this, the intensity in the cross-sections is normalized by the mean intensity
<figure><img src="image_10.png" /></figure>

FIG. 8. Final image of NAFZ. (a) Vertical North-South cross-section at 30.37°E. The North-South profile is oriented perpendicular to the fault traces. The location of the profile is shown in Fig. 1a. The locations of the southern (SNAF) and northern (NNAF), and the major crustal blocks (SZ: Sakarya zone, AA: Armutlu-Almacik and IZ: Istanbul zone) are labeled. The interpreted location of the fault at depth are indicated by a red line. The color scale refers to the scattering intensity. It is normalized by the maximum value of the scattering energy inside the volume. Our interpretation of the Moho's location is indicated by red dashed lines. (b) Depth slices retrieved from the 3D scattering volume at *z* = 15, 30 and 40 km with (c) their corresponding RPSFs.

calculated at each depth.

Three depth slices retrieved from the final 3D images at *z* = 15, 30 and 40 km are also represented in Fig. 8b with their corresponding RPSFs in Fig. 8c. Compared to the initial RPSFs (Figs. 4c), the resolution is significantly improved by a factor that increases from 7 at small depth (*z* < 15 km) to 9 beyond *z* = 40 km. The final matrix images in Fig. 8b can also be compared to the initial confocal ones in Fig. 4b. While the original images show random-like features, the images obtained after aberration correction reveal better defined features and a reflectivity that is mainly concentrated in the North. The same observation can be made by looking at the corrected cross-section in Fig. 8a. The differences between the corrected and raw cross-sections are pronounced. While in the raw image, no clear structures and layers are visible, the corrected image reveals sub-horizontal structures with a refined level of details, thanks to the drastic gain in resolution revealed by the RPSF.

Due to its significant seismic activity over the past 100 years, and to assess the ongoing hazard
posed by this activity, extensive research has been conducted on the NAFZ to image its structure and determine its mechanical characteristics. The scattering structure in Fig. 8a is interpreted with reference to prior studies conducted in the region.

The first thing to notice in the profiles is that the scattered energy is predominantly situated in the North, which corresponds with the location of the Northern branch. This observation may be associated with the greater seismic activity of the Northern strand compared to the seismic activity of the Southern strand. The scattering below this strand and to the North of it, that extends to at least 60 km, can be explained by the damage caused by the large deformation of this complex fault system with a cumulative slip of the order of 80 km [58, 59] during the last million years as well as the heterogeneities that have been inherited from the complex tectonic history of the region.

At the east of the Sea of Marmara, the Moho depth was reported to be between 30 and 35 km [60, 61]. A deepening of the Moho was identified (~40 km) in the IZ by [62], [63], [36], [37], [64] and [65]. In Fig. 8a, a high scattering zone is observed between 25 and 40 km depth corresponding to a heterogeneous lower crust. Its lower boundary indicates the presence of the Moho (red dashed line). The Moho depth varies from 35 km in the South to 42 km in the North. The reflectivity is disrupted around 40.75°N suggesting the presence of a step in the Moho below the Northern strand. The latter observation is in agreement with previous studies [64, 65]. Below the Moho, reflective structures are observed, mainly beneath AA and IZ, in agreement with [49]. These findings, supported with other studies [37, 49, 64, 65], suggest that the NNAF cuts through the entire crust and reaches the upper mantle [37, 49, 64, 65]. A signature of the NAFZ in the mantle has also been proposed by the long period analysis of [38].

The signature of the Northern strand at depth can be identified by the presence of discontinuities in the scattering distribution in the first 20 km of the crust (Fig. 8a) and also by the termination of sub-Moho structures below the Northern strand. The Southern strand, on the other hand, lacks significant scattering, indicating that it has a weaker signal compared to the Northern strand. This, along with the continuity of the Moho in the South, suggests that the SNAF is confined in the crust and does not extend to the upper mantle, Armutlu block being a crustal structure.

In this section, only one cross-section has been depicted to demonstrate the significant enhancements and the gain in resolution provided by the presented matrix approach. A more in-depth analysis of the scattering volume around the NAFZ will be provided in a future study.
## VI. Conclusion

Matrix imaging provides unprecedented view of the NAFZ. To that aim, we exploited seismic noise data from a dense deployment over the rupture region of the 1999 Izmit earthquake. Ambient noise cross-correlations enable the passive measurement of the reflection matrix associated with the dense array of geophones. The body wave component is then used to image the in-depth reflectivity of the NAFZ subsurface. Compared with our previous work that considered a sparse scattering medium [21], the NAFZ case is more general since it exhibits both specular reflectors such as Moho discontinuity and a random distribution of heterogeneities.

The strength of matrix imaging lies in the fact that it does not require an accurate velocity model. Here, a layered velocity model is employed but strong phase distortions subsist since lateral variations of the wave velocity are not taken into account. Nevertheless, such complex aberrations are compensated by two matrix methods previously developed in optical microscopy [18, 24] and ultrasound imaging [16, 25]. First, the CLASS algorithm exploits angular correlations and memory effect exhibited by the reflection matrix to compensate for spatially-invariant aberrations. Second, a local analysis of the distortion matrix enables a local compensation of spatially-distributed aberrations. Together, those two approaches provide a sharp estimate of the transmission matrix between the Earth surface and the subsurface, leading to a narrowing of the imaging PSF by a factor that goes from 7 to 9. Therefore, a diffraction-limited resolution is reached for any pixel of the image.

Thanks to matrix imaging, the scattering structure of the crust and upper mantle of the NAFZ continental strike slip fault is thus revealed. The 60 km depth profile, show terminations of crustal discontinuities mainly below the northern branch. The localized scattering around the NNAF is consistent with the fact that it is the most seismically active fault and that it ruptured during the last 7.6 Izmit earthquake. We identify a step in the Moho coinciding with the surface location of this branch in the East of DANA network. Moreover, the scattering extends to the upper mantle in the North. All these observations are consistent with previous studies and suggest that the NNAFZ is localized in the crust and extends to the upper mantle.

Even though the result are promising, several points remain that would allow improved images. First, potential conversion between S and P-waves is not considered by matrix imaging. The method could be improved in the future by considering both longitudinal and shear waves, as well as wave conversion between them. Second, only a broadband compensation of phase distortions
is performed. Yet, scattering phenomena or multiple reflections would require a procedure that
moves beyond the application of simple time delays to the impulse response between geophones.
Finally, a reflectivity image is only qualitative since it does not directly quantify the mechanical
properties of the subsurface. Yet matrix imaging offers the possibility of mapping the velocity
distribution inside the medium [16]. This will be the focus of a future study.

Acknowledgments

This project has received funding from the European Research Council (ERC) under the Eu-
ropean Union Horizon 2020 research and innovation program (grant agreement No 742335, F-
IMAGE and grant agreement No. 819261, REMINISCENCE).

Open Research

The data used for this study were recorded by the temporary Dense Array for North Anatolia (DANA) and can be found at [26]. The cross-correlation data and codes used to post-process the seismic data within this paper have been deposited at [66].
# Supplementary Information

This supplementary material includes further details about: (i) the reflection point spread function; (ii) the CLASS algorithm, (iii) the study of the reflection matrix in the plane wave basis to determine the nature of the scattering process.

## S1. Reflection point spread function

In this Supplementary Section, an analytical expression of the RPSF is derived in the specular and scattering regimes. To do so, a local isoplanatic assumption is made (Eq. 17 of the accompanying paper). Under this assumption, the wave-field along each antidiagonal of the focused R-matrix can be rewritten as follows (Eq. 12)

$$R(\rho_m - \Delta\rho/2, \rho_m + \Delta\rho/2, z) = \int d\rho' H_L(\rho' - \Delta\rho/2, \rho_m, z) \gamma(\rho' + \rho_m, z) H_L(\rho' + \Delta\rho/2, \rho_m, z). \quad (S1)$$

with $\rho' = \rho - \rho_m$. Injecting Eq. 12 into Eq. 16 of the accompanying paper leads to the following expression for the RPSF:

$$I(\Delta\rho, \rho_m, z) = \int \int d\rho'_1 d\rho'_2 H_L(\rho'_1 - \Delta\rho/2, \rho_m, z) H_L^*(\rho'_2 - \Delta\rho/2, \rho_m, z) \\ \times \gamma(\rho'_1 + \rho_m, z) \gamma^*(\rho'_2 + \rho_m, z) \\ \times H_L(\rho'_1 + \Delta\rho/2, \rho_m, z) H_L^*(\rho'_2 + \Delta\rho/2, \rho_m, z). \quad (S2)$$

To go beyond this general expression, one can consider two asymptotic regimes depending on the relative values between the characteristic length scale $\ell_\gamma$ of the reflectivity $\gamma(\rho, z)$ at the ballistic depth and the typical width $\delta\rho$ of the focal spots.

In the specular scattering regime, the characteristic size $\ell_\gamma$ of reflectors is larger than the width $\delta\rho$ of each focal spot. $\gamma(\rho, z)$ can thus be assumed as invariant over the PSF support $(\gamma(\rho_m + \rho', z) \simeq \gamma(\rho_m, z))$, such that:

$$I(\Delta\rho, \rho_m, z) = |\gamma(\rho_m, z)|^2 \left| \int d\rho' H_L(\rho' - \Delta\rho/2, \rho_m, z) H_L(\rho' + \Delta\rho/2, \rho_m, z) \right|^2 \quad (S3)$$

$$= |\gamma(\rho_m, z)|^2 \times \left| H_L \overset{\Delta\rho}{\otimes} H_L(\Delta\rho, \rho_m, z) \right|^2 \quad (S4)$$
In the scattering regime ($\ell_{\gamma} < \delta\rho$), the medium reflectivity can be considered, in first approximation, as random:

$$ \langle \gamma(\rho_1, z) \gamma^*(\rho_2, z) \rangle = \langle |\gamma|^2 \rangle \delta(\rho_2 - \rho_1), \qquad (S5) $$

where $\delta$ is the Dirac distribution. The ensemble average of the RPSF (Eq. S2) is given by:

$$ \begin{aligned} \langle I(\Delta\rho, \rho_m, z) \rangle = & \iint d\rho'_1 d\rho'_2 H_L(\rho'_1 - \Delta\rho/2, \rho_m, z) H_L^*(\rho'_2 - \Delta\rho/2, \rho_m, z) \\ & \times \langle \gamma(\rho'_1 + \rho_m, z) \gamma^*(\rho'_2 + \rho_m, z) \rangle \\ & \times H_L(\rho'_1 + \Delta\rho/2, \rho_m, z) H_L^*(\rho'_2 + \Delta\rho/2, \rho_m, z). \end{aligned} \qquad (S6) $$

By combining the previous equation with Eq. S5, the ensemble average of $I(\Delta\rho, \rho_m, z)$ can be expressed as follows:

$$ \langle I(\Delta\rho, \rho_m, z) \rangle = \langle |\gamma|^2 \rangle \times \left[ |H_L|^2 \stackrel{\Delta\rho}{\otimes} |H_L|^2 \right] (\Delta\rho, \rho_m, z). \qquad (S7) $$

## S2. CLASS algorithm

As stated in the accompanying paper, a full-field phase correction is first applied to the $R_{kk}$-matrix through the CLASS algorithm [18, 23].

In order to prove that $C(\mathbf{k}_+, z)$ is an estimator of the spatial frequency spectrum $\tilde{\gamma}(\mathbf{k}_+, z)$ of the medium reflectivity, one can inject the expression of $R(\mathbf{k}_{\text{out}}, \mathbf{k}_{\text{in}}, z)$ (Eq. 24 of the accompanying paper) into Eq. 26:

$$ C(\mathbf{k}_+, z) = \tilde{\gamma}(\mathbf{k}_+, z) \sum_{\mathbf{k}_{\text{out}}} \tilde{H}_I(\mathbf{k}_+ - \mathbf{k}_{\text{out}}, z) \tilde{H}_I(\mathbf{k}_{\text{out}}, z). \qquad (S8) $$

$C(\mathbf{k}_+, z)$ is thus equal to $\tilde{\gamma}(\mathbf{k}_+, z)$ modulated by the autocorrelation function of the aberration transmittance $\tilde{H}_I$. In first approximation, the latter quantity is real and the phase of $C(\mathbf{k}_+, z)$ is actually an estimator of the phase of $\tilde{\gamma}(\mathbf{k}_+, z)$.

In order to prove that $\phi_C(\mathbf{k}_{\parallel}, z)$ is actually an estimator of $\arg\{\tilde{H}_I(\mathbf{k}_{\parallel}, z)\}$ and to determine its bias, $R_{kk}$ can be replaced by its expression (Eq. 24 of the accompanying paper) in Eqs. 26-28
of the accompanying paper. It yields the following expression for $\phi_C$:

$$ \phi_C(\mathbf{k}_{\text{out}}, z) = \arg \left[ \tilde{H}_L(\mathbf{k}_{\text{out}}, z) \right] \\ + \arg \left[ \sum_{\mathbf{k}_{\text{in}}} |\tilde{\gamma}(\mathbf{k}_{\text{out}} + \mathbf{k}_{\text{in}}, z)|^2 \tilde{H}_L(\mathbf{k}_{\text{in}}, z) \sum_{\mathbf{k}'} \tilde{H}_L(\mathbf{k}', z) \tilde{H}_L^*(\mathbf{k}_{\text{out}} + \mathbf{k}_{\text{in}} - \mathbf{k}', z) \right]. \quad (S9) $$

The last expression shows that the estimator $\phi_C$ can be decomposed as a sum of its expectation $\arg[\tilde{H}_L(\mathbf{k}_{\text{out}}, z)]$ and its bias. For a medium of random reflectivity, the term $|\tilde{\gamma}(\mathbf{k}_{\text{out}} + \mathbf{k}_{\text{in}}, z)|^2$ can be replaced by its ensemble average, i.e a constant. It yields:

$$ \phi_C(\mathbf{k}_{\text{out}}, z) = \arg \left[ \tilde{H}_L(\mathbf{k}_{\text{out}}, z) \right] \\ + \arg \left[ H_L(\mathbf{0}, z) \sum_{\mathbf{k}'} \tilde{H}_L(\mathbf{k}', z) \tilde{H}_L^*(\mathbf{k}_{\text{out}} + \mathbf{k}_{\text{in}} - \mathbf{k}', z) \right]. \quad (S10) $$

The last expression shows that the bias directly depends on the autocorrelation of the aberration phase law. The more complex the aberration is, the more biased its estimator is. It is equivalent to the bias exhibited by standard adaptive focusing methods induced by the blurring of a virtual guide star induced by focusing [25].

To prove that the CLASS operation amounts to maximize the confocal intensity in the focused basis, one can express the diagonal coefficients of $R_{\rho\rho}^{(C)}$ as follows:

$$ \begin{aligned} R^{(C)}(\boldsymbol{\rho}_c, \boldsymbol{\rho}_c, z) &= \sum_{\mathbf{k}_{\text{out}}} \sum_{\mathbf{k}_{\text{in}}} R^{(C)}(\mathbf{k}_{\text{out}}, \mathbf{k}_{\text{in}}, z) e^{i(\mathbf{k}_{\text{in}} + \mathbf{k}_{\text{out}}) \cdot \boldsymbol{\rho}_c} \\ &= \sum_{\mathbf{k}_+} C'(\mathbf{k}_+) e^{i\mathbf{k}_+ \cdot \boldsymbol{\rho}_c} \end{aligned} \quad (S11) $$

with $C'(\mathbf{k}_+, z) = \sum_{\mathbf{k}_{\text{out}}} R^{(C)}(\mathbf{k}_{\text{out}}, \mathbf{k}_+ - \mathbf{k}_{\text{out}}, z)$, the sum of antidiagonal coefficients of $R_{\mathbf{k}\mathbf{k}}^{(C)}$. The confocal image is thus the Fourier transform of $C'$. By realigning the phase of $R_{\mathbf{k}\mathbf{k}}$'s coefficients located along the same antidiagonal, CLASS maximizes the intensity of $C'$ and thus of the confocal intensity, by virtue of Parseval identity:

$$ \sum_{\boldsymbol{\rho}_c} |R^{(C)}(\boldsymbol{\rho}_c, \boldsymbol{\rho}_c, z)|^2 \equiv \sum_{\mathbf{k}_+} |C'(\mathbf{k}_+)|^2. \quad (S12) $$
### S3. Nature of the scattering process

The aberration correction process depends on the scattering regime we are facing. To determine it, the plane wave basis is particularly adequate [17]. This section shows how the $R_{kk}$-matrix (Eq. 14 of the accompanying paper) can indicate the nature of the scattering processes taking place in NAFZ.

Indeed, assuming that the mismatch between the wave velocity model and reality only induces phase distortions between plane waves ($|\tilde{H}(k)|=1$), the norm-square of $R_{kk}$-coefficients, $R(k_{\text{out}}, k_{\text{in}}, z)$, is shown to be independent of aberrations [17]:

$$ |R(k_{\text{out}}, k_{\text{in}}, z)|^2 = |\tilde{\gamma}(k_{\text{out}} + k_{\text{in}}, z)|^2. \qquad (\text{S13}) $$

Each anti-diagonal of $R_{kk}$ ($k_{\text{out}} + k_{\text{in}} = \text{constant}$) encodes one spatial frequency of the medium's reflectivity. The spatial frequency spectrum of the medium's reflectivity can be estimated by averaging the intensity of the backscattered wave-field along each anti-diagonal of $R_{kk}$. The result is displayed in Figs. S1b and d at two different depths: $z = 25$ km and $z = 35$ km. The norm square of the spatial frequency spectrum $\tilde{\gamma}(k_{\parallel}, z)$ reveals the nature of the scattering process inside the medium. At $z = 25$ km, the frequency spectrum shows an almost flat spatial frequency spectrum (Fig. S1(b)) which is a manifestation of distributed heterogeneities (Fig. S1(a)). This regime is often referred to as a speckle wave-field in ultrasound imaging [17]. At depth $z = 35$ km, $\tilde{\gamma}(k_{\parallel}, z)$ still shows a flat background due to randomly distributed heterogeneities but it also exhibits an over-intensity in the vicinity of $k_{\parallel} = 0$ (Fig. S1(d)). This peak in the spatial frequency domain is characteristic of a specular reflector (Fig. S1(c)) at that depth that may be associated with an interface between lower crust layers. From these two examples, we can see that the subsurface of NAFZ consists of a mix between specular reflectors and distributed heterogeneities.
<figure><img src="image_11.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) Specular reflector</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_14.png" /></figure>

FIG. S1. Reflection matrix in the plane wave basis. (a) Sketch showing the angular decomposition of the reflected wave-field in the speckle regime for a plane wave illumination (blue). A set of plane waves (green) are reflected in all directions. (b) Spatial frequency spectrum of the reflectivity (Eq. S13) at z = 25 km. (c) Sketch showing the angular decomposition of the wave-field reflected by a planar interface. The incident plane wave (blue) is reflected with the same angle (green), such that $k_{\text{out}} + k_{\text{in}} = 0$; (d) Spatial frequency spectrum of the reflectivity (Eq. S13) at z = 35 km.
[1] A. Sengör, The North Anatolian transform fault: Its age, offset and tectonic significance, J Geol. Soc. London **136**, 269 (1979).

[2] A. Barka, The North Anatolian fault zone, in *Annales tectonicae*, Vol. 6 (1992) pp. 164–195.

[3] N. N. Ambraseys and C. Finkel, *The seismicity of Turkey and adjacent areas: A historical review*, 1500-1800 (MS Eren, 1995).

[4] R. S. Stein, A. A. Barka, and J. H. Dieterich, Progressive failure on the North Anatolian fault since 1939 by earthquake stress triggering, *Geophys. J. Int.* **128**, 594 (1997).

[5] A. Vauchez, A. Tommasi, and D. Mainprice, Faults (shear zones) in the Earth’s mantle, *Tectonophysics* **558**, 1 (2012).

[6] V. Lyakhovsky and Y. Ben-Zion, Evolving geometrical and material properties of fault zones in a damage rheology model, *Geochem. Geophys. Geosystems* **10** (2009).

[7] M. Campillo and P. Roux, Seismic imaging and monitoring with ambient noise correlations, *Treatise on Geophysics* **1**, 256 (2014).

[8] M. Campillo and A. Paul, Long-range correlations in the diffuse seismic coda, *Science* **299**, 547 (2003).

[9] N. M. Shapiro and M. Campillo, Emergence of broadband rayleigh waves from correlations of the ambient seismic noise, *Geophys. Res. Lett.* **31** (2004).

[10] D. Draganov, K. Wapenaar, W. Mulder, J. Singer, and A. Verdel, Retrieval of reflections from seismic background-noise measurements, *Geophys. Res. Lett.* **34** (2007).

[11] P. Poli, H. Pedersen, and M. Campillo, Emergence of body waves from cross-correlation of short period seismic noise, *Geophys. J. Int.* **188**, 549 (2012a).

[12] P. Poli, M. Campillo, H. Pedersen, and LAPNET Working Group, Body-wave imaging of Earth’s mantle discontinuities from ambient seismic noise, *Science* **338**, 1063 (2012b).

[13] L. Retailleau, P. Boué, L. Li, and M. Campillo, Ambient seismic noise imaging of the lowermost mantle beneath the North Atlantic Ocean, *Geophys. J. Int.* **222**, 1339 (2020).

[14] D. Zigone, Y. Ben-Zion, M. Lehujeur, M. Campillo, G. Hillers, and F. L. Vernon, Imaging subsurface structures in the San Jacinto fault zone with high-frequency noise recorded by dense linear arrays, *Geophys. J. Int.* **217**, 879 (2019).

[15] R. Qian and L. Liu, Imaging the active faults with ambient noise passive seismics and its application
to characterize the Huangzhuang-Gaoliying fault in Beijing Area, northern China, Eng. Geol. **268**,
105520 (2020).

[16] W. Lambert, L. A. Cobus, M. Couade, M. Fink, and A. Aubry, Reflection matrix approach for quantitative imaging of scattering media, Phys. Rev. X **10**, 021048 (2020a).

[17] W. Lambert, L. A. Cobus, T. Frappart, M. Fink, and A. Aubry, Distortion matrix approach for ultrasound imaging of random scattering media, Proc. Nat. Sci. Acad. **117**, 14645 (2020b).

[18] S. Kang, P. Kang, S. Jeong, Y. Kwon, T. D. Yang, J. H. Hong, M. Kim, K.-D. Song, J. H. Park, J. H. Lee, M. J. Kim, K. H. Kim, and W. Choi, High-resolution adaptive optical imaging within thick scattering media using closed-loop accumulation of single scattering, Nat. Commun. **8**, 1 (2017).

[19] A. Badon, V. Barolle, K. Irsch, A. C. Boccara, M. Fink, and A. Aubry, Distortion matrix concept for deep imaging in optical coherence microscopy, Sci. Adv. **6**, eaay7170 (2020).

[20] T. Blondel, J. Chaput, A. Derode, M. Campillo, and A. Aubry, Matrix approach of seismic imaging: Application to the Erebus volcano, Antarctica, J. Geophys. Res.: Solid Earth **123**, 10,936 (2018).

[21] R. Touma, T. Blondel, A. Derode, M. Campillo, and A. Aubry, A distortion matrix framework for high-resolution passive seismic 3-D imaging: Application to the San Jacinto fault zone, California, Geophys. J. Int. **226**, 780 (2021).

[22] E. Giraudat, A. Burtin, and A. Aubry, Passive seismic matrix imaging of La Soufrière of Guadeloupe volcano, in *EGU General Assembly Conference Abstracts* (2021) pp. EGU21-4516.

[23] C. Choi, K.-D. Song, S. Kang, J.-S. Park, and W. Choi, Optical imaging featuring both long working distance and high spatial resolution by correcting the aberration of a large aperture lens, Sci. Rep. **8**, 1 (2018).

[24] S. Yoon, H. Lee, J. H. Hong, Y.-S. Lim, and W. Choi, Laser scanning reflection-matrix microscopy for aberration-free imaging through intact mouse skull, Nat. Commun. **11**, 5721 (2020).

[25] W. Lambert, L. A. Cobus, J. Robin, M. Fink, and A. Aubry, Ultrasound matrix imaging – Part II: The distortion matrix for aberration correction over multiple isoplanatic patches, IEEE Trans. Med. Imag. **41**, 3921 (2022b).

[26] DANA, Dense Array for North Anatolia (DANA) [Online Dataset]. International Federation of Digital Seismograph Networks. (2012).

[27] A. Barka, H. Akyuz, E. Altunel, G. Sunal, Z. Cakir, A. Dikbas, B. Yerli, R. Armijo, B. Meyer, J. De Chabalier, et al., The surface rupture and slip distribution of the 17 august 1999 Izmit earthquake (M 7.4), North Anatolian fault, Bull. Seismol. Soc. Am **92**, 43 (2002).
[28] H. Akyuz, R. Hartleb, A. Barka, E. Altunel, G. Sunal, B. Meyer, and V. R. Armijo, Surface rupture and slip distribution of the 12 november 1999 Duzce earthquake (m 7.1), North Anatolian fault, Bolu, Turkey, Bull. Seismol. Soc. Am **92**, 61 (2002).

[29] R. Reilinger, S. McClusky, P. Vernant, S. Lawrence, S. Ergintav, R. Cakmak, H. Ozener, F. Kadirov, I. Guliev, R. Stepanyan, et al., Gps constraints on continental deformation in the Africa-Arabia-Eurasia continental collision zone and implications for the dynamics of plate interactions, J. Geophys. Res.: Solid Earth **111** (2006).

[30] N. Ambraseys, Seismic sea-waves in the Marmara Sea region during the last 20 centuries, Journal of Seismology **6**, 571 (2002).

[31] Y. Yılmaz, Ş. Genç, E. Yi?itbaş, M. Bozcu, and K. Yılmaz, Geological evolution of the late Mesozoic continental margin of Northwestern Anatolia, Tectonophysics **243**, 155 (1995).

[32] A. I. Okay and O. Tüysüz, Tethyan sutures of northern Turkey, Geological Society, London, Special Publications **156**, 475 (1999).

[33] F. Chen, W. Siebel, M. Satir, M. Terzioğlu, and K. Saka, Geochronology of the Karadere basement (NW Turkey) and implications for the geological evolution of the Istanbul zone, Int. J. Earth Sci. **91**, 469 (2002).

[34] M. K. Salah, S. Sahin, and M. Kaplan, Seismic velocity structure along the western segment of the North Anatolian fault zone imaged by seismic tomography, Bull. Earthq. Res. Inst. Univ. Tokyo **82**, 209 (2007).

[35] I. Koulakov, D. Bindi, S. Parolai, H. Grosser, and C. Milkereit, Distribution of seismic velocities and attenuation in the crust beneath the North Anatolian Fault (Turkey) from local earthquake tomography, Bull. Seismol. Soc. Am. **100**, 207 (2010).

[36] E. Papaleo, D. G. Cornwell, and N. Rawlinson, Seismic tomography of the North Anatolian Fault: New insights into structural heterogeneity along a continental strike-slip fault, Geophysi. Res. Lett. **44**, 2186 (2017).

[37] E. Papaleo, D. Cornwell, and N. Rawlinson, Constraints on North Anatolian Fault Zone width in the crust and upper mantle from S wave teleseismic tomography, J. Geophys. Res.: Solid Earth **123**, 2908 (2018).

[38] A. Fichtner, E. Saygin, T. Taymaz, P. Cupillard, Y. Capdeville, and J. Trampert, The deep structure of the North Anatolian fault zone, Earth Planet. Sci. Lett. **373**, 109 (2013).

[39] Y. Çubuk-Sabuncu, T. Taymaz, and A. Fichtner, 3-D crustal velocity structure of western Turkey:
Constraints from full-waveform tomography, Phys. Earth Planet. Inter. **270**, 90 (2017).

[40] G. Taylor, S. Rost, G. A. Houseman, and G. Hillers, Near-surface structure of the North Anatolian Fault zone from Rayleigh and Love wave tomography using ambient seismic noise, Solid Earth **10**, 363 (2019).

[41] G. Bensen, M. Ritzwoller, M. Barmin, A. L. Levshin, F. Lin, M. Moschetti, N. Shapiro, and Y. Yang, Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, Geophys. J. Int. **169**, 1239 (2007).

[42] M. S. Longuet-Higgins, A theory of the origin of microseisms, Philos. Trans. Royal Soc. A **243**, 1 (1950).

[43] K. Hasselmann, A statistical analysis of the generation of microseisms, Rev. Geophys. **1**, 177 (1963).

[44] L. Stehly, M. Campillo, and N. Shapiro, A study of the seismic noise from its long-range correlation properties, J. Geophys. Res.: Solid Earth **111** (2006).

[45] Ö. Emre, T. Y. Duman, S. Özalp, F. Şaroğlu, Ş. Olgun, H. Elmacı, and T. Çan, Active fault database of Turkey, Bull. Earthq. Eng. **16**, 3229 (2018).

[46] K. Akbayram, C. C. Sorlien, and A. I. Okay, Evidence for a minimum 52±1 km of total offset along the northern branch of the North Anatolian Fault in northwest Turkey, Tectonophysics **668**, 35 (2016).

[47] A. Berkhout, Wave field extrapolation techniques in seismic migration, a tutorial, Geophysics **46**, 1638 (1981).

[48] A. J. Berkhout and C. P. A. Wapenaar, A unified approach to acoustical reflection imaging. II: The inverse problem, J Acoust. Soc. Am. **93**, 2017 (1993).

[49] M. Kahraman, D. G. Cornwell, D. A. Thompson, S. Rost, G. A. Houseman, N. Türkelli, U. Teoman, S. A. Poyraz, M. Utkucu, and L. Gülen, Crustal-scale shear zones and heterogeneous structure beneath the north anatolian fault zone, turkey, revealed by a high-density seismometer array, Earth Planet. Sci. Lett. **430**, 129 (2015).

[50] A. E. Karahan, H. Berckhemer, and B. Baier, Crustal structure at the western end of the North Anatolian Fault Zone from deep seismic sounding, Ann. Geophys., **44** (2001).

[51] W. Lambert, J. Robin, L. A. Cobus, M. Fink, and A. Aubry, Ultrasound matrix imaging – Part I: The focused reflection matrix, the F-factor and the role of multiple scattering, IEEE Trans. Med. Imag. **41**, 3907 (2022a).

[52] M. Born and E. Wolf, *Principles of optics* (Seventh edition) (Cambridge University Press, Cambridge, 2003).
[53] A. Aubry and A. Derode, Detection and imaging in a random medium : A matrix method to overcome multiple scattering and aberration, J. Appl. Phys. **106**, 044903 (2009).

[54] S. Kang, S. Jeong, H. Choi, W. Ko, T. D. Yang, J. H. Joo, J.-S. Lee, Y.-S. Lim, Q.-H. Park, and W. Choi, Imaging deep within a scattering medium using collective accumulation of single-scattered waves, Nat. Photonics **9**, 1 (2015).

[55] I. Freund, M. Rosenbluh, and S. Feng, Memory effects in propagation of optical waves through disordered media, Phys. Rev. Lett. **61**, 2328 (1988).

[56] S. Shahjahan, A. Aubry, F. Rupin, B. Chassignole, and A. Derode, A random matrix approach to detect defects in a strongly scattering polycrystal: How the memory effect can help overcome multiple scattering, Appl. Phys. Lett. **104**, 234105 (2014).

[57] U. Najar, V. Barolle, P. Balondrade, M. Fink, A. C. Boccara, M. Fink, and A. Aubry, Non-invasive retrieval of the time-gated transmission matrix for optical imaging deep inside a multiple scattering medium, arXiv: 2303.06119 (2023).

[58] R. Armijo, B. Meyer, A. Hubert, and A. Barka, Westward propagation of the North Anatolian fault into the Northern Aegean: Timing and kinematics, Geology **27**, 267 (1999).

[59] M. Bohnhoff, P. Martínez-Garzón, F. Bulut, E. Stierle, and Y. Ben-Zion, Maximum earthquake magnitudes along different sections of the North Anatolian fault zone, Tectonophysics **674**, 147 (2016).

[60] E. Zor, E. Sandvol, C. Gürbüz, N. Türkelli, D. Seber, and M. Barazangi, The crustal structure of the East Anatolian plateau (Turkey) from receiver functions, Geophys. Res. Lett. **30** (2003).

[61] E. Vanacore, T. Taymaz, and E. Saygin, Moho structure of the Anatolian Plate from receiver function analysis, Geophys. J. Int. **193**, 329 (2013).

[62] A. Frederiksen, D. Thompson, S. Rost, D. Cornwell, L. Gülen, G. Houseman, M. Kahraman, S. Poyraz, U. Teoman, N. Türkelli, et al., Crustal thickness variations and isostatic disequilibrium across the North Anatolian Fault, western Turkey, Geophys. Res. Lett. **42**, 751 (2015).

[63] G. Taylor, S. Rost, and G. Houseman, Crustal imaging across the North Anatolian Fault Zone from the autocorrelation of ambient seismic noise, Geophys. Res. Lett. **43**, 2502 (2016).

[64] S. Rost, G. Houseman, A. Frederiksen, D. Cornwell, M. Kahraman, S. Altuncu Poyraz, U. Teoman, D. Thompson, N. Türkelli, L. Gülen, et al., Structure of the northwestern North Anatolian Fault Zone imaged via teleseismic scattering tomography, Geophys. J. Int. **227**, 922 (2021).

[65] J. Jenkins, S. N. Stephenson, P. Martínez-Garzón, M. Bohnhoff, and M. Nurlu, Crustal thickness variation across the Sea of Marmara region, NW Turkey: A reflection of modern and ancient tectonic
processes, Tectonics **39**, e2019TC005986 (2020).

[66] R. Touma, A. Le Ber, M. Campillo, and A. Aubry, Seismic matrix imaging [Software]. Figshare (2023).