Site effects can significantly increase the seismic hazard and risk locally. Unconsolidated deposits such as thick and soft sediments in sedimentary basins are prone to strongly amplify the ground motion. Site effects are caused, among others, by the seismic impedance between rock and sediments, the 1D, 2D and 3D resonances, and the edge-generated surface waves. In turn, the site response can vary significantly from one site to another (site-to-site variability, e.g., Bindi et al., 2009; Hollender et al., 2015; Bindi et al., 2017; Imtiaz et al., 2018; Perron et al., 2018) and from one earthquake to another (within-site variability, e.g., Thompson et al., 2012; Ktenidou et al., 2016; Ktenidou et al., 2017; Maufroy et al., 2017; Perron, 2017; Zhu et al., 2018; Zhu et al., 2022). At large ground motion levels, non-linear effects in specific soils will increase the site response uncertainty as well (Régnier et al., 2013; Régnier et al., 2016). Understanding and reducing the ground motion estimation uncertainty is important for Probabilistic Seismic Hazard Assessment, especially at a long return period (Bommer and Abrahamson, 2006). The site-to-site and within-site variabilities have practical implications for site-specific study and microzonation, for instance on the spatial resolution and required duration of the instrumentation.

The within-site variability is very small when estimated from 1D SH site response analysis because it is a strong simplification of the real phenomena. On the contrary, approaches based on direct observations from real earthquake recordings are appropriate for analyzing the variability of the site response. One of the most commonly used approaches to measure the empirical amplification function is the standard spectral ratio (SSR) introduced by Borcherdt (1970). It consists in performing the ratio in the Fourier domain between the signal recorded at one station on sediments and the signal obtained at another station located nearby on a stiffer site condition (i.e., a rock site) for the same earthquake. However, in noisy urban areas in regions of low-to-moderate seismicity, recording earthquakes with a good signal-to-noise ratio (SNR) can require several months, if not years. It is thus important to estimate the number of earthquakes that should be recorded at the sites to evaluate the empirical site amplification function based on the desired accuracy.

The main goal of this work is to define, for the specification in the Swiss building SIA 261/1 (SIA, 2020), the minimum number of earthquakes in empirical site effect assessment. We first evaluate the stability and validity of the mean amplification as a function of the number of earthquake recordings used to compute it. The variations of the mean amplification are expected to be directly related to the within-site variability at each site. To verify that, we estimate the SSR and surface-to-borehole spectral ratio (SBSR) amplification function for stations of the Swiss strong motion network and of the Japanese KiKnet network having recorded hundreds of earthquakes. We use this large amount of data to determine the statistical distribution of the site amplification. Based on the statistical distribution, we propose an analytic equation predicting the variation of the mean amplification according to the standard deviation and to the number of recorded events. We also determined the dependence on the mean amplification functions of the ground motion intensity, measured as the peak ground acceleration (PGA).

In Switzerland, we developed a waveform database covering the time period from January 1998 to September 2019. Waveforms at each Swiss site were selected according to a magnitude–distance filter. In Japan, the database is covering the time period from October 1997 to March 2016. The SSR is computed for each component individually or the mean of the two horizontal components and can be noted as follows: S S R i ( f ) = F A S S i ( f ) F A S R i ( f ) , (1) where S S R i ( f ) is the SSR for the i t h component as a function of frequency f and F A S S i and F A S R i are respectively the Fourier amplitude spectra (FAS) at the site and at the reference computed over the i t h component. The SSR approach is based on the assumption that the earthquake source and wave propagation along the path are the same between the site and the reference and thus canceled out when performing the spectral ratio between the two. This assumption is valid if the site-to-reference distance (R_STA) is much smaller than the hypocentral distance (Rh). In practice, adopting Rh > 10R_STA is considered to be enough, even though a certain part of the SSR variability can probably be explained by a remaining influence of the source and of the path (Borcherdt, 1970; Perron, 2017). The ground motion amplification at the reference station is assumed to be negligible, that is to say, equal to one at every frequency. In practice, it is never the case (Hollender et al., 2017; Hollender et al., 2018; Hobiger et al., 2021), so the SSR-based amplification factors are not absolute but are always relative to the considered reference. One of the main limitations of the SSR is of having a rock outcropping susceptible to be used for the reference site located not too far from the considered site of interest. An alternative to the classical SSR is to deploy one station at the earth’s surface on sediments and the second at the same location but in a borehole deep enough to reach the geophysical bedrock. This so-called SBSR approach has the advantage of solving the between-station distance limitation but introduces some new difficulties because of the seismic wave reflection at the earth’s surface. The upgoing and downgoing waves are indeed fully constructive at the earth’s surface, although they can be destructive at certain frequencies at depth (Cadet et al., 2012). However, in the context of analyzing only the variability of the site response, the downgoing wave interaction can reasonably be neglected (Cadet et al., 2012; Hollender et al., 2018). We followed the same procedure for every computation of the site response in Switzerland and Japan. This procedure is as follows:

