Peak Ground Motions and Site Response at Anza and Imperial Valley, California

Power spectra of shear-waves for eighteen earthquakes from the Anza-Imperial Valley region were inverted for source, mid-path Q, site attenuation and site response. The motivation was whether differences in site attenuation (parameterized as t*, r/cQ, where r is distance along ray path near the site, c is shear velocity and Q is the quality factor that parameterizes attenuation) and site response could be correlated with residuals in peak values of velocity or acceleration after removing the affect of distance-dependent attenuation. We decomposed spectra of S-waves from horizontal components of 18 earthquakes from 2010 to 2018 into a common source for each event with ω−2 spectral fall-off at high frequencies and then projected the residuals onto path and site terms following the methodology of Boatwright et al. (Bull Seismol Soc Am 81:1754–1782, 1991). The site terms were constrained to have an amplification at a particular frequency governed by VS30 at two of the sites which had downhole shear-wave logs. The 18 events, 3 < M < 4, had moments between approximately 1020 and 1022 dyne-cm, and stress drops between 1 and 100 bars. Average mid-crust attenuation had a Q of 844 reflecting the average path through the crystalline rock of the San Jacinto Mountains. t* for each station corresponded to the geologic environment such that stations on hard rock had low t* (e.g. stations KNW, PFO and RDM) a station in the San Jacinto fault zone (station SND) had a moderate t* of 0.035 s and stations in the Imperial Valley usually had higher t*s. Generally t* correlated with average amplification suggesting that sites characterized by low surface velocities and higher attenuation also have more amplification in the 1–6 Hz band. Residuals of peak values were determined by subtracting the prediction of Boore and Atkinson (2008). There is a correlation between average amplification and peak velocity, but not peak acceleration. Interestingly, there is less scatter at high values of amplification although there is also less data. Scatter in values of peak velocity and peak acceleration are higher at shorter compared to longer durations. When using a frequency-dependent form for Q, variances are higher, sometimes much higher; the dataset does not support frequency-dependent Q, which is not similar to results from the Imperial Valley and northeastern North America.


Introduction
Site response is known to contribute to the character of ground motion (Borcherdt and Gibbs 1976;Tucker and King 1984;Seed et al. 1987;Flores-Estrella et al. 2007;Hartzell et al. 1996;Molimari et al. 2015;Kaklamanos and Bradley 2018 for example). However, even though there is a general correspondence of higher ground motion at sites on soft soils, actually predicting quantitative values based on a given site conditions is difficult. In this paper, we try to determine several average parameters that characterize a site response and then compare them with peak velocity and accelerations for 18 earthquakes near Anza, CA (Table 1). This region is particularly well suited to this study because of the strong contrast in site conditions across the array from hard rock sites on exposed granite to sites on soft rock located in fault zones and deep sedimentary basins of the neighboring Imperial Valley. There is also a large catalog of waveforms from stations at Anza and Imperial Valley stations, which are part of the Southern California Seismic Network. We incorporate a method pioneered by Boatwright et al. (1991), that decomposes source, path and site contributions in the spectral domain. This algorithm evolves from Andrews (1986) study of Mammoth lakes aftershocks where he used a generalized inversion to determine a common source spectra and individual site spectra. In Boatwright et al. (1991) constraints are applied by using the expected amplification at sites where the velocity structure is well known as well as solving for an Brune, x -2 , source rather than components at a set of frequencies. The method has previously been used in the Marina district of San Francisco after the Loma Prieta earthquake ) and along the San Francisco Peninsula using aftershock data from the Loma Prieta earthquake . Many others have also used the constrained inversion to determine site response and source parameters (e.g. Tsuda et al. 2006;Yefei et al. 2013). At Anza, we find that there is a correlation of peak velocity with site conditions as determined using the method of Boatwright et al. (1991), but simple models of the site will not lead to predicted values of amplification except at the more extreme cases of hard or soft rock. We also explore whether duration has an effect on peak values and if there is evidence of a frequency-dependent Q.

