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

HAL
open science

# Detailed space-time variations of the seismic response of the shallow crust to small earthquakes from analysis of dense array data

Luis Fabian Bonilla, Yehuda Ben-Zion

► To cite this version:

Luis Fabian Bonilla, Yehuda Ben-Zion. Detailed space-time variations of the seismic response of the shallow crust to small earthquakes from analysis of dense array data. Geophysical Journal International, 2020, 225, pp.298 - 310. 10.1093/gji/ggaa544. hal-04402182

HAL Id: hal-04402182

https://univ-eiffel.hal.science/hal-04402182v1

Submitted on 18 Jan 2024

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.

Public Domain
<figure><img src="image_2.png" /></figure>

# Detailed space-time variations of the seismic response of the shallow crust to small earthquakes from analysis of dense array data

Luis Fabian Bonilla<fnref n="1" /> and Yehuda Ben-Zion<fnref n="2" />

Accepted 2020 November 11. Received 2020 November 9; in original form 2020 May 24

## SUMMARY

We compute high-resolution space-time variations of subsurface seismic properties from autocorrelation functions (ACF's) of noise and local earthquakes, recorded by the Sage Brush Flat dense array deployed around the Clark branch of the San Jacinto fault. The resolved temporal changes are referred to as apparent velocity changes because they reflect both non-linear response and variations of material properties such as cracking and damage. Apparent velocity changes are estimated at four frequency bands (10–15, 10–20, 15–30 and 20–40 Hz) for two local earthquake data sets. In one analysis, ACF's from *P*- and *S*-wave windows of 31 small events with magnitudes below 3.1 are used to compute the apparent velocity variations with respect to the mean ACF of each phase, and we also use the mean ACF of noise data as reference to estimate the changes. In a further analysis, the temporal evolution of properties is computed using moving time windows in continuous waveform over one-hour long data with noise and earthquake signals. The apparent velocity changes and recovery times are frequency dependent and present a strong spatial variability across the array. The resolved changes are larger and recovery time shorter with data associated with higher frequencies. At frequencies larger than 15 Hz, and using the mean ACF of noise data as a reference, the apparent average velocity changes across the array during the passage of the *P* and *S* waves from the small local events are 2.5 per cent and 6 per cent, respectively. The apparent velocity changes decrease by one order of magnitude when the earthquake data are used as a reference. The relatively large changes in response to very low ground motion have important implications on nonlinear processes involving degradation and healing of the subsurface material during common earthquake shaking.

**Key words:** Earthquake ground motions; Seismic interferometry; Seismic noise; Site effects.

## 1 INTRODUCTION

The properties and dynamics of material in the top few hundred metres of the crust are not well understood. Ongoing weathering processes and ground motion generated by earthquakes produce distributed cracking and other forms of rock damage that reduce the seismic velocities of the top material to values that can be 200–400 m s⁻<fnref n="1" /> or lower (e.g. Theodulidis et al. 1996; Bonilla et al. 2011; Zigone et al. 2019; Qin et al. 2020). The pervasive pre-existing damage, low seismic velocities that amplify ground motion, and low normal stress make the material in the top crust highly susceptible to additional failures. Damaged rocks with cracks and other defects also exhibit nonlinear elasticity manifested by shifts of the resonance frequency with increasing loading amplitude and related phenomena (e.g. Lyakhovsky et al. 2009; Remillieux et al. 2017). Nonlinear response to loadings and failure processes have been

documented in various laboratory experiments (e.g. Lockner et al. 1977; Johnson & Jia 2005; Stanchits et al. 2006; Pasqualini et al. 2007; Hamiel et al. 2009; Aben et al. 2019). *In situ* analyses of nonlinear rock behaviour and temporal changes of seismic properties are essential for improved understanding of wave propagation in the top crust, estimating seismic shaking hazard and other topics.

Various techniques have been developed to estimate temporal changes of seismic properties induced at the subsurface by earthquakes and other loadings. Since material failure is followed by rapid partial recovery of properties (e.g. Dieterich & Kilgore 1996; Nakatani & Scholz 2004; Pei et al. 2019; Wu et al. 2009), the results depend strongly on the analysis time step. Analyses based on the ambient seismic noise are wide spread since ambient noise data are available in essentially all regions and time intervals (e.g. Sens-Schönfelder & Wegler 2006; Brenguier et al. 2008; Froment et al. 2013; Brenguier et al. 2014; Hobiger et al. 2016; Brenguier

<footnote n="1">¹Geotechnical Engineering, Environment, Natural hazards and Earth sciences Department, Université Gustave Eiffel, 77447 Marne-la-Vallée Cedex 2, France. E-mail: luis-fabian.bonilla-hidalgo@univ-eiffel.fr</footnote>

<footnote n="2">²Southern California Earthquake Center, University of Southern California, Los Angeles CA 90089-0740, USA</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** A dense array with 1108 vertical-component geophones at the SGB site on the SJFZ. Red dots show the row and column intervals every 10 sensors. Black lines indicate fault surface traces including the MCF from Share *et al.* (2020). The row and column numbers of the dense array start from the SW corner and increase toward the NW and NE (blue arrows), respectively. The colour scale indicates the topography. Stations R0101 and R7010 are depicted by the white circles.</figcaption></figure>

et al. 2019). However, noise sources are relatively weak, so noise-based analysis requires stacking of relatively quiet periods, leading to analysis time steps of days or more. Waveforms generated by repeating earthquakes have also been used to derive temporal changes of seismic velocities (e.g., Poupinet et al. 1984; Schaff & Beroza 2004; Peng & Ben-Zion 2006), and also involve time steps that are typically days or more. Using spectral ratios of seismograms recorded by different stations (e.g. Rubinstein & Beroza 2004; Sawazaki et al. 2006; Wu et al. 2009, 2010; Nakata & Snieder 2012), cross-correlations of earthquake waveforms (Roux & Ben-Zion 2014), and autocorrelation in moving time windows of seismograms (Bonilla et al. 2019; Qin et al. 2020) can achieve time resolution of minutes to seconds and resolve finer details of temporal changes.

The Sage Brush Flat (SGB) site in the trifurcation area of the San Jacinto fault zone (SJFZ) southeast of Anza, California, had a dense array of 1108 vertical ZLand nodes recording continuous waveforms for ~30 d in spring of 2014 (Ben-Zion et al. 2015). The dense array covered ~600 m × 600 m region around the Clark branch of the SJFZ with a core grid consisting of 20 rows perpendicular to and centred on the fault trace (Fig. 1). Each row had 50 sensors at a nominal 10 m interstation spacing and the nominal separation between rows was 30 m. The remaining 108 sensors were deployed as extensions to multiple rows. The array recorded earthquake and noise data continuously at 500 samples s⁻¹. The site has a fault scarp producing elevation change of 70 m to the NE and local sedimentary basin to the SW. Several imaging studies provided high-resolution results on the velocity structure at the site (e.g. Hillers et al. 2016; Qin et al. 2018; Mordret et al. 2019; Share et al. 2020). The dense array data provide unique opportunities to improve significantly the *in situ* constraints on nonlinear effects and temporal changes of properties in the top crust, in response to various loadings and in relation to different shallow seismic velocities and other site conditions.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Velocity time-series recorded at station R0101; the red triangles indicate the onset of earthquake data, the largest one being at ~1800 s. (b) Spectrogram of the unfiltered signal. Horizontal white dashed lines represent the upper and lower bounds of usable frequency band.</figcaption></figure>

In the present study, we analyse changes induced by noise and earthquake data recorded at the SGB dense array. The resolved changes account for both nonlinear response and variations of seismic properties due to cracking and material damage, so are referred to as apparent velocity changes. Section 2 describes the data and methods used in the analyses. In Section 3, we first investigate apparent velocity changes at the subsurface by comparing delay times of mean autocorrelation functions (ACF's) from noise and earthquake data at different frequency bands. The resolved apparent velocity changes with respect to ambient noise at frequencies larger than 15 Hz are found to be around 2.5 per cent and 6 per cent for the *P*- and *S*-wave windows of earthquake records, respectively. These values are one order of magnitude smaller when earthquake data are used as a reference. Next, we evaluate apparent velocity changes in continuous noise recordings that contain some earthquake data. The results show similar apparent velocity changes, with respect to the background noise, along with some healing effects of the subsurface material. The derived values are frequency dependent and show larger changes for higher frequencies reflecting shallower material. The results are discussed and summarized in the final Section 4.

## 2 DATA AND METHODS

The rectangular array shown in Fig. 1 consists of 1108 vertical 10 Hz ZLand geophones and it recorded in a continuous mode with a sampling rate of 500 Hz between 2014 May 7 and June 13 (Ben-Zion et al. 2015). The high-resolution techniques developed and used in this work are illustrated using one hour of data recorded by the array beginning at 01:41 GMT on 2014 May 27. Fig. 2 shows the velocity time history at station R0101 (marked in Fig. 1) and corresponding spectrogram computed with the short-time Fourier transform. The waveform during this period of time includes four earthquakes that are marked by the red triangles in Fig. 2(a) (the vertical scale is augmented to show the small events). The signal is bandpass filtered between 2 and 50 Hz with an acausal Butterworth filter of four poles. Fig. 2(b) displays the corresponding spectrogram of the unfiltered time-series. The four events are easily identified by narrow time windows having energy between 5 and 40 Hz. Additional events with frequencies higher than 15 Hz are also observed
<table><caption>**Table 1.** Events used in this study. PGA and PGV are the mean of the network.</caption><thead><tr><th>Event no.</th><th>Date</th><th>Latitude</th><th>Longitude</th><th>Depth (km)</th><th>PGA (m s<sup>-2</sup>)</th><th>PGV (m s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1</td><td>09/05/2014 02:22</td><td>34.0080</td><td>-117.0223</td><td>21.4</td><td>4.09e-04</td><td>4.69e-06</td></tr><tr><td>2</td><td>09/05/2014 04:49</td><td>33.9272</td><td>-116.6928</td><td>25.0</td><td>1.35e-02</td><td>1.61e-04</td></tr><tr><td>3</td><td>09/05/2014 14:41</td><td>33.9272</td><td>-116.6954</td><td>23.6</td><td>8.13e-04</td><td>7.67e-06</td></tr><tr><td>4</td><td>09/05/2014 14:42</td><td>33.9261</td><td>-116.6992</td><td>23.9</td><td>3.97e-03</td><td>4.10e-05</td></tr><tr><td>5</td><td>09/05/2014 22:57</td><td>33.5527</td><td>-116.7558</td><td>7.2</td><td>1.20e-03</td><td>1.40e-05</td></tr><tr><td>6</td><td>10/05/2014 00:12</td><td>33.5521</td><td>-116.7557</td><td>7.4</td><td>1.16e-03</td><td>1.22e-05</td></tr><tr><td>7</td><td>10/05/2014 18:12</td><td>33.5508</td><td>-116.7562</td><td>7.3</td><td>4.64e-04</td><td>6.28e-06</td></tr><tr><td>8</td><td>13/05/2014 07:27</td><td>33.4891</td><td>-116.4541</td><td>5.1</td><td>1.10e-03</td><td>1.19e-05</td></tr><tr><td>9</td><td>15/05/2014 09:23</td><td>33.5137</td><td>-116.4521</td><td>8.1</td><td>1.67e-03</td><td>1.22e-05</td></tr><tr><td>10</td><td>15/05/2014 13:14</td><td>33.4715</td><td>-116.4506</td><td>7.9</td><td>7.93e-04</td><td>8.01e-06</td></tr><tr><td>11</td><td>16/05/2014 09:42</td><td>33.6228</td><td>-116.6830</td><td>15.7</td><td>3.27e-04</td><td>4.84e-06</td></tr><tr><td>12</td><td>16/05/2014 09:44</td><td>33.6223</td><td>-116.6834</td><td>15.5</td><td>7.39e-04</td><td>6.13e-06</td></tr><tr><td>13</td><td>16/05/2014 18:26</td><td>33.5035</td><td>-116.4776</td><td>5.9</td><td>9.57e-04</td><td>9.09e-06</td></tr><tr><td>14</td><td>17/05/2014 06:23</td><td>33.4703</td><td>-116.4560</td><td>7.6</td><td>1.78e-03</td><td>1.62e-05</td></tr><tr><td>15</td><td>17/05/2014 15:19</td><td>33.3540</td><td>-116.7422</td><td>15.1</td><td>3.49e-03</td><td>3.18e-05</td></tr><tr><td>16</td><td>20/05/2014 10:21</td><td>33.4993</td><td>-116.4781</td><td>6.3</td><td>1.85e-03</td><td>1.96e-05</td></tr><tr><td>17</td><td>20/05/2014 12:46</td><td>33.9222</td><td>-116.7040</td><td>22.0</td><td>5.53e-04</td><td>5.70e-06</td></tr><tr><td>18</td><td>22/05/2014 15:12</td><td>33.4157</td><td>-116.3515</td><td>13.1</td><td>4.23e-04</td><td>6.98e-06</td></tr><tr><td>19</td><td>22/05/2014 15:13</td><td>33.4141</td><td>-116.3548</td><td>14.2</td><td>4.24e-04</td><td>7.01e-06</td></tr><tr><td>20</td><td>23/05/2014 22:15</td><td>33.5135</td><td>-116.4541</td><td>8.2</td><td>5.35e-03</td><td>5.27e-05</td></tr><tr><td>21</td><td>25/05/2014 05:00</td><td>33.4756</td><td>-116.4962</td><td>11.8</td><td>1.45e-03</td><td>1.44e-05</td></tr><tr><td>22</td><td>25/05/2014 06:27</td><td>33.5131</td><td>-116.4541</td><td>8.0</td><td>9.94e-04</td><td>1.07e-05</td></tr><tr><td>23</td><td>27/05/2014 02:10</td><td>33.5856</td><td>-116.6183</td><td>14.5</td><td>1.48e-02</td><td>1.44e-04</td></tr><tr><td>24</td><td>27/05/2014 10:36</td><td>33.3840</td><td>-116.3837</td><td>7.4</td><td>3.22e-04</td><td>4.41e-06</td></tr><tr><td>25</td><td>28/05/2014 09:04</td><td>33.5123</td><td>-116.4548</td><td>8.5</td><td>1.99e-03</td><td>1.55e-05</td></tr><tr><td>26</td><td>28/05/2014 11:25</td><td>33.5535</td><td>-116.4082</td><td>11.1</td><td>7.90e-04</td><td>6.85e-06</td></tr><tr><td>27</td><td>01/06/2014 17:00</td><td>33.5020</td><td>-116.4557</td><td>7.0</td><td>2.50e-03</td><td>2.34e-05</td></tr><tr><td>28</td><td>02/06/2014 01:35</td><td>33.3838</td><td>-116.3840</td><td>7.3</td><td>7.76e-04</td><td>1.09e-05</td></tr><tr><td>29</td><td>03/06/2014 04:55</td><td>33.5048</td><td>-116.4769</td><td>6.2</td><td>5.85e-04</td><td>6.14e-06</td></tr><tr><td>30</td><td>03/06/2014 23:07</td><td>33.5035</td><td>-116.4773</td><td>6.1</td><td>1.11e-03</td><td>1.47e-05</td></tr><tr><td>31</td><td>04/06/2014 14:45</td><td>33.7303</td><td>-116.7206</td><td>17.1</td><td>6.16e-03</td><td>9.04e-05</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) Velocity time history filtered between 2 and 50 Hz for event 1 (Table 1) recorded at station R0101 with time windows of *P* and *S* waves marked. (b) Spectrogram of unfiltered signal with clear energy between 5 and 50 Hz for both the *P*- and *S*-wave windows.</figcaption></figure>

and correspond to other sources of ground motion such air-traffic and wind shaking obstacles above the ground (Ben-Zion *et al.* 2015; Meng & Ben-Zion 2018a,b; Johnson *et al.* 2019).

We select for analysis 31 events recorded by the array that have a signal-to-noise ratio (S/N) greater than three (Table 1). Fig. 3 shows an example of earthquake 1 in Table 1 recorded at station R0101. The time windows of the *P* and *S* waveforms used in the ACF

computation are 5 s long (Fig. 3a). The *P*- and *S*-wave windows are centred at [$t_P - 3, t_P + 2$] and [$t_S - 2, t_S + 3$], respectively, where $t_P$ and $t_S$ are the onsets of the *P*- and *S*-wave arrivals picked by the automated algorithm of Ross & Ben-Zion (2014) and manually corrected. Fig. 3(b) shows the spectrogram of the unfiltered signal with a clear energy between 5 and 50 Hz. In the following, the ACF are computed between 2 and 50 Hz for both noise and earthquake data.

In the next sections, we explain the signal processing methods to compute ACF's and the extraction of the time-frequency parameters needed to compute velocity changes.

## 2.1 Autocorrelation functions

Following Claerbout (1968), a reflection seismogram at a site can be obtained from the autocorrelation of the recorded earthquake waveforms. Wapenaar (2003) generalized this conjecture to three dimensions. Most studies using ACFs look for coherent phases at long delay times to map the subsurface geological structures (e.g. Ruigrok & Wapenaar 2012; Romero & Schimmel 2018). In this study, we are interested in velocity changes that occur at very shallow depths (e.g. the top 100 m). Bonilla *et al.* (2019) used the first positive peak of the ACF to study velocity changes at a site in Japan. However, it is not clear which peak should be used and how it is related to the location of a shallow interface at depth. To address these issues, we propagate a delta-like function having a peak ground acceleration (PGA) of $3 \text{ m s}^{-2}$, approximating an *S* wave, over a 1-D model of a soil column of 300 m depth. The model has five layers including the underlying bedrock (Table 2).
<table><caption>**Table 2.** Soil profile used in the numerical simulation. Layers 1–4 have a density of 1700 kg m⁻³, and Qs = 20. The bedrock (layer 5) has a density of 2000 kg m⁻³. For the nonlinear computation, layer 1 is assumed to be a cohesionless soil having a friction angle of 30°.</caption><thead><tr><td>Layer no.</td><td>Depth (m)</td><td>Vs (m s<sup>-1</sup>)</td></tr></thead><tbody><tr><td>1</td><td>0</td><td>200</td></tr><tr><td>2</td><td>10</td><td>400</td></tr><tr><td>3</td><td>50</td><td>600</td></tr><tr><td>4</td><td>100</td><td>800</td></tr><tr><td>5</td><td>300</td><td>1500</td></tr></tbody></table>

All layers have viscoelastic response and the top layer can also have nonlinear soil behaviour. The wave propagation calculation uses the code NOAH (Bonilla *et al.* 2005), which solves the SH equation of motion for a 1-D linear and nonlinear site response with a finite differences scheme. The nonlinear behaviour in layer 1 is approximated by a multishear mechanism using a hyperbolic soil model as described in Iai *et al.* (1990). The interested reader may find details of the nonlinear constitutive model in Iai *et al.* (1990) and some applications to modelling observed data in Bonilla *et al.* (2005). The maximum frequency resolution is 25 Hz for 40 points per wavelength, the grid space is 0.2 m, and pore pressure excess is not allowed.

Fig. 4(a) shows the resulting acceleration time history at the surface for the viscoelastic computation. Fig. 4(b) displays the phase correlation between solutions at each node at depth and the surface as a reference (Schimmel 1999). This technique has the advantage of using only the phase information of the analytical signal, so it is not necessary to perform any bit normalization nor pre-whitening of the data. Ventosa *et al.* (2019) showed that the method of Schimmel (1999) can be computed in the frequency domain, making its calculation very fast and inexpensive. The waveforms are filtered between 1 and 10 Hz before computing the correlation, and the results clearly show the up- and downgoing wavefields (e.g. Chandra *et al.* 2015). The incoming waves take around 0.49 s to travel the 300 m distance and 0.98 s to go back to depth. Fig. 4(c) shows the ACFs of the waveform at each node, which are also computed using the phase correlation method (Schimmel 1999). Fig. 4(d) plots both the ACF of the surface signal and the normalized acceleration from Fig. 4(a) for the S-wave shifted to the origin.

Note that the ACF's intercept the media interfaces at 0.1, 0.3, 0.48 and 0.98 s respectively (Figs 4c and d). These times (vertical dashed lines in Fig. 4d) indicate the two-way traveltime between the surface and reflected wavefield at the interfaces. This simple numerical experiment illustrates Claerbout's concept that the ACF provides information on the velocity structure underneath the station. Furthermore, since the input is a delta-like function, the resulting ground motion at the surface is close to the Green's function of the media. Indeed, both the acceleration at the surface (shifted in time to have the same origin) and its corresponding ACF are similar, confirming that the autocorrelation contains the Green's function of the system (Wapenaar 2003).

To examine the effects of nonlinear material behaviour on the computed ground motion, and in particular its ACF, the input motion of 3 m s⁻² is sufficiently strong to trigger nonlinear response at the first layer. Fig. 5(a) shows computed accelerograms at the surface for the linear and nonlinear cases (Table 2). Note the lower amplitude and slight delay of the nonlinear result compared to the linear case. Fig. 5(b) displays the spectral ratio between the surface and the input motion. As expected, the damping produced by hysteresis reduces

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Wave propagation of a delta-like function simulating an S-wave travelling in a 1-D soil column (Table 2). (a) Acceleration at the surface (GL-0m). (b) Phase correlation between the waveforms at each node and the one at the surface. The horizontal white dashed lines represent the layer interfaces. The acausal and causal parts represent the up- and downgoing wavefields, respectively. (c) ACFs for each time history at each node. (d) Autocorrelation of the surface accelerogram and the S-wave train of (a) shifted to the origin. The vertical dashed lines mark the two-way traveltimes between the surface and interfaces.</figcaption></figure>

the amplification above 4 Hz for the nonlinear case, producing a de-amplification of the ground motion. The ACF's for the two cases (Fig. 5c) also show that the nonlinear behaviour delays the autocorrelation, which may be interpreted as a velocity reduction. It is also seen that the autocorrelation for the nonlinear case does not have a clear correspondence with the two-way traveltimes shown in Fig. 4(d), indicating that even if the nonlinear behaviour is confined to the first 10 m (first layer), the total propagation can affect longer delay times. This suggests that picking the first negative peak may be best to monitor velocity changes in the top crust. If we assume that the delay time accumulates uniformly in the top layer, the velocity change can be computed as

$$ \frac{\Delta t}{t} = -\frac{\Delta v}{v} \qquad (1) $$

where Δt = t - t₀ with t₀ being the reference delay time (i.e. the linear response for the example described above). From Fig. 5(c), t = 0.127 s and t₀ = 0.105 s, leading to an estimate of 17 per cent velocity reduction.

During the computation of the nonlinear case, the stress and strain time histories at each node of the nonlinear layer are obtained. The instantaneous shear modulus μ(t) can be calculated as the slope of the stress–strain phase space (not shown here). Since we know the initial shear modulus at t = 0, we can obtain the normalized variation of shear modulus as μ/μ₀, where the time dependency is
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Acceleration time histories recorded at the surface for the linear and nonlinear wave propagation within the soil column of Table 2. (b) Spectral ratios between the surface and input motion for the linear and nonlinear results. (c) ACFs for the linear and nonlinear accelerograms at the surface. The delay times for the first minimum are also plotted and their values indicated in the legend. (d) Velocity changes at each node obtained from computed stress–strain time histories in gray lines (see the text for explanation). The thick red line represents the average of all nodes in first layer. The horizontal blue line marks the velocity changes using the values in Fig. 5(c) and eq. (1).</figcaption></figure>

dropped to simplify the notation. Since the density is the same, the square root of the normalized shear modulus represents the shear wave speed variability $v/v_0$, which can easily be converted to the form of eq. (1) as

$$ \frac{\Delta v}{v_0} = \frac{v - v_0}{v_0} = \frac{v}{v_0} - 1 \qquad (2) $$

Fig. 5(d) shows the velocity changes from eq. (2) at each node of the first layer (grey lines), along with the mean of all velocity changes (thick red line). Several aspects of the procedure and results should be noted. First, evaluating the instantaneous slope is not trivial and can produce noisy results. We use a Savitzky–Golay smoothing filter (Press et al. 1992) to compute a smoothed first derivative using five points and a second-order polynomial approximation. The resulting slope is low-passed filtered below 10 Hz with a non-causal Butterworth filter of four poles and then normalized to the initial shear modulus. The velocity changes show low values near the surface up to almost 30 per cent changes at the bottom of the layer. The largest reductions are near the strong ground motion pulses between 1.5 and 1.8 s in Fig. 5(a). Furthermore, the average velocity change is around 20 per cent, which is close to what we obtain from the ACF’s in Fig. 5(c). This suggests that the ACF’s capture the average behaviour of the layer. However, we cannot assume that the average will always represent the nonlinear behaviour of the whole layer. Indeed, the results will depend on the frequency band used in the analysis. Low frequencies will map averages associated with deeper velocity changes and high frequencies will provide results for the shallow structures. The constitutive model used in this computation does not have recovery effects, so the velocity changes remain negative after the passage of the seismic waves. However, there is an apparent healing due to the averaging of velocity changes at depth as shown by the thick red line. This simple numerical model provides useful insights of nonlinear site response and related velocity changes.

In the subsequent data analyses, we use delay times from ACF’s in moving time windows to estimate velocity changes. Since the

synthetic calculations done in this section do not model more complex mechanisms describing changes of material properties due to damage or cracking, the computed results from observed data are referred to as apparent velocity changes in the rest of the text.

## 2.2 Predominant Frequency

Fig. 5(b) implies that nonlinear material behaviour in the top layer also reduces the high-frequency content of the ground motion. This is consistent with observations associated with aftershocks of the 1994 Northridge earthquake in California (Field et al. 1997). Similar results were observed with spectral ratios computed in vertical arrays during weak and strong motion (Beresnev & Wen 1996; Bonilla et al. 2011; Régnier et al. 2013). In a soil layer over a half-space, the resonance frequency is proportional to the shear wave velocity, so velocity reduction during strong motion will be mapped into a frequency decrease. However, the results of Fig. 5 and eq. (1) indicate that an onset of nonlinear behaviour by ground motion larger than some threshold can be manifested as an apparent velocity reduction.

Following Cohen (1989), the predominant frequency $f_p$ and the spectral bandwidth $f_b$ are computed as

$$ f_p^2 = \frac{\int_0^\infty f^2 P(f) df}{\int_0^\infty P(f) df} \qquad (3a) $$

$$ f_b^2 = \frac{\int_0^\infty (f - f_c)^2 P(f) df}{\int_0^\infty P(f) df} \qquad (3b) $$

Here, $f_c$ is the central frequency and $P(f)$ is the power spectrum

$$ f_c = \frac{\int_0^\infty f P(f) df}{\int_0^\infty P(f) df} \qquad (4) $$

We could approximate the apparent velocity change using the predominant frequency as

$$ \frac{\Delta v}{v_0} = \frac{\Delta f}{f_{p0}} = \frac{f_p - f_{p0}}{f_{p0}} \qquad (5) $$
where $f_{p0}$ is the reference predominant frequency. In the case of Fig. 5(c), the corresponding predominant frequencies, between 1 and 10 Hz, of the ACF's for the linear and nonlinear cases are 4.42 and 3.82 Hz, respectively. This implies an apparent velocity change of 14 per cent, which is also close to what it has been computed before. However, using the frequency properties of ground motion exclusively may mix effects of structures located at different depths. In other words, if the frequency band is too large, the frequency changes are influenced not only by the most nonlinear layer, but also by the response of and interaction with other structures. This may produce inaccurate velocity changes, so we mainly use this information to confirm the computed velocity variations because they occur together. If there is no change of the predominant frequency and time delay, which should be inversely proportional, the estimated velocity change is not stable.

To avoid the problems discussed above, we analyse the data using narrow-band signals. In addition, studies of narrow-band audio signals have shown that the predominant period can be obtained from the time of the first positive peak of the ACF (e.g. Valero & Alias 2013). Narrow-band ACF's are computed by filtering the ACF to selected frequency intervals. In this study, we compute ACF between 2 and 50 Hz, where most of the energy of earthquake waveforms is. Since the sensors have a natural frequency of 10 Hz, the chosen frequency bands are between 10–15, 10–20, 15–30 and 20–40 Hz. Eqs (3) and (4) are applied to the Fourier spectra of the computed ACF only for frequencies within the chosen intervals.

# 3 RESULTS

## 3.1 Analysis of noise and earthquake data

As shown in previous section, apparent velocity changes can be computed from noise and earthquake records by tracking delay times in ACF's. We characterize the reference time $t_0$ using the mean ACF computed during one hour of noise data. Significant changes in media properties are assumed to be produced only by earthquake activity. Time delays generated by earthquakes are computed from mean ACF's for P- and S-wave time windows. Velocity changes can be computed also for each event, but since these events are quite small (see PGA and peak ground velocity, PGV, values in Table 1), we stack their corresponding ACF's to have more robust average estimations.

The computation of the mean ACF for the noise data is done with moving time windows of 10 s and 50 per cent overlap using the phase correlation technique (Schimmel 1999). The velocity time series is first converted to acceleration and then bandpass-filtered between 2 and 50 Hz before computing the autocorrelation. This is because nonlinear changes affect mostly high frequencies as shown in the synthetic example above; thus accelerograms show this better than velocity or displacement time histories. Earthquake and noise signals are processed with the same procedure. We compute ACF's for P and S waves, as depicted in Fig. 3(a), for 31 events having S/N larger than 3 (Table 1). The onset of the P and S arrival times were automatically picked (Ross & Ben-Zion 2014) and manually corrected if needed. Finally, the mean of all earthquake ACF's is obtained.

Fig. 6 presents example results of mean ACF's from noise and earthquake data at station R0101. The three panels show autocorrelations at different frequency bands. We note that for all cases, the ACF's based on the P wave and even more so the S wave lag behind those of the noise. The same procedure is applied at different

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Example of mean ACF from noise and earthquake data at station R0101. (a) ACF from bandpass filtered data between 2 and 50 Hz. (b) Narrow-band ACF between 10 and 20 Hz. (c) Narrow-band ACF between 20 and 40 Hz.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** PDF of the predominant frequency (left) and relative apparent velocity changes of P- and S-wave windows with respect to noise (right) for all sensors at the SBG site.</figcaption></figure>

frequency bands for all stations that have the selected earthquake data. Fig. 7 (left) shows the obtained probability density function (PDF) of the predominant frequency, computed as the inverse of the first positive peak from the ACF's (Valero & Alias 2013), for the noise and earthquake (P- and S-wave windows) data. The PDF of the apparent velocity changes of P and S waves from the delays of the ACFs with respect to noise are summarized in Fig. 7 (right).
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Spatial distribution of relative apparent temporal changes of seismic velocity for the P-wave windows with respect to noise for the four selected frequency bands.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Spatial distribution of relative apparent temporal changes of seismic velocity for the S-wave windows with respect to noise for the four selected frequency bands.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** PDF of the relative apparent velocity changes of *P*- and *S*-wave windows with respect to the mean ACF of earthquake data for all sensors at the SBG site.</figcaption></figure>

The predominant frequency decreases progressively from the noise data to *P*- and *S*-wave windows. The largest difference is between the noise and the *S* wave. This is expected because *S*-waves have lower frequencies than *P* waves, along with larger amplitude that can produce more nonlinear behaviour and changes of properties in the very shallow layers. We emphasize that the used events have low amplitude ground motion. Since the data are recorded on vertical components, the noise might be composed at low frequencies mainly of Rayleigh waves. However, the presented results are associated with high-frequency wavefield and the noise at these frequencies likely contains a mixture of signals as shown by Meng & Ben-Zion (2018b).

The apparent velocity changes for both the *P* and *S* waves have clear frequency dependency. The 10–15 and 10–20 Hz bands for the *P*-wave window have mainly positive changes with some small negative values. The two higher frequency bands have larger apparent velocity reduction with a mean of 2.5 per cent. Since the amplitude of *S* waves is larger than *P* waves, the apparent velocity changes are shifted to greater negative values with the majority around 6 per cent. Larger apparent velocity changes at higher frequencies suggest the presence of non-consolidated material at the subsurface. Given the low amplitude of the recorded wavefield, we cannot use lower frequencies and resolve deeper effects. Yet, this method involves the use of the first peak of the autocorrelation only and may allow resolving deeper structural effects with stronger sources producing larger low-frequency components.

Figs 8 and 9 display the spatial variability of the apparent velocity changes obtained for the *P*- and *S*-wave windows at the four selected frequency bands. The values of the apparent velocity variations at each sensor are represented by the colour scale. The obtained spatial variability, at a scale of around 10 m (which corresponds to the

mean separation between sensors in the ‘column’ direction), likely reflects the presence of strong heterogeneities in the shallow crust at the SGB site, as found in previous imaging and geological studies (e.g. Meng & Ben-Zion 2018a; Mordret et al. 2019; Share et al. 2020).

The part of the study area to the SW of the fault traces has a more heterogeneous distribution of apparent velocity changes for both *P* and *S* waves in all examined frequency bands. This region with transition in the topography is associated with heterogeneous amplification of the ongoing seismic wavefield (Ben-Zion et al. 2015, Fig. 6). We also observe that the eastern branch of the main Clark fault (MCF) delineates a clear separation between regions of low and large apparent velocity changes for both *P* and *S* waves (Figs 8 and 9). A similar zoning of seismic velocities was found by Qin et al. (2018) when studying residuals of *P* waves from teleseismic events recorded by the dense array and by Mordret et al. (2019) with detailed noise-based tomography of the shallow structure at the site.

The previous results use noise as a reference to compute the apparent velocity changes across the network. However, we can also estimate those changes using earthquake data only, involving a comparison with the same wavefield characteristics. This is done by computing the velocity change of each event relative to the mean ACF at each station. Fig. 10 shows the PDF, for the whole array, of the 50-percentile apparent velocity changes at each station. There are positive and negative values at all frequency bands, but for frequencies above 20 Hz the distribution is skewed toward negative values. In this case, the *P* and *S* waves produce average apparent velocity changes of 0.1 per cent and 0.3 per cent, respectively. Figs 11 and 12 show the spatial distribution of Fig. 10 for the *P* and *S* waves across the array, respectively. The frequency band of 20–40 Hz still indicates different material response for the sections separated by the MCF traces, especially for the *S*-wave response, confirming the previous results using noise as a reference.

## 3.2 Analysis of continuous waveforms

In the previous section, we use at each station the mean ACF from noise data as a reference to compute velocity changes produced by earthquakes. In this section, we calculate velocity changes from continuous recordings. We follow Bonilla et al. (2019) where only earthquake signals were used and the reference was defined as the mean ACF for accelerations lower than $0.01 \text{ m s}^{-2}$ for time windows before the PGA. We estimate ACF’s using moving windows of 2 s length with 25 per cent overlap. The time resolution of the obtained ACF is 0.5 s. The mean ACF is used as a reference for computing the evolution of apparent velocity changes. As in the previous section, we compute the autocorrelations using the method of Schimmel (1999) at the four frequency bands 10–15, 10–20, 15–30 and 20–40 Hz. In the following, we present results for the waveform containing the largest event shown in Fig. 2, between 1790 and 1910 s, at sensors of row 10 that cross four branches of the fault at columns 24, 30, 32 and 47 approximately (Fig. 1). The event has a magnitude of 2.3 and epicentral distance from the array of about 10 km.

Fig. 13 shows results for the frequency band of 10–20 Hz. The acceleration time history at station R7010 is displayed to show the ground motion duration and the *P* and *S* waves. The mean peak ground velocity for all stations computed at each moving time window is also plotted in logarithmic scale. The duration of the coda and surface waves lasts at least five times longer than the *P*
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Spatial distribution of relative apparent temporal changes of seismic velocity for the *P*-wave windows with respect to the mean ACF of earthquake data for the four selected frequency bands.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Spatial distribution of relative apparent temporal changes of seismic velocity for the *S*-wave windows with respect to the mean ACF of earthquake data for the four selected frequency bands.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** (a) Example of the waveform recorded at station R7010 for the frequency band of 10–20 Hz. The solid line shows the mean peak ground velocity of all sensors in row 10 in logarithmic scale. (b) Evolution of apparent velocity changes at each station in row 10. The horizontal dashed lines show the location of the fault branches in Fig. 1. (c) Mean velocity changes of all sensors in row 10.</figcaption></figure>

and *S* waves (Fig. 13a). The time evolution of the apparent velocity changes for all stations in row 10 is shown in Fig. 13(b) together with the location of the surface fault traces (dashed lines). Clear apparent velocity drops are observed in the data of all stations, with spatial and temporal variability across the array. The apparent velocity changes at the times of the *P* and *S* waves are coherent along lines of sensors. The mean of the stacked results (Fig. 13c) have values of 19 per cent and 8 per cent during the times of the *P* and *S* waves, respectively, with a slow recovery that lasts approximately until 1880 s. Positive apparent velocity changes slowly begin around 1870 s for the sensors in columns between 32 and 70. On the other hand, for stations between columns 10 and 32, there is a strong change from negative to positive apparent velocity changes. For this particular frequency band, the largest velocity reductions occur between columns 24 and 30 in the low velocity zone between the two branches of the Clark fault (e.g. Mordret *et al.* 2019; Share *et al.* 2020). Corresponding results for the other three other frequency bands are presented in the Supporting Information (Figs S1–S3).

The mean velocity changes in Fig. 13(c) shows a strong velocity drop at the beginning of the earthquake signal (i.e. *P*-wave window). This is also seen in Figs S1–S3 in the Supporting Information and reflects the transition from ambient noise to the *P*-wave window. This confirms that the results depend on the reference chosen to compute the velocity changes, which is another reason why we use the term apparent velocity changes. Furthermore, the relatively large values in Figs 8 and 9 can be explained because of this effect; therefore, the apparent velocity changes in Figs 10–12 are probably closer to the real changes in the upper crust because there is a common reference using earthquake data only. The larger variations at the onset of the earthquake signal could be used to identify events recorded over a large frequency band. This is not the purpose of the current study, but the method is fast enough to implement this for continuous recordings.

## 4 DISCUSSION AND CONCLUSIONS

We use ACFs of noise and earthquake data recorded by a dense seismic array to compute maps of apparent velocity changes across the SGB site on the SJFZ (Figs 8–12). Numerical simulations indicate that the derived values reflect a combination of changes of seismic velocities and nonlinear material response of the top material (Fig. 5). The simulations also show that changes in the top 10 m of the crust can produce time delays for long propagation distances (Fig. 5c), in agreement with studies on effects of shallow velocity changes on large-scale wavefields (Yang *et al.* 2018; Juarez & Ben-Zion 2020). The observational analyses on the dense array data reveal significant apparent velocity changes during the passage of low amplitude *P* and *S* waves from small local earthquakes with magnitudes below 3.1. The observed changes in response to low amplitude incoming wavefield demonstrate that nonlinear processes in the shallow crust are more pervasive than is typically assumed.

The spatial variations of the apparent velocity changes are correlated with the shallow velocity structure at the SGB site. The NE part of the site has some clear zones of large negative values of apparent velocity changes for both *P*- and *S*-wave windows, especially for the frequency bands of 15–30 and 20–40 Hz (Figs 8–12). One zone of significant apparent velocity reductions extends parallel to the eastern branch of the MCF (between columns 30–40) up to the fault junction (row 10 and column 30). This zone corresponds to a seismic trapping structure reflecting damaged fault zone rocks (Ben-Zion *et al.* 2015; Qin *et al.* 2018), and it continues beyond that fault junction about 10 m to the SE between columns 40 and 50. The region between columns 40 and 60 associated with relatively high topography (Fig. 1), low seismic velocities (Mordret *et al.* 2019) and a recently mapped structure (Share *et al.* 2020) also has significant apparent velocity reductions. The general correlations between the derived changes and damaged low velocity zones at the site are consistent with expectations based on damage rheology models (e.g. Lyakhovsky *et al.* 1997; Sens-Schönfelder *et al.* 2018).

The obtained apparent velocity changes and associated recovery are clearly frequency dependent (Fig. 13 and Figs S1–S3, Supporting Information), with larger changes and faster recovery times for higher frequencies. This suggest the presence of very shallow non-consolidated material that is highly affected during weak ground motion. The healing phase lasts much longer than the duration of the direct *P* and *S* waves and continues during the occurrence of the coda and surface waves. At lower frequencies sampling deeper material, the small amplitude signals produce very small apparent velocity changes as seen in the frequency band of 10–15 Hz (Fig. S1, Supporting Information). Stronger incoming wavefield associated with larger sources can produce significant changes at frequencies lower than 10 Hz, as observed during the Tohoku-Oki Mw9 2011 earthquake in Japan at a site where the PGA reached 0.6g (Bonilla *et al.* 2019).

Fig. 14 displays relative changes $v/v_0$ as a function of the mean peak ground velocity shown in Fig. 13 and Figs S1–S3 in the Supporting Information for the four frequency bands. The results use $v/v_0 = 1 + \Delta v/v_0$ and have been binned to have at least three points per decade. The use of the peak ground velocity allows estimating the dynamic strain from the ratio between the PGV and wave speed, if the local velocity is known (Chandra *et al.* 2015, 2016). The observations exhibit frequency dependency of the relative apparent velocity changes, with increasing values as the PGV becomes larger. The frequency band 10–15 Hz shows very small changes, whereas the 20–40 Hz band indicate relative changes up to 10 per cent and the frequency bands 10–20 and 15–30 Hz have intermediate values.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Relative apparent velocity changes *v*/*v*₀ with respect to noise for all stations in row 10 as function of mean peak ground velocity for the four used frequency bands.</figcaption></figure>

The derived values of changes depend on the reference used to compute Δv/v₀ and the analysis time step. In this study, when using a reference based on the average of ACF's for 1 hr of noise data, the derived relative apparent velocity changes are surprisingly high given the low PGV values and associated dynamic strains, 2.5 per cent and 6 per cent for the *P* and *S* waves, respectively. The use of noise data as a reference to compute apparent velocity changes in earthquake waveforms may seem inappropriate since the wavefield composition is different. In this study, we are interested in relative velocity changes and the use of noise provides a reasonable stable reference. We note, however, that the high frequency noise field at the site includes significant energy propagating vertically from below that is likely generated by ongoing small-scale earthquakes in the shallow crust (e.g. Hillers et al. 2016; Gradon et al. 2020). That said, when using the average of earthquake data only, the related changes are greatly reduced by at least one order of magnitude, for the *P*- and *S*-wave time windows, to 0.1 per cent and 0.3 per cent, respectively. These results are likely more coherent with the dynamic strain levels induced by the recorded events at the site. Continuous velocity changes are followed by rapid recoveries of properties as has also observed in other studies (e.g. Dieterich & Kilgore 1996; Nakatani & Scholz 2004; Wu et al. 2009). However, the material recovery is frequency dependent, faster at higher frequencies, and it is probably related to the mobilized volume of material changes at different wavelengths. The derived results are likely dominated by changes in the top tens of metres of the crust and may involve local opening and closing of fractures during the passage of earthquake waves (e.g., Lyakhovsky et al. 2009; Ostrovsky et al. 2019; Qin et al. 2020). We note that dynamic reduction of elastic moduli of material under elastic strain is expected to produce isotropic radiation that can aid the failure process (Ben-Zion & Ampuero 2009; Ben-Zion & Lyakhovsky 2019). For shallow material under low initial strain

the effect is small, but so is the confining pressure that governs the material strength.

The application of the ACF-based technique to earthquake waveform data recorded at high sampling rate by a dense seismic array allows us to derive high-resolution results on temporal changes of properties at the top crust. The results demonstrate significant changes of properties, up to one order of magnitude higher, may occur in response to weak ground motion with PGV around 10⁻⁷ m s⁻¹ (Fig. 14), compared to traditional studies using only noise correlation data. The onset of nonlinear behaviour and changes of seismic properties at such low levels of dynamic strain are consistent with results of laboratory experiments (e.g. Pasqualini et al. 2007; Remillieux et al. 2017). The use of the first negative peak value of the ACF limits the analysis to small time delays corresponding to the first reflector sensitive to a given frequency band. Using longer time delays may provide more information about the structure and changes of properties at greater depth (e.g. Pham & Tkaick 2017; Romero & Schimmel 2018).

These results show that nonlinear material behaviour and temporal changes of properties are far more pervasive that commonly thought. Relative apparent velocity reductions, with respect to the mean ACF of earthquake data, up to 2 per cent are observed to be produced by very small ground motion generated by frequent low magnitude earthquakes (Fig. 10). The changes are obtained using relatively high-frequency bands, 15–40 Hz, and are hence limited to very shallow depths. However, deeper changes may be resolved for larger sources producing wavefields with larger amplitudes and longer wavelengths. The increasing availability of high-quality data recorded by dense seismic arrays in areas with significant seismicity provide important opportunities to extend the results of this study to analyses of property changes produced by stronger sources. Analysis of some such data will be attempted in a follow-up study.
## ACKNOWLEDGEMENTS

We thank Lei Qin for help with pre-processing of the data used in the study and thoughtful discussions. The manuscript benefitted from constructive comments by Marine Denolle and an anonymous reviewer. We use Matplotlib, Librosa, Pandas and Obspy, Python libraries for processing and plotting the data and results. The study was supported by the Southern California Earthquake Center (based on NSF cooperative agreement EAR-1600087 and USGS cooperative agreement G17AC00047), the U.S. Department of Energy (award number DE-SC0016520) and the collaborative agreement 2019–2020 between the University Gustave Eiffel and the Ministry of Ecology (Direction of Disaster Prevention, DPRG) of France on earthquake risk and hazard studies.

## REFERENCES

Aben, F.M., Brantut, N., Mitchell, T.M. & David, E.C., 2019. Rupture energetics in crustal rock from laboratory-scale seismic tomography, *Geophys. Res. Lett.*, **46**, 7337–7344.

Ben-Zion, Y. & Ampuero, J.-P., 2009. Seismic radiation from regions sustaining material damage, *Geophys. J. Int.*, **178**, 1351–1356.

Ben-Zion, Y., et al., 2015. Basic data features and results from a spatially-dense seismic array on the San Jacinto fault zone, *Geophys. J. Int.*, **202**, 370–380.

Ben-Zion, Y. & Lyakhovsky, V., 2019. Representation of seismic sources sustaining changes of elastic moduli, *Geophys. J. Int.*, **217**, 135–139.

Ben-Zion, Y. et al., 2015. Basic data features and results from a spatially-dense seismic array on the San Jacinto fault zone, *Geophys. J. Int.*, **202**, 370–380.

Beresnev, I.A. & Wen, K.L., 1996. Nonlinear site response—a reality? *Bull. seism. Soc. Am.*, **86**, 1964–1978.

Bonilla, L.F., Archuleta, R.J. & Lavallée, D., 2005. Hysteretic and dilatant behavior of cohesionless soils and their effects on nonlinear site response: Field data observations and modeling, *Bull. seism. Soc. Am.*, **95**, 2373–2395.

Bonilla, L.F., Guéguen, P. & Ben-Zion, Y., 2019. Monitoring coseismic temporal changes of shallow material during strong ground motion with interferometry and autocorrelation, *Bull. seism. Soc. Am.*, **109**(1), 187–198.

Bonilla, L.F., Tsuda, K., Pulido, N., Régnier, J. & Laurendeau, A., 2011. Nonlinear site response evidence of K-NET and KiK-net records from the 2011 off the Pacific coast of Tohoku earthquake, *Earth Planets Space*, **63**, 785–789.

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**, 80–82.

Brenguier, F., Shapiro, N.M., 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. et al., 2019. Turning vehicle traffic into a powerful seismic source for monitoring active faults, *Geophys. Res. Lett.*, **46**, 9529–9536.

Chandra, J., Guéguen, P. & Bonilla, L.F., 2016. PGA-PGV/V~S~ considered as a stress–strain proxy for predicting nonlinear soil response, *Soil Dynam. Earthq. Eng.*, **85**, 146–160.

Chandra, J., Guéguen, P., Steidl, J.H. & Bonilla, L.F., 2015. In-situ assessment of the G-γ curve for characterizing the nonlinear response of soil: application to the Garner Valley Downhole Array (GVDA) and the Wildlife Liquefaction Array (WLA), *Bull. seism. Soc. Am.*, **105**(2A), 993–1010.

Claerbout, J.F., 1968. Synthesis of a layered medium from its acoustic transmission response, *Geophysics*, **33**(2), 264–269.

Cohen, L., 1989. Time-frequency distributions—a review, *Proc. IEEE*, **77**, 941–981.

Dieterich, J.H. & Kilgore, B.D., 1996. Imaging surface contacts; power law contact distributions and contact stresses in quarts, calcite, glass, and acrylic plastic, *Tectonophysics*, **256**, 219–239.

Field, E.H., Johnson, P.A., Beresnev, I.A. & Zeng, Y., 1997. Nonlinear ground-motion amplification by sediments during the 1994 Northridge earthquake, *Nature*, **390**, 599–602.

Froment, B.M., Campillo, M., Chen, J.H. & Liu, Q.Y., 2013. Deformation at depth associated with the 12 May 2008 MW 7.9 Wenchuan earthquake from seismic ambient noise monitoring, *Geophys. Res. Lett.*, **40**, 78–82.

Gradon, C., Roux, P., Moreau, L., Lecointre, A. & Ben-Zion, Y., 2020. Characterization with dense array data of seismic sources in the shallow part of the San Jacinto fault zone, *Geophys. J. Int.*, doi:10.1093/gji/ggaa411.

Hamiel, Y., Lyakhovsky, V., Stanchits, S., Dresen, G. & Ben-Zion, Y., 2009. Brittle deformation and damage-induced seismic wave anisotropy in rocks, *Geophys. J. Int.*, **178**, 901–909.

Hillers, G., Roux, P., Campillo, M. & Ben-Zion, Y., 2016. Focal spot imaging based on zero lag cross correlation amplitude fields: Application to dense array data at the San Jacinto fault zone, *J. geophys. Res.*, **121**, 8048–8067.

Hobiger, M., Wegler, U., Shiomi, K. & Nakahara, H., 2016. Coseismic and post-seismic velocity changes detected by Passive Image Interferometry: comparison of one great and five strong earthquakes in Japan, *Geophys. J. Int.*, **205**(2), 1053–1073.

Iai, S., Matsunaga, Y. & Kameoka, T., 1990. Strain space plasticity model for cyclic mobility, *Rep. Port Harbour Res. Inst.*, **29**, 27–56.

Johnson, C.W., Meng, H., Vernon, F. & Ben-Zion, Y., 2019. Characteristics of ground motion generated by wind interaction with trees, structures, and other surface obstacles, *J. geophys. Res.*, **124**, 8519–8539.

Johnson, P.A. & Jia, X., 2005. Nonlinear dynamics, granular media and dynamic earthquake triggering, *Nature*, **473**, 871–874.

Juarez, A. & Ben-Zion, Y., 2020. Effects of shallow velocity reductions on three-dimensional propagation of seismic waves, *Seismol. Res. Lett.*, **91**, 3313–3322.

Lockner, D., Walsh, J. & Byerlee, J., 1977. Changes in seismic velocity and attenuation during deformation of granite, *J. geophys. Res.*, **82**, 5374–5378.

Lyakhovsky, V., Ben-Zion, Y. & Agnon, A., 1997. Distributed damage, faulting, and friction, *J. geophys. Res.*, **102**, 27635–27649.

Lyakhovsky, V., Hamiel, Y., Ampuero, J. & Ben-Zion, Y., 2009. Non-linear damage rheology and wave resonance in rocks, *Geophys. J. Int.*, **178**, 910–920.

Meng, H. & Ben-Zion, Y., 2018a. Detection of small earthquakes with dense array data: example from the San Jacinto fault zone, southern California, *Geophys. J. Int.*, **212**, 442–457.

Meng, H. & Ben-Zion, Y., 2018b. Characteristics of airplanes and helicopters recorded by a dense seismic array near Anza California, *J. geophys. Res.*, **123**, 4783–4797.

Mordret, A., Roux, P., Boué, P. & Ben-Zion, Y., 2019. Shallow 3-D structure of the San Jacinto Fault zone revealed from ambient noise imaging with a dense seismic array, *Geophys. J. Int.* doi:0.1093/gji/ggy464.

Nakata, N. & Snieder, R., 2012. Estimating near-surface shear wave velocities in Japan by applying seismic interferometry to KiK-net data, *J. geophys. Res.*, **117**, doi:10.1029/2011JB008595.

Nakatani, M. & Scholz, C. H., 2004. Frictional healing of quartz gouge under hydrothermal conditions: 1. Experimental evidence for solution transfer healing mechanism, *J. geophys. Res.*, **109**, B07201, doi:10.1029/2001JB001522.

Ostrovsky, L., Lebedev, A., Riviere, P., Shokouhi, C., Wu, C., Stuber Geesey, M.A. & Johnson, P.A., 2019. Long-time relaxation induced by dynamic forcing in geomaterials, *J. geophys. Res.*, **124**, 5003–5013.

Pasqualini, D., Heitmann, K., TenCate, J.A., Habib, S., Higdon, D. & Johnson, P., 2007. Nonequilibrium and nonlinear dynamics in Berea and Fontainebleau sandstones: low-strain regime, *J. geophys. Res.*, **112**, B01204, doi:10.1029/2006JB004264.

Pei, S., Niu, F., Ben-Zion, Y., Sun, Q., Liu, Y., Xue, X. & Shao, Z., 2019. Seismic velocity reduction and accelerated recovery due to earthquakes on the Longmenshan fault, *Nat. Geosci.*, **12**(5), 387–392.

Peng, Z. & Ben-Zion, Y., 2006. Temporal changes of shallow seismic velocity around the Karadere-Duzce branch of the north Anatolian fault and strong ground motion, *Pure appl. Geophys.*, **163**, 567–599.
Pham, T.-S. & Tkalčić, H., 2017. On the feasibility and use of teleseismic P wave coda autocorrelation for mapping shallow seismic discontinuities, *J. geophys. Res. Solid Earth*, **122**, 3776–3791.

Poupinet, G., Ellworth, W.L. & Frechet, J., 1984. Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras Fault, California, *J. geophys. Res.* **89**, 5719–5731.

Press, W.H., Teukolsky, S.A., Vetterling, W.T. & Flannery, B.P., 1992. *Numerical Recipes in FORTRAN*, 2nd ed., Cambridge: University Press.

Qin, L., Ben-Zion, Y., Bonilla, L.F. & Steidl, J.H., 2020. Imaging and monitoring temporal changes of shallow seismic velocities at the Gardner Valley near Anza, California, following the M7.2 2010 El Mayor-Cucapah earthquake, *J. geophys. Res.: Solid Earth*, **125**, e2019JB018070, doi:10.1029/2019JB018070.

Qin, L., Ben-Zion, Y., Qiu, H., Share, P.E., Ross, Z.E. & Vernon, F.L., 2018. Internal structure of the San Jacinto fault zone in the trifurcation area southeast of Anza, California, from data of dense seismic arrays, *Geophys. J. Int.*, **213**, 98–114.

Remillieux, M.C., Ulrich, T.J., Goodman, H.E. & Ten Cate, J.A., 2017. Propagation of a finite-amplitude elastic pulse in a bar of Berea sandstone: a detailed look at the mechanisms of classical nonlinearity, hysteresis, and non equilibrium dynamics, *J. geophys. Res.: Solid Earth*, **122**, 8892–8909.

Romero, P. & Schimmel, M., 2018. Mapping the basement of the Ebro Basin in Spain with seismic ambient noise and autocorrelations, *J. geophys. Res.: Solid Earth*, **123**, doi:10.1029/2018JB015498.

Ross, Z. & Ben-Zion, Y., 2014. An earthquake detection algorithm with pseudo probabilities of multiple indicators, *Geophys. J. Int.*, **197**, 458–463.

Roux, P. & Ben-Zion, Y., 2014. Monitoring fault zone environments with correlations of earthquake waveforms, *Geophys. J. Int.*, **196**, 1073–1081.

Rubinstein, J.L. & Beroza, G.C., 2004. Evidence for widespread nonlinear strong ground motion in the Mw6.9 Loma Prieta Earthquake, *Bull. seism. Soc. Am.*, **94**, 1595–1608.

Ruigrok, E. & Wapenaar, K., 2012. Global-phase seismic interferometry unveils P-wave reflectivity below the Himalayas and Tibet, *Geophys. Res. Lett.*, **39**, L11301, doi:10.1029/2012GL051672.

Régnier, J., Cadet, H., Bonilla, L.F., Bertand, E. & Semblat, J., 2013. Assessing nonlinear behavior of soil in seismic site response: Statistical analysis on KiK-net strong motion data, *Bull. seism. Soc. Am.*, **103**(3), 1750–1770.

Sawazaki, K., Sato, H., Nakahara, H. & Nishimura, T., 2006. Temporal change in site response caused by earthquake strong motion as revealed from coda spectral ratio measurement, *Geophys. Res. Lett.*, **33**, L21303, doi:21310.21029/22006GL027938.

Schaff, D.P. & Beroza, G.C., 2004. Coseismic and postseismic velocity changes measured by repeating earthquakes, *J. geophys. Res.*, **109**, B10302, doi:10.1029/2004JB003011.

Schimmel, M., 1999. Phase cross-correlations: design, comparisons, and applications, *bull. seism. soc. am.*, **89**(5), 1366–1378.

Sens-Schönfelder, C., Snieder, R. & Li, X., 2018. A model for nonlinear elasticity in rocks based on friction of internal interfaces and contact aging, *Geophys. J. Int.*, **216**, 319–331.

Sens-Schönfelder, C. & Wegler, U., 2006. Passive image interferometry and seasonal variations of seismic velocities at Merapi volcano, Indonesia, *Geophys. Res. Lett.*, **33**, L21302, https://doi.org/10.1029/2006GL027797.

Share, P.-E. et al., 2020. Characterizing the uppermost 100 m structure of the San Jacinto fault zone southeast of Anza, California, through joint analysis of geologic, topographic, seismic and resistivity data, *Geophys. J. Int.*, **222**, 781–794.

Stanchits, S., Vinciguerra, S. & Dresen, G., 2006. Ultrasonic velocities, acoustic emission characteristics and crack damage of basalt and granite, *Pure appl. Geophys.*, **163**, 974–993.

Theodulidis, N., Bard, P.Y., Archuleta, R. & Bouchon, M., 1996. Horizontal-to-vertical spectral ratio and geological conditions: the case of Garner Valley Downhole Array in Southern California, *Bull. seism. Soc. Am.*, **86**(2), 306–319.

Valero, X. & Alias, F., 2013. Narrow-band autocorrelation function features for the automatic recognition of acoustic environments, *J. acoust. Soc. Am.*, **134**, 880–890.

Ventosa, S., Schimmel, M. & Stutzmann, E., 2019. Towards the processing of large data volumes with phase cross-correlation, *Seismol. Res. Lett.*, **90**, 1663–1669.

Wapenaar, K., 2003. Synthesis of an inhomogeneous medium from its acoustic transmission response, *Geophysics*, **68**(5), 1756–1759.

Wu, C., Peng, Z. & Ben-Zion, Y., 2009. Non-linearity and temporal changes of fault zone site response associated with strong ground motion, *Geophys. J. Int.*, **176**, 265–278.

Wu, C., Peng, Z. & Ben-Zion, Y., 2010. Refined thresholds for nonlinear ground motion and temporal changes of site response associated with medium size earthquakes, *Geophys. J. Int.*, **182**, 1567–1576.

Yang, W., Wang, B., Yuan, S. & Ge, H., 2018. Temporal variation of seismic-wave velocity associated with groundwater level observed by a downhole airgun near the Xiaojiang Fault Zone, *Seismol. Res. Lett.*, **89**(3), 1014–1022.

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

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Figure S1.** Same as Fig. 13 for the frequency band of 10–15 Hz.

**Figure S2.** Same as Fig. 13 for the frequency band of 15–30 Hz.

**Figure S3.** Same as Fig. 13 for the frequency band of 20–40 Hz.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.