1) Automatic quality checks of earthquake recordings and automatic picking of the P and S wave arrival (T_P, T_S) through a time–frequency analysis;

Figure 1 shows an example of the SBSR computation in Japan.

In total, SSR is estimated from three pairs of stations where approximately 100 good-quality earthquakes have been recorded in Switzerland, and SBSR is computed from 60 pairs of surface-to-borehole stations with up to 2000 good-quality earthquakes in Japan. Figure 2 and Figure 3 respectively show the distribution of the SBSR for 60 pairs of surface-to-borehole stations in Japan and the SSR for the three pairs of surface stations in Switzerland. Figure 4 provides a summary of the number of good-quality earthquake recordings, geometric mean, and geometric standard deviation as a function of frequency in Japan (gray curves) and Switzerland (red curves).

Figure 2, Figure 3, and Figure 4 clearly show that the amplification functions are different from one site to another, both in terms of mean and standard deviation. It reflects the differences in the geological conditions of the sites, which determine, among others, the fundamental resonance frequency of the site (f₀), here corresponding to the first peak on the amplification function. The within-site standard deviation can also vary drastically from one site to another and depending on the frequency. In Japan, we can separate the amplification functions into two groups: the first group with f₀ > 0.5 Hz, with an amplification function equal to one and low standard deviation (close to 1.05) for frequency below f₀; a second group with f₀ below the minimum frequency of the analysis here (0.1 Hz), and having significant amplification (above one) and high variability at low frequency. It is also clear that the variability of the site response is on average higher in Switzerland than that in Japan. For the Swiss sites, this is probably because of the SSR method imposing relatively high site-to-reference distances and non-negligible site effects at the surface reference station. In Japan, we can observe some anomalies (eye shapes departing from the log-normal distribution) in the amplification function at high frequency (e.g., for stations: KiK-IBRH13; KiK-IBRH17; KiK-TCGH16). It is not possible to clearly determine its origin, but from our experience, this is very probably an artificial artifact because of coupling issues of the borehole instrumentation or because of a modification on the instrumentation at some point due to maintenance of the station for instance.

As we have seen in the previous section, both the mean and standard deviation of the amplification function as a function of frequency are dependent to the geological characteristics of the site itself. However, the nature of the site response distribution is the same independently to the site or to the frequency and has been shown to be well modeled by a log-normal distribution (Ktenidou et al., 2011). In other words, the distribution of the logarithm of the relative amplification of the ground motion between two sites is Gaussian. To qualitatively verify the log-normal distribution of the site response at every frequency, the quantile–quantile (Q–Q) plot and the histogram are represented at frequencies 0.5, 1.0, 2.5, 5.1, 9.9, and 20.6 Hz in Figure 5. The shape of the histograms of the logarithm of the amplification factors represents a Gaussian and Q–Q curves of every site at every frequency are well aligned along the 1/1 line, in particular in the interval ± 2 σ to the mean. These indicate that the site response is very well approximated by log-normal distribution at least up to ± 2 σ . Beyond 2 σ , the few non-natural outliers and the limited number of samples increase the scatter of the Q–Q curves, meaning that the log-normal distribution is still valid but interpretations made out of it are less reliable.