Data and Method
This method inverts spectra of shear-waves from horizontal components and then decomposes the contribution from source and site. The seismograms being inverted are from the Anza array (Berger et al. 1984) and neighboring stations in the Imperial Valley from the Southern California Seismic Network (SCSN). All stations have three component sensors and are broadband with natural periods typically of 120 s. They are usually digitized at 100 samples per second (sps) although some use a sample rate of 200 sps with a 24-bit analog-to-digital converter. All seismograms are interpolated to a sample rate of 200 sps if necessary. Seismograms from the Anza and SCSN networks were instrument-corrected and then displayed to manually pick windows surrounding the S-wave arrival on the horizontal components. Typically the window was started a few tenths of a second before the S-wave arrival to prevent tapering the main signal and ended close to where the envelope of the major S-wave arrivals were moving into coda. A short window between the P and S wave is used as an estimate of noise. To emphasize S-waves we only use data from horizontal components. Power spectrum of the windowed velocity traces are determined and is the dataset used in the inversion. Spectral amplitudes were normalized to a standard window length of 10 s.
The inversion Fletcher and Boatwright 1991) is briefly summarized here. The method expands on the original idea of Andrews (1986) that a spectrum recorded at a particular site can be thought of as the multiplication (in Fourier space) of a source spectrum and site response. Consequently an inversion scheme is constructed such that a set of shear-wave spectra are inverted for a common source and a response spectrum for each site. Following Eq. 2 of Andrews (1986) taking logarithms of the source, site, and observed spectra where k is the index for each record R, i is the index for the station, and j is the index for each earthquake, E are the source spectra and S are the site spectra. Fourier spectra for each record are inverted for site and source spectra at the same Fourier components. This method suffers by having a degree of ambiguity in that a factor of x could be taken out of the source and put in the site spectra. This was a much better way to interpret site-specific spectra rather than just averaging S-wave displacement spectra because it explicitly tried to separate site from source effects. Boatwright et al. (1991) expands on this method by assuming a Brune source with x -2 fall-off above a corner frequency and a long period level X and then parameterizing the attenuation in the following way where g is the attenuation factor, r is the hypocentral distance, f is frequency, c is the exponent to r that defines the distance component such as geometrical spreading, t* is a site attenuation factor (r/cQ), c is body wave velocity, T is the total travel time and Q is the quality factor (for attenuation). This assumes that a constant mid-crustal Q is an adequate model for this region. However, with such a dramatic difference between the Anza region and the Imperial Valley this is not likely to be the case. Any biases between the two regions will probably be reflected in a bias in values for t*. We use t* to describe site attenuation (Solomon 1973) because t* can be assigned to different layers and is additive. This splits the total attenuation into a regional value common to all paths and site-specific terms. Geometrical spreading is defined to be r -1 for distances less than 50 km and r -1/2 at greater distances. After taking natural logs of the model and attenuation terms the equation to be inverted for becomes where f n are the frequency components, B 4 is a term from expanding the source term in a Taylor series (see Boatwright et al. 1991) and the u 0 are the long period levels of the source spectrum . The inversion runs in two steps; the first solves for source and wave propagation effects and the second projects the residuals from the first step onto site response . Site responses are further constrained by using estimates of V S30 and basement shear-wave velocity to determine an amplification from impedance at a particular frequency using the quarter wavelength rule. The amplification from impedance is where S is the amplification and q is the density and v is the velocity. The numerator includes the values at the surface compared to values of the basement in the denominator (Joyner et al. 1976). The frequency at which the amplification is constrained is where v l is the velocity of the top layer and d l is the thickness of the layer. We used the velocity of a shallow layer near the surface at stations PFO and KNW using bore logs from Fletcher et al. (1990). All responses are then rescaled in a least squares sense so that the responses at stations KNW and PFO have the appropriate amplification at the desired frequency from Eqs. (4) and (5). Velocities for stations PFO and KNW are given in Table 2 with the half space having a density of 2.7 gm/cm 3 and a shear wave velocity of 3.6 km/s. Figure 1 shows the location of the stations used in the study, which are the Anza stations and the nearby broadband stations of the SCSN in the Imperial Valley. Epicenters from the SCSN catalog are shown for all 18 earthquakes used in the study (see Table 1). Seismicity that trends northwestsoutheast fall along fault lines of the San Jacinto and Elsinore fault systems. There are cross cutting trends both in the Imperial Valley and Anza. The stations are bunched in two groups; one near Anza and the other in the Imperial Valley south of the Salton Sea. This configuration reflects historical interest in the San Jacinto fault zone and the Imperial and Brawley fault zones, respectively. The maximum number of

