# Effects of Nonlinear Soil Behavior on Kappa ($\kappa$): Observations from the KiK-Net Database

Chunyang Ji, Ashly Cabas, Luis Fabian Bonilla, Céline Gelis

► To cite this version:

Chunyang Ji, Ashly Cabas, Luis Fabian Bonilla, Céline Gelis. Effects of Nonlinear Soil Behavior on Kappa ($\kappa$): Observations from the KiK-Net Database. Bulletin of the Seismological Society of America, 2021, 111 (4), pp.2138 - 2157. 10.1785/0120200286. hal-04402168v2

HAL Id: hal-04402168

https://univ-eiffel.hal.science/hal-04402168v2

Submitted on 30 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
1 **Effects of Nonlinear Soil Behavior on Kappa (κ): Observations from the KiK-net**

2 **database**

3 Ji, Chunyang, Cabas, Ashly, Bonilla, Luis Fabian, and Gelis, Céline

4

5 **Declaration of competing interests**

6 The authors acknowledge there are no conflicts of interest recorded.

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21
Soil nonlinear behavior is often triggered in soft sedimentary deposits subjected to strong ground shaking and has led to catastrophic damage to civil infrastructure in many past earthquakes.

Nonlinear behavior in soils is associated with large shear strains, increased material damping ratio and reduced stiffness. However, most investigations of the high-frequency spectral decay parameter $\kappa$, which captures attenuation, have focused on low-intensity ground motions inducing only small shear strains. Because studies of the applicability of the $\kappa$-model when larger

deformations are induced are limited, this paper investigates the behavior of $\kappa$ (both $\kappa_r$ per record and site-specific $\kappa_0$ estimates) beyond the linear-elastic regime. Twenty stations from the KiK-net database, with time-average shear wave velocities in the upper 30 m between 213 and 626 m/s, are

used in this study. We find that the classification scheme used to identify ground motions that trigger soil nonlinear behavior biases estimates of $\kappa_0$ in the linear and nonlinear regimes. A hybrid method to overcome such bias is proposed considering proxies for in situ deformation (via the shear strain index) and ground shaking intensity (via peak ground acceleration). Our findings show

that soil nonlinearity affects $\kappa_r$ and $\kappa_0$ estimates, but this influence is station-dependent. Most $\kappa_0$ at our sites had a 5-20% increase at the onset of soil nonlinear behavior. Velocity gradients and

impedance contrasts influence the degree of soil nonlinearity and its effects on $\kappa_r$ and $\kappa_0$. Moreover, we observe that other complexities in the wave propagation phenomenon (e.g., scattering and amplifications in the high-frequency range) impose challenges to the application of the $\kappa_0$-model,

including the estimation of negative values of $\kappa_r$.
43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

The anelastic attenuation of seismic waves as they travel through sedimentary deposits is a function of the deformations induced, which in turn depend on the material properties (e.g., plasticity of the soil) and the intensity of the ground shaking. Material damping ratio, ξ, is commonly used in geotechnical earthquake engineering to quantify viscous and anelastic energy dissipation in soils subjected to dynamic loading. Empirical models of ξ often have a constant minimum value (known as the minimum shear-strain damping, ξmin) for small shear strains considered in the linear-elastic regime (e.g., Darendeli, 2001). Yet, values of ξ increase as larger shear strains are induced in soil deposits by stronger ground excitations (Idriss et al. 1978, Seed et al. 1986, Darendeli 2001). The characterization of ξ across a wide range of strains is essential to model the effects of local soil conditions on earthquake ground motions.

The high-frequency spectral decay parameter κ was introduced by Anderson and Hough (1984) based on the Fourier spectrum characteristics of a ground motion's shear-wave window recorded directly in the field, which makes it an observable parameter that quantifies total attenuation (e.g., energy dissipation caused by scattering and anelasticity). Estimates of κ have proven useful in multiple applications, from stochastic modeling of ground motions (Boore 2003) to the development of host-to-target adjustments of ground motion models (e.g., Campbell 2003, Al Atik et al. 2014). The site-specific, distance-independent component of κ, known as κ₀ (Anderson 1991), is defined as a site parameter that captures the attenuation due to the propagation of seismic waves through near-surface materials. The relationship between κ₀ and ξmin has been investigated in previous studies (e.g., Cabas et al. 2017; Ktenidou et al. 2015; Xu et al. 2019) for weak motion data, but the quantification of κ and κ₀ beyond the linear-elastic regime remains unsolved.
66
67 Most studies on individual estimates of κ (i.e., the value measured from the observed Fourier
68
69 Amplitude Spectra (FAS) per record following the Anderson and Hough (1984) approach),
70 hereafter referred to as κr, and its site-specific component κ₀ have used ground motion records that
71 do not trigger nonlinear behavior at the sites of interest or that are not considered influenced by
72 the site's nonlinear response (e.g., Ktenidou et al. 2013, Van Houtte et al. 2011, Laurendeau et al.
73 2013, Edwards et al. 2015, Perron et al. 2017, Xu et al. 2019). However, nonlinear soil behavior
74 has often been responsible for increasing the damage potential of strong ground motions in past
75 earthquakes (e.g., Darragh and Shakal, 1991, Trifunac and Todorovska, 1994, Bonilla et al., 2011,
76 Rong et al., 2016). Understanding near-surface attenuation effects in the nonlinear regime is then
77 necessary for a thorough assessment of seismic hazards and risks imposed to civil infrastructure
78 (Anderson and Hough 1984). Hence, this paper investigates the relationships among κ, shear
79 strains and ground motion intensity to understand the behavior of κ at the onset of nonlinear soil
80 behavior.

81

82
83
84
85
86
87
88
The first paper that attempted to connect soils' nonlinear response and Kr was conducted by Yu et
al. (1992), where the authors compared two simulated records: one from a linear site response
analysis, and the other from a time-domain nonlinear site response analysis. Yu et al. (1992)
observed that the value of Kr estimated with the Anderson and Hough (1984) approach and
corresponding to the motion affected by soil nonlinearity was smaller. However, later studies
found a positive correlation between Kr (and K0), strain amplitudes and the intensity of ground
shaking (e.g., Durward et al 1996; Lacave-Lachet et al 2000; Dimitriu et al. 2001).
89

90

91

92

93

94

95

96

97

98

99

100

101

102

103

104

105

106

107

108

109

110

111

Durward et al. (1996) found that $κ_r$ values were a function of peak ground velocity (PGV, varying from 0.01 to 1 m/s), which was used as a proxy for deformation. Values of $κ_r$ were computed for more than 60 records observed at 23 sites in the Imperial Valley, California based on the acceleration spectrum approach (Anderson and Hough, 1984). Durward et al. (1996) hypothesized that soil nonlinearity had affected $κ_r$ because higher $κ_r$ values correlated well with higher PGVs.

Moreover, Lacave-Lachet et al. (2000) analyzed ground motions from the 1995 Kobe earthquake in Japan (i.e., the main shock and aftershocks), and found that $κ_r$ values increased with increasing peak ground acceleration (PGA). Hence, Lacave-Lachet et al. (2000) proposed to use $κ_r$ to detect the onset of soil nonlinearity. Dimitriu et al. (2001) investigated the dependency between site-specific $κ_0$ and ground shaking intensity. Values of $κ_r$ for 23 ground motions (i.e., 46 horizontal components with values of $κ_r$ reported for each individual component) were computed at two adjacent sites in Lefkas, western Greece, based on the acceleration spectrum method. Dimitriu et al. (2001) provided evidence that $κ_0$ was a proxy for soil nonlinearity based on the observed dependency between $κ_0$ and ground shaking amplitudes, which were represented by mean horizontal acceleration in the S-wave window (MGA), PGA, and PGV. Positive correlations were found between the 46 $κ_0$ values and MGA, PGA, and PGV in log-scale, while a negative correlation was observed between $κ_0$ and the site dominant-resonance frequency. However, Van Houtte et al. (2014) identified an opposite correlation between estimates of $κ_0$ (computed as the individual measured $κ_r$ with epicentral distance less than 30 m) and PGA at hard sites (i.e., with $V_{s30}$ varying from 422 to 1073 m/s) using ground motions from the 2010-2011 Canterbury earthquake sequence in New Zealand. The authors suggested further investigations to understand the associated physical mechanism.
There are still few and contradicting observations of the effects of nonlinearity on $κ_r$ and $κ_0$ estimates (Ktenidou et al. 2015). Previous studies only considered a very limited database of ground motions. This paper takes advantage of the unique Japanese Kiban-Kyoshin network (KiK-net), which is rich in high-quality ground motions, to further investigate the effects of soil nonlinearity on $κ_r$ and $κ_0$. More specifically, we explore the dependence of $κ_r$ and $κ_0$ on ground shaking intensity (i.e., weak, moderate or strong ground motions as parameterized by PGA), and on the level of shear strains induced in near-surface materials at 20 KiK-net stations. First, we present the conceptual basis for the relationship between $κ$, shear strains and ground motion intensity. Then, we describe our database and methods, starting with the identification of an appropriate classification scheme for linear and nonlinear ground motions. The analysis of the effects of ground shaking amplitudes on $κ_r$ at each study site follows. Lastly, we compare the ratio of nonlinear and linear site-specific $κ_0$ across all selected stations to assess the variation of near-surface attenuation estimates from the linear-elastic to the nonlinear regime.

Conceptual basis for the interpretation of $κ$ beyond the linear-elastic regime

The induced strain level in a given soil layer is a function of the material properties, and the amplitude and frequency content of the incoming wavefield at the site. Stronger ground shaking results in larger-strain responses, which produce an increase in material damping ratio (in combination with a reduction of shear modulus). The short wavelength of high frequency waves allows for multiple cycles of shearing in near-surface sedimentary layers, which makes them more sensitive to the effects of a higher material damping ratio. Thus, we hypothesize that stronger
134 ground shaking inducing larger deformations in sedimentary deposits will affect estimates of the
135 high-frequency spectral decay parameter κ.

136
137 Figure 1 serves as an example to illustrate this hypothesis. The acceleration FAS and empirical
138 transfer function (ETF) corresponding to a pair of ground motions recorded at depth and at the
139 ground surface at FKSH14 are shown. FKSH14 is one of the 20 KiK-net station analyzed in this
140 work (see Database description). One of the ground motion pairs has a low ground shaking
141 intensity (with a surface PGA of less than 0.1 m/s² for both horizontal components), while the
142 other one has a higher ground shaking intensity (with a surface PGA of 0.71 m/s² and 0.25 m/s²
143 for the H₁ and H₂ components, respectively). Values of κr per record are also estimated for both
144 pairs. To minimize the bias from path effects and isolate local site effects on κr values per record,
145 the selected weak and strong ground motions in Figure 1 correspond to events with similar focal
146 depths, azimuths and epicentral distances (i.e., the focal depths, azimuths and epicentral distances
147 are 5 km, 319°, and 15.04 km, respectively for the low-intensity event, and 5 km, 323°, and 15.79
148 km, respectively for the high-intensity event). The moment magnitudes of the events associated
149 with the weak and strong ground motions in Figure 1 are 4.0 and 5.1, respectively. To reduce the
150 variability associated with the calculation of κr from the acceleration spectrum method, we use the
151 same frequency window for all ground motions (Edwards et al. 2011).

152
153 Larger values of κr at the ground surface (e.g., 0.071 s for the horizontal component H₁) are
154 obtained for the high-intensity ground motion compared to the corresponding κr values for the
155 low-intensity ground motion (e.g., 0.057 s for H₁). Likewise, the high-intensity ground motion
156 results in a larger κTF (Drouet et al. 2010), which is computed on the decaying portion of the
empirical transfer function at high frequencies and is equivalent to Δκ (i.e., κ<sub>r_sur</sub> - κ<sub>r_bor</sub>) estimates when using the same frequency band for calculation purposes. The value of κ<sub>TF</sub> corresponding to the high-intensity and low-intensity motions are 0.026 s and 0.019 s, respectively. Meanwhile, the ETF corresponding to the high-intensity event shows lower amplifications at higher frequencies (e.g., amplification ranges from 1 to 2 between 15 and 30 Hz for H<sub>1</sub> approximately) than its counterpart for the low-intensity ground motion (i.e., amplification ranges from 2 to 3 between 15 and 30 Hz for H<sub>1</sub> approximately), which reflects the stronger influence of increased material damping ratio on high frequencies.