Proving the log-normal distribution of the amplification function is important because then peculiar statistical properties apply. For example, if a variable x is normally distributed then the distribution of sample means ( x n ¯ ) computed from subsets of n samples also are normally distributed. One major output of that is the confidence interval ( I c ). Given that a sample mean ( x n ¯ ) and unbiased standard deviation ( s n ) have been estimated from a finite number of samples ( n ), the confidence interval is the interval inside which the population mean ( μ ) for an infinite number of samples has a certain confidence level to be included in. It is defined as follows: I c 1 − α % = [ x n ¯ − Z α / 2 s n n ;     x n ¯ + Z α / 2 s n n ] , (2) where Z α / 2 is the critical value that defines the confidence level ( 1 − α ). For a normal distribution and a confidence level of 95%, Z 0.025 is equal to 1.96. However, because the number of samples can be sometimes very limited (i.e., only a few earthquakes have been recorded), it is preferable to use the Student distribution, also called t-distribution. This distribution correctly accounts for a small number of samples and tends to be a normal distribution as the number of samples increases. For a Student distribution, the formulation of I c 1 − α % is the same (Eq. 2), but the estimation of Z α / 2 is different, as it now also depends on n . The evolution of Z α / 2 , n as a function of n and for the confidence levels 68, 95, 99, and 99.9% is given in Figure 6, left panel. In the following, we will keep using the notation x n ¯ and s n for the measured sample geometric mean and standard deviation, whereas μ and σ represent the population geometric mean and standard deviation of the distribution. For an infinite number of samples, the two notations become equivalent: x ∞ ¯ = μ and s ∞ = σ . Moreover, we will only focus on the confidence level of 95%, because the 95% confidence interval corresponds approximately to the interval comprised between [ − 1.96 σ   1.96 σ ] , which in turn corresponds to the portion where the Q–Q plot best fit the 1/1 line (Figure 5). As the distribution is not normal but log-normal, we accordingly modified the confidence interval formulation. The 95% confidence interval for a log-Student distribution is finally: I c 95 % = [ x n ¯ ∗ 1 exp ( Z 0.025 ,   n l n ( s n ) n ) ;   x n ¯ ∗ exp ( Z 0.025 ,   n l n ( s n ) n ) ] , (3) with x n ¯ and s n respectively the sample geometric mean and standard deviation computed as x n ¯ = exp ( 1 n ∑ i = 1 n ln ( x i ) ) , (4) s n = exp ( 1 ( n − 1 ) ∑ i = 1 n ( ln ( x i ) − ln ( x ¯ n ) ) 2 ) , (5)

Figure 6 (right panel) shows the evolution of I c 68 % , I c 95 % , I c 99 % , and I c 99.9 % for a standard normal and standard Student distribution ( μ = 0 ;   σ = 1 ). It illustrates the very rapid reduction of the confidence interval as the number of samples increases, from more than 10 σ when n < 10 to less than 1 σ when n > 10 .

After demonstrating the validity of the log-normal assumption, we verified the validity of the prediction of I c 95 % for a Student distribution as a function of the number of earthquakes n by comparing I c 95 % with the observations in Switzerland and Japan. First, we defined two different confidence intervals: I c 95 N ( n ) = [ x N ¯ ∗ 1 exp ( Z 0.025 ,   N l n ( s N ) n ) ;   x N ¯ ∗ exp ( Z 0.025 ,   N l n ( s N ) n ) ] (6) I c 95 n ( n ) = [ x n ¯ ∗ 1 exp ( Z 0.025 ,   n l n ( s n ) n ) ;   x n ¯ ∗ exp ( Z 0.025 ,   n l n ( s n ) n ) ] (7) where x N ¯ and s N are respectively the total geometric mean and standard deviation computed over the entire dataset of N events. x n ¯ and s n are respectively the local geometric mean and standard deviation computed over a subset of n randomly selected events. I c 95 N ( n ) is the total 95% confidence interval used to predict the variation of any local mean x n ¯ computed from n events. I c 95 n ( n ) is the local 95% confidence interval used to predict the interval of variation of the total mean x N ¯ . This assumes that x N ¯ = μ and s N = σ , which is reasonably correct here since N is most of the time much higher than 100 earthquakes.