Results
The inversion for the 18 events and 21 sites (3 sites in 21 had no metadata or poor signal to noise ratio, SMTC, SSW, CLI2) was run with a final variance of 3.75% over a bandwidth of 0.4-25 Hz. The regional Q was 844 and the difference from 1 or 1/2 for c was -0.0113 or essentially the assumed model over a distance range of about 1-120 km. The average t* for all sites was 0.013 s before the V S30 constraint was applied.
The regional path Q value of 844 is marginally higher than the value of 830 for the San Francisco peninsula (Fletcher et al. 2003), which is characterized by the Franciscan terrane on the east side of the San Andreas fault and Salinian rocks on the west side. After the V S30 constraint is applied individual t* values range from 0.0 to 0.047 s and are given in Table 3. The lowest values are for stations KNW, RDM, and PFO with a value of 0.0 s to a high of 0.047 s at station WES. Generally these values correspond to the geology of the site. Stations KNW, PFO, and RDM are in the San Jacinto mountains (see Fig. 1) with exposed granites or intrusives at the site or close by (Vernon 1989). Station PFO has a layer of grus (decomposed granite) at the surface, but is closely underlain by hard rock. Sites in the Imperial Valley have high attenuation (high t*) reflecting the soft, low-velocity alluvium in the Imperial Valley. Figure 1 shows the geology of the region and most Anza stations are on pinkish colors, which are the Mesozoic granitic rocks of the San Jacinto or Santa Rosa Mountains. Stations in the Imperial Valley are mostly on river terraces or alluvium.
Values of t* can be compared to j  defined by Anderson and Hough (1984) to be a parameter that describes high-frequency roll-off when using a linear frequency axis that masks the corner frequency (Brune 1970;Hanks and Wyss 1972). The formula for attenuation is where tÃ ¼ r Qc Solomon (1973). From Anderson and Hough (1984) and Hough et al. (1988) aðf Consequently, t* is exactly j especially when j is interpreted as a parameter that describes attenuation. Kilb et al. (2012) determined j for many of the Anza stations using a large number of events in the Anza region. While they determined j in several ways they found many of the same features we also see in that stations PFO and KNW have small js, stations such as SND have larger js. For example, they found values of j at station SND and TRO that are about 25 ms greater than j values at stations KNW and PFO. These differences are similar to that found here of about 26-36 ms. Hough et al. (1988) also provide estimates of j and they separate an initial value j 0 from a distance dependent term. Their values at least have a constant bias compared with the results here. Stations TRO, SND and PFO have a j 0 value of 2.5, 18.1, and 3.5 ms, respectively, from Hough et al. (1988). While some of the trends in their results are similar to the results here such as a low value for station KNW, other values do not agree such as the higher value at station PFO (3.5 ms) compared to that at station TRO (2.5 ms). The main difference between these studies and the one used here is we try to decompose path from source and site response in the main inverse system rather than infer from systematics in the residuals.
Seismic response functions, after re-normalization, for four sites are shown in Fig. 2a-c. Each shows the response curve from 0.4 to 25.0 Hz. These curves are the remainder after removing the regional attenuation and the local t* appropriate for that site (given in Table 3). The error bars show the variance for each estimate. The response curves tend to have strong amplification in the 1-6 Hz range at sites on alluvium or in the San Jacinto fault zone. Sites that are on hard rock tend to have low amplification in the 1-6 Hz range, but have resonances above 10 Hz (e.g. station KNW). Apparent in Fig. 2 is a stark contrast between station ERR and KNW. The average 1-6 Hz amplification at ERR is over 10 but close to 2 at KNW. On the other hand there appears to be a resonance at stations KNW above 10 Hz. These responses are expected based on a simple model of soft/hard (high attenuation, high amplification at low frequencies versus low amplification and low attenuation, respectively) rock geology of the site. Station ERR is on the deep sediments of the Imperial Valley if not the NE-SW trending fault zones (e.g. Elmore Ranch fault) near the southwest corner of the Salton Sea. Station SND is also located on a fault zone, the San Jacinto fault zone. The response at station TRO does not easily conform to the soft/hard rock idea as it is on Toro Mountain, which is mostly granite or exposed intrusives (Lancaster et al. 2012). The largest attenuation (t* = 0.047 s) is at station WES, which is at the southern edge of Fig. 1 in the Imperial Valley.   Sahakian et al. (2018) determine sites terms from the residuals of applying a ground motion prediction equation (GMPE) to peak values from earthquakes in the Anza region. We compared these site terms to our estimates of 1-6 Hz amplification in Fig. 3a. There is no apparent correlation with high values of amplification at both ends of the range in site terms. Apparently the decomposition in Sahakian et al. (2018) does not decompose site terms in a way similar to Boatwright et al. (1991). On the other hand Klimasewski et al. (2019) uses a similar procedure to Andrews (1986), which is a less constrained decomposition compared to Boatwright et al. (1991). However, they normalize their results by constraining the source spectra to have a Brune (1970) like spectra with x -2 at high frequencies. Their results are quite similar to ours as shown in Fig. 3b, c with exception of their value of amplification for station ERR.
There is a general correspondence of t* with amplification as shown in Fig. 4. This is not surprising as sites with the most amplification typically are underlain by soft soils that would attenuate more. The lower the near surface velocity the more amplification and attenuation the site is likely to endure.

