# Monitoring of seismic anisotropy at the time of the 2008 Iwate-Miyagi (Japan) earthquake

M. Saade, J. P. Montagner, P. Roux, K. Shiomi, B. Enescu, F. Brenguier

► To cite this version:

M. Saade, J. P. Montagner, P. Roux, K. Shiomi, B. Enescu, et al. Monitoring of seismic anisotropy at the time of the 2008 Iwate-Miyagi (Japan) earthquake. Geophysical Journal International, 2017, 211, pp.483-497. 10.1093/gji/ggx321. insu-03596048

HAL Id: insu-03596048

https://insu.hal.science/insu-03596048

Submitted on 3 Mar 2022

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
<figure><img src="image_2.png" /></figure>

# Monitoring of seismic anisotropy at the time of the 2008 Iwate-Miyagi (Japan) earthquake

M. Saade,<fnref n="1" /> J.P. Montagner,<fnref n="1" /> P. Roux,<fnref n="2" /> K. Shiomi,<fnref n="3" /> B. Enescu<fnref n="4" /> and F. Brenguier<fnref n="2" />

Accepted 2017 July 28. Received 2017 July 27; in original form 2017 March 1

## SUMMARY

We investigate temporal variations of the polarization anomaly of surface waves (Rayleigh and Love) relative to an isotropic medium, before and after the 2008 June 13, Iwate-Miyagi Nairiku earthquake ($M_w = 6.9$) in Japan, using passive image interferometry. We use continuous data records of 20 high-sensitivity three-component seismic stations (National Institute for Earth Sciences and Disaster Resilience—NIED Hi-net) from mid-2007 to the end of 2009. We compute cross-correlations of ambient seismic noise between each pair of stations, in the frequency range of 0.08–0.22 Hz. The temporal variations of the surface wave polarization show slow and weak variations due to seasonal changes to which rapid and strong variations are superposed at the time of the earthquake, and more surprisingly, also a month and a half before the earthquake. Strong changes in the polarization of surface waves might be related to changes in the distribution of crack orientation affecting seismic anisotropy.

**Key words:** Seismic anisotropy; Seismic interferometry; Seismic noise; Surface waves and free oscillations.

## 1 INTRODUCTION