To estimate the reliability of the confidence interval more quantitatively, we bootstrapped the amplification factors at each frequency over 1000 random selections of n events, with n = [ 2   3   4   6   8   10   14   18   24   32 ] . We evaluated the proportion of local means included inside the total confidence interval ( P 1 =   x n ¯ ⊂ I c 95 N ( n ) ), and the proportion of total means included inside the local confidence interval ( P 2 =   x N ¯ ⊂ I c 95 n ( n ) ). Following Eq. 3, P 1 and P 2 can be written: P 1 ( f , n ) = 1 1000 ∑ k = 1 1000 ( x N ¯ ( f ) ∗ 1 exp ( Z 0.025 ,   N ln ( s N ( f ) ) n ) ≤ x n ¯ k ( f ) ≤ x N ¯ ( f ) ∗ exp ( Z 0.025 , N ln ( s N ( f ) ) n ) ) (8) P 2 ( f , n ) = 1 1000 ∑ k = 1 1000 ( x n ¯ k ( f ) ∗ 1 exp ( Z 0.025 , n ln ( s n k ( f ) ) n )     ≤     x N ¯ ( f )     ≤     x n ¯ k ( f ) ∗ exp ( Z 0.025 , n ln ( s n k ( f ) ) n ) ) (9) with the inequation equal to 1 when it is true and 0 otherwise. x n ¯ k and s n k are respectively the k t h local geometric mean and standard deviation computed over a subset of n   randomly selected events. If the distribution is perfectly normal, then both P 1 and P 2 are equal to 95%. However, we do not expect the site response distribution to be perfectly normal at every site and for all frequencies, so a certain convergence to 95% should be observed as the number of events n increases.

Figure 7 shows the bootstrap estimation of P 1 and P 2 from the amplification function of SIOO/SIOV in Switzerland. First, it is clear that the variability between the 1000 x n ¯ k decreases (blue points) as n increases (from top-left to bottom-right panel). This decay seems well predicted by I c 95 N (orange lines). This observation is also supported by P 1 which is relatively close to the value of 95% at all frequency and for any n > 2 . For n = 2 we can observe that P 1 is slightly higher than 95% between 10 and 20 Hz. In contrast, P 2 shows some significant low values for any n < 10 . However, P 2 shows a better agreement with the 95% value as n increases. This observation confirms the good approximation of using the log-normal distribution to model the site amplification variability. I c 95 % makes a relatively good prediction of the observed variability of x n ¯ , even when the number of samples is low.

Now, we follow the same procedure for every three SSR in Switzerland and SBSR in Japan. The corresponding results are given in Figure 8. We can make a similar observation as in Figure 7, P 1 is the average equal to the 95% value at all frequency and for every number of events. For the Swiss SSR, we can, however, observe a stronger scatter when the numbers of events are minimum ( n < 4 ). Again, we observe a stronger deviation of 95% in P 2 both in Switzerland and in Japan. In Switzerland, the discrepancy of P 2 is higher, especially close to 1 Hz and for n < 6 . P 2 is an average lower than 95% but tends to it as n increases. A good agreement is found for n > 6 and a complete stabilization is observed above 14 events. In Japan, we observed a different behavior, with P 2 being too low when n = 2 , and then too high when 2 < n < 8 mainly at low frequency (f < 2 Hz). For n > 8 , we observed a good stabilization of P 2 with mean values slightly below 95%.

The confidence interval computed from a large site response dataset is a good estimator of what is going to be the behavior of the mean computed from much smaller subsets of even only three earthquakes and for any frequency. However, it is clear that using 10 recordings of earthquakes or above greatly improves the quality of the prediction and the significance of the results. In conclusion, at least 10 events should be considered to have a good statistical significance and to make good use of the confidence interval predicting power.

Some questions which arise when evaluating the amplification function at a specific site are as follows: Is the number of earthquake recordings sufficient to accurately estimate the amplification function? Which minimal number of earthquakes ( n m i n ) should be used to evaluate the site response? Based on the confidence interval definition (Eq. 3), it is clear that the variability of x n ¯ depends both on s n and n . Because s n is site- and frequency-dependent (Figure 4), n m i n is by consequence also site- and frequency-dependent. In other words, there is no unique value of n m i n which can be considered for every site response analysis in the world. On the other hand, the property of the site response to be log-normally distributed can be supposed as universal. It is then possible to determine n m i n for any site response analysis, based on the log-normal distribution assumption and the use of the confidence interval definition.