Effect of t* and Average Amplification on Peak Values
Peak values are used to determine ground motion prediction equations for the engineering community (Joyner and Boore 1989). These peak values, however, can scatter considerably about some mean such that estimates of error can be large. While low velocity soils have an obvious effect, can they be used to explain differences in peak values? We removed the effect of a general attenuation, which would include the effect of geometrical spreading and a general mid-crustal attenuation by applying the ground motion model (GMM) from Boore and Atkinson (2008) for each station knowing the distance and magnitude and assuming a strike slip focal mechanism (which they mostly are). Admittedly we use an older ground motion model but the emphasis here is in trends versus various site parameters and biases between models would not matter.
To show the effect of amplification we plot all residuals of peak velocities with respect to average amplification (1-6 Hz) to show that peak values generally increase at stations with more amplification (Fig. 5). There is a slight negative bias at sites with average amplification less than about 7, where as at higher levels of amplification there is a positive bias. This result is similar to that from Choi and Stewart (2005) and Borcherdt (1994), who showed amplification factors by site class. One characteristic of the figure is the apparent tendency for high scatter in the peak velocity values at low levels of amplification and clustered peaks at high values of amplification. Figure 5b shows the standard deviation of the peak values at each amplification level to quantify the low levels of scatter at high levels of amplification. Complicating this observation, however, is that there are less data at the higher levels of amplification. Figure 6a, b shows the peak acceleration with just a slight tendency for peak values to correlate with amplification, but the same tendency for high levels of amplification to correlate with low scatter in peak values. In conclusion, higher attenuation seems to be masked by higher amplification, typically at sites with low velocity soil layers such as at sites in the Imperial Valley. As t* is weakly correlated with average amplification in the 1-6 Hz band, there will be a correlation between t* and peak velocity. This suggests that the loss of scatter at high values of amplification seen in previous figures may be caused by the higher attenuation with the assumption that the scatter is associated with the high frequency content of the seismogram.
The amplitude residuals for all peak values are shown in Fig. 7 with different distance metrics.  Fig. 7a could also include radiation pattern effects such as directivity. Figure 7b shows the effect of using hypocentral distance, which in the case of the most negative value in Fig. 7a increased the distance from around 4.5 km to almost 12 km. Clearly the most negative values in Fig. 7b have been vastly reduced. Using very small epicentral distances for deep (11.5 km) earthquakes was not realistic. Figure 7a also had fairly symmetric distribution in residuals, which has become much more biased to large values in Fig. 7b, which means the epicentral-distance based prediction is under estimating the observed data, which are mostly at close distances. Figure 7c shows all the residuals for peak acceleration from differentiating the velocity signal. These values are much larger than for velocity. Similar to Fig. 7b the values are biased to positive values, but have few negative values. This analysis ignores the effect of topography, which would have higher relief and stronger gradients in elevation near the town of Anza close to the San Jacinto mountains compared to the Imperial Valley which is mostly flat. Spudich et al. (1996) studied the origins of a very high peak acceleration at the Tarzana site in Los Angeles, California from the Northridge earthquake (M = 6.7, Jan. 17, 1994). They deployed 21 stations around the accelerograph site to study the directional dependence of the amplification and found a strong directional resonance close to the accelerograph site but also that it died off within tens of meters. Neuberg and Pointer (2000) studied how particle motion changes on a theoretical side of volcano and then applied it to Montserrat volcano in the Caribbean. They find their  Oct. 18, 1989). They found the patterns of amplification to be very complicated and probably related to the 3-D structure of the hills, but a suggestion that ray paths parallel to the axis the hill produced higher ground motion. The analysis here will average out much of these effects because of the variation in azimuths at most stations. However, to the extent there is a bias in the range of azimuths relative to the dominant trend of the San Jacinto mountains, then some topographic effect may remain that would provide a bias between the stations near the town of Anza and the Imperial Valley.