Over the last decade, the use of ambient seismic noise has been a key aspect of advanced studies of the crustal structure. Seismic interferometry has provided the possibility to image the subsurface and monitor physical parameters in the crust during seismic cycles (Sens-Schönfelder & Wegler 2006; Wegler & Sens-Schönfelder 2007; Brenguier et al. 2008). Indeed, the cross-correlation tensor (CCT) of seismic ambient noise is an estimate of the transfer function (Green's function) between two receivers. For a homogeneous random distribution of noise sources, the CCT represents the impulse response of the Earth between these two receivers (Weaver & Lobkis 2001; Wapenaar 2004; Gouédard et al. 2008). In practice, the CCT of ambient noise is dominated by surface waves, which permits the measurement of the horizontal seismic properties of the shallow crust (Shapiro & Campillo 2004). For instance, the computation of the relative traveltime shift between a perturbed and a reference ZZ component for each station pair (Sens-Schönfelder & Wegler 2006; Brenguier et al. 2008) enables the monitoring of seismic velocity changes in the subsurface structures that might be related to stress changes in seismogenic zones.

For example, Brenguier et al. (2008) investigated the seismic velocity changes for the 2004 Parkfield earthquake. They reported a velocity change of 0.08 per cent at the time of the $M_w 6.0$ earth-

quake, which was interpreted as co-seismic stress changes and damage near the San Andreas fault.

Another approach for monitoring seismogenic zones is to monitor the off-diagonal terms of the CCT: vertical-transverse (ZT), transverse-vertical (TZ), radial-transverse (RT) and transverse-radial (TR). These non-diagonal terms of the Green's function are null when the medium is homogeneous and isotropic, and when the distribution of sources is homogeneous. For anisotropic or strongly heterogeneous media (Crampin et al. 1980), the polarization plane of surface waves will be affected and will deviate with respect to the isotropic homogeneous case, leading to quasi-Rayleigh and quasi-Rayleigh waves. The surface wave polarizations are no longer parallel or perpendicular to the direction of propagation. An optimal rotation algorithm [ORA, Roux (2009)] is then applied to the CCT, and hence the deviation anomaly angles (horizontal $\psi_p$ and vertical $\delta_p$) at both receivers (seismic stations A and B) can be retrieved. These anomaly angles provide information about the anisotropy in elastic properties of the medium or the inhomogeneity of the sources (Durand et al. 2011). Temporal changes in the polarization anomaly of surface waves can be induced by changes in crack-induced anisotropy; hence an overall rotation of the fast axis caused by rotation of the stress field that can close or open some cracks.

Measurements of the changes in surface wave polarization were applied to the 2004 Parkfield earthquake by Durand et al. (2011),

<footnote n="1">¹Seismology Laboratory, Institut de Physique du Globe de Paris, France. E-mail: maria.saade90@gmail.com</footnote>

<footnote n="2">²Institut des Sciences de la Terre, University Joseph Fourier, Grenoble, France</footnote>

<footnote n="3">³National Research Institute for Earth Science and Disaster Resilience, Tsukuba, Japan</footnote>

<footnote n="4">⁴Department of Geophysics, Kyoto University, Japan</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Map of the topography of the Iwate-Miyagi area. Black triangles: Hi-net stations; brown star: epicentre of the 2008 Iwate Miyagi earthquake; red lines: upper edges of the fault proposed by Takada *et al.* (2009); grey dots: aftershocks with magnitudes >3.5 within 7 d of the main earthquake and red triangles: active volcanoes. The beachball plot describes the source mechanism associated with the earthquake. Inset: position of the Iwate-Miyagi area within Japan.</figcaption></figure>

where the fault displacement was a right-lateral strike-slip. Strong and rapid variations of $\psi_p$ were observed for a specific area near the fault zone, with a maximum of 30. Variations in $\delta_p$ are negligible. Durand *et al.* (2011) explained that these rapid changes in the polarization might be due to changes in the orientation of the anisotropy following co-seismic rotation of the stress field. However, shear wave splitting measurements did not show any precursor, co-seismic or post-seismic temporal changes in the region near the rupture of the $M_w = 6.0$ Parkfield earthquake (Liu *et al.* 2007). The shear wave splitting technique used for measuring anisotropy has many limitations due to the vertical propagation of S waves, such as the difficulty to detect changes if the anisotropy is weak or if the anisotropic layer is thin and shallow.

Furthermore, numerical investigations have been carried out to define the effects of anisotropy on the polarization of surface waves, and to estimate the order of magnitude of the changes to this parameter (Saade *et al.* 2015). These data showed that the polarization anomaly $\psi_p$ of surface waves in an anisotropic medium can be very large (i.e. $\psi_p$ can reach as much as 20), and that it mainly contains $2\psi$ and $4\psi$ azimuthal dependence.

The objective of this study is to investigate the relation of co-seismic changes in polarization of surface waves to stress changes before and after an earthquake, using seismic noise interferometry, in a different geodynamic context than Parkfield, such as a compressive environment. Therefore, we selected the 2008 Iwate-Miyagi earthquake ($M_w = 6.9$), the epicentre of which was in the middle of the dense Hi-net seismic network operated by the National Research Institute for Earth Science and Disaster Resilience (NIED, Japan). The focal depth of the epicentre was about 7.8 km. The deduced focal mechanism corresponds to a reverse fault with a compression axis in WNW-ESE direction (Takada *et al.* 2009).

# 2 METHOD

## 2.1 Data

Data were available from 20 three-component, short-period, velocity-type seismometers installed at Hi-net borehole stations located within 50 km of the epicentre of the 2008 Iwate-Miyagi earthquake and exhibiting an approximate interstation distance of 20 km (Fig. 1). We processed continuous data recorded from 2007 July to the end of 2009, that is, almost 1 yr before and 1.5 yr after the main earthquake.

## 2.2 Pre- and post-processing

The different steps of the data processing are presented in the Appendix (see Fig. A1). The pre-data processing is partly similar to the one presented in Bensen *et al.* (2007). First, we applied spectral whitening to the continuous data within the frequency band of 0.08–0.22 Hz (i.e. corresponding to wavelengths of ~10–45 km). The choice of the frequency range is fundamental, since the polarization measurement depends on the amplitude of the noise signal. The chosen frequency domain covers the dominant noise source that comes from secondary microseisms. In this low-frequency domain, we obtain the most stable waveform of the seismic noise CCTs that are sensitive to the deep subsurface structure. Hence, CCTs computed in this frequency range are not particularly influenced by shallow co-seismic damage, especially by the very strong ground displacement that was associated with this 2008 Iwate-Miyagi earthquake (Takada *et al.* 2009). Then, we applied clipping to these data, to cut the remaining event signals. The clipping is smooth as the energy of the signal is reduced by a maximum of 7 per cent during
calm periods (no events) and by a maximum of 48 per cent during the Iwate-Miyagi earthquake.

The CCTs were then computed between the pairs of three-component seismograms. The CCTs are nine-component tensors $C_{AB}^G(t)$ computed in the geographical coordinate system G (vertical [Z], east [E] and north [N]) for each receiver pair (A and B). For the cross-correlation formula, each component of the CCT is normalized with respect to the total energy, according to:

$$[C_{AB}^G(t)]_{ij} = \frac{\int_0^{\Delta T} S_{A,i}(\tau)S_{B,j}(t+\tau)d\tau}{\sqrt{\int_0^{\Delta T} \sum_i S_{A,i}^2(\tau)d\tau \int_0^{\Delta T} \sum_j S_{B,j}^2(\tau)d\tau}}, \quad (1)$$

where $S_{A,i}(t)$ and $S_{B,j}(t)$ are the vertical and horizontal components of the signal at receivers A and B, 'i' is the (Z, E and N) component, and choosing a recording time window $\Delta T = 1$ hr, over which the correlation was performed. The nine components of the CCT as defined by the (Z, E and N) geographical coordinate system are ZZ, ZE, ZN, EZ, EE, EN, NZ, NE and NN, which are rotated in the (Z, R and T) coordinate system, where the radial direction [R] is defined by the orientation of the pair of receivers. We will denote the azimuth of the path between stations A and B by $\psi_{AB}$. The CCT then reads:

$$C_{AB} = \begin{pmatrix} ZZ & ZR & ZT \\ RZ & RR & RT \\ TZ & TR & TT \end{pmatrix}. \quad (2)$$

Unlike velocity measurements that depend on the phase change and are measured using the coda part of the ZZ correlation function, polarization measurements depend on the amplitude change of the ballistic part of all the components of the CCT.

The main part of the procedure is the ORA (Roux 2009) that rotates the CCTs, minimizes the off-diagonal components, and then determines the set of angles $\delta_p$ and $\psi_p$ at each station (Fig. A2), $\delta_{pA}$, $\delta_{pB}$ and $\psi_{pA}$, $\psi_{pB}$ that minimize the components ZT, RT, TZ and TR of the original CCTs (see following details). The ORA computes a final misfit that is the sum of the energy on ZT, RT, TZ and TR of the rotated CCTs divided by the total energy of the tensors. We used the final misfit as a criterion to determine the performance of the ORA, and we sought to minimize the final misfit down to at least 10 per cent of the total energy of the CCT. We did not need to find the Green's tensor; nonetheless, we needed stable CCTs and a good signal-to-noise ratio, to obtain good ORA performance. To do this, we stacked the CCTs, but first we applied a singular value decomposition filter to the CCTs of each pair of receivers, to improve the signal-to-noise ratio and the temporal resolution by optimizing the stack.

Singular value decomposition (Hegde & Ye 2008) decomposes matrix C, which contains one component of a CCT for one pair of receivers over the entire time period of the study. This produces the product of three matrices, $C = U*S*V^T$, where $S(M \times N)$ is the diagonal matrix of the singular values, $U(M \times M)$ and $V(N \times N)$ are orthogonal matrices that consist of singular vectors and T indicates the transpose of the matrix. N is then the number of CCTs for the entire period of the study, and M is the time length of the cross-correlation (200 s). Consequently, for every pair of receivers and each component of the CCT, matrix C was reconstructed by keeping the first four effective singular values that can significantly reproduce the image component. The remaining singular values represent noise. This procedure helped to reduce the number of stacks of CCTs down to 7 d which is necessary in order to improve to signal-to-noise ratio.

Next, we investigated the off-diagonal terms of the CCTs: ZT, TZ, RT and TR. Indeed, for a random distribution of seismic sources, these non-diagonal terms are expected to be zero on average in an isotropic homogeneous medium. However, we observed that the ZT, TZ, RT and TR components of the CCTs were not zero and could be very large (Fig. A3).

Then, we use the ORA, which rotates the CCTs and minimizes the off-diagonal components, to find the tensor of the Rayleigh and Love waves. The ORA determines the set of angles $\delta_p$ and $\psi_p$ (Fig. A2) that minimizes the components ZT, RT, TZ and TR of the original CCTs. $\delta_p$ and $\psi_p$ are the vertical and horizontal angles, respectively, computed for each receiver of a pair. These angles determine the deviation of the plane of surface waves with respect to $\psi_{AB}$ and provide information about the anisotropic properties of the medium and the inhomogeneity of the sources (Durand et al. 2011).

In the general case, the value of the ORA azimuthal rotation angle that minimizes the components ZT, RT, TZ and TR of the CCTs is the result of different effects, as the azimuth of the receiver pairs, the noise-source incidence (if the sources are not isotropically distributed), and an additional angle that can arise from deviations of the Rayleigh and Love waves in the heterogeneous and anisotropic medium. $\psi_{pA}$ and $\psi_{pB}$ can be different if the source incidence is not the same at both of the receivers. The results of this study show that $\psi_{pA}$ was almost equal to $\psi_{pB}$ for all of the receiver pairs that were processed. We then deduced the value of the azimuth of the pairs of receivers from the horizontal angle computed by ORA, and hence the resulting angle is the horizontal polarization anomaly angle (HPA = $\psi_p$).

The temporal variations of $\psi_p$ can be due first to seasonal variations, and secondly to changes in the structure of crack distributions or heterogeneities. To separate these different contributions to the HPA, we applied a singular value decomposition filter to the matrix $P(M' \times N')$ that contained all of the processed $\psi_p(t)$ of all receiver pairs ($P = U*S*V^T$). In this case, $U(M' \times M)$ is the matrix of the station eigenvectors and $V(N' \times N)$ is the matrix of the temporal eigenvectors. $N'$ is the number of pair of receivers processed, and $M'$ is the time length of the period of study. Looking at the main temporal eigenvectors (Fig. 2) helps determine the dominant effects and remove eigenvectors that are dominated by seasonal changes, while keeping those dominated by strong and rapid changes at the time of the 2008 Iwate-Miyagi earthquake. However, in this case, the separation is difficult, knowing that the main first two eigenvectors contain both seasonal variations and rapid variations at the time of the earthquake. Nevertheless, we reconstructed P by keeping the first two eigenvectors, and thus 35 per cent of the energy of $\psi_p(t)$.

## 3 RESULTS AND INTERPRETATION

We computed the temporal changes in the HPA of surface waves for the different pairs of receivers (Fig. A4), from mid-2007 to the end 2009. Generally, the HPA can be different at the two receivers (A and B) of a specific pair; that is, $\psi_{pA}$ is different from $\psi_{pB}$. The near equality of the two angles for a specific pair of receivers indicates that the interstation distances considered (with a maximum of 30 km) are small relative to the distance between the noise source and the pair of receivers, which makes the noise incidence almost equal at the two receivers of the pair. The static value of the HPA that was determined during calm periods, as for example between 2007 October and 2008 March, was slightly different from zero because
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Temporal eigenvectors obtained by singular value decomposition of the horizontal polarization anomalies computed by the optimal rotation algorithm, for all pairs of receivers. In each panel, the value given at the top is the normalized eigenvalue that corresponds to the temporal eigenvector.</figcaption></figure>

it was influenced by the strong heterogeneity under the study area (Cheng *et al.* 2011), and/or by the anisotropy induced by the crack distribution or mineralogical composition in the area.

## 3.1 Seasonal variations

The temporal changes in the HPA of the surface waves reveal their slow, but large, seasonal behaviour. We explored different possible causes for these seasonal variations: the water level changes in the subsurface related to the precipitation changes (rain and snow), and the seasonal change associated with the incidence of ambient seismic noise provided by interactions between the ocean waters and the coast.

These seasonal changes are long-period changes, as the micro-seismic noise incidence is mostly from the Northern Hemisphere in the northern winter, and from the Southern Hemisphere in the southern winter (Stutzmann *et al.* 2012). According to Stutzmann *et al.* (2012), relative to Japan, the dominant noise incidence is mostly directive and the noise amplitude is maximum between November and February. Hobiger *et al.* (2012) also reported strong seasonal variations in the velocity measurements in 2008 and 2009 that had a period of about 1 yr, with a maximum velocity change in mid-June. Fig. 3 shows that the average HPA had seasonal variations with a period of almost 1 yr, and a maximum around August in 2008, and June in 2009. These repetitive slow variations are consistent with seasonal changes observed for the seismic velocity measurements of Hobiger *et al.* (2012), and the period of seasonal changes of the incidence of the seismic noise.

Furthermore, we computed the seasonally averaged azimuthal distributions of the normalized energy flow relative to the network under study (Fig. A5). The measurements of noise incidence revealed changes in the amplitudes of the seismic noise, which was stronger in winter for Japan and came mainly from the Pacific Ocean (azimuth around 90), and also from the Sea of Japan (azimuth between 250 and 300). In summer, the noise sources were distributed mainly in the Pacific Ocean, and also in the Southern Hemisphere (azimuth around 150). The seasonal changes in the noise incidence might have had an influence on the seasonal changes in the HPA, as their annual period was almost the same. The minimum HPA change in winter (as for the velocity change) might be partly due to gravita-

tional compaction caused by the snow coverage and associated with water circulations in aquifers. We note the high negative correlation between HPA change and snow (Fig. 4). In summer, the changes in the HPA increase might have been a result of the snow melt and the rainfall (especially during the rainy season, as the derivative of the precipitation shows in Fig. 4) that reduced the effective pressure in the crust.

## 3.2 Anisotropy variations

Temporal changes in the HPA can also be influenced by changes in the orientation of the crack distribution during the seismic cycle and more particularly the co-seismic phase, which induces changes in the apparent anisotropy. Indeed, changes in the seismic anisotropy can strongly affect the polarization of surface waves (Saade *et al.* 2015) and cause strong and rapid changes in the HPA. These changes can be distinguished from seasonal changes as they are rapid (i.e. over a short timescale) and large.

The results of the HPA measured around the area of Iwate-Miyagi (Fig. 3) also show relatively rapid changes for different receiver pairs at the time of the 2008 Iwate-Miyagi earthquake, and smaller, but more surprisingly, also fast changes before the earthquake, in 2008 April. We computed the derivative of the HPA (Fig. 3), in order to use the rapidity of the change as a criterion to separate the slow variations primarily seasonal and the rapid changes due to the earthquake. Note that a similar analysis was done in Brenguier *et al.* (2008a) to separate the slow seasonal change from the rapid traveltime change of seismic velocity prior to and eruption at the Piton de la Fournaise volcano. Only two peaks in the derivative of the HPA exceed the value of four times the variance of the signal. A maximum peak for the derivative of the HPA occurs at the time of the earthquake (Fig. 3). A second large peak occurred at the end of April, and it was not compatible with any seasonal changes. This second peak might be related to a deep (160 km) earthquake of magnitude 5.7 that occurred on 2008 April 17, in the area covered by the seismic network used in this study (Fig. A6). Changes in the crack-induced anisotropy following the co-seismic rotation of the stress field in the shallow depth are a possible explanation for these rapid changes.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Top panel: average of the horizontal polarization anomaly from mid-2007 to the end of 2009. Bottom panel: the first derivative of the average horizontal polarization anomaly (red signal). The flat values at the end of 2009 May are due to a gap in the data. The variance σ of the average derivative of HPA has been computed for the whole period of study.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Top panel: average change in the horizontal polarization anomaly. Middle panel: seasonal changes in snow depth (from the Japan Meteorological Agency) in the area of study for 2008 and 2009. Bottom panel: average rainfall per week (from the Japan Meteorological Agency).</figcaption></figure>

### 3.2.1 Regionalization of the polarization anomaly

To obtain an estimate of the spatial distribution of the polarization anomalies at the time of the Iwate-Miyagi earthquake, a regionalization technique was used that is based on the forward problem for the effect of slight anisotropy on the amplitudes of Rayleigh and Love waves of Tanimoto (2004), only valid for fundamental

modes. We consider the case of first-order perturbation theory in a homogeneous half-space where the fourth-order elastic tensor has 21 independent elastic coefficients. The tangent of the angle $ψ_p$ formed by the direction of polarization of the quasi-Rayleigh wave and the direction of propagation $ψ$ can be written approximately as a function with azimuthal dependence terms of $2ψ$ and
4ψ, as the elastic tensor is a fourth-order tensor. To get these different azimuthal terms, we modified the regionalization code of Montagner & Nataf (1986).

As shown by Tanimoto (2004), in a weakly anisotropic medium, polarization of the fundamental mode of surface waves contain 2ψ and 4ψ azimuthal dependence. This has also been found in numerical experiments (Saade et al. 2015), with similar azimuthal dependence as the phase velocity variations (Smith & Dahlen 1973). For small values of ψₚ, the polarization anomaly, ψₚ, can then be written as a Fourier series, as follows:

$$ \psi_p = a_1 \cos(2\psi) + b_1 \sin(2\psi) + a_2 \cos(4\psi) + b_2 \sin(4\psi), \quad (3) $$

where ψₚ is the static HPA, ψ is the azimuth of the pair of receivers relative to the north, aₙ = ∫<sup>+π</sup><sub>-π</sub> ψₚ cos(nψ)dψ and bₙ = ∫<sup>+π</sup><sub>-π</sub> ψₚ sin(nψ)dψ. The main steps of the theory used for getting this expression are reminded in Appendix A2.

Considering that the 2ψ terms are dominant for Rayleigh waves, ψₚ can be written as follows:

$$ \psi_p = \alpha \sin[2(\psi - \psi_\alpha)], \quad (4) $$

where ψ_α is the orientation of the anisotropy and α is the amplitude coefficient of ψₚ that depends on the amplitude of anisotropy in the medium. a₁ = -2αsin(2ψ_α) and b₁ = 2αcos(2ψ_α). Hence, ψ_α = 1/2αtan[-a₁/b₁].

The temporal variation of the horizontal polarization anomaly Δψₚ is obtained by deriving eq. (4), knowing that the time-dependent parameters are ψ_α the orientation of anisotropy and α that depends on the amplitude of anisotropy. Δψₚ is then written as follows:

$$ \Delta\psi_p = -2\alpha\Delta\psi_\alpha\cos[2(\psi - \psi_\alpha)] + \Delta\alpha\sin[2(\psi - \psi_\alpha)]. \quad (5) $$

Using the Fourier decomposition as in eq. (3), we obtain:

$$ \Delta\psi_p = [-2\alpha\Delta\psi_\alpha\cos(2\psi_\alpha) - \Delta\alpha\sin(2\psi_\alpha)]\cos(2\psi) + [-2\alpha\Delta\psi_\alpha\sin(2\psi_\alpha) + \Delta\alpha\cos(2\psi_\alpha)]\sin(2\psi), \quad (6) $$

with a = -2αΔψ_αcos(2ψ_α) - Δαsin(2ψ_α) and b = -2αΔψ_αsin(2ψ_α) + Δαcos(2ψ_α).

Consequently, the three possible cases that can give rise to a temporal variation of the HPA are:

(i) Δα = 0 and Δψ_α ≠ 0: in this case, the temporal variation of the HPA Δψₚ is due to a change in the orientation of anisotropy, associated with a change in the stress field causing a variation in the cracks alignment. The orientation of anisotropy in this case is then ψ_α = 1/2αtan[b/a].

(ii) Δα ≠ 0 and Δψ_α = 0: here the value of Δψₚ is due to a change in the amplitude of anisotropy that can also change the polarization anomaly of surface waves. Normally, the polarization anomaly increases with the amplitude of anisotropy in the medium. In this case, ψ_α = 1/2αtan[-a/b].

(iii) Δα ≠ 0 and Δψ_α ≠ 0: the third case is the combination of both effects, ψₚ is due to a change in the amplitude and the orientation of anisotropy. The difficulty here is to separate both effects.

Before moving to inversion, it is important to note that the sensitivity kernels used in the regionalization of anisotropy change is quite complicated. The complexity arises from the fact that the polarization anomaly of surface waves is not cumulative along the path between one source and one receiver, unlike traveltime measurements. Hence, the polarization anomaly should not depend on

the interstation distance and consequently, the latter is not considered in the inversion process at each gridpoint.

We then consider that the sensitivity kernels of the polarization anomaly measured after noise cross-correlation are mostly sensitive to the structure around the two receivers (Fig. A7, Appendix A2). Consequently, for every grid cell of (0.05° × 0.05°) and for all paths passing through the cell, we compute the average of the HPAs assigned with the corresponding weight, assuming that this weight of the polarization measurement is empirically defined as a regionalization distance that decreases away from the receivers. Using this average value, we determine the orientation of the 2ψ anisotropy.

In the inversion process, we limit ourselves to small inter-station distances (around 20 km). Consequently, we get a good coherence between the two receivers. Besides, the polarization measurement, which is not cumulative, will be influenced by the regionalization distance. This way, we measure local effects and we avoid averaging in each grid cell very different measurements influenced by different parts of the crust.

The top panel in Fig. 5 shows the map of the mid-orientation of anisotropy (before, between and after the change) and the amplitude of the HPA change at the time of the Iwate-Miyagi earthquake. The amplitude of the HPA was computed for a time window of 1 week, and its changes were not homogeneous in the region. It was strong for specific areas (Fig. 5, top panel, darker blue) where it reached a maximum of 30. These strong changes are relatively large and surprising, but they were only observed in specific areas at the time of the seismic event, similar to the case of the large HPA change that was measured at the time of the 2004 Parkfield earthquake (Durand et al. 2011). According to eq. (6) and the regionalization results, the HPA change at the time of the Iwate-Miyagi earthquake is mostly influenced by a change of the orientation of anisotropy (Δψ_α ≠ 0). But the amplitude of anisotropy could have also locally changed and affected the HPA, even thought this seems to be a secondary effect. For instance, Nakata & Snieder (2012) found that the 2011 Tohoku-Oki earthquake (M<sub>w</sub> = 9.0) increased the difference between fast and slow shear wave velocities arising from shear wave splitting in most parts of northeastern Japan, but it did not significantly change fast shear wave polarization directions in the near surface.

The largest area of strong HPA changes was near the epicentre and around the fault zone, which consisted of five fault patches as proposed by Takada et al. (2009). The geometry of the fault system is complex and included at least one conjugate faulting. The patches in the southern (SSW) part are characterized by a slip that is larger than the northern (NNE) part, and they extended from the middle portion of the focal region.

We observe a consistent pattern of the orientation of the fast axis of anisotropy. This is because the HPA is sensitive to deep subsurface structure, which is quite simple in this area. This is not the case at other places in Japan, where the shallow subsurface structure is very complex, which influences shear wave splitting measurements (Iidaka et al. 2013). We also note that the presence of active volcanoes in the region indicates the abundance of volcanic fissures with inclusions of highly pressurized fluids.

For instance, in the eastern part of the area under study (Fig. 5, top panel), the fast direction of anisotropy appears to be oriented mainly ESE-WNW. This shows that the orientation of the cracks is mostly influenced by the regional stress field in the area associated with the westward subduction of the Pacific Plate. In the northern and western parts of the area, the fast axis of anisotropy seems to be mainly oriented north-south, more precisely parallel to the quaternary active faults. In the central dark blue zone (Fig. 5, top panel), the orientation of the fast axis was mostly parallel to
<figure><img src="image_7.png" /><figcaption>Polarization change at the end of April</figcaption></figure>

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

**Figure 5.** Map of the $2\psi$ anisotropy regionalization. Red lines, mid-orientation of the fast axis of anisotropy between the orientation before and the orientation after the tectonic change; green dots, seismic stations; black lines, upper edges of the fault proposed by Takada *et al.* (2009). Top panel: the lengths of the lines define the amplitude of the polarization anomaly change at the time of the Iwate-Miyagi earthquake; thin grey lines, quaternary active faults; yellow star, epicentre of the 2008 Iwate-Miyagi earthquake. Bottom panel: the lengths of the lines define the amplitude of the polarization anomaly change at the end of 2008 April; left-hand star, epicentre of the earthquake that occurred on 2008 April 17.

the fault (on the western part of the surface trace of the fault). This indicates that the concentration of the cracks in this zone consists of intense microcracking, probably related to the presence of active volcanoes, with an alignment adjacent to the fault. Near the epicentre and the volcanoes, the pattern of the orientation of the fast axis of anisotropy is more complex. This might be due to the mix of different phenomena that contribute to the closing and opening of microcracks, which have different orientations, such as (1) the microcracking along the fault, (2) the WNW-ESE compres-

sion associated with the earthquake, (3) the ESE-WNW regional compression in the background and (4) the local stress field associated with the loading of volcanoes and controlling the magma pathway. We compared the short-term changes in the HPA at the time of the earthquake (Fig. 5, top panel), the short-term changes in the HPA at the end of 2008 April (Fig. 5, bottom panel), and the long-term changes in the HPA associated with the slow seasonal variations (Fig. 6, bottom panels). We note the similarity of the patterns of the spatial distributions of the HPA changes in these three cases. The dark blue patches where the changes are high (Figs 5 and 6) appear to be correlated with zones where the orientation of anisotropy is coherent at large scales. In these zones, the strong amplitude of anisotropy and its coherent orientation on larger zones might be due to the better alignment of cracks. The latter can be the consequence of higher concentration of microcracks, mainly volcanic fissures, thus higher sensitivity to stress fields due to the presence of highly pressurized volcanic fluids (Pritchard *et al.* 2013; Takada & Fukushima 2013; Brenguier *et al.* 2014), also of the influence of a dominant stress field in these dark blue patches inducing a privilege alignment of the microcracks. This makes these zones more susceptible to changes in the crust. Seasonal changes are mainly located in these sensitive zones but extend more in the area under study (Fig. 6). This might be due to the different phenomena inducing these changes which are not only related to stress field but most probably to the water level. This is why at low topography, where the water level is higher, the area is less sensitive to seismic changes than to seasonal changes.

The fast changes in the HPA that occurred at the end of April, 2008 (Fig. 5, bottom panel), and at the time of the earthquake (Fig. 5, top panel) are both short-term changes that extended over 1 week at most. Their spatial distributions are very similar, with larger amplitudes at the time of the earthquake suggesting a common origin. We also considered other fast changes of HPA. For example, short-term HPA change distribution that is computed for a random peak in 2007 December is presented in the top panel in Fig. 6. It shows a different spatial distribution and much smaller amplitudes than for the HPA changes in April and June. Moreover, the measurements of HPA change at the end of April correspond to the case where $\Delta\psi_{\alpha} \neq 0$ (change in the orientation of anisotropy), unlike seasonal changes that correspond to the case where $\Delta\alpha \neq 0$ and $\Delta\psi_{\alpha} = 0$ (change in the amplitude of anisotropy). This means that the phenomenon that occurred at the end of 2008 April, and that induced fast changes in the HPA is similar to the effects of the earthquake on the HPA. The possible interpretation is that the changes in the stress field in the dark blue patches induced changes in the orientation of the cracks, and hence changes in the crack-induced anisotropy.

The orientation of the fast axis of anisotropy was almost the same throughout the period under study. This is understandable for the area where the crack orientation is mainly ESE-WNW because it is probably caused by the regional stress field associated with the subduction of the Pacific Plate. However, for the area near the fault zone, the pattern of the microcracks was similar even in 2007, almost a year before the earthquake. It is possible that the microcracks were already organized along the zone where the rupture finally occurred. In this case, the pattern of the microcracks would represent a tool to delimit the fault and demonstrate its presence.

# 4 DISCUSSION

The aim of this study is to advance our understanding of seismic zones through an exploration of the entire CCT of seismic
Seasonal polarization change in December 2007
(over 7 days)

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

Total seasonal change in 2007

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

Total seasonal change in 2009

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Map of the 2ψ anisotropy regionalization. Red lines, average orientation of the fast axis of anisotropy, where the lengths of the lines define the amplitude of the polarization anomaly change; green dots, seismic stations; black lines, upper edges of the fault proposed by Takada *et al.* (2009). Top panel: regionalization of the seasonal horizontal polarization anomaly for the HPA changes over 7 d in 2007 December. Bottom panels: regionalization of the total seasonal changes in the HPA in 2007 (left) and 2009 (right).</figcaption></figure>

noise, and especially of the off-diagonal components. Using the CCT components ZT, TZ, RT and TR, we extract an observable, the polarization anomaly change from surface waves that could be related to seismic anisotropy (orientation $\psi_{\alpha}$) before and after earthquakes. In this study, we measure the horizontal polarization change ($\psi_{p}$) in the surface waves in a different tectonic context to that of the strike-slip environment in Parkfield (Durand *et al.* 2011). We monitor the HPA of the surface waves in the area of Iwate-Miyagi, before and after the 2008 earthquake ($M_{w} = 6.9$), where a thrust fault earthquake occurred.

Previously, Hobiger *et al.* (2012, 2016) studied co-seismic and post-seismic shear wave velocity changes caused by the 2008 Iwate-Miyagi earthquake using the combined nine components of CCTs computed between pairs of seismic stations. The co-seismic velocity drops (of maximum 0.63 per cent) were mostly observed

in a frequency range of 0.5–1 Hz, and were concentrated in the southern part of the fault zone. In the lower frequency range of 0.125–0.25 Hz, the maximum co-seismic velocity drop observed was about 0.10 per cent. Hobiger *et al.* (2012) interpreted these velocity drops as a result of near-surface damage underground.

Compared to velocity measurements, the order of magnitude of the polarization anomaly angle (of the order of several degrees and can reach as much as 30) makes it an interesting physical parameter to monitor. Although its interpretation is complex, as this parameter might be induced by heterogeneities in the medium, anisotropy, and/or heterogeneity of the noise-source distribution (non-alignment of the noise incidence and the pair of receivers). The temporal variation of the polarization anomaly angle is also difficult to interpret, as it might be caused by stress-induced anisotropy change (amplitude change or change in the orientation of the fast
axis) or by seasonal changes that could be related to water level, noise incidence, thermoelastic strain induced by atmospheric temperature variations (Ben-Zion & Allam 2013), etc. The separation of these different contributions is not simple but we could use as a criterion the timescale of the change in order to separate the seasonal changes (slower, extend on few months) from the co-seismic changes (faster, extend on few days).

For instance, the HPA changes in the area of Iwate-Miyagi show a superposition of two phenomena. The first aspect of these variations in the HPA occurred on an annual basis and was relatively slow, and is correlated with the seasonal changes, such as the rainfall and snow in the area of study.

Another aspect of these variations in the HPA was the rapid and large variations (maximum of about 30) that occurred at the time of the 2008 Iwate-Miyagi earthquake. Indeed, the data processing for retrieving $\psi_p$ is designed in a general way with no a priori assumptions onto the nature of the medium. And the multitude of parameters involved in the creation of $\psi_p$ makes it difficult to find a unique interpretation of the variation of $\psi_p$. But when the HPA is interpreted we believe that the simplest explanation is the change in the stress field that induces a change in the local seismic anisotropy. In the light of previously reported numerical experiments (Saade et al. 2015), a possible interpretation of such large and fast changes in surface wave polarization is the temporal change in the anisotropy of the medium induced by a stress perturbation that changes the crack distribution (Crampin 1981).

We also observe a new phenomenon: one and a half months before the earthquake, large and fast changes in the HPA also occurred near the fault zone. This phenomenon is not correlated with any seasonal changes, and it was very similar to the changes during the earthquake, in terms of the spatial distribution and the timescale, unlike the seasonal changes. So far we do not have any satisfactory explanation of this phenomenon but we suppose that it could be related to a deep earthquake ($M_w = 5.7$) that occurred on 2008 April 17. This earthquake is the only one with a magnitude larger that 5 that occurred in the area of study, besides the main Iwate-Miyagi earthquake. As for the distribution of anisotropy it is also complex, as it depends on the different concentrations of the cracks in the area, on mineralogical anisotropy, and on the varying phenomena that contribute to the alignment of these cracks. In the area of Iwate-Miyagi, the main effects that seem to be influencing the alignment of the crack distribution appear to be the regional stress fields where the maximum stress axis is ESE-WNW, and the fault-related microcracking. In both cases, the alignment of the crack distribution was mainly horizontal, hence it induces anisotropy with horizontal fast axis. This might explain the deviation of the polarization of the surface waves, which was mainly horizontal, and thus also the small values of $\delta_p$, the vertical polarization anomaly, that has an average of about $\pm 5$ and did not show significant temporal changes.

## 5 CONCLUSIONS

Using passive noise interferometry we extract an observable, the polarization anomaly of surface waves, that we monitor before and after the 2008 Iwate-Miyagi earthquake ($M_w = 6.9$). The polarization anomaly change show a superposition of two phenomena: slow seasonal variations and rapid and large co-seismic variations (that can reach 30). A possible interpretation of the co-seismic change is the stress perturbation that induces a change in the crack distribution therefore a change in the local seismic anisotropy, in

specific zones of the area that are more sensitive than others. In fact, a stress perturbation induces a change in the distribution of fissures which may cause a seismic velocity drop and/or a deviation of seismic anisotropy. But given the order of magnitude of the HPA change, independently of its interpretation, it seems that it is a reliable parameter that can be used to monitor seismogenic zones especially during a seismic cycle. And if we were to go further in such a study, a denser network and smaller interstation distances (of about 5 km) are required for improved determination of the orientation of the fast axis and the relative amplitude of anisotropy.

## ACKNOWLEDGEMENTS

The authors would like to thank the National Research Institute for Earth Science and Disaster Resilience (NIED) for providing the data. IPGP is part of the Labex USPC. ISTerre is part of Labex OSUG@2020. Jean Paul Montagner acknowledges the support of IUF (Institut Universitaire de France, CNRS, UMR 7154). Most of the computations were performed using the S-CAPAD platform of Institut de Physique du Globe de Paris (IPGP).

## REFERENCES

* Ben-Zion, Y. & Allam, A.A., 2013. Seasonal thermoelastic strain and post-seismic effects in Parkfield borehole dilatometers, *Earth planet. Sci. Lett.*, **379**, 120–126.

* Bensen, G.D., Ritzwoller, M.H., Barmin, M.P., Levshin, A.L., Lin, F., Moschetti, M.P., Shapiro, N.M. & Yang, Y., 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements, *Geophys. J. Int.*, **169**, 1239–1260.

* Brenguier, F., Campillo, M., Hadzioannou, C., Shapiro, N.M., Nadeau, R.M. & Larose, E., 2008. Postseismic relaxation along the San Andreas Fault at Parkfield from continuous seismological observations, *Science*, **321**(5895), 1478–1481.

* Brenguier, F., Shapiro, N., Campillo, M., Ferrazzini, V., Duputel, Z., Coutant, O. & Nercessian, A., 2008. Towards forecasting volcanic eruptions using seismic noise, *Nat. Geosci.*, **1**, 126–130.

* Brenguier, F., Campillo, M., Takeda, T., Aoki, Y., Shapiro, N.M., Briand, X., Emoto, K. & Miyake, H., 2014. Mapping pressurized volcanic fluids from induced crustal seismic velocity drops, *Science*, **345**(6192), 80–82.

* Cheng, B., Zhao, D. & Zhang, G., 2001. Seismic tomography and anisotropy in the source area of the 2008 Iwate-Miyagi earthquake (M 7.2), *Phys. Earth planet. Inter.*, **184**(3–4), 172–185.

* Crampin, S., 1981. A review of wave motion in anisotropic and cracked elastic-media, *Wave Motion*, **3**, 343–391.

* Crampin, S., Evans, R., Ucer, B., Doyle, M., Davis, J.P., Yegorkina, G.V. & Miller, A., 1980. Observations of dilatancy-induced polarization anomalies and earthquake prediction, *Nature*, **286**(5776), 874–877.

* Durand, S., Montagner, J.P., Roux, F., Brenguier, F., Nadeau, R.M. & Ricard, Y., 2011. Passive monitoring of anisotropy change associated with the Parkfield 2004 earthquake, *Geophys. Res. Lett.*, **38**(13), doi:10.1029/2011GL047875.

* Gouédard, P. et al., 2008. Cross-correlation of random fields: Mathematical approach and applications, *Geophys. Prospect.*, **56**(3), 375–393.

* Hegde, G.M. & Ye, C., 2008. Singular value decomposition filter: an effective method to enhance the SwissRanger SR-3000- range images, *Int. J. Intell. Control Syst.*, **13**(4), 242–250.

* Hobiger, M., Wegler, U., Shiomi, K. & Nakahara, H., 2012. Coseismic and postseismic elastic wave velocity variations caused by the 2008 Iwate-Miyagi Nairiku earthquake, Japan, *J. geophys. Res.*, **117**, B09313, doi:10.1029/2012JB009402.

* Hobiger, M., Wegler, U., Shiomi, K. & Nakahara, H., 2016. Coseismic and postseismic velocity changes detected by passive image interferometry of
one great and five strong earthquakes in Japan, Geophys. J. Int., 205(2), 1053-1073.

Iidaka, T. & Obara, K., 2013. Shear-wave splitting in a region with newly-activated seismicity after the 2011 Tohoku earthquake, Earth Planets Space, 65, 1059-1064.

Liu, Y., Zhang, H., Thurber, C. & Roecker, S., 2007. Shear wave anisotropy in the crust around the San Andreas fault near Parkfield: spatial and temporal analysis, Geophys. J. Int., 172, 957-970.

Montagner, J.P. & Nataf, H.-C., 1986. A simple method for inverting the azimuthal anisotropy of surface Waves, J. geophys. Res., 91(B1), 511-520.

Nakata, N. & Snieder, R., 2012. Time-lapse change in anisotropy in Japan's near surface after the 2011 Tohoku-Oki earthquake, Geophys. Res. Lett., 39(11), doi:10.1029/2012GL051979.

Pritchard, M.E., Jay, J.A., Aron, F., Henderson, S.T. & Lara, L.E., 2013. Subsidence at southern Andes volcanoes induced by the 2010 Maule, Chile earthquake, Nat. Geosci., 6, 632-636.

Roux, P., 2009. Passive seismic imaging with directive ambient noise: application to surface waves and the San Andreas Fault in Parkfield, CA, Geophys. J. Int., 179(1), 367-373.

Saade, M., Montagner, J.M., Roux, P., Cupillard, P., Durand, S. & Brenguier, F., 2015. Influence of seismic anisotropy on the cross correlation tensor: numerical investigations, Geophys. J. Int., 201, 595-604.

Sens-Schönfelder, C. & Wegler, U., 2006. Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia, Geophys. Res. Lett., 33(21), doi:10.1029/2006GL027797.

Shapiro, N. & Campillo, M., 2004. Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise, Geophys. Res. Lett., 31, L07614, doi:10.1029/2004GL019491.

Smith, M. & Dahlen, F., 1973. The azimuthal dependence of Love and Rayleigh wave propagation in a slightly anisotropic medium, J. geophys. Res., 78(17), 3321-3333.

Stutzmann, E., Arduhuin, F., Schimmel, M., Mangeney, A. & Patau, G., 2012. Modelling long-term seismic noise in various environments, Geophys. J. Int., 191(2), 707-722.

Stehly, L., Campillo, M. & Shapiro, N. M., 2006. A study of the seismic noise from its long-range correlation properties, J. geophys. Res., 111, B10306, doi:10.1029/2005JB004237.

Takada, Y. & Fukushima, Y., 2013. Volcanic subsidence triggered by the 2011 Tohoku earthquake in Japan, Nat. Geosci., 6, 637-641.

Takada, Y., Kobayashi, T., Furuya, M. & Murakami, M., 2009. Coseismic displacement due to the 2008 Iwate-Miyagi Nairiku earthquake detected by ALOS/PALSAR: preliminary results, Earth Planets Space, 61, e9-e12.

Tanimoto, T., 2004. The azimuthal dependence of surface wave polarization in a slightly anisotropic medium, Geophys. J. Int., 156(1), 73-78.

Wapenaar, K., 2004. Retrieving the elastodynamic Green's function of an arbitrary inhomogeneous medium by cross correlation, Phys. Rev. Lett., 93(25), 254301-1-254301-4.

Weaver, R. & Lobkis, O., 2001. Ultrasonics without a source: thermal fluctuation correlations at MHz frequencies, Phys. Rev. Lett., 87(13), 134301-1-134301-4.

Wegler, U. & Sens-Schönfelder, C., 2007. Fault zone monitoring with passive image interferometry, Geophys. J. Int., 168(3), 1029-1033.

# APPENDIX

## A1 Supporting material for the data processing and interpretation of the results

Additional supporting information are presented in this Appendix. Figs A1-A4 are supporting materials for the data processing: Fig. A1: A scheme of the data processing followed in this study.

Fig. A2: Illustration of the angles computed by the optimal rotation algorithm.

Fig. A3: Example of the cross-correlation tensors computed for the entire period of the study, showing the signals on the off-diagonal components ZT, RT, TZ and TR.

Fig. A4: Examples of the temporal changes in the horizontal polarization anomalies and the locations of the processed pairs of receivers.

As for Figs A5 and A6, they are supporting materials for the interpretation part of the results.

Fig. A5 is the polar plot of the azimuth of incidence of the noise in summer and winter relative to the area and time of the study, computed using the method described in Stehly et al. (2006). For each pair of receivers, we measure the normalized amplitude of the Rayleigh waves reconstructed from the cross-correlation tensors, for positive and negative correlation times. The relative amplitude for the positive and negative parts of the cross-correlation tensors can provide information about the incidence of the noise. Indeed, we attribute the normalized amplitudes of the positive part to the azimuth of the pair of receivers, and the normalized amplitudes of the negative part to the azimuth of the pair of receivers plus 180.

Fig. A6 represents the magnitude and timing of all the earthquakes that occurred in the area of this study with magnitudes larger than 4.

And finally, Table A1 represents the reduction of variance as a function of the correlation length used in the regionalization of the polarization anomalies of surface waves. We note that the covariance function of the polarization anomaly between two stations A and B is defined by $C_{\psi}(A, B) = \sigma_{\psi}(A)\sigma_{\psi}(B)e^{-\Delta^2(A,B)/2\lambda^2}$, where $\sigma_{\psi}(A)$ and $\sigma_{\psi}(B)$ are a priori errors on the polarization anomalies at A and B, $\Delta(A, B)$ is the interstation distance and $\lambda$ is the wavelength.

## A2 Theoretical solution of the polarization anomaly

The calculation of sensitivity kernels of polarization anomalies is not simple, contrarily to sensitivity kernels of surface wave phase velocity. The latter can be theoretically formulated in an equation where time is additive and increases with epicentral distance. For the polarization of surface waves, it is acquired as soon as the wave enters the anisotropic medium.

We will consider two simple cases in order to explain how we interpret the polarization anomalies and why we chose sensitivity kernel maxima centred at the two receivers:

(1) a horizontal transverse isotropy (HTI) medium with uniform azimuthal anisotropy, of amplitude $\alpha_0$ and orientation $\Psi_0$,

(2) and a medium with two different anisotropic media with anisotropy amplitudes $\alpha_1$ and $\alpha_2$ and anisotropy orientations $\Psi_{\alpha1}$ and $\Psi_{\alpha2}$ (Fig. A7).

We calculate the three components of displacements in both cases at two stations A and B, located in different anisotropic regions in the second case.

We follow the approach of Tanimoto (2004), at first order of perturbation, for calculating the effect of a general anisotropy for surface waves. This approach is only valid for the Rayleigh wave fundamental mode, which is the case since seismic noise is dominated by shallow sources.

The unperturbed Rayleigh wave displacement is of the form:

$$ \mathbf{u}_R(\mathbf{r}, t) = \begin{pmatrix} V(z) \cos \Psi \\ V(z) \sin \Psi \\ iU(z) \end{pmatrix} e^{i[k(x \cos \Psi + y \sin \Psi) - \omega t]}, \quad (A1) $$
<figure><img src="image_12.png" /><figcaption>**Figure A1.** Scheme of the data processing followed in this study.</figcaption></figure>

where $V(z)$ and $U(z)$ are the scalar depth eigenfunctions for Rayleigh waves, $k$ is the horizontal wavenumber and $\Psi$ is the azimuth of the wavenumber $k$ measured clockwise from the north.

The unperturbed Love wave displacement is of the form:

$$ \mathbf{u}_{\mathrm{L}}(\mathbf{r}, t) = \begin{pmatrix} -W(z) \sin \Psi \\ W(z) \cos \Psi \\ 0 \end{pmatrix} e^{i[k(x \cos \Psi + y \sin \Psi) - \omega t]}, \quad (\mathrm{A}2) $$

where $W(z)$ is the scalar depth eigenfunction for Love waves.

We use the coordinate system $x = \text{south-north}$, $y = \text{west-east}$ and $z = \text{vertical}$.

In the perturbed anisotropic medium, we search for a solution such as the eigenfunctions of anisotropic medium are:

$$ \mathbf{u} = a_{\mathrm{L}}\mathbf{u}_{\mathrm{L}} + a_{\mathrm{R}}\mathbf{u}_{\mathrm{R}}. \quad (\mathrm{A3}) $$
<figure><img src="image_13.png" /><figcaption>Figure A2. Azimuths $\psi_{pA}$ and $\psi_{pB}$, and tilts $\delta_{pA}$ and $\delta_{pB}$, computed by the optimal rotation algorithm.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure A3. Example of the nine components of the cross-correlation tensors for the NRKH-TAJH receiver pair computed between mid-2007 and the end of 2009. The asymmetry of the signal is the result of the non-uniform distribution of the seismic noise. Moreover, the signals on the off-diagonal components of ZT, RT, TZ and TR are large.</figcaption></figure>

The problem then reduces to a set of two equations defined as:

$$ \begin{pmatrix} T_{LL} & T_{LR} \\ T_{RL} & T_{RR} \end{pmatrix} \begin{pmatrix} a_L \\ a_R \end{pmatrix} = \omega^2 \begin{pmatrix} a_L \\ a_R \end{pmatrix}. \quad (A4) $$

The terms $T_{LL}$, $T_{LR}$, $T_{LR}$ and $T_{RR}$ are azimuthally dependent terms (Tanimoto 2004). The coupling between Rayleigh and Love waves comes up in $T_{RL} = T_{LR}$. It is possible to find solutions in terms of quasi-Rayleigh and quasi-Love waves.
<figure><img src="image_15.png" /><figcaption>**Figure A4.** Top panel: example of temporal change in the horizontal polarization anomalies computed by the optimal rotation algorithm, for pairs of receivers for station HNRH. Bottom left panel: distribution of the HPA computed for all pairs of receivers as a function of the azimuth of the pairs. Bottom right panel: locations of the processed pairs of receivers.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure A5.** Polar plot of the azimuth of incidence of the noise in summer and winter relative to the area and time of the study, computed using the method described in Stehly *et al.* (2006). In the grey colour scale of the polar plots, white corresponds to the maximum amplitude of noise and black corresponds to the absence of noise.</figcaption></figure>

The quasi-Rayleigh waves propagate with a phase velocity $V_R$ as given to first order, by Smith & Dahlen (1973) and Montagner & Nataf (1986) from the term $T_{RR}$ and now has three components (where we omit the propagation term):

$$ \mathbf{u}_{qR} = \begin{pmatrix} \cos\Psi V - \frac{T_{RL}}{T_{LL}-T_{RR}}\sin\Psi W \\ \sin\Psi V + \frac{T_{RL}}{T_{LL}-T_{RR}}\cos\Psi W \\ iU \end{pmatrix}. \quad (A5) $$

The phase velocity $V_L$ of Love waves can be found from $T_{LL}$ and the three components of quasi-Love waves are given by:

$$ \mathbf{u}_{qL} = \begin{pmatrix} -\sin\Psi W + \frac{T_{RL}}{T_{LL}-T_{RR}}\cos\Psi V \\ \cos\Psi W + \frac{T_{RL}}{T_{LL}-T_{RR}}\sin\Psi V, \\ i\frac{T_{RL}}{T_{LL}-T_{RR}}U \end{pmatrix}. \quad (A6) $$
<figure><img src="image_17.png" /><figcaption>**Figure A6.** Top panel: first derivative of the average of the horizontal polarization anomaly (from mid-2007 to the end of 2009). Bottom panel: representation of the magnitude and timing of earthquakes in the study area with magnitude >4. First red line: earthquake with $M_w = 5.7$ at a depth of 160 km and second red line: main Iwate-Miyagi earthquake.</figcaption></figure>

<table><caption>**Table A1.** Variance reduction relative to the initial model obtained for different correlation lengths (km) used in the regionalization of the polarization anomalies of surface waves. The chosen correlation length is 8 km.</caption><thead><tr><th>Correlation length (km)</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th><th>13</th><th>14</th><th>15</th></tr></thead><tbody><tr><th>Variance reduction</th><td>45.8</td><td>68.5</td><td>81.4</td><td>87.0</td><td>88.2</td><td>89.4</td><td>89.7</td><td>89.8</td><td>88.9</td><td>88.4</td><td>87.7</td><td>87.0</td><td>86.1</td><td>85.0</td><td>83.7</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Figure A7.** Map view of the two anisotropic regions around stations A and B. $\Psi_{P,A}$ is the polarization anomaly of the wave arriving at station A and $\Psi_{P,AB}$ is the final polarization anomaly of the wave propagating from A to B.</figcaption></figure>

We assume that Rayleigh and Love waves have sufficiently different phase velocities ($T_{LL} \neq T_{RR}$) to limit the coupling between Rayleigh and Love waves at the transition between two anisotropic media. This means that Rayleigh waves will be transformed solely into quasi-Rayleigh waves at this transition, neglecting in first approximation the energy component that is modified into Love waves. The same reasoning applies to Love waves.

## A2.1 One anisotropic region

Without loss of generality, we consider the simplest case of an incoming wave along AB such that $\cos \Psi = 0$ and $\sin \Psi = 1$. The quasi-Rayleigh wave components are given by:

$$ \mathbf{u}_{qR} = \begin{pmatrix} -\frac{T_{RL}}{T_{LL}-T_{RR}} & W \\ \frac{V}{iU} & 0 \end{pmatrix} e^{i\omega(\Delta/V_{R}-t)}, \quad (A7) $$

and the quasi-Love wave components are given by:

$$ \mathbf{u}_{qL} = \begin{pmatrix} -W \\ \frac{T_{RL}}{T_{LL}-T_{RR}}V \\ i\frac{T_{RL}}{T_{LL}-T_{RR}}U \end{pmatrix} e^{i\omega(\Delta/V_{L}-t)}. \quad (A8) $$

Let us define the coupling parameter $\tau = \frac{T_{RL}}{T_{LL}-T_{RR}}$ and $\tau \ll 1$.

We define by $\Psi_P$ the polarization angle of the quasi-Rayleigh waves by using the horizontal components given in eq. (A7):

$$ \sin \Psi_P = \frac{-\tau W}{\sqrt{V^2 + (-\tau W)^2}} \approx \left(-\tau \frac{W}{V}\right) - \frac{1}{2} \left(-\tau \frac{W}{V}\right)^2. \quad (A9) $$

To first order, by using the first-order expansion of eq. (A9), we get:

$$ \Psi_P \approx \alpha_1 \cos(2\Psi) + \alpha_2 \sin(2\Psi), \quad (A10) $$

where we neglect the $4 - \Psi$ terms. To first order, the component along the propagation vector $k$ (x-axis) is not perturbed and is still equal to eigenfunction $V$.

If we consider a wave originating from the left, at both stations, the wave will display a polarization anomaly $\Psi_P$. For a wave originating from the right, the same polarization anomaly will be observed at both stations.

## A2.2 Two anisotropic regions

Let us consider a wave originating from station A, located in the first anisotropic region, where the anisotropy direction is $\Psi_{\alpha1}$, which will induce a polarization anomaly $\Psi_{P1}$.

Somewhere before reaching station B, the wave enters a second anisotropic region with a different anisotropy direction $\Psi_{\alpha2}$.

Let us consider the propagation of the quasi-Rayleigh wave into the second anisotropic region and consider separately its different horizontal components from eq. (A7).

The radial component of the quasi-Rayleigh arriving from the region 1, will be perturbed again in the region 2 as follows:

$$ (0, V) e^{i\omega\Delta_1/V_{R1}} \rightarrow (-\tau_2 W, V) e^{i\omega(\Delta_1/V_{R1} + \Delta_2/V_{R2})}, \quad (A11) $$
where $V_{R1}$ is the Rayleigh wave velocity in region 1 and $V_{R2}$ is the Rayleigh wave velocity in region 2. The terms $\Delta_1$ and $\tau_1$ are related to the region 1, and $\Delta_2$ and $\tau_2$ are related to the region 2.

The transverse component of the quasi-Rayleigh becomes then:

$$(-\tau_1 W, 0)e^{i\omega(\Delta_1/V_{R1})} \rightarrow (-\tau_1 W, \tau_1 \tau_2 V)e^{i\omega(\Delta_1/V_{R1} + \Delta_2/V_{R2})}. \quad (A12)$$

The final equation for the displacement of the quasi-Rayleigh wave arriving at the second station B is then:

$$\mathbf{u}_{qR,AB} = \begin{pmatrix} -(\tau_2 + \tau_1)W \\ (\tau_2 \tau_1 + 1)V \\ iU \end{pmatrix} e^{i\omega(\Delta_1/V_{R1} + \Delta_2/V_{R2})}. \quad (A13)$$

The final horizontal polarization anomaly $\Psi_{P,AB}$ is then:

$$\begin{align}
\Psi_{P,AB} &\approx \sin\Psi_{P,AB} = \frac{(\tau_2 + \tau_1)W}{\sqrt{(\tau_2 + \tau_1)^2 W^2 + (\tau_2 \tau_1 + 1)^2 V^2}} \nonumber \\
&\approx (\tau_2 + \tau_1) \frac{W}{V}. \tag{A14}
\end{align}$$

The polarization anomaly $\Psi_{P,AB}$ of a wave propagating through two anisotropic regions characterized by the location of the two stations A and B, seems to be influenced by $\tau_1$ and $\tau_2$, the coupling parameters of the two regions. It means that the observation of the polarization anomaly of a wave propagating from station A to B is primarily sensitive to the anisotropy distribution around the stations. Furthermore, the result of $\Psi_{P,AB}$ shows that the reciprocity is valid.

In fact, the average distance between pairs of receivers is around 20 km. With a correlation length about 8 km, we implicitly assume that there might be two anisotropic regions with different orientations and amplitudes of anisotropy. And since the polarization anomaly is independent of the extension of the anisotropic region, we give it the maximum sensitivity at the receiver.

Our model is the simplest model that can explain data. The distribution of anisotropy might be more complicated but it is difficult to find additional complexity. In addition, since the variance reduction is really large (~90 per cent) (Table A1), we can a posteriori consider that it is a reasonable assumption.