Provided that the geometric mean x n ¯ and standard deviation s n of the site response has been measured at a particular site over a certain number of earthquakes n , it is possible to determine in which confidence interval the population mean for an infinite number of events μ has a certain confidence level (here 95%) to be included in. It is also possible to predict what will be the reduction of this interval if the number of earthquake observations increases. In the same way, it is possible to determine the number of earthquakes required to limit to a certain level the width of the interval where μ has a 95% probability to be found within. The width of the interval is independent to the x n ¯ and can be defined from Eq. 3 by C 95 % = exp ( Z 0.025 ,   n l n ( s n ) n ) , (10) C 95 % is the coefficient of variation between μ and x n ¯ such as x n ¯ C 95 % ≤ μ ≤ C 95 % x n ¯ with a 95% probability. It is now possible to estimate the minimum number of earthquakes required to limit the variation between μ and x n ¯ below a certain coefficient C 95 % as n m i n = ( Z 0.025 ,   n l n ( s n ) l n ( C 95 % ) ) 2 (11)

For example, if the amplification at 1 Hz has been measured from n = 10 earthquakes with a geometric standard deviation of s 10 = 1.5 , we can estimate the minimum number of earthquake n m i n to have C 95 % = 1.2 (20% of variation) with a probability of 95% as n m i n = ( Z 0.025 , 10 l n ( s 10 ) l n ( C 95 % ) ) 2 = ( 2.26 l n ( 1.50 ) l n ( 1.20 ) ) 2 = 25.31 → 26   e a r t h q u a k e s

It is important to note that for a Student distribution, Z 0.025 ,   n is the function of n . Z 0.025 ,   n will decrease very rapidly as the number of measured earthquakes increases (Figure 6). Using Eq. 11 and measured s n (Figure 4), n m i n is computed for every site in Switzerland and Japan, and at every frequency. The results are reported in Figure 9. As already discussed, n m i n is dependent on s n , so it is variable for the different sites and frequency. Swiss SSRs have the highest uncertainty and logically required the highest number of earthquakes for a given coefficient of variation C 95 % . Table 1 summarizes the minimum number of earthquakes which is valid for 99, 95, and 84% of our sites and frequencies as a function of C 95 % . For 10 earthquakes recorded, the estimation of the mean is only 40% accurate approximately ( C 95 % = 1.4 ). It is possible to reduce this uncertainty to 25% by recording 20 events ( C 95 % = 1.25 ). Depending on the desired limit for the coefficient of variation of the mean, one can make own estimations of the minimum number of earthquakes using Eq. 7.

It has to be highlighted that s n is the key parameter for the estimation of n min . If s n is wrongly determined, so will be n m i n . One difficulty to have a representative determination of s n is how to deal with the outliers. Including erratic outliers will artificially increase s n , while removing natural outliers from rare events will truncate the true distribution and reduce s n . Another difficulty is that looking only at the value of n m i n might not be enough for all sites. One could claim that because the site response has been measured from 30 earthquakes, the statistical significance of the result is good and the coefficient of variation of the mean is low. However, if all the events present the same characteristic and location because they belong to the same cluster of events, then the significance of the results is not good and the true variability of the site response might be strongly underestimated. For instance, Perron (2017) showed that approximately 50% of the within-site variability in 2D and 3D basins comes from the lighting effect, which strongly depends on the source location. This implies that both the number of events and their spatial distribution around the site should be considered in site response analysis.

The dependence of the site response on the intensity of the ground motion is a complex research topic that interests the community for several decades (e.g., Sánchez-sesma, 1987; Aki, 1993). The non-linear behavior of unconsolidated soil to strong ground motion solicitations is of major interest in engineering seismology. Non-linearity tends to reduce the fundamental resonance frequency of the site, leading to an increase of the hazard at low frequency and a decrease at high frequency (Régnier et al., 2016). In extreme cases, it can also lead to liquefaction phenomena.