Duration
Duration as well as amplitude is important when designing nuclear critical structures and is part of the design criteria for nuclear power plants in France (Hernandez and Cotton 2000). Seismograms from earthquakes at Anza show a wide range of impulsiveness from fairly drawn out S-wave/surface wave/coda wave trains to spikey single pulse arrivals as shown in Fig. 8. Vertical components tend to have much more scattered energy then the horizontal components possibly from the addition of conversions or scattering associated with SV. One Correlation of t* with amplification in the 1-6 Hz band. Both variables are determined from decomposing site from source. Apparent in the figure is a strong correlation between the two. This is expected from the model of sites characterized by low velocity sediments near the surface underlain by high velocity basement that will amplify in the 1-6 Hz range but will also attenuate because of the slow velocity and weak soils. Site DRE has a large amplification and is located in the Imperial Valley We computed the 5-95% time for the cumulative sum (or integral) of the squared velocity (shown in Fig. 9a). As v 2 is kinetic energy this procedure determines the evolution of energy in the seismogram. For example, Frankel and Wennerberg (1987) developed a model of energy flux to explain the evolution of coda with respect to direct arrivals. We also applied this to acceleration (Fig. 9b), but there is more physical meaning to v 2 . There is more scatter in the residual peak values at short duration and conversely, the residual peak values at long duration have little scatter. This suggests that the anomalous peak values all occur from the source and cannot be attributed to scattered energy (e.g. coda) or surface waves both of which would occur later in the record. Residual peak velocity and acceleration versus 2 measures of distance. All of the peak velocities are shown in a against epicentral distance to emulate Joyner-Boore distance. Replotted in b using hypocentral distance to more accurately measure the distance from the source of the high frequency radiation. Note the loss of most negative residuals. c The residuals for peak acceleration versus hypocentral distance. Note that the scale of acceleration is two orders of magnitude larger than for velocity 2764 J. B. Fletcher and J. Boatwright Pure Appl. Geophys.