Values of κ<sub>r</sub> and Δκ (or κ<sub>TF</sub> in Figure 1) are identical for the high- and low-intensity records for the H<sub>2</sub> component. It must be noted that the difference in PGA between the low- and high-intensity surface ground motions is less significant for the H<sub>2</sub> components than that for the H<sub>1</sub> components, while the PGAs at depth are very similar (i.e., the borehole PGAs for the high-intensity motion are 0.06 and 0.07 m/s<sup>2</sup> for H<sub>1</sub> and H<sub>2</sub> components, respectively; the borehole PGAs for the low-intensity motion are 0.01 m/s<sup>2</sup> for both H<sub>1</sub> and H<sub>2</sub> components). This may explain why Δκ values are the same for the low and high intensity motions in the H<sub>2</sub> component direction. Additionally, this observation also hints that the near-surface attenuation and site effects may be affected by the ground motion directionality as shown by Ji et al. (2020). Finally, soil nonlinearity is commonly observed at shallower layers (Régnier et al. 2013), and their effects on κ<sub>r</sub> at borehole are smaller than at surface. Thus, the κ<sub>r</sub> at borehole should be less affected by soil nonlinearity because it nonlinear soil behavior is less likely to be triggered at depth (e.g., the borehole V<sub>s</sub> at FKSH14 is 1210 m/s with borehole sensor depth of 147 m). In Figure 1, values of κ<sub>r</sub> at borehole are similar for the high- and low-intensity motions. The changes in κ<sub>r</sub> shown in Figure 1 exemplify the need
180 to further investigate the influence of the onset of soil nonlinearity on κ<sub>r</sub> as well as the potential
181 implications on κ<sub>0</sub> values at a specific site.

182
183
Database description

184 In this study, we use ground motions from the KiK-net database, which provides high quality
185 strong ground motions recorded at more than 600 stations installed uniformly across Japan. Each
186 station possesses a pair of sensors, one at the surface and another one at depth that is typically
187 between 100 to 200 m deep. The sampling frequency of the observed acceleration series is either
188 100 or 200 Hz. The P- and S-wave velocity profiles are measured by downhole PS logging and
189 available at the KiK-net website (see Resources and Data). The earthquake information, including
190 the seismic moment magnitude M<sub>w</sub>, focal depth and epicenter location are provided by the
191 broadband seismography network (F-net) catalog. The dataset used in this work was processed by
192 Bahrampouri et al. (2020) with an automated protocol, which corrects the baseline and removes
193 the background noise with a bandpass/high-pass acausal filter. The low-frequency filter corner
194 frequency is determined based on a required value of the (signal+noise)/noise equal to 3,
195 corresponding to a signal-to-noise ratio (SNR) of 2.0 (Boore and Bommer, 2005); this frequency
196 ranges between predetermined minimum and maximum values of 0.05 and 0.5 Hz, respectively.
197 The minimum high-frequency filter corner frequency is determined for a SNR of 1.0 (Douglas and
198 Boore 2011). The minimum bandpass width is 60% of the range defined from zero to the Nyquist
199 frequency. Further detailed descriptions on the ground motion processing are available in
200 Bahrampouri et al. (2020).

201
We use surface and borehole records (only horizontal components) in this paper. The criteria applied to select records are as follows: (1) epicentral distance is less than 150 km, (2) the SNR ratio is larger than 3.0 at each frequency from 1.0 to 30 Hz, (3) focal depth is less than 35 km (Ji et al. 2020), and (4) the seismic wave path does not cross the Japanese volcanic belt (Nakano et al. 2015). Thus, twenty stations with more than five nonlinear ground motions (the definition of nonlinear ground motions is described next in the Methods section) are used in this work (with 18 stations having more than 10 nonlinear records and 8 stations having more than 15 nonlinear records, see Table 1). Table 1 provides local soil conditions and the number of ground motions at each selected site. The locations of selected stations are shown in Figure 2a, while the magnitude and distance distribution of selected ground motions is provided in Figure 2b.

## Method

### Identification of the onset of nonlinearity

Identifying ground motions that trigger soil nonlinear behavior is key to evaluate $\kappa$ estimates beyond the linear-elastic regime empirically. The shear-strain index ($I_y$ = PGV/Vs), which is a proxy for the in-situ deformation, and PGA, which describes the peak amplitude of the ground motion, are commonly used to differentiate linear from nonlinear ground motions (e.g., Xu et al. 2019, Cabas et al. 2017, Wang et al. 2019). Moreover, the correlation between PGA and $I_y$ has been shown to be an effective proxy to capture in-situ stress-strain relationships. This correlation has been characterized via the classic hyperbolic model, which fits empirical observations (Chandra et al. 2014, 2015, Guéguen et al. 2018). However, there is lack of consensus regarding the sufficiency and efficiency of existing proxies associated with the onset of nonlinear behavior. For example, Xu et al. (2019) assumed that records with $I_y$ less than 0.01% are linear ground
motions, while Cabas et al. (2017) adopted 0.1%. On the other hand, Ktenidou et al. (2013) chose a PGA of 0.1 m/s² as the threshold for linear ground motions. Régnier et al. (2013) conducted a statistical analysis on the KiK-net dataset to understand nonlinear site response at their stations. They characterized linear soil behavior as that associated with motions with a PGA at depth less than 0.1 m/s².

We develop appropriate criteria to identify nonlinear ground motions based on PGA and $I_γ$. In this paper, the shear-strain index at the surface ($I_{γ,0}$) is defined as follows:

$$I_{γ,0} = PGV_{rotD50} / V_{s,0} \qquad (1)$$

where $V_{s,0}$ is the shear-wave velocity at the ground surface, and $PGV_{rotD50}$ is the median PGV for all rotated surface ground motions following the approach of Boore (2010). The use of $PGV_{rotD50}$ rather than the PGV from recorded ground motion horizontal components can minimize directionality effects. The use of $I_{γ,0}$ (Equation 1) as an indicator of the onset of soil nonlinearity has some limitations though. First, the selection of single $V_s$ and PGV values to capture the depth-dependent deformation in the profile may underestimate or overestimate the level of nonlinearity experienced across the whole column. Thus, $I_{γ,0}$, as defined in Equation 1, serves simply as a proxy for a representative shear strain to take place in the sedimentary deposit of interest. Notably, there is no consensus regarding the most appropriate choice of $V_s$ for $I_γ$ estimates. For example, $V_{s30}$ is a commonly used site proxy and hence often selected to compute $I_γ$ (e.g., Kim et al. 2016, Guéguen et al. 2018). The equivalent $V_s$ measured between two successive sensors with seismic interferometry by deconvolution has also been adopted for $I_γ$ estimates (e.g., Chandra et al. 2015, 2016, Wang et al. 2019). Second, values of $I_{γ,0}$ are not directly comparable to shear strains measured in the laboratory, not only because $I_γ$ is a proxy and not a measured value, but also
because the dissipation of seismic energy as captured in the laboratory may not fully represent the attenuation mechanisms taking place in the field (Cabas et al. 2017).

By applying the classic hyperbolic model to describe the correlation between PGA<sub>rotD50</sub> and I<sub>γ,0</sub>, we find that there is no unique threshold to identify nonlinear ground motions across all study sites. Figure 3 provides examples of theoretical hyperbolic fitting curves at four stations (i.e., IBRH16, IBRH17, IBRH20 and IWTH21) with varying V<sub>s30</sub> values (from 244 to 626 m/s) to demonstrate the limitations associated with using a single parameter to identify nonlinear ground motions at multiple sites. Scatter points represent PGA<sub>rotD50</sub> and I<sub>γ,0</sub> pairs from recorded ground motions at the sites of interest, while the lines correspond to the fitting curves from the hyperbolic model. It is observed that the same deformation at various sites would be triggered by different levels of ground shaking (e.g., I<sub>γ,0</sub> of 0.05% will be caused by a PGA<sub>rotD50</sub> around 1 m/s<sup>2</sup> at a NEHRP D site, such as IBRH20 with V<sub>s30</sub> of 244 m/s, and by a PGA<sub>rotD50</sub> of 2 m/s<sup>2</sup> at a NEHRP C site, such as IBRH16 with V<sub>s30</sub> of 626 m/s). On the other hand, if nonlinearity is assumed to be triggered when the PGA<sub>rotD50</sub> is larger than a predetermined threshold, different levels of I<sub>γ,0</sub> will be associated with the onset of soil nonlinearity. Hence, in this work, we propose a hybrid method (further described in the next section) based on both, the intensity of the excitation and in-situ deformation to classify ground motions.

**Linear, transitional, and nonlinear datasets**