One question often arises when speaking about empirical site effects assessment which is: is the measured amplification function from weak ground motion representative of site response to strong ground motion? To address this question, we compute the equivalent of the standard normal distribution ( μ = 0 ,   σ = 1 ) for every individual amplification function at all sites in Switzerland and Japan as Z i = l n ( x i ) − l n ( x n ¯ ) l n ( s n ) , (12)

This common standard normal distribution formulation allows using the site response of every site together. Z i ( f ) represent the ith normalized amplification function normally distributed with x n ¯ = 0 and s n = 1 . Together, it represents about 28,000 normalized amplification functions obtained from thousands of earthquakes recorded at 63 pairs of stations (three Swiss sites and 60 Japanese sites). For each normalized amplification function, we computed on the corresponding waveforms the horizontal mean PGA.

Figure 10 shows the number of events per frequency, the distribution of the PGA and the normalized amplification function for four PGA bins [(0.001 0.01), (0.01 0.1), (0.1 1), and (1 10) m/s²]. First, it should be mentioned that the number of events varies strongly from one PGA bin to another. This explains the apparent differences when looking at the normalized amplification function (black curves) of the different bin. We observe that the normalized amplification function for every PGA bin can be explained by the standard normal distribution, which indicates that no non-linear behavior is observed here. The mean is fairly equal to 0 and the standard deviation is equal to 1 for every frequency of every bin. That demonstrates, first, that the linear behavior characterizes the vast majority of the sites, and second, that the linear site response is independent to the ground motion intensity. Therefore, if we consider a specific site having a linear behavior, the amplification function observed from the weak motion of a small magnitude earthquake will be the same as the one for the strong motion of a large magnitude earthquake, all other things being the same. This highlights the importance and the validity of using the recording of low-to-moderate earthquakes to assess the anelastic amplification functions for larger earthquakes as long as there is no significant non-linear site response at the site of interest.

Site effect is a major contributor to the seismic hazard, and its evaluation at specific sites of interest generally requires the recording of several earthquakes. We address here the question of the site response variability and of the minimum necessary number of earthquakes to be recorded.

To address this question, we carefully compute empirical amplification functions at 60 KiKnet sites from several hundred earthquakes and three Swiss sites from several tens of earthquakes. We performed statistical analysis on the amplification function to estimate the geometric mean and standard deviation, and more importantly to determine the distribution law of the amplification factor at each frequency. Independent to the site and to the frequency, we found that the log-normal distribution is a very good approximation for the site response. Based on that we developed a strategy to estimate the minimum number of earthquakes from the confidence interval definition. We first demonstrate the validity of the use of the confidence interval to model the uncertainty of the geometric mean estimation. We found that between 8 and 14 earthquakes are necessary to have a good prediction by the confidence interval, that is to say, a good statistical significance. For most of the sites, 10 samples seem to be the best compromise between minimizing the number of recordings and having a good statistical significance of the results. Based on the confidence interval, we provide the analytic formula to estimate the minimum number of earthquakes to be recorded, as a function of the within-site standard deviation (Eq. 11). We used it on the Swiss and Japanese amplification function and determine, among others, that with a 95% probability: the mean varies by less than 40% for 10 earthquakes, and less than 25% for 20 events.

It is very important to point out that satisfying the minimal number of earthquakes by itself is not sufficient. The selected earthquakes should be uncorrelated and as much evenly distributed around the site as possible to cover the entire variability of the site response. Therefore, one should not use only earthquakes belonging to a single cluster of events. In our dataset, the linear site response is observed to be independent to the intensity of the ground motion. In other words, assessing the site response from the recording of low PGA and low magnitude earthquakes, provides the same amplification functions as from recording of high PGA and large magnitude earthquakes, as far as the soil behaves linearly.

As a general rule, a minimum of 10 uncorrelated earthquakes should be considered, but the higher the number of earthquakes, the lower the uncertainty on the geometric mean site response assessment. Based on our results, the specification in the Swiss building SIA 261/1 recommends taking a minimum of 10 uncorrelated earthquakes to perform site-specific studies.