Frequency Dependence of Q
A number of studies raise the question of whether Q, a measure of attenuation, is dependent on frequency or not. The most common form for attenuation (Kovach and Anderson 1964) which is Eq. (6) rewritten, where r is distance, c is velocity, T is period (1/f), and Q is called the quality factor and controls attenuation. This function reduces the amplitude a certain fraction per wavelength. Attenuation is dependent on Example of scattered (a) and b impulsive seismograms for the earthquake on 2016 08/31. While a has an impulsive P-wave, but more of an envelope for the S-wave, b has a much smaller coda implying more attenuation (scattered energy is more attenuated) or there is less scattering to build the coda. The red vertical bars show the window used to calculate the S-wave power spectra frequency if Q is frequency independent. Most early studies (Thatcher and Hanks 1973;Fletcher et al. 1987) use frequency-independent Q, compared to a number of more recent studies that use frequencydependent Q (Singh et al. 1982;Boatwright and Seekins 2011). Early lab tests suggested a frequencyindependent Q (Liu and Peselnick 1983), but raises the question of whether these experiments represent actual conditions in the earth at appropriate frequencies. We can address this question in two ways. The first is to perform the previous inversion using a frequency-dependent Q model and see if the fit is better. The second is to investigate the residual from the fit to a frequency-independent Q for any frequency-dependent effects. We assume the functional form for Q is where we vary n. Boatwright and Seekins (2011) used a value of 0.5 for n. Figure 10 shows the variance of the final iteration (measure of fit) using various values of n. Clearly within the model assumptions of the inverse procedure, the system does not support any frequency dependence as any frequency dependence added to the model yields a higher (much higher in some cases) final variance. The inversion program calculates the residual from the fit of model to the data . This residual is frequency-dependent and is the residual after subtracting the site responses and sources over the distance range. Figure 11 shows the amplitude of the residual plotted by distance and frequency. Apparent in the Figure is that there is a strong change near 90 km, which is close to the border of the San Jacinto Mountains/ Santa Rosa Mountains and the low-lying basin area near Borrego Springs. This change is from being more positive to more negative (less attenuation). The figure also shows that there is little frequency dependence in the contours. Ridges and valleys mostly go straight across the figure (the ordinate axis is frequency). Some ridges show some bending such as at a ridge from 80 to 65 km, but this has a low amplitude and does not represent the general trend.

Conclusions
Power spectra for 18 earthquakes (3 \ M \ 4) have been decomposed to determine source properties, regional path effects, site response, and site attenuation. Peak values were normalized using the GMM of Boore and Atkinson (2008). Site response (1-6 Hz average amplification) determined from the inverse procedure correlates with higher attenuation at most sites. Station amplification and attenuation mostly correlate with site geology at least at the more obvious locales such as those on hard rock (e.g., stations PFO, KNW, and RDM). Peak velocity correlates with the 1-6 Hz amplification, but not peak acceleration. Site attenuation correlates with site amplification as expected because low velocities that amplify also attenuate more than hard rock sites. Scatter in peak values is markedly reduced at sites with high amplification for both velocity and attenuation. Large peak values occurred at seismic envelopes with short durations but there was high scatter in peak values at these short durations. The dataset does not support frequency-dependent values of Q. are indebted to Lawrence Baker for developing the waveform server that allows us to analyze seismograms from both the Anza array and the Southern California Seismic Network easily. Internal reviews by Alex Grant and Alan Yong as well as journal reviews by Debi Kilb and Pieter-Ewald Share were very helpful and markedly improved the paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4. 0/.

Funding
Work funded by National Earthquake Hazards Reduction Program.