Surface and borehole ground motions are considered separately herein. Régnier et al. (2013) considered that records with borehole PGA less than 0.1 m/s² cannot trigger nonlinear site response at the ground surface. They also observed that soil nonlinearity is mainly triggered at superficial
layers. In this study, there are only 3% records with borehole PGA larger than 0.1 m/s². Therefore the borehole records are assumed to remain in the linear-elastic regime (i.e., they do not trigger nonlinear behavior at depth). Surface records are separated into three sub-datasets, namely linear, transitional (i.e., soil's behavior is between the linear-elastic and nonlinear regimes), and nonlinear ground motions. First, we define a threshold based on Iγ,₀ to differentiate linear from transitional records, which is hereafter referred to as Iγ,₀,t. Likewise, a transitional threshold, Iγ,₀,t, is defined to separate transitional and nonlinear ground motions. The linear Iγ,₀,t threshold is defined as the onset of soil nonlinearity by visual inspections of the corresponding PGArotD50 versus Iγ,₀ curve, and corresponds to the point where PGArotD50 values begin to increase at a higher rate with increasing Iγ,₀. The transitional Iγ,₀,t threshold captures when the soil nonlinearity becomes more apparent, which corresponds to the point where the second change in slope of the PGArotD50 versus Iγ,₀ curve takes place. Figure 4 provides an example of the selection of Iγ,₀,t and Iγ,₀,l at station MYGH10. The threshold separating the linear and transitional ground motions is Iγ,₀,l = 0.001%, while the threshold separating transitional and nonlinear ground motions is Iγ,₀,t = 0.007%.

A maximum PGA<sub>rotD50</sub> of 0.25 m/s<sup>2</sup>, which is the value adopted by Régnier et al. (2016) to define low-amplitude motions, is chosen as an additional constraint to avoid linear ground motions being erroneously included into the nonlinear dataset. Thus, linear, transitional, and nonlinear datasets are defined as follows:

290 ○ Linear ground motions: records with I<sub>γ,0</sub> less than the I<sub>γ,0,1</sub> threshold.

291 ○ Nonlinear ground motions: records with (a) I<sub>γ,0</sub> larger than the I<sub>γ,0,t</sub> threshold and (b)
292 PGA<sub>rotD50</sub> larger than 0.25 m/s<sup>2</sup>.

293 ○ Transitional ground motions: records not classified as either linear or nonlinear.
The validity of the proposed linear, transitional, and nonlinear datasets is tested by examining the behavior of the shear modulus, G against $I_{\gamma,0}$ at the study sites. The reduction of G for empirical ground motions is estimated as follows (after Guéguen et al., 2019):

$$ \frac{G}{G_{\max}} = \frac{PGA_{rotD50}}{PGV_{rotD50}/V_{s,0}} \left( \frac{PGA_{rotD50}}{PGV_{rotD50}/V_{s,0}} \right)_{\max} \quad (2) $$

The term $\left( \frac{PGA_{rotD50}}{PGV_{rotD50}/V_{s,0}} \right)_{\max}$ is computed from the corresponding average ratio of records with $I_{\gamma,0}$ less than 0.001%, which is the predetermined threshold of $I_{\gamma,0,1}$ for the linear-elastic deformation limit in this work. Figure 5 shows the $G/G_{\max}$ versus $I_{\gamma,0}$ curves at all study sites. Even though Figure 5 does not directly correspond to laboratory-based curves, it serves as a first order verification of the distinct behavior corresponding to the linear, transitional, and nonlinear ground motions identified with the categorization scheme proposed herein. One challenge when interpreting these data is the characterization of the soil volume being sampled when using $I_{\gamma}$, which is related to the frequency band that PGV is acting on. Identified linear ground motions mainly have $G/G_{\max}$ values around 1 ($G/G_{\max}$ values higher than 1 result from using mean G values as a proxy for $G_{\max}$), while the ratios corresponding to the nonlinear dataset are generally less than 1 due to the onset of soil nonlinearity. Notably, values of $G/G_{\max}$ associated with the transitional dataset are between the linear and nonlinear datasets. It is not clear whether the site response associated with records identified as transitional could be equivalent to a linear-elastic or a nonlinear response, because associated $G/G_{\max}$ values vary within a single station and across stations. Hence, the characterization as transitional is deemed appropriate herein.
314

315 $\kappa_{r\_AS}$ estimates

316 We use the acceleration spectrum approach (Anderson and Hough, 1984) to estimate $\kappa_{r\_AS}$. To
317 minimize the variability introduced by the selection of the S-wave window, the entire time series
318 is used. Because compatibility with engineering analysis such as geotechnical site response
319 analysis and ground motion models is also desired, and such applications use complete time series,
320 calculating $\kappa_{r\_AS}$ using the entire time series is convenient. Moreover, the differences of $\kappa_{r\_AS}$
321 values measured from S-wave window and the entire time series are not significant in most cases
322 (Ji et al. 2020). In this study, a subset of randomly selected ground motions is used to further assess
323 potential discrepancies between using the S-wave window relative to the entire time series, and no
324 significant differences are observed. Hence, the entire time histories are used.
325

325

326 The variability in estimates of $\kappa_{r\_AS}$ is a function of the selection of the frequency band (Edwards
327 et al. 2015, Perron et al. 2017) among other factors (Ji et al. 2020, Ktenidou et al. 2013). Moreover,
328 soil nonlinearity affects low and high frequencies differently. The onset of nonlinear soil behavior
329 can influence high frequencies first (Bonilla et al. 2011, Bonilla and Ben-Zion, 2020) because
330 larger shear strains are induced in softer, thinner layers located at shallower depths (i.e., in a profile
331 with increasing stiffness with depth). Hence, we compute $\kappa_{r\_AS}$ based on a pre-determined fixed
332 frequency window ($[f_1, f_2]$). The pre-determined $f_1$ corresponds to the maximum value between 1.4
333 $f_c$ (where $f_c$ is the earthquake source corner frequency of each record) and $1.4 f_0$ (where $f_0$ is the
334 site's predominant frequency). Indeed, if $f_1$ is lower than $f_0$, the value of $\kappa_{r\_AS}$ will be biased by the
335 site amplification in the high-frequency range (Parolai and Bindi 2004). On the other hand, the $f_c$
336 requirement is added to reduce the effects of the earthquake source. The value of $f_2$ is set to be 25
337 Hz due to consideration of KiK-net instrument's response (Aoi et al., 2004, Fujiwara et al. 2004,
338 Oth et al. 2011, Laurendeau et al. 2013). These limits ensure a broad frequency bandwidth for κ
339 calculations of at least 10 Hz per record, which reduces potential bias from local amplification
340 effects (Parolai and Bindi, 2004; Ktenidou et al. 2016). The arithmetic average of the resulting
341 κ<sub>r_AS</sub> estimated from two orthogonal horizontal components is used because it is not affected by
342 the record azimuth and its implementation can reduce the variability in κ<sub>r_AS</sub> caused by ground
343 motion directionality (Ji et al. 2020). However, a prescribed, fixed-frequency band does not
344 guarantee the most appropriate linear regression for the high-frequency spectral decay in all cases.
345 Thus, we further investigate the performance of the fixed-frequency window and the effects of the
346 frequency band selection for weak and strong ground motion records by comparing individual
347 κ<sub>r_AS</sub> values estimated from a pre-determined fixed frequency window with their counterparts,
348 κ<sub>r_auto</sub>, resulting from an automated algorithm which does account for the most appropriate linear
349 regression.

350
351 The automated procedure used in this paper follows a similar protocol as those presented in
352 Sonnemann, and Halldorsson (2017) and Pilz et al. (2019), which focus on finding an appropriate
353 frequency band ([$f_1, f_2$]) to describe the linear decay in the high frequency range over a relatively
354 broad frequency window. As part of the automated protocol, the minimum $f_1$ is selected as the
355 maximum value between $1.4f_0$, and $1.4f_c$. To ensure a minimum frequency bandwidth of 10 Hz,
356 the maximum value of $f_1$ is 15 Hz and the minimum $f_2$ corresponds to ($f_1+10$) Hz. With 0.5 Hz
357 increments in $f_1$ and $f_2$, $f_1$ is varied from the maximum value between $1.4f_0$ and $1.4f_c$ to 15 Hz, while
358 $f_2$ changes from ($f_1+10$) Hz to 25 Hz. Going through all the possible combinations of $f_1$ and $f_2$, the
359 frequency range with the minimum root mean square error over the frequency bandwidth is set as
360 the optimal frequency band. The errors are computed with the following equation after Pilz et al.

361 (2019):

$$P = \frac{RMS}{\sqrt{\Delta f}} \qquad (3)$$

362 Where $\Delta f$ is the frequency bandwidth, and RMS is the root mean square error between the fitting line and smoothed FAS. The FAS is smoothed with the Konno-Ohmachi filter with a coefficient of 40 (Konno and Ohmachi, 1998). It should be noted that this automated procedure returns $\kappa_{r,auto}$ values associated with an appropriate regression for a broad frequency band. However, the changes of the FAS shape in high frequencies caused by site effects or soil nonlinearity (e.g., bumps or multiple linear trends) may not be properly accounted for by the automated algorithm.

369 Figure 6 compares $\kappa_{r\_AS}$ and $\kappa_{r, auto}$ for all selected ground motions at FKSH14, where overall similar $\kappa_{r\_AS}$ and $\kappa_{r, auto}$ estimates are observed and discrepancies are more significant at the ground surface than at depth. The remaining stations also show an acceptable agreement between the two methods at the surface and at depth. However, there are few records that show significant differences between $\kappa_{r\_AS}$ and $\kappa_{r, auto}$ (e.g., one surface record has $\kappa_{r\_AS}$ of about 0.025 s, while $\kappa_{r, auto}$, auto is almost 0.05 s). This reflects the variability of $\kappa_{r}$ as a function of the frequency band selected. For example, there are 13% of records at AICH17 showing clear variations with the frequency band selection based on the coefficient of variation (COV) obtained for all frequency windows evaluated by the automatic procedure (i.e., COV larger than 0.15). Differences between $\kappa_{r\_AS}$ and $\kappa_{r, auto}$ are mainly caused by some empirical FAS shapes, for example, when multiple linear decaying trends are present in the high-frequency range. The latter cannot be captured by the single linear decay assumption within the $\kappa$-model. More research is needed to study how more
complex wave propagation effects in the high frequency range can be captured. Testing the appropriateness of the κ-model as introduced by Anderson and Hough (1984) for these cases is outside the scope of this study, but it constitutes an area of relevant future research.

Negative $\kappa_{r\_AS}$ values are excluded from this work. Overall, there are 14 of 20 stations with less than 20% negative $\kappa_{r\_AS}$ values (i.e., arithmetic means of two horizontal components), and 1 of 20 stations with more than 50% negative $\kappa_{r\_AS}$ estimates. To further understand the observed negative $\kappa_{r\_AS}$ values, the corresponding FAS are reinspected. These FAS generally show bumps or multiple linear decays in high frequencies. The pre-determined frequency band used in this study is then not able to capture the linear decay appropriately and leads to negative $\kappa_{r\_AS}$ values. Moreover, amplifications at higher frequencies may also affect estimates of $\kappa_{r\_AS}$ and lead to negative values.

The original $\kappa_r$ model (i.e., the linear decay of FAS in log-linear scale for high frequency ranges per record proposed by Anderson and Hough, 1984), requires the site response at the site of interest to be almost flat in the high frequency range. Complex in-situ site conditions that lead to high frequency amplification (e.g., heterogeneities in the near-surface, and/or shallow impedance contrasts) can challenge this simple linear decay model. Our observations suggest the need to further explore the limitations and simplified assumptions of the Anderson and Hough (1984) $\kappa_r$ model for it to be applicable or extended to more complex environments and stronger ground shaking. Accounting for the discrepancies between actual field conditions and assumptions suggested by the Anderson and Hough (1984) $\kappa$ model could reduce the large variability observed in $\kappa_{r\_AS}$ estimates. Further research should focus on evaluating the limitations of the Anderson and Hough (1984) model and its potential modification to capture more complex wave propagation patterns in heterogeneous media, especially near the surface. In this work, such an investigation is
405 not included, but future efforts of the authors envision the assessment of negative $\kappa_{r\_AS}$ estimates
406 observed in this work as a first step to improve the existing $\kappa$ model.

407
408 **$\kappa_0$-model**
409 $\kappa_{r\_AS}$ is generally modeled with contributions from a site-specific component ($\kappa_0$), a path
410 component ($\kappa_R$), and a source component ($\kappa_s$). The source component, $\kappa_s$, is often assumed to be
411 negligible and its contribution is reduced by using a dataset with sufficient records (Van Houtte et
412 al. 2011, Ktenidou et al. 2014). A linear distance-dependency model is commonly applied to
413 capture the path component $\kappa_R$, which represents source-to-site effects or regional attenuation
414 (Hough et al. 1988, Anderson, 1991, Ktenidou et al. 2013, Boore and Campbell, 2017). Thus, the
415 most commonly accepted model is described below:

$$416 \qquad \kappa_{r\_AS} = \kappa_0 + \kappa_R \times R_e \qquad (4)$$

417 where $\kappa_0$ and $\kappa_{r\_AS}$ are in units of time (s), $\kappa_R$ is in units of second per kilometer (s/km) and $R_e$ refers to epicentral distance in km. This model is valid when a unique source-to-site path is assumed for each record along with a homogeneous and frequency-independent seismic quality factor $Q$ (Knopoff 1964). In this paper, the assumption of a unique source-to-site path is supported by using ground motions with $R_e$ less than 150 km (e.g., Palmer and Atkinson 2020, Cabas et al. 2017, Ktenidou et al. 2013). In addition, ground motions whose travel path crosses Japan's volcanic belt are not included in our database to minimize the likelihood of seismic waves propagating through regions with varying $Q$ values (Pei et al. 2009, Nakano et al. 2015).

425
426 The model described by Equation (4) is straightforward to apply when only surface linear ground
427 motion datasets are used. However, the incorporation of nonlinear and borehole ground motions
428 adds complexity to the estimation of regional attenuation as captured by $κ_R$. In this paper, we
429 assume that soil nonlinearity is triggered near the surface, which is consistent with previous studies
430 showing that nonlinear behavior occurs mostly in the superficial soil layers (i.e., Régnier et al.,
431 2013; Bonilla et al., 2019; Qin et al., 2020). Thus, nonlinear soil behavior is treated as a site
432 contribution rather than a path contribution, and we assume the regional attenuation to be identical
433 for linear and nonlinear ground motions recorded at the ground surface and at depth.
434

435 Analogous to the formulation suggested by Douglas et al. (2010) for soil and rock sites, we propose
436 a model based on Equation (4), which includes linear and nonlinear surface and borehole records:

$$ \kappa_{r\_AS} = N_1\kappa_{0\_depth} + N_2\kappa_{0\_lin\_sur} + N_3\kappa_{0\_nl\_sur} + N_4\kappa_{0\_tran\_sur} + \kappa_R \times R_e \quad (5) $$

437 where $κ_{0_depth}$ is the site-specific $κ_0$ at depth (i.e., depth of borehole sensor), and $κ_{0_bin\_sur}$, $κ_{0_trans\_sur}$ and $κ_{0_nls\_sur}$ are the site-specific linear, transitional and nonlinear $κ_0$ at the surface. The coefficients
438 $N_1$, $N_2$, $N_3$, and $N_4$ are defined as follows:

$$ N_1 = \begin{cases} 1 & \text{for dataset at depth} \\ 0 & \text{otherwise} \end{cases} $$

$$ N_2 = \begin{cases} 1 & \text{for linear dataset at surface} \\ 0 & \text{otherwise} \end{cases} $$

$$ N_3 = \begin{cases} 1 & \text{for nonlinear dataset at surface} \\ 0 & \text{otherwise} \end{cases} $$

$$ N_4 = \begin{cases} 1 & \text{for transitional dataset at surface} \\ 0 & \text{otherwise} \end{cases} $$

441 The parameter $N_4$ in Equation (5) only takes a value of 1 when ground motions classified as
442 transitional are considered independently and not as part of the linear or nonlinear datasets. That
443 means that when transitional ground motions are excluded from analysis (i.e., AP1, see Table 2)
444 or included into either the linear (i.e., AP2, see Table 2) or nonlinear (i.e., AP3, see Table 2)
445 datasets, the coefficient $N_4$ will be equal to zero.
446
450

451

## Results and Discussion

452

### Effects of soil nonlinearity on empirical $\kappa_{r\_AS}$

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

First, we study the influence of soil nonlinearity on $\kappa_{r\_AS}$ estimates per record at each site. It should be noted that soil nonlinearity is commonly triggered at shallower soil layers (Régnier et al. 2013), so we only focus on surface records in this section. Figure 7 depicts calculated $\kappa_{r\_AS}$ values at the surface against the corresponding $I_{\gamma,0}$ values at FKSH14 ($V_{s30} = 237$ m/s) and MYGH10 ($V_{s30} = 348$ m/s). As described in Equation (4), $\kappa_{r\_AS}$ is affected by both local site conditions and path effects in the context of a linear-elastic deformation analysis. Hence, the colorbar in Figure 7 represents varying epicentral distances, and the sizes of markers represent the corresponding $PGA_{rotD50}$. An overall increasing trend of $\kappa_{r\_AS}$ with increasing intensity of ground shaking (either evidenced by increased PGA or $I_{\gamma,0}$ values) is observed at FKSH14 for events that share similar epicentral distances. A slightly decreasing trend is observed for short-distance records with $R_e$ less than about 50 km and high $PGA_{rotD50}$. However, $\kappa_{r\_AS}$ values corresponding to those shorter distance and higher $PGA_{rotD50}$ events (i.e., the largest circles in Figure 7a) are larger than their counterparts for low-intensity motions (i.e., the smallest circles in Figure 7a) regardless of the $R_e$. To the best knowledge of the authors, this observation has not been reported before and could be associated with the depth of influence of $\kappa_r$. Variations in $\kappa$ due to strong nonlinear effects may be a function of a more significant contribution of the site to the overall attenuation, which may not be necessarily the case for smaller amplitude events. There might be several mechanisms of attenuation combined, and their contributions as captured by $\kappa$ need to be further investigated. Similar trends are observed at other seven sites with $V_{s30}$ less than 400 m/s, which include AICH17,
472 CHBH13, FKSH11, IBRH20, IWTH26, MYGH07, and TCGH16, and at KMMH12 with $V_{s30}$
473 greater than 400 m/s.

474
475 The increasing $\kappa_{r\_AS}$ trend with increasing $PGA_{rotD50}$ and $I_{\gamma,0}$ is not as significant at MYGH10
476 (Figure 7b), which has relatively stiffer site conditions than FKSH14. Either no correlation or a
477 slightly decreasing trend is found at other stiff sites with $V_{s30}$ greater than 400 m/s (i.e., FKSH21,
478 NIGH12, NGNH29, NIGH07, KMMH01, and IBRH16), and at four softer sites with $V_{s30}$ between
479 300 and 400 m/s (i.e., IWTH21, FKSH18, FKSH19, and IBRH17). We note that the number of
480 available nonlinear records for the $R_e$ ranges at the sites where the decreasing trend is observed is
481 rather limited. Additional nonlinear records at those sites are necessary (i.e., stronger intensity
482 ground motions) to further evaluate the contributing factors to a potential decreasing trend in $\kappa_{r\_AS}$
483 values. However, in general, we observe that positive correlations between $\kappa_{r\_AS}$ and the intensity
484 of ground shaking are more significant at softer sites (e.g., TCGH16 with $V_{s30}$ of 213 m/s) than at
485 stiffer sites (e.g., KMMH12 with $V_{s30}$ of 408 m/s). These data support that the onset of soil
486 nonlinearity can affect $\kappa_{r\_AS}$ estimates, but such influence is station-dependent. The level of soil
487 nonlinearity can be unique at each site (for a similar intensity of ground shaking) because of the
488 characteristics of shallow geologic structures (e.g., differences in velocity gradients and seismic
489 impedance contrasts) and the location of low shear-wave velocity layers. Thus, subsurface
490 conditions can play a key role on the effects of nonlinearity on $\kappa_{r\_AS}$. We observe the same patterns
491 shown in Figure 7 when using our results from the automated procedure to compute $\kappa_{r,auto}$.

492
493 **Effects of soil nonlinearity on the empirical $\kappa_0$-model**

494 Linear, transitional, and nonlinear ground motion datasets are used in this section to evaluate the
495 $\kappa_0$-model beyond the linear-elastic regime. We explore four approaches (i.e., AP1 to AP4) to
496 incorporate records within the transitional dataset into the $\kappa_0$-model presented in Equation (5).
497 Table 2 summarizes how the identified linear, transitional, and nonlinear datasets are used to
498 estimate $\kappa_{0\_lin\_sur}$ and $\kappa_{0\_nl\_sur}$.

499
500 Figures 8 and 9 present the resulting $\kappa_0$ models from each approach at FKSH14 and MYGH10,
501 respectively. Figure 8 shows that $\kappa_{r\_AS}$ and $\kappa_0$ values corresponding to the nonlinear ground
502 motions (regardless of the selected approach to construct the nonlinear dataset) are larger than their
503 linear counterparts at FKSH14. However, results at the stiffer station presented in Figure 9 show
504 little disagreement between $\kappa_{r\_AS}$ and $\kappa_0$ values corresponding to the linear and nonlinear motions
505 (regardless of the approach to construct each dataset).

506
507 Variations in $\kappa_{0\_lin\_sur}$ estimates are observed as a function of the approach considered to construct
508 the linear datasets (see specific values in Table 3 for FKSH14 and MYGH10, the results for other
509 selected stations are available in the electronic supplement). Similarly, variations in $\kappa_{0\_nl\_sur}$ values
510 are also found across the different approaches to define the nonlinear datasets. At FKSH14 (Figure
511 8), $\kappa_{0\_lin\_sur}$ estimates are more variable as a function of the dataset definitions (with a maximum
512 difference of 15.83% across approaches AP1 to AP4), compared to $\kappa_{0\_nl\_sur}$ values (with a
513 maximum difference of 5.10%). In addition, Figure 8 (d) shows that data points corresponding to
514 the transitional dataset are more compatible with their counterparts within the nonlinear dataset,
515 which may indicate that at FKSH14, the level of nonlinearity induced by the transitional dataset is
closer to that induced by the ground motions in the nonlinear dataset. Other study sites such as AICH17 ($V_{s30}$ = 314 m/s) and IWTH21 ($V_{s30}$ = 521 m/s) also show that $\kappa_{0\_lin\_sur}$ estimates are more sensitive to dataset selections. In contrast, variations of $\kappa_{0\_lin\_sur}$ and $\kappa_{0\_nl\_sur}$ across datasets at MYGH10 (Figure 9) are small, with maximum differences of only 1.56% and 2.53%, respectively. Large differences in $\kappa_{0\_nl\_sur}$ estimates across datasets are observed at eight sites, but the limited number of nonlinear records at some of those sites may be the main contributing factor (e.g., there are only six nonlinear records at IWTH21, which results in a maximum difference of 47.05% for $\kappa_{0\_lin\_sur}$ and 11.17% for $\kappa_{0\_nl\_sur}$). Adding transitional records to either the linear or the nonlinear dataset at such sites can bias the regression model. In general, differences in the $\kappa_0$-model as a function of the selected dataset are observed in 45% of our study sites (with differences in $\kappa_{0\_lin\_sur}$ or $\kappa_{0\_nl\_sur}$ values greater than 10%). This is a relevant observation because it demonstrates the importance of selecting appropriate ground motions even for typical $\kappa_0$ estimations (i.e., in the linear-elastic regime) at a given site.

Our findings suggest that the development of a $\kappa_0$-model beyond the linear-elastic regime requires an evaluation of the definition of what constitutes linear and nonlinear ground motion datasets. The identification of transitional ground motion datasets in this study allows us to assess which records provide estimates of $\kappa_{r\_AS}$ that are closer to either the linear or the nonlinear behavior at different sites. Differences in behavior triggered by the records within the transitional database may be caused by unique local site conditions (i.e., the level of soil nonlinearity developed at each site) or by limitations of the simplified definition used herein to classify transitional records (i.e., as a function of PGA and $I_\gamma$). Identifying appropriate linear and nonlinear datasets for $\kappa_{r\_AS}$ estimations requires further research to provide consistent models of near-surface attenuation that
can more effectively be implemented from small to large shear strains. However, the site-specific
response at a site of interest may impose challenges in determining appropriate dataset
classifications based on a simple, generalized criterion.

Figure 10 provides ratios of $\kappa_{0\_nl\_sur}/\kappa_{0\_lin\_sur}$ at the 20 study sites against the corresponding time-average $V_s$ value in the top 5 m ($V_{s5}$), 10 m ($V_{s10}$), and 30 m ($V_{s30}$). The ratios are computed based on the AP3 and AP4 approaches to construct linear and nonlinear datasets (Table 2). Larger ratios are observed at softer sites regardless of the dataset chosen (i.e., AP3 and AP4) for the $\kappa_0$-model.

Differences between $\kappa_0$ values in the linear and nonlinear regimes seem to be reconciled at sites with higher $V_{s5}$ (> 300 m/s), $V_{s10}$ (> 300 m/s) and $V_{s30}$ (> 400 m/s) values, where the ratios fluctuate more closely around unity particularly when using AP3. The trend of increasing ratios of $\kappa_{0\_nl\_sur}/\kappa_{0\_lin\_sur}$ with softer site conditions is better captured by $V_{s5}$ and $V_{s10}$ than by $V_{s30}$, because soil nonlinearity is more likely triggered at shallower and softer layers (Régnier et al. 2013). Hence, lower $V_s$ layers may dominate soil nonlinearity effects on $\kappa_0$. Thus, site proxies that can characterize such near-surface layers may be more informative when evaluating nonlinear soil effects on $\kappa_0$.

When grouping transitional and nonlinear ground motions (i.e., AP3), most stations result in ratios of $\kappa_{0\_nl\_sur}/\kappa_{0\_lin\_sur}$ larger than one, which can be interpreted as the signature of soil nonlinearity on the near-surface attenuation estimates (i.e., near-surface attenuation increases with increasing deformations as soil nonlinearity is triggered). These findings are consistent with the behavior of material damping ratio observed in dynamic laboratory testing of soils (i.e., increased damping ratio with increasing shear strain; Darendeli 2001; Menq 2003; Ishibashi and Zhang 1993). When
treating linear, transitional, and nonlinear datasets independently (i.e., AP4), there are 12 sites with ratios larger than one. The instances where ratios of $\kappa_{0\_nl\_sur}/\kappa_{0\_lin\_sur}$ are lower than one may result from the limited nonlinear records available at those sites coupled with the uncertainties associated with $\kappa_{0\_lin\_sur}$ (e.g., Ji et al., 2020).

Overall, the variations observed in the $\kappa_{0\_nl\_sur}/\kappa_{0\_lin\_sur}$ ratio support our hypothesis that soil nonlinearity plays a role on the estimates of near-surface attenuation from recorded ground motions. This effect is station-dependent, and further research is needed to identify the most appropriate parameter or vector of parameters capable of capturing the influence of nonlinear soil behavior on near-surface attenuation. Moreover, the relatively weaker correlation between $V_{s30}$ and the $\kappa_{0\_nl\_sur}/\kappa_{0\_lin\_sur}$ ratio evidences the challenges in connecting site conditions and soil nonlinearity via a single site parameter. Multiple parameters that can describe attenuation and impedance effects from the shallow and deep geologic structures should be investigated. The $\kappa_{0\_nl\_sur}/\kappa_{0\_lin\_sur}$ ratio corresponding to IWTH21 ($V_{s30} = 521$ m/s) is not shown in Figure 10 because it is very large (i.e., approximately 1.8). This observation may result from uncertainties associated with $\kappa_{r\_AS}$ values propagating to estimates of $\kappa_0$ when the fixed frequency band approach is applied for all records without consideration of the optimal linear decay trend. In fact, the corresponding $\kappa_{0\_nl\_sur}/\kappa_{0\_lin\_sur}$ ratio when implementing the automated procedure is approximately 0.90 for this station. Finally, Figure 10 shows less scatter in $\kappa_{0\_nl\_sur}/\kappa_{0\_lin\_sur}$ ratios when using datasets defined by AP3. In addition, using AP3 results in ratios either larger than one or approaching one for most stations (i.e., only FKSH19 and KMMH01 results in a ratio lower than unity), which is consistent with our conceptual basis for increased attenuation with the onset of nonlinear soil behavior.
Therefore, we adopt the AP3 approach (which includes transitional records into nonlinear dataset) to evaluate predictions of near-surface attenuation in the next section of this paper.

Effects of soil nonlinearity on predicted near-surface attenuation

Site-specific $\kappa_{0\_lin\_sur}$ or $\kappa_{0\_nl\_sur}$ values from Equation (5) allow for the comparison of empirical estimates of near-surface attenuation, but these two parameters represent the average attenuation of all records in linear and nonlinear regimes. Thus, in this section, we introduce the predicted near-surface attenuation at zero-distance ($\kappa_{0\_pred}$), which is expected to capture the attenuation contributed by the superficial soil layers per event by removing the path contributions from $\kappa_{r\_sur}$.

$\kappa_{0\_pred}$ is modeled as:

$$ \kappa_{0\_pred} = \kappa_{r\_AS\_sur} - R_e \cdot \kappa_R \quad (6) $$

where $\kappa_{r\_AS\_sur}$ refers to the individual $\kappa_{r\_AS}$ value for a surface ground motion, and the path-component, $\kappa_R$, corresponds to the values derived with Equation (5) at each site of interest. We assume that by removing the effect of the path-component $\kappa_R$ from $\kappa_{r\_AS\_sur}$ values per record, the remaining $\kappa_{0\_pred}$ becomes an approximation to the attenuation contributed by the shallower sedimentary deposits per event. Thus, we can explore how the near-surface attenuation changes with the various input ground motion amplitudes at the site of interest.

Figure 11 provides comparisons between $\kappa_{0\_pred}$, ground shaking intensity, and deformation as captured by PGA$_{rotD50}$ and $I_{\gamma,0}$ at FKSH14 and MYGH10. Both colors and sizes of markers represent the PGA$_{rotD50}$ values per record. The red dashed-lines result from a local regression model characterizing the $\kappa_{0\_pred}$ versus $I_{\gamma,0}$ data. Triangles and circles represent the linear and nonlinear
ground motions (identified with AP3), respectively. Values of $\kappa_{0\_pred}$ first increase and then decrease with increasing $PGA_{rotD50}$ and $I_{\gamma,0}$ at FKSH14. This behavior is also observed at other 7 sites (i.e., AICH17, CHBH13, FKSH11, IWTH26, KMMH12, MYGH07, TCGH16). Even with the decreasing trend for large $I_{\gamma,0}$, $\kappa_{0\_pred}$ is still generally higher than its counterpart in the linear regime (i.e., the means of linear and nonlinear $\kappa_{0\_pred}$ are 0.050 sec and 0.0605 sec, respectively). Overall, $\kappa_{0\_pred}$ values at FKSH14 corresponding to larger deformations and higher $PGA_{rotD50}$ are larger than those corresponding to weaker ground motions. In contrast, only a weak correlation to the intensity of ground shaking and deformation in situ is observed at MYGH10. These results are consistent with our estimations of $\kappa_r$ shown in Figure 7. Soil nonlinear behavior can influence near-surface attenuation as captured by $\kappa_r$ and $\kappa_0$, and local site conditions may play a key role in this process. The remaining approaches explored in this study (i.e., AP1, AP2, and AP4) provide similar results as those shown in Figure 11.

Figure 12 compares the probability distribution of $\kappa_{0\_pred}$ values from the linear and nonlinear datasets (AP3 case) at FKSH14 and MYGH10. The resulting $\kappa_{0\_pred}$ values are fitted with a Gaussian distribution and the corresponding probability density functions (PDFs) are represented by red lines. A shift to the right (i.e., toward larger $\kappa_{0\_pred}$ values) of the theoretical PDF is observed at FKSH14 as ground motions from the linear and nonlinear datasets are considered. The mean $\kappa_{0\_pred}$ estimates change from 0.05 s for the linear dataset to 0.0605 s for the nonlinear dataset at FKSH14 (i.e., a difference of 21%). In contrast, the variation of mean $\kappa_{0\_pred}$ between linear and nonlinear datasets at MYGH10 is 4%. Most of our study sites have either a significant increase in their mean $\kappa_{0\_pred}$ when using the nonlinear dataset (i.e., an increase of more than 20%) or only a slight increase. There are only 4 stations that show a decrease in their mean $\kappa_{0\_pred}$ values with
respect to the linear dataset when using the nonlinear one (i.e., FKSH19, KMMH01, NIGH07, and NIGH12). Statistical hypothesis tests (i.e., t-test) are conducted to analyze whether there is a statistically significant difference between the means of linear and nonlinear $\kappa_{0\_pred}$ distributions. Considering a critical value of 5%, p-values at each station are shown in Table 4. Nine out of 20 stations display statistically significant differences between their mean $\kappa_{0\_pred}$ corresponding to the linear and nonlinear datasets. Table 4, Figures 11 and 12 show that soil nonlinear behavior can affect $\kappa_{0\_pred}$ at the sites selected in this study, although this influence is station-dependent. At the stations that display apparent effects of nonlinearity on $\kappa$, an increasing trend in predicted near-surface attenuation with increasing ground shaking intensity and/or increasing deformation is observed.

## Conclusions

In this work, we investigated the influence of soil nonlinear behavior on $\kappa_{r\_AS}$ values per record and site-specific $\kappa_0$ estimates at 20 stations selected from the KiK-net database. To avoid potential bias on our results due to the calculation process, we also examined the effects of the frequency band selection on $\kappa_{r\_AS}$ estimates, and the differences between using the S-wave window or the entire time series FAS. We compared results from a predetermined fixed-frequency window approach with an automated procedure that considers multiple frequency windows. The latter is capable of finding the optimal frequency band per record for all records at each site. The selection of a common, fixed and broader frequency band for $\kappa_{r\_AS}$ estimations reduced the scatter and bias in the data, while providing reasonable estimations of $\kappa_{r\_AS}$. On the other hand, values of $\kappa_{r\_AS}$ computed from the S-wave window FAS were reasonably similar to their counterparts based on the entire time series FAS. Hence, the analyses presented in this paper were conducted with $\kappa_r$
values estimated by the fixed-frequency band approach and the FAS corresponding to the entire time series.

A consistent identification of ground motions that trigger nonlinear behavior in sedimentary deposits is also necessary to quantify near-surface attenuation beyond the linear-elastic regime. Based on the examination of an in-situ stress-strain proxy, namely the correlation between PGA<sub>rotD50</sub> and I<sub>γ,0</sub>, we found that the variation of shear strains with ground shaking intensity at the onset of nonlinear soil behavior is site-specific. A unique threshold for a single parameter, whether it is PGA<sub>rotD50</sub> or I<sub>γ,0</sub>, was not able to capture the onset of soil nonlinearity at our study sites in a consistent manner across all sites. Hence, we proposed a hybrid method to classify linear and nonlinear ground motions considering both, PGA<sub>rotD50</sub> and I<sub>γ,0</sub>, which resulted in linear, transitional, and nonlinear datasets at each site.

Increasing κ<sub>r_AS</sub> values with increasing PGA<sub>rotD50</sub> or I<sub>γ,0</sub> for ground motions with similar epicentral distances were observed at about half of our study sites. This trend was more consistently observed at softer sites. Additionally, we found that κ<sub>0</sub>-models could be biased by the definition of linear and nonlinear ground motion datasets. Hence, we studied the effects of ground motion categorization and proposed a hybrid classification scheme for linear and nonlinear records. We defined transitional ground motions as those associated with soil behavior between the linear-elastic and nonlinear regimes. Even though more research is necessary to define robust classification schemes for linear and nonlinear ground motions, we observed that including the transitional motions into the nonlinear dataset reduced the variability associated with κ estimations at our study sites.
Our results also revealed differences between $\kappa_{0\_lin\_sur}$ (i.e., $\kappa_0$ corresponding to the linear-elastic regime) and $\kappa_{0\_nl\_sur}$ (i.e., $\kappa_0$ for the nonlinear regime) at most sites when implementing the $\kappa_0$-model using ground motions classified by AP3 (which includes transitional records into the nonlinear dataset). Such differences were more prevalent among softer sites. Site parameters such as $V_{s5}$, $V_{s10}$, and $V_{s30}$ were used in this study to investigate the influence of soil conditions on the effects of nonlinearity on $\kappa_0$. Considering that high frequencies have short wavelengths, and that nonlinear soil behavior is triggered in low velocity layers more often located at a shallow depth, site proxies such as $V_{s5}$ and $V_{s10}$ may be more informative than $V_{s30}$ when assessing effects of nonlinearity on $\kappa_0$. For instance, large $V_{s30}$ values do not imply that all near-surface layers have a large $V_s$. The ratio of $\kappa_{0\_nl\_sur}$ and $\kappa_{0\_lin\_sur}$ decreases and approaches one for increasing $V_{s5}$, $V_{s10}$, and $V_{s30}$, when using the AP3 method to define nonlinear ground motion datasets.

The hypothesis posed and tested in this paper focused on the effects of ground shaking intensity on induced shear strains in sedimentary deposits and associated consequences on the attenuation experienced by seismic waves (particularly in the high frequency range). In general, we find that soil nonlinear behavior can affect estimates of $\kappa_{r\_AS}$ and $\kappa_0$, but our results show that this influence is station-dependent. This is reasonable because the wave propagation of short wavelength waves is highly affected by heterogeneities in the soil or rock, local geologic structures, and topography. Moreover, the level of soil nonlinearity can be distinct at a given site (even when site parameterizations such as $V_{s30}$ are similar and the considered intensity of ground shaking is also similar) because of the complexities of the in situ subsurface conditions (e.g., differences in velocity gradients and seismic impedance contrasts). We note that 2D/3D site effects may affect
697 ground motions recorded at six of our stations. The influence of soil nonlinearity on κ values
698 computed at these stations (i.e., AICH17, CHBH13, FKSH21, IBRH20, KMMH01, KMMH12;
699 based on the classification of Thompson et al. 2012 and Pilz and Cotton 2019) may be masked by
700 the combined effects of wave scattering and topographic effects. Further research is necessary to
701 evaluate the contributions of the aforementioned mechanisms on κ estimates at stations subjected
702 to 2D/3D site effects. Likewise, future work should focus on collecting and analyzing additional
703 strong ground motion data to identify local site conditions more conducive to generate significant
704 changes in near-surface attenuation as captured by κ₀ when nonlinear soil behavior is triggered.
705

706 Complexities in the wave propagation phenomenon driven by scattering effects and amplification
707 in the high-frequency range can result in negative estimates of κr. In this study, we obtained
708 negative κr_AS estimates when multiple linear decaying trends, bumps, and high frequency
709 amplifications affected the corresponding FAS spectral shape. The identification of multiple linear
decays in the high-frequency range supports previous work on the bias in κr_AS associated with the
710 selection of the frequency band. The bumps and amplifications in the high frequencies present in
711 the FAS of some of the ground motions in our database hint that the site response may not be
712 approximately flat within the frequency range of interest for κr_AS calculation. Considering that a
713 flat site response is one of the assumptions of the Anderson and Hough (1984) κ-model, further
714 research is needed to overcome this limitation at sites where this is not the case. This work not
715 only provides evidence of the need to understand and quantify κ in both, the linear and nonlinear
716 regimes, but it also presents the limitations of the current κ-model when it comes to characterizing
717 attenuation when conditions deviate from the original assumptions embedded in the Anderson and
718 Hough (1984) model.
719
## Data and Resources

Accelerograms and geotechnical data are downloaded from the KiK-net network at http://www.kyoshin.bosai.go.jp (last accessed May 2020). The earthquake information is available from F-net network at http://www.fnet.bosai.go.jp/top.php (last accessed May 2020). The supplemental material to this article includes two tables and three figures. The tables provide the $\kappa_0$-model results when using the datasets defined by AP1, AP2, AP3, and AP4, and the $\kappa_{0\_pred}$ results estimated with datasets defined by AP3. The three sets of figures presented depict the results corresponding to our 20 stations as follows:

* Surface $\text{PGA}_{\text{rotD50}}$ against $I_{\gamma,0}$ (i.e., results analogous to those presented in Figure 4 for MYGH10).

* Surface $\kappa_{r\_AS}$ estimates and their corresponding $\text{PGA}_{\text{rotD50}}$, $I_{\gamma,0}$ and $R_e$ values for selected ground motions at each study site (i.e., results analogous to those presented in Figure 7 for FKSH14 and MYGH10).

* Estimated surface $\kappa_{0\_pred}$ and their corresponding ground shaking intensity and in situ deformation characterized by $\text{PGA}_{\text{rotD50}}$ against $I_{\gamma,0}$, respectively (i.e., results analogous to those presented in Figure 11 for FKSH14 and MYGH10).

## Acknowledgements

The authors would like to appreciate that the National Research Institute for Earth and Disaster Prevention (NIED) provides the ground motion and site information. The authors acknowledge Dr. Adrian Rodriguez-Marek and Mahdi Bahrampouri for sharing their processed dataset. We also are thankful to Associate Editor Stefano Parolai and two anonymous reviewers for the useful comments. This work was sponsored by the U.S. Geological Survey under Grant Number
G19AP00058 and the NCSU Internationalization Seed Grant. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the opinions or policies of the U.S. Geological Survey.

### Reference

Anderson, J. G., & Hough, S. E. (1984). A model for the shape of the Fourier amplitude spectrum of acceleration at high frequencies. *Bull Seismol Soc Am.* **74**, no. 5, 1969-1993.

Anderson, J.G. (1991). A preliminary descriptive model for the distance dependence of the spectral decay parameter in southern California. *Bull Seismol Soc Am.* **81**, no. 6, 2186-2193.

Al Atik, L., Kottke, A., Abrahamson, N. and Hollenback, J., (2014). Kappa (κ) scaling of ground-motion prediction equations using an inverse random vibration theory approach. Bulletin of the Seismological Society of America, 104, no. 1, 336-346.

Aoi, S., T. Kunugi, and H. Fujiwara (2004). Strong-motion seismograph network operated by NIED: K-NET and KiK-net, J. Japan Assoc. Earthq.Eng. 4, no. 3, 65-74.

Bahrampouri, M., Rodriguez-Marek, A., Shahi, S., & Dawood, H. (2020). An update database for ground motion parameters for KiK-net records. Earthquake Spectra.

Bonilla, L.F., Tsuda, K., Pulido, N., Régnier, J. and 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 and Space, 63(7), p.50.

Bonilla, L.F., Guéguen, P. and Ben-Zion, Y., 2019. Monitoring Coseismic Temporal Changes of Shallow Material during Strong Ground Motion with Interferometry and Autocorrelation. Monitoring Coseismic Temporal Changes of Shallow Material during
765 Strong Ground Motion. Bulletin of the Seismological Society of America, 109(1), pp.187-198.

767 Bonilla, L.F. and Ben-Zion, Y., 2020. Detailed space-time variations of the seismic response of the shallow crust to small earthquakes from analysis of dense array data. Geophysical Journal International.

770 Boore, D.M., 2003. Simulation of ground motion using the stochastic method. *Pure and applied geophysics*, 160(3-4), pp.635-676.

772 Boore, D.M. and Bommer, J.J., 2005. Processing of Strong-Motion Accelerograms: Needs, Options and Consequences. Soil Dynamics and Earthquake Engineering, 25(2), 93-115.

774 Boore, D. M. (2010). Orientation-independent, nongeometric-mean measures of seismic intensity from two horizontal components of motion. *Bull Seismol Soc Am.* 100, no. 4, 1830-1835.

776 Boore, D.M. and Campbell, K.W., 2017. Adjusting central and eastern North America ground-motion intensity measures between sites with different reference-rock site conditions. Bulletin of the Seismological Society of America, 107(1), pp.132-148.

779 Campbell, K.W., 2003. Prediction of strong ground motion using the hybrid empirical method and its use in the development of ground-motion (attenuation) relations in eastern North America. Bulletin of the Seismological Society of America, 93(3), pp.1012-1033.

782 Cabas, A., Rodriguez-Marek, A., & Bonilla, L. F. (2017). Estimation of Site-Specific Kappa (κ0)-Consistent Damping Values at KiK-Net Sites to Assess the Discrepancy between Laboratory-Based Damping Models and Observed Attenuation (of Seismic Waves) in the Field. *Bull Seismol Soc Am.* **107**, no. 5, 2258-2271.
786 Chandra, J., Guéguen, P. and Bonilla, L.F., 2014. Application of PGV/VS proxy to assess
787 nonlinear soil response-from dynamic centrifuge testing to Japanese K-NET and KiK-net
788 data. In Second European Conference on Earthquake Engineering and Seismology,
789 Istanbul Aug (pp. 25-29).

790 Chandra, J., Guéguen, P., Steidl, J.H. and Bonilla, L.F., 2015. In situ assessment of the G-γ curve
791 for characterizing the nonlinear response of soil: Application to the Garner Valley
792 downhole array and the wildlife liquefaction array. Bulletin of the Seismological Society
793 of America, 105(2A), pp.993-1010.

794 Chandra, J., Gueguen, P. and Bonilla, L.F., 2016. PGA-PGV/Vs considered as a stress-strain
795 proxy for predicting nonlinear soil response. Soil dynamics and earthquake engineering, 85,
796 pp.146-160.

797 Darragh, R.B. and Shakal, A.F., 1991. The site response of two rock and soil station pairs to strong
798 and weak ground motion. Bulletin of the Seismological Society of America, 81(5),
799 pp.1885-1899.

800 Darendeli M. B. (2001). Development of a new family of normalized modulus reduction and
801 material damping curves, Ph.D. Thesis, University of Texas at Austin, Austin, Texas.

802 Dimitriu, P., Theodulidis, N., Hatzidimitriou, P., & Anastasiadis, A. (2001). Sediment non-linearity and attenuation of seismic waves: a study of accelerograms from Lefkas, western Greece. *Soil Dynam Earthquake Eng.* **21**, no.1, 63-73.

805 Durward, J.A., Boore, D.M., Joyner, W.B., Durward, J.A., Boore, D.M. and Joyner, W.B., 1996.
806 The amplitude dependence of high-frequency spectral decay: constraint on soil
807 nonlinearity. In Proc, of the International Workshop on Site Response Subjected to Strong Earthquake Motions (pp. 82-103).
808

809 Drouet, S., F. Cotton, and P. Gueguen (2010). VS30, κ, regional attenuation and Mw from accelerograms: Application to magnitude 3–5 French earthquakes, Geophys. J. Int. 182, 880–898.
810 811

812 Douglas, J., Gehl, P., Bonilla, L.F. and Gélis, C., (2010). A κ model for mainland France. Pure Appl Geophys. 167, no.11, 1303-1315.
813

814 Douglas, J. and Boore, D.M., 2011. High-frequency filtering of strong-motion records. Bulletin of Earthquake Engineering 9(2), 395-409.
815

816 Edwards, B., Ktenidou, O.J., Cotton, F., Abrahamson, N., Van Houtte, C. and Fäh, D., (2015). Epistemic uncertainty and limitations of the κ 0 model for near-surface attenuation at hard rock sites. Geophys J Int. 202, no. 3, 1627-1645.
817 818

819 Fujiwara, H., S. Aoi, T. Kunugi, and S. Adachi (2004). Strong-motion observation networks of NIED: K-NETand KiK-net, National Research Institute for Earth Science and Disaster Prevention.Garcia, D., Wald, D. J., & Hearne, M. G. (2012). A global earthquake discrimination scheme to optimize ground-motion prediction equation selection. Bull Seismol Soc Am. 102, no.1, 185-203.
820 821 822 823

824 Guéguen, P., Bonilla, L.F. and Douglas, J., 2019. Comparison of Soil Nonlinearity (In Situ Stress–Strain Relation and G/Gmax Reduction) Observed in Strong-Motion Databases and Modeled in Ground-Motion Prediction EquationsComparison of Soil Nonlinearity Observed in Strong-Motion Databases and Modeled in GMPEs. Bulletin of the Seismological Society of America, 109(1), pp.178-186.
825 826 827 828
829 Hough, S.E., Anderson, J.G., Brune, J., Vernon III, F., Berger, J., Fletcher, J., Haar, L., Hanks, L. and Baker, L., 1988. Attenuation near Anza, California. Bulletin of the Seismological Society of America, 78(2), pp.672-691.

832 Idriss, I.M., Dobry, R.U. and Sing, R.D., 1978. Nonlinear behavior of soft clays during cyclic loading. *Journal of geotechnical and geoenvironmental engineering*, 104.

834 Ishibashi, I. and Zhang, X., 1993. Unified dynamic shear moduli and damping ratios of sand and clay. *Soils and foundations*, 33(1), pp.182-191.

836 Ji, C., Cabas, A., Cotton, F., Pilz, M. and Bindi, D., 2020. Within-Station Variability in Kappa: Evidence of Directionality Effects. *Bull Seismol Soc Am.*, 110(3), pp.1247-1259.

838 Knopoff, L. (1964). *Q, Rev. Geophys. 2*, no. 4, 625-660.

839 Ktenidou, O. J., Gélis, C., & Bonilla, L. F. (2013). A study on the variability of kappa (κ) in a borehole: Implications of the computation process. *Bull Seismol Soc Am.* **103**, no.2A, 1048-1068.

842 Ktenidou, O. J., Cotton, F., Abrahamson, N. A., & Anderson, J. G. (2014). Taxonomy of κ: A review of definitions and estimation approaches targeted to applications. *Seismol Res Lett.* **85**, no.1, 135-146.

845 Ktenidou, O.J., Abrahamson, N.A., Drouet, S. and Cotton, F., (2015). Understanding the physics of kappa (κ): Insights from a downhole array. *Geophys J Int.* 203, no. 1, 678-691.

847 Ktenidou, O.J., Abrahamson, N., Darragh, R. and Silva, W., 2016. A methodology for the estimation of kappa (κ) from large datasets, example application to rock sites in the NGA-East database, and implications on design motions. *PEER Report 2016*, 1.
850 Konno, K. and Ohmachi, T., 1998. Ground-motion characteristics estimated from spectral ratio between horizontal and vertical components of microtremor. Bulletin of the Seismological Society of America, 88(1), pp.228-241.

853 Laurendeau, A., Cotton, F., Ktenidou, O. J., Bonilla, L. F., & Hollender, F. (2013). Rock and stiff-soil site amplification: Dependency on VS 30 and kappa (κ 0). Bull Seismol Soc Am. 103, no. 6, 3131-3148.

856 Lacave-Lachet, C., Bard, P.Y., Gariel, J.C. and Irikura, K., 2000. Straightforward methods to detect non-linear response of the soil. Application to the recordings of the Kobe earthquake (Japan, 1995). Journal of seismology, 4(2), pp.161-173.

859 Menq, F.Y., (2003). Dynamic properties of sandy and gravelly soils, Ph.D. Thesis, University of Texas at Austin, Austin, Texas.

861 Nakano, K., Matsushima, S., & Kawase, H. (2015). Statistical properties of strong ground motions from the generalized spectral inversion of data observed by K-NET, KiK-net, and the JMA Shindokei network in Japan. *Bull Seismol Soc Am.* **105**, no. 5, 2662-2680.

864 Oth, A., Parolai, S. and Bindi, D., (2011). Spectral analysis of K-NET and KiK-net data in Japan, Part I: Database compilation and peculiarities. Bull Seismol Soc Am. 101, no.2, 652-666.

866 Parolai, S. and Bindi, D., (2004). Influence of soil-layer properties on k evaluation. Bulletin of the Seismological Society of America, 94(1), pp.349-356.

868 Pilz, M. and Cotton, F., 2019. Does the One-Dimensional Assumption Hold for Site Response Analysis? A Study of Seismic Site Responses and Implication for Ground Motion Assessment Using KiK-Net Strong-Motion Data. Earthquake Spectra, 35(2), pp.883-905.
871 Pilz, M., Cotton, F., Zaccarelli, R. and Bindi, D., (2019). Capturing Regional Variations of Hard-Rock Attenuation in Europe. *Bull Seismol Soc Am.*

873 Perron, V., Hollender, F., Bard, P. Y., Gélis, C., Guyonnet-Benaize, C., Hernandez, B., & Ktenidou, O. J. (2017). Robustness of kappa (κ) measurement in low-to-moderate seismicity areas: Insight from a site-specific study in Provence, France. *Bull Seismol Soc Am.* 107, no.5, 2272-2292.

877 Pei, S., Cui, Z., Sun, Y., Toksöz, M.N., Rowe, C.A., Gao, X., Zhao, J., Liu, H., He, J. and Morgan, F.D., 2009. Structure of the upper crust in Japan from S-wave attenuation tomography. Bulletin of the Seismological Society of America, 99(1), pp.428-434.

880 Qin, L., Ben-Zion, Y., Bonilla, L.F. and Steidl, J.H., 2020. Imaging and Monitoring Temporal Changes of Shallow Seismic Velocities at the Garner Valley Near Anza, California, Following the M7.2 2010 El Mayor-Cucapah Earthquake. *Journal of Geophysical Research: Solid Earth*, 125(1), p.e2019JB018070.

884 Régnier, J., Cadet, H., Bonilla, L.F., Bertrand, E. and Semblat, J.F., 2013. Assessing nonlinear behavior of soils in seismic site response: Statistical analysis on KiK-net strong-motion data. *Bulletin of the Seismological Society of America*, 103(3), pp.1750-1770.

887 Régnier, J., Cadet, H. and Bard, P.Y., 2016. Empirical Quantification of the Impact of Nonlinear Soil Behavior on Site Response. *Empirical Quantification of the Impact of Nonlinear Soil Behavior on Site Response*. Bulletin of the Seismological Society of America, 106(4), pp.1710-1719.
891 Rong, M., Wang, Z., Woolery, E.W., Lyu, Y., Li, X. and Li, S., 2016. Nonlinear site response
892 from the strong ground-motion recordings in western China. Soil Dynamics and
893 Earthquake Engineering, 82, pp.99-110.

894 Seed, H.B., Wong, R.T., Idriss, I.M. and Tokimatsu, K., 1986. Moduli and damping factors for
895 dynamic analyses of cohesionless soils. *Journal of geotechnical engineering*, 112(11),
896 pp.1016-1032.

897 Sonnemann, T. and Halldorsson, B., 2017, June. Towards an Automated Kappa Measurement
898 Procedure. In International Conference on Earthquake Engineering and Structural
899 Dynamics (pp. 39-52). Springer, Cham.

900 Thompson, E.M., Baise, L.G., Tanaka, Y. and Kayen, R.E., 2012. A taxonomy of site response
901 complexity. Soil Dynamics and Earthquake Engineering, 41, pp.32-43.

902 Trifunac, M.D. and Todorovska, M.I., 1996. Nonlinear soil response—1994 Northridge,
903 California, earthquake. Journal of geotechnical engineering, 122(9), pp.725-735.

904 Van Houtte, C., Drouet, S., & Cotton, F. (2011). Analysis of the origins of κ (kappa) to compute
905 hard rock to rock adjustment factors for GMPEs. *Bull Seismol Soc Am.101*, no. 6, 2926-
906 2941.

907 Van Houtte, C., Ktenidou, O.J., Larkin, T. and Holden, C., (2014). Hard-site κ 0 (kappa)
908 calculations for Christchurch, New Zealand, and comparison with local ground-motion
909 prediction models. Bull Seismol Soc Am.104, no. 4, 1899-1913.

910 Xu, B., Rathje, E.M., Hashash, Y., Stewart, J., Campbell, K. and Silva, W.J., 2020. κ 0 for soil
911 sites: Observations from Kik-net sites and their use in constraining small-strain damping
912 profiles for site response analysis. Earthquake Spectra, 36(1), pp.111-137.
913 Yu, G., Anderson, J.G. and Siddharthan, R.A.J., 1993. On the characteristics of nonlinear soil
914 response. Bulletin of the Seismological Society of America, 83(1), pp.218-244.

915 Zalachoris, G. and Rathje, E.M., 2015. Evaluation of one-dimensional site response techniques
916 using borehole arrays. Journal of Geotechnical and Geoenvironmental Engineering,
917 141(12), p.04015053.
918 Chunyang Ji

919 Department of Civil, Construction, and Environmental Engineering,

920 3101 Fitts-Woolard Hall, NC State University

921 Raleigh, NC 27695, the United States

922 Email: cji3@ncsu.edu

923

924 Ashly Cabas

925 Department of Civil, Construction, and Environmental Engineering,

926 3175 Fitts-Woolard Hall, NC State University

927 Raleigh, NC 27695, the United States

928 Email: amcabasm@ncsu.edu

929

930 Luis Fabian Bonilla

931 GERS-SRO, Univ. Gustave Eiffel,

932 F-77447 Marne-la-Vallée, France

933 Email: luis-fabian.bonilla-hidalgo@univ-eiffel.fr

934

935 Céline Gelis

936 IRSN

937 Paris Fontenay-aux-Roses, France

938 Email: celine.gelis@irsn.fr
<table><caption>**Table 1.** Local soil conditions, number of ground motions per dataset, predetermined fixed-frequency bandwidth and thresholds for shear strain index, $I_y$ at all study sites.</caption><thead><tr><th>Station</th><th>V<sub>s30</sub><sup>1</sup><br>(m/s)</th><th>V<sub>s0</sub><sup>2</sup><br>(m/s)</th><th>V<sub>s,depth</sub><sup>3</sup><br>(m/s)</th><th>Hole<br>Depth<br>(m)</th><th>Number<br>of<br>linear<br>records</th><th>Number of<br>transitional<br>records</th><th>Number<br>of<br>nonlinear<br>records</th><th>1<sup>st</sup><br>resonant<br>frequency<br>(Hz)</th><th>Predominant<br>frequency<br>(Hz)</th><th>f<sub>1</sub><sup>4</sup><br>(Hz)</th><th>f<sub>2</sub><sup>5</sup><br>(Hz)</th><th>I<sub>γ,0,1</sub><sup>6</sup><br>(%)</th><th>I<sub>γ,0,7</sub><sup>7</sup><br>(%)</th></tr></thead><tbody><tr><td>AICH17</td><td>314</td><td>150</td><td>2200</td><td>101</td><td>23</td><td>27</td><td>12</td><td>4.07</td><td>4.07</td><td>12.65</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>CHBH13</td><td>235</td><td>220</td><td>2920</td><td>1300</td><td>139</td><td>49</td><td>11</td><td>1.78</td><td>1.78</td><td>8.95</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>FKSH11</td><td>240</td><td>110</td><td>700</td><td>115</td><td>148</td><td>140</td><td>26</td><td>1.51</td><td>9.98</td><td>13.97</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>FKSH14</td><td>237</td><td>120</td><td>1210</td><td>147</td><td>114</td><td>221</td><td>28</td><td>1.12</td><td>4.15</td><td>10.05</td><td>25.00</td><td>0.001</td><td>0.007</td></tr><tr><td>FKSH18</td><td>307</td><td>140</td><td>2250</td><td>100</td><td>158</td><td>103</td><td>16</td><td>2.59</td><td>5.69</td><td>8.95</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>FKSH19</td><td>338</td><td>170</td><td>3060</td><td>100</td><td>185</td><td>95</td><td>21</td><td>3.27</td><td>3.27</td><td>10.05</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>FKSH21</td><td>365</td><td>200</td><td>1600</td><td>200</td><td>60</td><td>17</td><td>8</td><td>3.90</td><td>3.90</td><td>12.65</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>IBRH16</td><td>626</td><td>140</td><td>2050</td><td>300</td><td>137</td><td>81</td><td>15</td><td>1.71</td><td>7.08</td><td>10.05</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>IBRH17</td><td>301</td><td>90</td><td>2300</td><td>510</td><td>117</td><td>177</td><td>18</td><td>0.93</td><td>9.30</td><td>13.01</td><td>25.00</td><td>0.001</td><td>0.007</td></tr><tr><td>IBRH20</td><td>244</td><td>180</td><td>1200</td><td>923</td><td>133</td><td>86</td><td>11</td><td>0.27</td><td>0.27</td><td>8.95</td><td>25.00</td><td>0.001</td><td>0.007</td></tr><tr><td>IWTH21</td><td>521</td><td>150</td><td>2460</td><td>100</td><td>39</td><td>24</td><td>6</td><td>5.27</td><td>5.27</td><td>7.38</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>IWTH26</td><td>371</td><td>130</td><td>680</td><td>108</td><td>79</td><td>32</td><td>11</td><td>2.12</td><td>10.17</td><td>14.24</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>KMMH01</td><td>575</td><td>150</td><td>1900</td><td>100</td><td>94</td><td>24</td><td>15</td><td>4.03</td><td>9.98</td><td>13.97</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>KMMH12</td><td>410</td><td>210</td><td>1000</td><td>123</td><td>134</td><td>34</td><td>11</td><td>3.27</td><td>8.17</td><td>11.44</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>MYGH07</td><td>366</td><td>130</td><td>740</td><td>142</td><td>59</td><td>39</td><td>11</td><td>0.93</td><td>8.61</td><td>12.06</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>MYGH10</td><td>348</td><td>110</td><td>770</td><td>205</td><td>229</td><td>132</td><td>16</td><td>0.95</td><td>10.66</td><td>14.93</td><td>25.00</td><td>0.001</td><td>0.007</td></tr><tr><td>NGNH29</td><td>465</td><td>150</td><td>1040</td><td>110</td><td>81</td><td>38</td><td>16</td><td>1.95</td><td>6.93</td><td>10.05</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>NIGH07</td><td>528</td><td>200</td><td>1600</td><td>106</td><td>29</td><td>10</td><td>11</td><td>4.12</td><td>7.08</td><td>10.05</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>NIGH12</td><td>553</td><td>240</td><td>780</td><td>110</td><td>29</td><td>9</td><td>11</td><td>2.00</td><td>5.00</td><td>12.65</td><td>25.00</td><td>0.001</td><td>0.003</td></tr><tr><td>TCGH16</td><td>213</td><td>80</td><td>680</td><td>112</td><td>112</td><td>334</td><td>35</td><td>1.27</td><td>4.81</td><td>11.27</td><td>25.00</td><td>0.001</td><td>0.007</td></tr></tbody></table>

941 ¹V<sub>s30</sub>: time averaged shear-wave velocity in the top 30 m of the soil profile

942 ²V<sub>s0</sub>: shear-wave at the ground surface

943 ³V<sub>s,depth</sub>: shear-wave velocity at the depth of the borehole sensor
944 $^{4}f_i$: the lower frequency limit to estimate individual $\kappa_i$

945 $^{5}f_2$: the upper frequency limit to estimate individual $\kappa_r$

946 $^{6}I_{\gamma,0,l}$: the shear-strain index threshold to separate linear and transitional datasets

947 $^{7}I_{\gamma,0,t}$: the shear-strain index threshold to separate transitional and nonlinear datasets
<table><caption>**Table 2.** Ground motion datasets constructed via alternative approaches (AP1 to AP4) explored in this study to implement the $κ_0$-model.</caption><thead><tr><th>Approach</th><th>κ<sub>0</sub>_lin_<sub>sur</sub></th><th>κ<sub>0</sub>_nl_<sub>sur</sub></th><th>κ<sub>0</sub>_tran_<sub>sur</sub></th><th>κ<sub>0</sub>_depth</th><th>κ<sub>0</sub>-model</th></tr></thead><tbody><tr><td>AP1</td><td>Linear dataset</td><td>Nonlinear dataset</td><td>--</td><td>Borehole dataset</td><td>Equation (5)</td></tr><tr><td>AP2</td><td>Linear and transitional datasets</td><td>Nonlinear dataset</td><td>--</td><td>Borehole dataset</td><td>Equation (5)</td></tr><tr><td>AP3</td><td>Linear dataset</td><td>Nonlinear and transitional datasets</td><td>--</td><td>Borehole dataset</td><td>Equation (5)</td></tr><tr><td>AP4</td><td>Linear dataset</td><td>Nonlinear dataset</td><td>Transitional datasets</td><td>Borehole dataset</td><td>Equation (5)</td></tr></tbody></table>

<table><caption>**Table 3.** Site-specific $κ_0$ values obtained from different dataset definitions at stations FKSH14 and MYGH10.</caption><thead><tr><th>Approach</th><th>κ<sub>0</sub>_lin_<sub>sur</sub> (s)</th><th>κ<sub>0</sub>_nl_<sub>sur</sub> (s)</th><th>κ<sub>0</sub>_tran_<sub>sur</sub> (s)</th></tr></thead><tbody><tr><td>AP1</td><td>0.0488</td><td>0.0633</td><td>--</td></tr><tr><td>AP2</td><td>0.0565</td><td>0.0638</td><td>--</td></tr><tr><td>AP3</td><td>0.0500</td><td>0.0607</td><td>--</td></tr><tr><td>AP4</td><td>0.0499</td><td>0.0638</td><td>0.0602</td></tr><tr><td>Maximum difference</td><td>15.83%</td><td>5.10%</td><td></td></tr><tr><td>AP1</td><td>0.0567</td><td>0.0563</td><td>--</td></tr><tr><td>AP2</td><td>0.0568</td><td>0.0560</td><td>--</td></tr><tr><td>AP3</td><td>0.0559</td><td>0.0572</td><td>--</td></tr><tr><td>AP4</td><td>0.0560</td><td>0.0558</td><td>0.0575</td></tr><tr><td>Maximum difference</td><td>1.56%</td><td>2.53%</td><td></td></tr></tbody></table>
<table><caption>**Table 4.** T-test results for $\kappa_{0\_pred}$ values estimated with linear and nonlinear datasets defined by AP3 (which includes records categorized as transitional into the nonlinear dataset).</caption><thead><tr><th rowspan="2">Station</th><th rowspan="2">V<sub>s30</sub> (m/s)</th><th>p-value</th><th rowspan="2">Statistically different</th></tr><tr><th>Linear κ<sub>0</sub><sub>pred</sub> vs. nonlinear κ<sub>0</sub><sub>pred</sub>*</th></tr></thead><tbody><tr><td>AICH17</td><td>314</td><td>10.37%</td><td>--</td></tr><tr><td>CHBH13</td><td>235</td><td>&lt;0.01%</td><td>Yes</td></tr><tr><td>FKSH11</td><td>240</td><td>1.07%</td><td>Yes</td></tr><tr><td>FKSH14</td><td>237</td><td>&lt;0.01%</td><td>Yes</td></tr><tr><td>FKSH18</td><td>307</td><td>2.36%</td><td>Yes</td></tr><tr><td>FKSH19</td><td>338</td><td>57.66%</td><td>--</td></tr><tr><td>FKSH21</td><td>365</td><td>47.12%</td><td>--</td></tr><tr><td>IBRH16</td><td>626</td><td>13.80%</td><td>--</td></tr><tr><td>IBRH17</td><td>301</td><td>55.92%</td><td>--</td></tr><tr><td>IBRH20</td><td>244</td><td>5.97%</td><td>--</td></tr><tr><td>IWTH21</td><td>521</td><td>9.23%</td><td>--</td></tr><tr><td>IWTH26</td><td>371</td><td>0.62%</td><td>Yes</td></tr><tr><td>KMMH01</td><td>575</td><td>0.32%</td><td>Yes</td></tr><tr><td>KMMH12</td><td>410</td><td>3.38%</td><td>Yes</td></tr><tr><td>MYGH07</td><td>366</td><td>0.04%</td><td>Yes</td></tr><tr><td>MYGH10</td><td>348</td><td>7.45%</td><td>--</td></tr><tr><td>NGNH29</td><td>465</td><td>72.74%</td><td>--</td></tr><tr><td>NIGH07</td><td>528</td><td>88.97%</td><td>--</td></tr><tr><td>NIGH12</td><td>553</td><td>66.61%</td><td>--</td></tr><tr><td>TCGH16</td><td>213</td><td>&lt;0.01%</td><td>Yes</td></tr></tbody></table>

<footnote>*The transitional dataset is included into the nonlinear dataset (i.e., datasets following the criteria of AP3).</footnote>
963 **List of Figure Captions:**

964 **Figure 1.** Comparisons between a weak and a strong ground motion recorded at FKSH14 ($V_{s30}$ = 237 m/s). The $M_w$ and $R_c$ are 4 and 15 km for the low-intensity ground motion, while 5.1 and 15 km for the high-intensity ground motion. The frequency window ([10.05, 30] Hz) applied in this plot is picked manually. The left and right columns correspond to analyses conducted on the horizontal components $H_1$ and $H_2$, respectively. The color version of this figure is available only in the electronic version of this article.

970 **Figure 2.** (a) Locations of selected Japanese recording stations in this study, and (b) magnitude and distance distribution of selected ground motions.

972 **Figure 3.** Hyperbolic models fitted to observed $PGA_{rotD50}$ and $I_{\gamma,0}$ data at four study sites. The $V_{s30}$ for IBRH16, IBRH17, IBRH20, and IWTH21 are 626, 301, 244, and 521 m/s, respectively.

974 **Figure 4.** Surface $PGA_{rotD50}$ against $I_{\gamma,0}$ at MYGH10 ($V_{s30}$ = 348 m/s). The red dot-dashed lines present the linear ($I_{\gamma,0,1}$) and transitional ($I_{\gamma,0,t}$) thresholds for $I_{\gamma,0}$. The color version of this figure is available only in the electronic version of this article.

977 **Figure 5.** $G/G_{max}$ versus $I_{\gamma,0}$ at study sites. The $G_{max}$ is computed from the average values of $(\frac{PGA_{rotD50}}{PGV_{rotD50}/V_{s,0}})$ for records with $I_{\gamma,0}$ less than 0.001%. The colors represent $PGA_{rotD50}$ values. The color version of this figure is only available in the electronic version of this article. The color version of this figure is available only in the electronic version of this article.

981 **Figure 6.** Comparisons of individual $\kappa_r$ estimates from our automated algorithm, $\kappa_{r,auto}$, and the fixed-frequency band method, $\kappa_{r\_AS}$ at FKSH14 ($V_{s30}$ = 237 m/s) for surface (left) and borehole (right) records. The color version of this figure is available only in the electronic version of this article.
**Figure 7.** Surface $κ_{r_AS}$ estimates and their corresponding PGA$_{rotD50}$, $I_{γ,0}$ and $R_e$ values for selected ground motions recorded at (a) FKSH14 ($V_{s30}$ = 237 m/s) and (b) MYGH10 ($V_{s30}$ = 348 m/s). Different colors represent varying epicentral distances per record, and the size of markers indicate the corresponding PGA$_{rotD50}$. The color version of this figure is only available in the electronic version of this article. The color version of this figure is available only in the electronic version of this article.

**Figure 8.** $κ_0$-model at FKSH14 ($V_{s30}$ = 237 m/s) from datasets defined by (a) AP1, which only considers the linear and nonlinear datasets, (b) AP2, where transitional records are included as part of the linear dataset, (c) AP3, where transitional records are included as part of the nonlinear dataset, and (d) AP4, where the linear, transitional, and nonlinear datasets are considered separately. The color version of this figure is available only in the electronic version of this article.

**Figure 9.** $κ_0$-model at MYGH10 ($V_{s30}$ = 348 m/s) with datasets defined by (a) AP1, which only considers the linear and nonlinear datasets, (b) AP2, where transitional records are included as part of the linear dataset, (c) AP3, where transitional records are included as part of the nonlinear dataset, and (d) AP4, where the linear, transitional, and nonlinear datasets are considered separately. The color version of this figure is only available in the electronic version of this article. The color version of this figure is available only in the electronic version of this article.

**Figure 10.** Ratio of $κ_{0NL\_sur}/κ_{0LIN\_sur}$ at study sites estimated using the dataset definitions based on AP3 (left panel) and AP4 (right panel) against to $V_{s5}$, $V_{s10}$, and $V_{s30}$. The color version of this figure is available only in the electronic version of this article.

**Figure 11.** Estimated surface $κ_{0\_pred}$ and their corresponding ground shaking intensity and in situ deformation characterized by PGA$_{rotD50}$ and $I_{γ,0}$, respectively at (a) FKSH14 ($V_{s30}$ = 237 m/s) and (b) MYGH10 ($V_{s30}$ = 348 m/s). Both, the color and the size of markers represent varying PGA$_{rotD50}$.
1008 values. Triangles and circles represent the linear and nonlinear datasets defined by AP3. The red dashed lines depict the local regression model based on the $\kappa_{0\_pred}$ and $I_{\gamma,0}$ data. The color version of this figure is only available in the electronic version of this article. The color version of this figure is available only in the electronic version of this article.

1012 **Figure 12.** Observed distribution of $\kappa_{0\_pred}$ at (a) FKSH14 ($V_{s30} = 237$ m/s) and (b) MYGH10 ($V_{s30} = 348$ m/s) for the linear and nonlinear datasets. The red lines depict the theoretical probability density function (PDF) fitted with a Gaussian distribution. The color version of this figure is available only in the electronic version of this article.

1016
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Comparisons between a weak and a strong ground motion recorded at FKSH14 ($V_{s30}$ = 237 m/s). The M<sub>w</sub> and R<sub>e</sub> are 4 and 15 km for the low-intensity ground motion,</figcaption></figure>

1018

1019

1020
while 5.1 and 15 km for the high-intensity ground motion. The frequency window ([10.05, 1021 30] Hz) applied in this plot is picked manually. The left and right columns correspond to 1022 analyses conducted on the horizontal components $H_1$ and $H_2$, respectively. The color 1023 version of this figure is available only in the electronic version of this article.

1025

1026

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

1027
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Locations of selected Japanese recording stations in this study, and (b) magnitude and distance distribution of selected ground motions.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Hyperbolic models fitted to observed PGA<sub>rotD50</sub> and I<sub>γ,0</sub> data at four study sites. The V<sub>s30</sub> for IBRH16, IBRH17, IBRH20, and IWTH21 are 626, 301, 244, and 521 m/s, respectively.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Surface PGA<sub>rotD50</sub> against I<sub>γ,0</sub> at MYGH10 (V<sub>s30</sub> = 348 m/s). The red dot-dashed lines present the linear (I<sub>γ,0,1</sub>) and transitional (I<sub>γ,0,t</sub>) thresholds for I<sub>γ,0</sub>. The color version of this figure is available only in the electronic version of this article.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** G/G<sub>max</sub> versus I<sub>γ,0</sub> at study sites. The G<sub>max</sub> is computed from the average values of (PGA<sub>rotD50</sub>/PGV<sub>rotD50</sub>/V<sub>s,0</sub>) for records with I<sub>γ,0</sub> less than 0.001%. The colors represent PGA<sub>rotD50</sub> values (in units of m/s<sup>2</sup>). The color version of this figure is only available in the electronic version of this article. The color version of this figure is available only in the electronic version of this article.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparisons of individual $\kappa_r$ estimates from our automated algorithm, $\kappa_{r,auto}$, and the fixed frequency band method, $\kappa_{r,AS}$ at FKSH14 ($V_{s30}$ = 237 m/s) for surface (left) and borehole (right) records. The color version of this figure is available only in the electronic version of this article.</figcaption></figure>

1048

1049

1050

1051

1052

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

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Surface $κ_{r\_AS}$ estimates and their corresponding $PGA_{rotD50}$, $I_{γ,0}$ and $R_e$ values for selected ground motions recorded at (a) FKSH14 ($V_{s30}$ = 237 m/s) and (b) MYGH10 ($V_{s30}$ = 348 m/s). Different colors represent varying epicentral distances per record, and the size of markers indicate the corresponding $PGA_{rotD50}$. The color version of this figure is only available in the electronic version of this article. The color version of this figure is available only in the electronic version of this article.</figcaption></figure>
1063

1064

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

1065

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

1066

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

1067

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

1068
**Figure 8.** $\kappa_0$-model at FKSH14 ($V_{s30}$ = 237 m/s) from datasets defined by (a) AP1, which only considers the linear and nonlinear datasets, (b) AP2, where transitional records are included as part of the linear dataset, (c) AP3, where transitional records are included as part of the nonlinear dataset, and (d) AP4, where the linear, transitional, and nonlinear datasets are considered separately. The color version of this figure is only available in the electronic version of this article.

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

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>Figure 9. $\kappa_0$-model at MYGH10 ($V_{s30}$ = 348 m/s) with datasets defined by: (a) AP1, which only considers the linear and nonlinear datasets, (b) AP2, where transitional records are included as part of the linear dataset, (c) AP3, where transitional records are included as part of the nonlinear dataset, and (d) AP4, where the linear, transitional, and nonlinear datasets are considered separately. The color version of this figure is available only in the electronic version of this article.</figcaption></figure>

1079

1080

1081

1082

1083

1084

1085

1086

1087
<figure><img src="image_18.png" /><figcaption>**Figure 10** Ratio of $κ_{0\_nl\_sur}/κ_{0\_lin\_sur}$ at study sites estimated using the dataset definitions based on AP3 (left panel) and AP4 (right panel) against to $V_{s5}$, $V_{s10}$, and $V_{s30}$. The color version of this figure is available only in the electronic version of this article.</figcaption></figure>
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 11.** Estimated surface $\kappa_{0\_pred}$ and their corresponding ground shaking intensity and in situ deformation characterized by $\text{PGA}_{\text{rotD50}}$ and $I_{\gamma,0}$, respectively at (a) FKSH14 ($V_{s30} = 237$ m/s) and (b) MYGH10 ($V_{s30} = 348$ m/s). Both, the color and the size of markers represent varying $\text{PGA}_{\text{rotD50}}$ values. Triangles and circles represent the linear and nonlinear datasets defined by AP3. The red dashed lines depict the local regression model based on the $\kappa_{0\_pred}$ and $I_{\gamma,0}$ data. The color version of this figure is only available in the electronic version of this article. The color version of this figure is available only in the electronic version of this article.</figcaption></figure>
<figure><img src="image_21.png" /></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 12.** Observed distribution of $\kappa_{0\_pred}$ at (a) FKSH14 ($V_{s30}$ = 237 m/s) and (b) MYGH10 ($V_{s30}$ = 348 m/s) for the linear and nonlinear datasets. The red lines depict the theoretical probability density function (PDF) fitted with a Gaussian distribution. The color version of this figure is available only in the electronic version of this article.</figcaption></figure>