Physically Based and Empirical Ground Motion Prediction Equations for Multiple Intensity Measures (PGA, PGV, I a , FIV3, CII, and Maximum Fourier Acceleration Spectra) on Sakhalin Island

: In the present study, empirical attenuation relations for multiple ground motion intensity measures (PGA, PGV, I a , FIV3, CII, and MFAS) were developed for Sakhalin Island (in the far east of Russia). A recorded strong motion dataset was used, making GMPEs applicable in active crustal regions with an earthquake magnitude range of 4–6 and a distance range of up to 150 km. The hypocentral distance was used as a basic distance metric. For the ﬁrst time in the research, an analytical representation of Arias intensity ( I a ) was obtained in the framework of a multi-asperity source model. Asperities are considered as sub-sources of high-frequency incoherent radiation. The physical representation of the attenuation model in our study was based on a stress drop on the asperities and the ratio of the total rupture area to the combined area of asperities. The average stress drop on asperities for the examined earthquakes was approximately 13.4 MPa, and the ratio of the total rupture area to the asperity area was 0.22, which is generally close to similar estimates for crustal earthquakes. The coefﬁcients and statistical scattering of the attenuation models were also analyzed. Moreover, a magnitude scale based on a modiﬁed Arias intensity is proposed in the present study. The new magnitude scale has an explicit physical meaning and is characterized by its simplicity of measurement. It is associated with the acceleration source spectrum level and can be successfully used in early warning systems.


Introduction
The investigation of physically based ground shaking measures and their use to conduct quantitative assessments of seismic hazard and risk factors is an urgent issue of earthquake engineering. The number of scientific articles devoted to the development of ground motion prediction equations increases every year (http://www.gmpe.org.uk, accessed on 1 June 2022); seismic intensity scales are also becoming more diverse in this area of research.
In addition to traditional ground shaking measures that can include simple peak measurements, such as PGA, PGV, or PGD, specialized metrics are also being used, which account for the amplitude and duration of a seismic signal. For example, the authors of [1] developed an attenuation model based on the maximum incremental ground velocity (MIV). The MIV is calculated for each horizontal component by integrating the acceleration record between two consecutive zero crossings. Essentially, this transformation provides an estimate of the area under the seismic pulse, and, therefore, accounts for the amplitude and duration of an individual pulse. The maximization of the new time series consisting of successive values of the areas (integrals) of each pulse presents the desired value of the maximum incremental velocity.
Cumulative absolute ground velocity (CAV) has been proposed in the literature as an instrumental indicator of structural damage (e.g., [2][3][4][5]). The CAV value is determined by the integration of the absolute acceleration over the whole seismic record. In [6], strong ground motion attenuation models were developed to obtain the CAV, its modification CAV 5 (when acceleration is above the 5 cm/s 2 threshold), and the standardized absolute velocity (CAVSTD). These models are applicable to shallow crustal and subduction (intraand interplate) earthquakes with magnitudes ranging from 4 to 9 and distances of up to 400 km.
Another example of the practical consideration of incremental velocity is mentioned in [7]. Effective incremental ground velocity (EIGV) employed as a new ground shaking measure was proposed in this article. It was developed to predict the response of a rigid sliding block subjected to oscillations from a half-sinusoidal pulse. EIGV is a function of pulse duration, peak amplitude, and the coefficient of friction occurring between the basement and the sliding block.
The filtered incremental ground velocity (FIV3) is a recently proposed ground shaking measure [8]. An undeniable advantage of FIV3 over single-peak metrics is the measurement of several unipolar peak amplitudes, which should generally contribute to a more accurate prediction of the damage state of a risk element. In [9], the coefficients of the ground motion prediction equations used to perform multiple intensity measurements were analyzed, and it was concluded in the research that the FIV3 attenuation models were characterized by a relatively smaller standard deviation value compared to similar characteristics based on the measurements of a single-peak value.
In addition to engineering-consisted ground shaking metrics, a category of instrumental measures with physical meaning exists in the literature. Arias intensity (I a ) can be used to refer to such measures. It is determined by the integral of the quadratic acceleration. The spectral representation of the source time function according to the ω 2 -Brune model provides an analytical expression for the Arias intensity as a functional dependence on the source-to-site distance, earthquake magnitude, material parameters, and stress parameter [10,11]. The I a metric has been successfully applied in research to assess earthquake-induced landslide hazards (e.g., [12,13]).
Generally, the idea of using Fourier acceleration spectra in engineering seismology has been suggested in the literature for quite a long period of time (e.g., [14]). A new instrumental intensity scale expressed through the average value of the high-frequency acceleration spectra (frequency band: 0.4-13 Hz) has been proposed [15]. The same study showed that macroseismic intensity, according to the modified Mercalli (MM) scale, correlated well with the parameters of the high-frequency spectra, and practically did not depend on regional features.
The development and physical interpretation of Fourier spectral attenuation models is an important part of recent seismological research [16,17]. It should be noted that the empirical characteristics of high-frequency (incoherent) radiation allow us to determine the physical parameters of the source. From the recent studies conducted on earthquake sources and rupture processes, it is known in the literature that high-frequency incoherent radiation is related to slip heterogeneity. Asperities are recognized as regions located on the fault that have large slips relative to the average slip of the rupture area. The asperity areas, as well as the entire rupture areas, scale with the total seismic moment. The size and stress drop of asperity are significant parameters for predicting strong ground motion parameters. These parameters are used as the input data in the technique of the stochastic simulation of seismic records (e.g., [18]). Therefore, the use of Fourier acceleration spectra simultaneously as engineering and physical measures is a practical method.
In addition to a variety of ground shaking measures utilized in the field, there are also various characteristics of the distance between the source and the site. Hypocentral distance (Rhyp) is quite suitable for point-source models, while other measures, such as Rrup (the closest distance between the site and the three-dimensional rupture plane) and Rjb (the shortest distance from a site to the surface projection of the rupture plane), are most commonly used for finite sources. The role of high-frequency radiation is also reflected in modified distance metrics (e.g., [19]).
Therefore, different strong ground motion metrics are used in various applied scenarios. For example, PGA/PGV are most commonly used to perform seismic hazard assessments. Arias intensity (I a ) is most commonly used to conduct seismic risk assessments, including geotechnical and structural applications. The use of the FIV3 metric has recently been proposed in the literature. It is used consistently with the potential-damageassessment method. Although this metric is not as popular as PGA, PGV, or I a , in our study we use FIV3 for empirical-testing purposes and for the comparison of ground motion equations with the most commonly used strong ground motion metrics. Additionally, we examine high-frequency spectral measurements, such as high-frequency Arias intensity and maximum Fourier acceleration spectra.
In this paper, regional attenuation equations based on the measures of PGA, PGV, I a , FIV3, CII, and maximum Fourier acceleration spectra (MFAS) are developed. A recorded strong motion dataset and felt reports are also used. The study area is considered an active crustal region. We analyzed the statistical scattering of ground motion prediction equations. A physical interpretation of the coefficients is presented within the framework of the multi-asperity source model. Attenuation equations for the multiple intensity measures are developed for the first time for the target area.

Study Area and Tectonic Settings
The study area is located on Sakhalin Island and the surrounding shelf (Russia). The island is located in the northwestern part of the Pacific Ring of Fire, which is a continent-toocean transition zone with a powerful system of meridional faults ( Figure 1). The boundary of the Okhotsk (or North American) and Amur (or Eurasian) plates is identified by these faults. The plate convergence velocity is 10-13 mm/year [20].
The interaction of the Amur and Okhotsk plates generates crustal seismicity, both at the boundaries and the periphery ( Figure 1). Strong earthquakes occur quite frequently at this location. The largest earthquake (Mw 7.3) occurred in 1971 in the south of the island, far from settlements. However, the 1995 Neftegorsk earthquake (Mw 7) was the most catastrophic event and killed more than 2000 people (in the north of the island).
The study area is located in the active crustal region with shallow seismicity (source depth up to 30 km), according to the classification proposed in [21]. The focal mechanisms are most commonly reverse or thrust faults, and, more rarely, strike-slip [22]. Geographical location of the study area, major seismic events, tectonic boundaries, active faults, instrumental sites, and epicenters of the earthquakes considered in the study, AM-Amurian Plate, OK-Okhotsk Plate, black arrows indicate the principal axes of compression.

Strong Motion Database
The strong motion recordings were collected as a result of long-term seismic monitoring and routine earthquake data-processing techniques. The duration of the recordings did not exceed three minutes. All the recordings were instrumentally corrected to an accelerogram measured in nm/s 2 .
Crustal earthquakes with source depths of up to 30 km were considered in the present study (Figure 2). The magnitude ranged from ML 4 to 6.1. We selected records with epicentral distances up to 150 km in order to parameterize the geometric spreading of body waves by a term close to 1/ . Figure 1. Geographical location of the study area, major seismic events, tectonic boundaries, active faults, instrumental sites, and epicenters of the earthquakes considered in the study, AM-Amurian Plate, OK-Okhotsk Plate, black arrows indicate the principal axes of compression.

Strong Motion Database
The strong motion recordings were collected as a result of long-term seismic monitoring and routine earthquake data-processing techniques. The duration of the recordings did not exceed three minutes. All the recordings were instrumentally corrected to an accelerogram measured in nm/s 2 .
Crustal earthquakes with source depths of up to 30 km were considered in the present study ( Figure 2). The magnitude ranged from ML 4 to 6.1. We selected records with epicentral distances up to 150 km in order to parameterize the geometric spreading of body waves by a term close to 1/r.
All recordings were preprocessed. The corresponding intensity measures were selected from each station component. All selected values were placed into the database (see Data Availability Statement). Figure 2 presents the distribution of the selected data on a magnitude-distance plot. The total numbers of recordings, earthquakes, and stations in the dataset were 285, 28, and 22, respectively. Local and moment magnitudes were tested in the attenuation relations. The calibration equation used to calculate the local magnitude is presented in [23]. Moment magnitude values were imported from earthquake data centers: USGS (https://www.usgs.gov/, accessed on 1 June 2022) and GFZ (https://www.gfz-potsdam.de/en/home/, accessed on 1 June 2022). If the moment magnitude was unavailable in the earthquake list, the magnitude conversation to the Mw scale was applied using global intermagnitude relations [24]. This was achieved by the polynomial approximation of the cubic-type equation: Almost all of the instrumental sites were classified as soil category D according to the NEHRP classification [25]. The average measured S-wave velocity in the upper 30 m ground layer (VS30) was 300-400 m/s (see Data Availability Statement). In our study, we proposed using a site term from Japanese GMPEs that was characterized by the similar reference VS30.

Felt Reports
In addition to the instrumental ground shaking measures utilized in the research, the community Internet intensity (CII) derived from felt reports was considered in the present study. For this purpose, questionnaires in Russian were used in accordance with the regionalized technique DYFI? [26,27]. Felt reports were collected through the webpage of the regional seismological service eqalert.ru or its mobile applications (see Data Availability Statement).
Between 1 August 2016 and 30 April 2022 771 questionnaires were completed, with over 600 responses reported by Sakhalin respondents. The magnitude versus distance All recordings were preprocessed. The corresponding intensity measures were selected from each station component. All selected values were placed into the database (see Data Availability Statement). Figure 2 presents the distribution of the selected data on a magnitude-distance plot. The total numbers of recordings, earthquakes, and stations in the dataset were 285, 28, and 22, respectively.
Local and moment magnitudes were tested in the attenuation relations. The calibration equation used to calculate the local magnitude is presented in [23]. Moment magnitude values were imported from earthquake data centers: USGS (https://www.usgs.gov/, accessed on 1 June 2022) and GFZ (https://www.gfz-potsdam.de/en/home/, accessed on 1 June 2022). If the moment magnitude was unavailable in the earthquake list, the magnitude conversation to the Mw scale was applied using global intermagnitude relations [24]. This was achieved by the polynomial approximation of the cubic-type equation: Almost all of the instrumental sites were classified as soil category D according to the NEHRP classification [25]. The average measured S-wave velocity in the upper 30 m ground layer (V S30 ) was 300-400 m/s (see Data Availability Statement). In our study, we proposed using a site term from Japanese GMPEs that was characterized by the similar reference V S30 .

Felt Reports
In addition to the instrumental ground shaking measures utilized in the research, the community Internet intensity (CII) derived from felt reports was considered in the present study. For this purpose, questionnaires in Russian were used in accordance with the regionalized technique DYFI? [26,27]. Felt reports were collected through the webpage of the regional seismological service eqalert.ru or its mobile applications (see Data Availability Statement).
Between 1 August 2016 and 30 April 2022 771 questionnaires were completed, with over 600 responses reported by Sakhalin respondents. The magnitude versus distance distribution of the processed macroseismic intensity is presented in Figure 2. The figure indicates that large earthquakes with ML 5-6 magnitudes are felt at distances up to 200 km from the epicenter. The total amount of CII measures is 131. The number of CIIs with two or more felt reports is 83.

Intensity Measures
For the first time in the research, we considered the maximum Fourier acceleration spectra (MFAS) as an intensity measure. We also introduced a modified Arias intensity (see Table 1). The physical meaning of these metrics is discussed below.
We calculated the Fourier amplitude spectrum of the horizontal-component acceleration combined with the quadratic mean value. Such spectral representations are closely related to the amplitudes of SH-wave. The maximum value of the amplitude spectrum was referred to as the MFAS metric. Arias intensity was defined as the sum of all the squared acceleration values obtained from the horizontal records. For further details, see Table 1.
For the other measurements, the maximum values of the two horizontal (NS and EW) station components were selected. The filtered incremental ground velocity (FIV3) was computed for several natural oscillation periods: 0.01, 0.2, 1, and 3 s. Table 1 summarizes the intensity measures, for which regional attenuation relationships were developed, and the corresponding transformations over signals of the original accelerogram a(t), where: T is the duration of record; T 0 is the natural period of oscillation; α is the parameter (0.7, according to ( [9]); Ω is the spectral domain; g is the acceleration of gravity (9.81 m/s 2 ); LP1 is the notation of a second-order low-pass Butterworth filter with a cut-off frequency of 1 Hz; HP1 is the notation of a second-order high-pass Butterworth filter with a cut-off frequency f H = 1 Hz; HP3 is the notation of a second-order high-pass Butterworth filter with a cut-off frequency f H = 3 Hz; NS and EW are the north-south and east-west directions of the horizontal recordings, respectively.
As an example of ground shaking intensity measurements, we considered a moderate earthquake with an ML 4.5 magnitude and a source depth of 8.6 km, which occurred on 18 April 2021 on Sakhalin Island. Several buildings closest to the epicenter settlement (Rhyp = 9.5 km) were lightly damaged. The macroseismic Internet intensity level derived from the felt reports was approximately CII = 6.3. The accelerograph located in this settlement recorded strong ground motions. The original records and transformations over signals are presented in Figure 3. The corresponding quadratic acceleration spectrum is shown in Figure 4. As an example of ground shaking intensity measurements, we considered a moderate earthquake with an ML 4.5 magnitude and a source depth of 8.6 km, which occurred on 18 April 2021 on Sakhalin Island. Several buildings closest to the epicenter settlement (Rhyp = 9.5 km) were lightly damaged. The macroseismic Internet intensity level derived from the felt reports was approximately CII = 6.3. The accelerograph located in this settlement recorded strong ground motions. The original records and transformations over signals are presented in Figure 3. The corresponding quadratic acceleration spectrum is shown in Figure 4.  The felt reports in each locality were reduced to the community weighted sum (CWS), taking into account various indicators of ground motion activity in accordance with [26]: where Felt is the "felt" index (values from 0 to 1); Mo is the "motion" index (values from 0 to 5); Re is the "reaction" index (values from 0 to 5); St is the "stand" index (values from 0 to 1); Sh is the "shelf" index (values from 0 to 3); Pic is the "picture" index (values from 0 to 1); Fur is the "furniture" index (values from 0 to 1); and Dam is the "damage" index (values from 0 to 3). The felt reports in each locality were reduced to the community weighted sum (CWS), taking into account various indicators of ground motion activity in accordance with [26]: where Felt is the "felt" index (values from 0 to 1); Mo is the "motion" index (values from 0 to 5); Re is the "reaction" index (values from 0 to 5); St is the "stand" index (values from 0 to 1); Sh is the "shelf" index (values from 0 to 3); Pic is the "picture" index (values from 0 to 1); Fur is the "furniture" index (values from 0 to 1); and Dam is the "damage" index (values from 0 to 3).
Then, the decimal value of the community Internet intensity (CII) is calculated and rounded to the first decimal place (if CWS ≥ 6.53): If CWS < 6.53, but the earthquake was felt (Felt > 0), then CII = 2 is obtained; if it was not felt (Felt = 0), then CII = 1.

Attenuation Model
The study area was comparable to a neighboring seismically active region (Japan) in terms of the tectonics and soil conditions. The reference value of the S-wave velocity in a 30 m soil layer was 350 m/s [28,29]. This was closely related to the soil conditions in Sakhalin Island.
The point-source model is most commonly used in the field for seismic sources where Mw < 5.5. We considered earthquakes with magnitude values up to Mw 5.8. The nearest instrumental site for an Mw 5.8 earthquake is located at the hypocentral distance of 21.5 km. Therefore, we proposed that the finite-fault effects for distances greater than 21.5 km were negligible. The point-source approximation seemed to be reasonable for a given magnitude and the distance ranges. Therefore, the following attenuation model was considered: Then, the decimal value of the community Internet intensity (CII) is calculated and rounded to the first decimal place (if CWS ≥ 6.53): If CWS < 6.53, but the earthquake was felt (Felt > 0), then CII = 2 is obtained; if it was not felt (Felt = 0), then CII = 1.

Attenuation Model
The study area was comparable to a neighboring seismically active region (Japan) in terms of the tectonics and soil conditions. The reference value of the S-wave velocity in a 30 m soil layer was 350 m/s [28,29]. This was closely related to the soil conditions in Sakhalin Island.
The point-source model is most commonly used in the field for seismic sources where Mw < 5.5. We considered earthquakes with magnitude values up to Mw 5.8. The nearest instrumental site for an Mw 5.8 earthquake is located at the hypocentral distance of 21.5 km. Therefore, we proposed that the finite-fault effects for distances greater than 21.5 km were negligible. The point-source approximation seemed to be reasonable for a given magnitude and the distance ranges. Therefore, the following attenuation model was considered: where I M ij is the intensity measure (macroseismic intensity or logarithm of instrumental measure) of the i-th earthquake at the j-th seismic station; a, k, b, and c are the regression parameters to be estimated; σ is the standard deviation of the empirical attenuation model; M i is the magnitude of the i-th earthquake; R ij is the hypocentral distance from the i-th source to the j-th seismic station; and ε ij is a N(0, 1) normally distributed random variable with zero mean and unit variance. For the last 25 years mixed-effects regressions have been the standard approach in the research used for GMPEs to account for correlations in the data. The alternative approach that is widely used in Japan is the two-step stratified regression analysis method (e.g., [30,31]). This method was proposed by Joyner and Boore [32] to avoid the interaction between the coefficients a and k. In other words, errors produced when measuring earth-quake magnitudes would affect the distance coefficient obtained from the ordinary least squares regression. In recent decades, the accuracy of moment magnitude determination has significantly increased. The ordinary least squares regression method seemed to be appropriate for this dataset. Therefore, the best coefficients of the model (4) were fitted by multiple linear regression. We analyzed the interaction between the magnitude and distance terms by fixing the coefficients a and/or k at their physical values (see Section 3.1). The standard deviation (σ) was considered as the basic metric of the statistical scatter. Additionally, the coefficient of determination (R 2 ) was used for each model representation. The statistical parameter R 2 is a metric of the quality of fit between the predicted and observed data, and it is expressed in relative units.

Physical Representation of High-Frequency Measures
The idea of developing physically based GMPEs was proposed in [33]. Hanks and McGuire used Parseval's theorem, which states that the integral of the square of a function is equal to the integral of the square of its Fourier transform.
In [10], a source-related acceleration spectrum was modeled by the simplified bandlimited function. It obtains constant values between a lower frequency, related to the source size, and an upper bounding frequency, related to the upper-crustal attenuation. The source spectrum was proposed to be zero outside the mentioned frequency range. Geometrical spreading was considered as a major attenuation factor. Anelastic attenuation was neglected. As a result, semi-empirical relations between the Arias intensity, moment magnitude, source-to-site distance, and static stress drop were derived.
In their study, Stafford et al. [11] considered a more general case where the earthquake acceleration spectrum was given by the multiplication of the attenuation and source functions. The authors proposed a ω 2 -Brune spectrum in the framework of a single-crack source model. The attenuation was described by geometrical spreading and the frequencyindependent Q-factor.
In our study, we introduced an improved analytical predictive equation to obtain the Arias intensity. We collected all the attenuation functions, such as frequency-dependent anelastic attenuation, upper-crustal attenuation, the upper-crustal amplification factor, and geometrical spreading attenuation. The novelty of the presented approach was that high-frequency radiation was described in terms of the multi-asperity source model.
Considering the fragmented structure of the earthquake source, many researchers (e.g., [18,34,35] and others) concluded that there were spots located on the fault plane characterized by increased values of the strength. Such spots are referred to asperities. Asperities are relatively small and occupy 10-30% of the total rupture area. They are subsources of incoherent high-frequency radiation on small space-time scales. The average stress drop on asperities is 10 (sometimes more) times greater than the average stress drop over the fault, which is approximately 3 MPa.
Once the largest asperity has been determined, it is possible to estimate the stress drop of the same asperity (∆τ) using the flat level of the acceleration spectra (A 0 ) given by the following equation [34]: where ρ and β are the density and shear-wave velocity in the vicinity of the source, respectively, and S a is the area of asperity. Generally, the amplitude acceleration spectrum of the S-wave can be written as where f is the frequency, As( f ) is the source-related acceleration spectrum (the shape of the spectrum), G(R) is the attenuation attributed to the geometrical spreading, Ae( f , R) is the anelastic whole-path attenuation factor, Am( f ) is the upper-crustal amplification factor, and An( f ) is the upper-crustal attenuation factor.
Recent studies conducted on earthquake spectra corrected for the attenuation functions showed that their shape was complicated by several corner frequencies [36]. However, the approximation of the observed spectra by the ω 2 model seems to be quite reasonable (e.g., [37,38]). Therefore, to simplify the theoretical calculations, the shape of the sourcerelated acceleration spectrum is presented by the following equation: where f c is the corner frequency associated with the source dimension.
For distances that do not exceed 150 km, the geometrical spreading function G(R) can be written as where R 0 is a reference distance, usually set as being equal to 1 km, and ζ is the indicator of the geometrical spreading.
The anelastic whole-path attenuation factor is presented by the following equation: where Q( f ) is the quality factor as a function of frequency f . In order to determine the necessary attenuation characteristics, the study area was examined in the framework of the single back-scattering model [39]. The quality factor Q( f ) was determined by where Q 0 is the constant. In general, a functional dependence in (10) was typical for active crustal regions. For example, a similar dependence was obtained in Japan for the crustal type of seismicity (e.g., [40]).
The amplification factor was calculated by solutions that accounted for the seismic energy conversation in the upper crust. According to [37], Am( f ) is determined by the square root of the impedance ratio between the source and surface, that is, where ρ 1 and β 1 are the time-weighted average values from the surface to a depth equivalent to a quarter wavelength density and shear-wave velocity, respectively. The upper-crustal attenuation factor, which represents the attenuation of waves in the upper 4 km of the Earth's crust, is defined by a generalization of Equation (9): where κ is the "kappa" parameter. By substituting (10) into (9) and then (5), (7)-(9), and (12) into (6), we obtained an analytical representation of the amplitude spectrum of the accelerogram containing the S-wave: where F is the free surface factor (usually set as being equal to 2) and V is the factor partitioning energy into two orthogonal directions (usually set as being equal to 1/ √ 2).
Equation (13) accounts for the filling factor k f given by the ratio of the total rupture area (S) to the combined area of asperities (S a ): According to Parseval's theorem, The seismic record is defined by the time interval [0, T], where the signal at each moment in time is produced by a real value. Formally, we considered that, outside the specified time interval, the signal was zero. Then, accounting for the conditions, and we can rewrite (15) as Accounting for (18), the Arias intensity can be written as Furthermore, by substituting (13) into (19), following the elementary transformations, we obtained the subsequent expression: where We considered the frequency-independent amplification factor denoted by Am.
In their study, Stafford et al. [11] showed that the parameter Ψ was close to 1 and possible deviations from 1 did not exceed 20%. We numerically integrated the expression presented in (21) and entered the results in a graph depending on λ (see Figure 5). At relatively small values of the variable λ, parameter Ψ was close to 1. For example, for κ = 0.03 and f c = 1, parameter Ψ was equal to 0.698 (λ = 0.188). Therefore, at f c < 1 (which corresponds to the considered magnitude range), this approximation seemed to be quite reasonable.
Empirical relationships between earthquake magnitude and source parameters were used in this study [18,35,41,42]: where M 0 is the seismic moment, Mw is the moment magnitude, and ∆σ is the stress drop of the circular source; the parameters k f and ∆τ are most commonly referred to as inner fault parameters, while M 0 and ∆σ are referred to as the outer characteristics of the source. Empirical relationships between earthquake magnitude and source parameters were used in this study [18,35,41,42]: lg 0 = 1.5 lg + lg ∆ + lg 16 where 0 is the seismic moment, Mw is the moment magnitude, and ∆ is the stress drop of the circular source; the parameters and ∆ are most commonly referred to as inner fault parameters, while 0 and ∆ are referred to as the outer characteristics of the source.
By taking the logarithm of (20), accounting for (23)-(25), we obtained an analytical attenuation equation for the Ia metric: where By comparing (4) and (26), we concluded that the produced attenuation model was identical to the physical representation of the intensity Ia; the coefficients of Model (4) corresponded to the following physical constants obtained from (26)  By taking the logarithm of (20), accounting for (23)-(25), we obtained an analytical attenuation equation for the I a metric: where By comparing (4) and (26), we concluded that the produced attenuation model was identical to the physical representation of the intensity I a ; the coefficients of Model (4) corresponded to the following physical constants obtained from (26) In [43], the empirical relation between V S30 and κ was developed using instrumental measures obtained from the dense seismic network KiK-net in Japan. For the range 155 < V S30 < 2000 m/s, the "kappa" parameter monotonically decreases with the increasing V S30 value. Therefore, the global (or regional) V S30 maps in combination with the abovementioned empirical relation can be used to parameterize the "kappa" attenuation within the physical representation of the attenuation equation.
For the modified metrics I a (1) and I a (3), according to their definitions presented in Table 1, the magnitude-distance relations will be similar to (26); the constant c 0 in this case will be produced by the following expression: where f H = 1 Hz or f H = 3 Hz. Expression (30) was obtained under the assumption ( f H > f c ), where the source acceleration spectrum (7) was presented by the flat level.
The spectral metric MFAS has a physical meaning similar to that of I a . Gusev [35] considered a multi-asperity source model in his study, assuming that the ratio of the total rupture area to the combined area of asperities was constant. He obtained an analytical expression where the maximum value of the acceleration spectrum, corrected for the attenuation functions and geometrical spreading, was produced by the functional dependence on the source radius and the stress drop on asperities; that is, Therefore, the high-frequency metrics I a and MFAS were dependent on inner and outer fault parameters. These parameters were more closely related to the generation of strong motions occurring in the near-source zone. Once these parameters were determined in the target region, it allowed us to follow a method proposed by Irikura [18] in order to successfully obtain stochastically simulated ground motion parameters. Table 2 presents a list of material parameters and physical constants that were additionally used in the calculations; they were mainly obtained from [39].

Performance Metrics of Attenuation Models
For the analysis of the performance metrics of GMPEs we used several analytical representations of the attenuation model (4): (I) statistical representation-simultaneous determination of the best model parameters; (II) physical representation-coefficient a was fixed (=1 for I a and =0.5 for all other ground shaking measures); and (III) geometric representation-coefficient k was fixed (=2 for I a and =1 for all other ground shaking measures). The coefficients of the GMPE for seismic intensity (CII) were estimated only in the framework of the statistical representation. The coefficients of model (4) were calculated for different types of magnitude (Mw and ML). If coefficient b took negative values, it was fixed to zero and the remaining coefficients were recalculated.
The results of the statistical fitting are illustrated in Figure 6. The original data and best attenuation curves according to Equation (4) within the statistical framework are presented in the figure. The coefficients for each model representation and the performance metrics for different types of magnitudes are presented in Table 3a-c. Statistical scattering was selected as the basic performance metric.      The high-frequency metric (MFAS) was characterized by the lowest scattering of GMPE in the framework of the statistical representation (see Table 3a-c). The attenuation equations based on the FIV3 metric were also characterized by a relatively reduced scatter pattern, which generally reproduces the conclusions obtained from [9]. This characteristic appeared to be the best predictor of seismic risks. We suggested that an even better predictive power of damage state from the physical perspective would be presented by a metric that not only considered the unipolar amplitudes of ground shaking in the low-frequency band, but also high-frequency radiation, represented in the time domain by short pulses.
The scattering of the physical models was only 20-30% higher than the same performance metric for the statistical representations of the attenuation model. This meant that physically based GMPEs could be applied in regions with relatively low instrumental measurements.
The small increase in the standard deviation for constraining either the magnitude (a) or geometrical spreading (k) coefficients shows that these two factors are strongly correlated in the dataset. If one factor is constrained, there is a tradeoff with the other. To test the physical model, both a and k were fixed at their physical values. From Table 3d, it can be observed that the statistical scattering changes insignificantly in this case. Value k is slightly dependent on value a, and if coefficient a decreases, then k also decreases, and vice versa. This indicates the non-zero correlation between coefficients a and k.
The attenuation of community Internet intensity (CII) derived from the felt reports is characterized by relatively high scattering (see Table 3a). The statistical scattering was significantly reduced by using the CII values derived from two or more felt reports.

Average Inner-Fault Parameters and the Physical Interpretation
We used Equation (29) and the empirical constants for the I a metric to estimate the average stress drop on asperities in the examined area. By substituting the constants obtained from Table 2 and coefficient c (see Table 3b) in Equation (29), we obtained the estimates of stress drop on asperity (∆τ) as a function of the filling factor k f . The given functional dependence is presented in Figure 7. The figure shows a graph of the function ∆σ k f by a simple calculation performed according to Equation (25). Estimates of the stress drop (∆σ) of local earthquakes in the study area were performed in [39]; the median value of the stress parameter was approximately 3 MPa. This reference value was used to estimate the inner-fault parameters k f and ∆τ (see Figure 7). Modified Ia metrics, proposed in this study, were based on the prefiltering of acceleration records with high-pass filters. By analyzing the empirical attenuation equations based on the Ia metric and its modified measures, we found that the coefficient associated with the magnitude term in GMPE decreased and tended to become unity when the cutoff frequency of the filter increased (see Table 3b). According to the theoretical model (see Equation (26)), this coefficient was strictly unit. We attributed the observed deviations from the theory to a violation of the simplified source spectral model with one corner frequency. Recent studies conducted on earthquake source spectra indicate that the shape of acceleration spectra is presented by several corner frequencies [36].
An increase in the statistical scattering behavior can be observed when the cut-off frequency of the HP filter increases. A possible explanation for the observed phenomena lies in the uncertainty of the parameter . Varying Equation (26) by variable produces Varying the attenuation equation based on the modified Ia metric leads to the following expression (see Equation (30)): The analysis of Expressions (32) and (33) shows that the absolute error | lg | increases as the cut-off frequency of the high-pass filter increases, which is in good agreement with the empirical data.

Acceleration Source Spectral Level
The modified Ia(3) metric is defined by the quadratic of the acceleration spectral level (see Equation (19)). Moreover, the coefficient associated with the magnitude term in the statistical representation of the Ia(3) empirical model was close to 1. It follows that the modified Arias intensity can be used to calibrate the new magnitude scale.
The amplitude spectrum of the accelerogram after applying a high-pass filter with a cut-off frequency of = 3 Hz can be written in the following form: A review of this problem regarding the stress drop on asperities showed that the most common estimates were provided for inland crustal earthquakes occurring in Japan. According to [44], the average stress drops were 21.2 ± 9.2, 13.3 ± 5.3, and 18.0 ± 8.6 MPa for reverse, strike-slip, and all types of faults, respectively. Taking into account the fact that the source mechanisms of Sakhalin earthquakes are predominantly reverse and strike-slip faults [22], the value of 13.3 MPa seems to be reasonable for comparison with the same value presented in the present study (13.4 MPa). The comparison clearly showed that the stress drop on asperities for the examined earthquakes was generally close to similar estimates for crustal strike-slip earthquakes occurring in Japan.
Therefore, the inland area of Japan may be proposed as a master region for selecting seismic records and further developing the regional ground-motion-prediction equations at a wide range of magnitudes.
Modified I a metrics, proposed in this study, were based on the prefiltering of acceleration records with high-pass filters. By analyzing the empirical attenuation equations based on the I a metric and its modified measures, we found that the coefficient associated with the magnitude term in GMPE decreased and tended to become unity when the cut-off frequency of the filter increased (see Table 3b). According to the theoretical model (see Equation (26)), this coefficient was strictly unit. We attributed the observed deviations from the theory to a violation of the simplified source spectral model with one corner frequency. Recent studies conducted on earthquake source spectra indicate that the shape of acceleration spectra is presented by several corner frequencies [36].
An increase in the statistical scattering behavior can be observed when the cut-off frequency of the HP filter increases. A possible explanation for the observed phenomena lies in the uncertainty of the parameter κ. Varying Equation (26) by variable κ produces Varying the attenuation equation based on the modified I a metric leads to the following expression (see Equation (30)): The analysis of Expressions (32) and (33) shows that the absolute error δlgI a f H increases as the cut-off frequency f H of the high-pass filter increases, which is in good agreement with the empirical data.

Acceleration Source Spectral Level
The modified I a (3) metric is defined by the quadratic of the acceleration spectral level (see Equation (19)). Moreover, the coefficient associated with the magnitude term in the statistical representation of the I a (3) empirical model was close to 1. It follows that the modified Arias intensity can be used to calibrate the new magnitude scale.
The amplitude spectrum of the accelerogram after applying a high-pass filter with a cut-off frequency of f H = 3 Hz can be written in the following form: Then, the modified Arias intensity according to (19) is given by the following expression: Equation (35) accounts for the term 4πρβ 3 A 0 , which is denoted as A. Then, the corresponding attenuation equation can be written as By comparing Expression (36) and the equation presented above, we obtain the following expression: Therefore, according to Equation (38), the I a (3) metric can be used to estimate the highfrequency acceleration level based on the empirical coefficients of the attenuation relations.
In recent decades, the high-frequency acceleration spectra have been investigated in different seismically active regions [45,46]. We summarized all the available data and plotted the empirical relationship between the moment magnitude (Mw) and acceleration source spectral level (A) for different regions and tectonic types (Figure 8). The Relation (38) after the substitution of the constants (see Figure 8) was in good agreement with the empirical data and theoretical predictions (e.g., [18]). coverage and other technical details. Therefore, the rapid and accurate estimation of an earthquake's magnitude is still an urgent issue in the field.
In seismic hazard assessments, the moment magnitude is more suitable among the other magnitude scales because it directly addresses the rupture size. The alternative approach is to use Mw as a basic measure of the asperity size, while the stress drop parameter is used to calibrate high-frequency amplitudes. This is one of the areas that requires additional research.

Applicability of Spectral Metrics
The results show that the attenuation relations based on the MFAS spectral metric are characterized by the smallest scattering pattern among the considered ground-motion measures. The distribution of the maximum value of a set of random numbers presented a smaller standard deviation than the standard deviation of the random values themselves. The lower standard deviation seemed to be an expected result.
The MFAS metric has one significant disadvantage. The peak amplitude of the spectrum can not only be caused by the asperity radiation, but also by soil amplification. In practice, it is rather problematic to discriminate the contribution of these phenomena. It requires geotechnical studies and/or the investigation of H/V spectral ratios at the site. Figure 9 presents the distribution of the frequencies corresponding to the maximum FAS. The median value is approximately 1.6 Hz. This is closely related to the value of 2.4 Hz proposed in [35]. The variability in the frequency dependance is caused by source and The figure indicates the linear dependence between the logarithm of the high-frequency level of the acceleration source spectrum and the moment magnitude in a wide range of magnitudes, ranging from Mw 4 to 8. This demonstrated the promise of using the proposed I a (3) metric for calibrating the corresponding magnitude, which we denoted as MIa (3). By replacing Mw with MIa(3) in (37) and substituting the corresponding coefficients obtained from Table 3b, we obtained the calibration equation for a new magnitude scale: The advantages of a new earthquake magnitude are as follows: (1) it addresses both the asperity size and stress drop on asperity; (2) the size and stress drop of asperity are significant parameters for predicting strong ground motion parameters; (3) the asperity areas, as well as the entire rupture areas, scale with the total seismic moment (e.g., [18]); (4) from the log-linear dependence between the acceleration source spectral level and moment magnitude, it follows that the spectral-based magnitude does not saturate in the high-value range; (5) the proposed magnitude scale is calibrated in the near-source distance range; and (6) the calibration equation, in terms of point-source distance metrics, is considered to be a better metric, not least because studies have shown that the hypocenter is often located close to areas of relatively large slips, which are considered as sub-sources of high-frequency incoherent radiation (e.g., [47]).
Efficient magnitude determination immediately following the occurrence of significant earthquakes is extremely important for early warning issues. The most commonly used regional and teleseismic magnitudes are saturated. This results in the significant underestimation of earthquake size, potential damage, and tsunami hazard. USGS reports reliable moment magnitude estimations based on the W-phase inversion, 30-60 min after the origin time of a large earthquake. Some approaches propose a shorter time of order of 10 min (e.g., [48]); meanwhile, the efficiency greatly depends on the available network coverage and other technical details. Therefore, the rapid and accurate estimation of an earthquake's magnitude is still an urgent issue in the field.
In seismic hazard assessments, the moment magnitude is more suitable among the other magnitude scales because it directly addresses the rupture size. The alternative approach is to use Mw as a basic measure of the asperity size, while the stress drop parameter is used to calibrate high-frequency amplitudes. This is one of the areas that requires additional research.

Applicability of Spectral Metrics
The results show that the attenuation relations based on the MFAS spectral metric are characterized by the smallest scattering pattern among the considered ground-motion measures. The distribution of the maximum value of a set of random numbers presented a smaller standard deviation than the standard deviation of the random values themselves. The lower standard deviation seemed to be an expected result.
The MFAS metric has one significant disadvantage. The peak amplitude of the spectrum can not only be caused by the asperity radiation, but also by soil amplification. In practice, it is rather problematic to discriminate the contribution of these phenomena. It requires geotechnical studies and/or the investigation of H/V spectral ratios at the site. Figure 9 presents the distribution of the frequencies corresponding to the maximum FAS. The median value is approximately 1.6 Hz. This is closely related to the value of 2.4 Hz proposed in [35]. The variability in the frequency dependance is caused by source and attenuation factors. Therefore, upper-crustal attenuation and the amplification factor should be determined prior to spectral fitting.  At the same time, the modified Arias intensity is an integral characteristic (averaged over the corresponding frequency range in the spectral domain). In our study, we expected that the soil amplification factor would not significantly affect the Ia-based GMPEs. Figure 8 shows that regional differences in the average high-frequency amplitude are present. For example, the Japan intraplate events present systematically higher values than the other regions. A good use of the local data would be to estimate the region-specific high-frequency amplitudes to develop an improved model over the use of a global average model for the stress drop.

Comparison of GMPEs and Finite-Fault Effects
In our study, we compared the resulting equation for PGA with several global and regional models for the crustal seismicity (see Figure 10). The four GMPEs were: AS1997 [49], ASB2013 [47], JSGGA2022 [50], and MF2013_1 [29]. The comparisons are presented on each panel corresponding to the earthquake magnitudes M 4, M 5.0, M 6, and M 7. Attenuation equations were reduced to VS30 = 350 m/s and the reverse faulting mechanism. Models JSGGA2022 and ASB2013 were developed using the Rhyp metric, whereas other GMPEs used Rrup and Rjb. Although GMPEs are guided with different uses in distance metrics, this does not significantly affect the overall results due to the point-source approximation for the considered magnitude range. Figure 10 presents the model-to-model variability caused by regional differences. As can be observed in these figures, the simple point-source model developed in the current At the same time, the modified Arias intensity is an integral characteristic (averaged over the corresponding frequency range in the spectral domain). In our study, we expected that the soil amplification factor would not significantly affect the I a -based GMPEs. Figure 8 shows that regional differences in the average high-frequency amplitude are present. For example, the Japan intraplate events present systematically higher values than the other regions. A good use of the local data would be to estimate the region-specific high-frequency amplitudes to develop an improved model over the use of a global average model for the stress drop.

Comparison of GMPEs and Finite-Fault Effects
In our study, we compared the resulting equation for PGA with several global and regional models for the crustal seismicity (see Figure 10). The four GMPEs were: AS1997 [49], ASB2013 [47], JSGGA2022 [50], and MF2013_1 [29]. The comparisons are presented on each panel corresponding to the earthquake magnitudes M 4, M 5.0, M 6, and M 7. Attenuation equations were reduced to VS30 = 350 m/s and the reverse faulting mechanism. Models JSGGA2022 and ASB2013 were developed using the Rhyp metric, whereas other GMPEs used Rrup and Rjb. Although GMPEs are guided with different uses in distance metrics, this does not significantly affect the overall results due to the point-source approximation for the considered magnitude range. The lineament-domain source model was used in the seismic hazard assessments conducted on Sakhalin Island. The maximum magnitude for area sources did not exceed Mw 6. For such seismic sources, one can use the regional attenuation Equation (40), while for extended line sources the global models should be applied. Figure 11 presents a comparison of the attenuation curves with the data normalized to Mw 7.0. The observed data were normalized as follows: where ′ is the normalized data, is the observed data, (M) is the predicted value by Equation (40) for identical parameters with observed data, and (Mw = 7) is the predicted value for Mw 7.
The PGA values at the near-source distances were close to the flat level at approximately 800 cm/s 2 ( Figure 11). The PGA values converted from macroseismic intensity (https://earthquake.usgs.gov/earthquakes/eventpage/usp0006y50/shakemap/stations?source=atlas&code=usp0006y50, accessed on 7 June 2022) for the catastrophic 1995 Neftegorsk earthquake (Mw 7.1) fit the prediction model well. The hypocentral distance was replaced by the rupture distance.  Figure 10 presents the model-to-model variability caused by regional differences. As can be observed in these figures, the simple point-source model developed in the current study did not extrapolate well to magnitudes of M ≥ 5 at short distances that often control seismic hazard outcomes. To be useful for seismic hazard applications, the GMPE needs to perform extrapolations in a reasonable manner. The finite-fault term added to the distance metric in Equation (4) was one way to force the point-source model to perform extrapolations properly. We imported this term from the model MF2013_1. The MF2013_1 model used for crustal seismicity in Japan was developed in the form of a geometric representation of GMPE. This means that the geometrical spreading term presented in Equation (4) was taken as k = 1. Therefore, we added an additional term to Equation (4). The hybrid model is presently formulated as where d and e are the constants d = 0.006875 and e = 0.5. The parameters of the developed attenuation model in the geometric representation are presented in Table 3c. The corresponding illustration is presented in Figure 10.
The lineament-domain source model was used in the seismic hazard assessments conducted on Sakhalin Island. The maximum magnitude for area sources did not exceed Mw 6. For such seismic sources, one can use the regional attenuation Equation (40), while for extended line sources the global models should be applied. Figure 11 presents a comparison of the attenuation curves with the data normalized to Mw 7.0. The observed data were normalized as follows: where obs is the normalized data, obs is the observed data, pre(M) is the predicted value by Equation (40) for identical parameters with observed data, and pre(Mw = 7) is the predicted value for Mw 7.

Selection of Point-Source Distance Metrics
The choice of distance metrics in ground motion attenuation models is frequently discussed in the scientific community (e.g., https://usgs.github.io/shakemap/man-ual3_5/tg_processing.html, accessed on 1 June 2022). The distance metrics Rrup and Rjb are most commonly used as the characteristic distance between the finite source and the site for seismic hazard assessments. That is why developers of GMPEs often select the abovementioned distance metrics. However, to rapidly generate shaking maps, hypocentral (or epicentral) distances are preferred in the research, since the finite-source characteristics are reported by seismological data centers more slowly than the hypocenter parameters.
An analysis of regional and global empirical attenuation relationships (http://www.gmpe.org.uk, accessed on 1 June 2022) indicates the lack of Rhyp-based GMPEs in active crustal regions. Apparently, the only exception to this practice is the ASB2013 model [47], developed by the authors specifically for European countries, which can be classified as a region with a stable continental crust. In this model, the Rhyp distance metric, as well as Rjb and Rrup, were used. Therefore, the model developed in the current study concerning the empirical attenuation relationships based on the Rhyp distance metric can be used in early warning systems (e.g., [51]) to perform an accurate prediction of the ground shaking level in the study area.

Conclusions
1. Empirical attenuation relations for multiple ground motion intensity measures (PGA, PGV, Ia, FIV3, CII, and MFAS) were developed in the current study. Some of these factors referred to low-frequency metrics (e.g., FIV3), while the Arias (Ia) and spectral The PGA values at the near-source distances were close to the flat level at approximately 800 cm/s 2 ( Figure 11). The PGA values converted from macroseismic intensity (https://earthquake.usgs.gov/earthquakes/eventpage/usp0006y50/shakemap/stations? source=atlas&code=usp0006y50, accessed on 7 June 2022) for the catastrophic 1995 Neftegorsk earthquake (Mw 7.1) fit the prediction model well. The hypocentral distance was replaced by the rupture distance.

Selection of Point-Source Distance Metrics
The choice of distance metrics in ground motion attenuation models is frequently discussed in the scientific community (e.g., https://usgs.github.io/shakemap/manual3 _5/tg_processing.html, accessed on 1 June 2022). The distance metrics Rrup and Rjb are most commonly used as the characteristic distance between the finite source and the site for seismic hazard assessments. That is why developers of GMPEs often select the abovementioned distance metrics. However, to rapidly generate shaking maps, hypocentral (or epicentral) distances are preferred in the research, since the finite-source characteristics are reported by seismological data centers more slowly than the hypocenter parameters.
An analysis of regional and global empirical attenuation relationships (http://www. gmpe.org.uk, accessed on 1 June 2022) indicates the lack of Rhyp-based GMPEs in active crustal regions. Apparently, the only exception to this practice is the ASB2013 model [47], developed by the authors specifically for European countries, which can be classified as a region with a stable continental crust. In this model, the Rhyp distance metric, as well as Rjb and Rrup, were used. Therefore, the model developed in the current study concerning the empirical attenuation relationships based on the Rhyp distance metric can be used in early warning systems (e.g., [51]) to perform an accurate prediction of the ground shaking level in the study area.

Conclusions
1. Empirical attenuation relations for multiple ground motion intensity measures (PGA, PGV, I a , FIV3, CII, and MFAS) were developed in the current study. Some of these factors referred to low-frequency metrics (e.g., FIV3), while the Arias (I a ) and spectral (MFAS) intensities were characteristic of high-frequency incoherent radiation, according to the multi-asperity source model. A recorded strong motion dataset was used in our study, making the GMPEs applicable to an earthquake magnitude range of 4-6 and a distance of up to 150 km. The hypocentral distance was used as a basic distance metric. The target area was the transition zone from the Pacific Ocean to the continental regions of Russia.
2. We proposed a GMPE corresponding to the simple point-source model with the physical value of geometrical spreading for the prediction of PGA in Sakhalin Island and the surrounding shelf. The applicability of the presented model was strictly bound to the magnitude range of 4-6. The reference V S30 took a value of 350 m/s. Site and finite-fault terms describing the PGA model in Japan were preferred for use in the present GMPE.
3. The attenuation relations based on the MFAS spectral metric, as well as a filtered incremental velocity of FIV3(0.01), were characterized by the smallest scattering pattern among the considered ground motion measures. 4. For the first time, an analytical representation of the Arias intensity (I a ) was obtained in the framework of a fragmented source model with asperities, which are sub-sources of high-frequency radiation. This allowed us to use this metric for estimating the average stress drop on asperities occurring in the target region.
5. The average stress drop on asperities for the examined earthquakes was approximately 13.4 MPa, and the ratio of the total rupture area to the combined area of asperity was 0.22, which is generally similar to the estimates for crustal earthquakes occurring in Japan. 6. The modified I a metrics proposed in the current study were based on the prefiltering of acceleration records with high-pass filters. It was shown that the coefficient associated with the magnitude term in the GMPE decreased and tended towards unity when the cut-off frequency of the filter increased. We attributed the observed deviations from the theory to a violation of the source model, according to which the amplitude spectra decayed with the frequency according to the ω 2 -law.
7. When the cut-off frequency of the HP filter increased, an increase in the statistical scattering in the attenuation equation for the modified I a metric was evident. This was explained in terms of the kappa attenuation parameter. 8. A modified I a (3) metric was proposed to calibrate the new magnitude scale associated with the acceleration source spectrum level in the high-frequency range.