Strong motion estimation in Costa Rica at non-record sites using spectral inversion method

We performed the spectral inversion to separate the source, path, and site effects using strong ground motion records from essential facilities such as powerhouses and hydroelectric dams, distributed all over Costa Rica. The spectral inversion technique was adapted to a sparsely observed strong motion data set. We used the results to estimate strong ground motions at sites where there were no observation records in and around the intensity anomaly that appeared during the 2012 Sámara earthquake of Mw 7.6. To ensure the quality, we incorporated a broadband station to estimate levels of acceleration source spectra for six events 5.3 ≤ Mw ≤ 6.5, using them as the reference events for the spectral inversion. In this paper, we introduced two consistency checks: confirmation of the correctness of the input data, and verification of the output using synthetic spectra. We used them as an iterative implementation of sequences of elimination of suspicious data and inversions. We regulated the source spectra of the reference events for the final inversion, considering the estimate of site effect for one-dimensional SH wave propagation at the site of the lowest site amplification. We obtained the frequency-dependent quality factor Q = 179f0.5598 for the northern and central parts of Costa Rica. From the output of the spectral inversion, we reproduced the acceleration spectra of earthquakes in the sites where the event was not recorded. We applied this formulation using the earthquake mentioned above. The resulting synthetic spectra are consistent with the anomalous intensity distribution reported at that time.


Introduction
The subduction of the Cocos and the Nazca plates under the Caribbean plate controls the tectonic frame of Costa Rica. The bathymetric features of the Cocos plate have significant variations from Northwest to Southeast. These features have different origins, such as geomorphological structures and thermal variations (Arroyo et al. 2013) and have influenced the seismogenic process (Morell 2016). The generated stress caused crustal deformation, and formed critical tectonic structures, such as the Central Costa Rica Deformed Belt (CCRDB)which crosses from the Pacific coast to the central part of Costa Rica and reaches the North Panamá Deformed Belt (NPDB) in the Caribbean coast (Fig. 1). It is considered the diffused limit of the Panamá microplate within the Caribbean plate (Montero 2001). The CCRDB also defines the south and east limit of the Great Metropolitan Area, making these fault systems a significant seismic threat in addition to the subduction zone.
The volcanic deposits that cover a significant part of the country are the predominant geological features in Costa Rica. These deposits due to volcanism, along with the deposits of sedimentary basins, can amplify seismic ground motion during a major earthquake. The influence of site amplification was observed during the 2012 Sámara earthquake of 7.6 Mw (National Seismological Network of Costa Rica (RSN)), which is the second biggest earthquake in the instrumental history of Costa Rica. The hypocenter of this event was located 10 km off the Nicoya peninsula, at 13 km depth (Yue et al. 2013). The reported damages were concentrated around the Nicoya peninsula, near the epicenter, and around Naranjo township in the central part of the country. The surface geology of the Naranjo township corresponds to weathered pyroclastic flows, that create significant soil deposits in mountainous areas, a common feature in the country (Red Sismológica Nacional (RSN) 2012). The ground motion caused partial damage to 1990 houses, 38 bridges, 56 schools, and 33 buildings for health services in these places, and caused many landslides that correspond to a maximum intensity of VII on the MMI scale (Red Sismológica Nacional (RSN) 2012).
Therefore, it is required to analyze this type of anomalous ground motion for an appropriate seismic hazard assessment. Thus, it is necessary to understand and quantify the study area's source, path, and site effects separately to quantitatively estimate the strong ground motion for future major events (e.g., Oth et al. 2011).
We based our investigation on the spectral inversion method of Moya and Irikura (2003), which uses reference events instead of a reference site.
We expect this information, especially the site effects, to serve for pre-and post-seismic evaluation of buildings before and after an earthquake and improve the maintenance and reinforcement. In the case of the 2012 Sámara earthquake, we performed a synthetic  Montero (2001). The arrows represent the velocity of subduction, while the white triangles represent the active volcanoes reproduction at sites that were installed in 2015, and do not have any record of the event.

Spectral inversion
Spectral inversion (hereafter SI) has been applied before to separate the source spectra, the attenuation during propagation, and the site amplification (Iwata and Irikura 1988). They applied this method, assuming that the earthquake source is a point source and that the quality factor Q depends only on the frequency and travel distance. Then, in the frequency domain, the spectra of the observed record O ij (f ) can be formulated using the source spectra S i (f ), the site effect G j (f ) that includes the free surface effect, and the hypocentral distance from the i-th event at the j-th site R ij (Iwata and Irikura 1988).
When the distances between the sites are smaller than hypocentral distances, this problem is generally solved by selecting a reference site (Moya and Irikura 2003). However, the reference sites themselves can amplify seismic waves and are not always available depending on the geological condition of the target site. Moya and Irikura (2003) also presented a formulation to solve this equation without needing a reference site. Instead, they proposed using reference events that give constraints based on the ω −2 model of source spectra (Aki 1967). The requirements to implement this method are that the seismic moment and the corner frequency of the reference events should be accurate and that the assumption of ω −2 model is applicable.
One of the controversial topics related to the shape of observed acceleration spectra is their departure from the ω −2 model, specifically their decay, at high frequencies. Anderson and Hough (1984) introduced the parameter κ (kappa) factor, which describes their decaying behavior with a constant slope in the plot of the logarithm of spectra against the frequency on the linear axis (Ktenidou et al. 2014). The cause of this κ-decay, termed exp(−πκf ), has been discussed as a product of a source factor (Papageorgiou and Aki 1983), a site effect (Boore 2003;Oth et al. 2011), or their combination (Tsai and Chen 2000;Purvance and Anderson 2003;Van Houtte et al. 2011;Kilb et al. 2012). Recent comprehensive reviews were given by Ktenidou et al. (2014) and Tsurugi et al. (2020). Also, Oth et al. (2011) clearly showed that, although SI allocates it into the source spectra, it should be related to attenuation within the near-surface weathered layers The path effect includes the attenuation in terms of the quality factor (Q factor), which depends on the frequency and the geometrical spreading effect (Oth et al. 2011). For the former, we employed the often-used formulation Q(f ) = Q 0 f N to describe its behavior, Q 0 is the Q value at 1 Hz, and N represents the dependence on the frequency.
For the source spectra, we also performed a correction of the radiation pattern to simplify the data processing and the analysis. We used the empirical frequency-dependent radiation pattern proposed by Satoh (2002) for the correction. Namely, we used theoretically calculated coefficients for SV and SH components for the frequencies below 2 Hz, the root mean square (RMS) of the radiation coefficient over the sphere √ (2/5) = 0.632 (Andrews 1986) above 5 Hz, and between them their linear interpolation. For example, Takemura et al. (2009) used the threshold frequencies between 2 and 5 Hz for the seismic records in Japan, but no such study exists for Costa Rica, although they should be determined area by area. Since SV/SH values for frequencies above 5 Hz are close to unity (Satoh 2002) (Takemura et al. 2009), the power is evenly allocated to these components. The theoretical radiation coefficients used below 2 Hz can take very small values along their nodal planes and may cause an abrupt behavior of the corrected records. Therefore, we replaced them with the RMS values mentioned above when they become smaller than this value. Hereafter, "the observed spectra denotes the spectra including for the radiation pattern mentioned above unless there is a special comment. All related quantities are based on them.  applied SI to a large-good quality data set of seismic records from small to moderate size events in central Italy to check the evidence of forwarding directivity effects for the Mw 6.1 main shock of the 2009. L'Aquila seismic sequence. For further analysis,  showed that the directivity effect above the source corner frequency diminishes with increasing frequency for small earthquakes. In the present study, we analyze moderate-sized, with a working frequency range that covers from 0.5 Hz to high frequencies. Therefore an ambiguity remains on the influence of source directivity effect on the result of SI. However, the events used in the present study do not contain information of the rupture direction of the seismic source, with the exception of the 2012 Sámara earthquake, 7.6 Mw. Thus, any correction that requires this information could not be conducted. Therefore, the resulting source spectra from SI is averaged in terms of azimuth and take-off angle.
We set the target frequency range of analysis from 1 Hz to 10 Hz for the engineering interest, but the data processing is done in the working frequency range from 0.5 Hz to 10 Hz to expect better results in a little broader frequency range.

Strong motion synthesis in the frequency domain
Once SI results obtained, it is possible to estimate the acceleration spectra of the target record O I J (f ), (i.e., the record at the J-th site during the I-th event),even if there was no ground motion records for the I-th. Two-element records are required together with the results of SI. One is the observed record of the J-th site O iJ (f ). The other, is a record of the I-th event at the j-th site O Ij (f ). From SI, the source effect of the i-th event S i (f ) and the site effect of the j-th site G j (f ) are estimated.
From O Ij (f ), the source spectra of the I-th event S I (f ) can be calculated using the site effect of the j-th site G j (f ) and the Q factor obtained in SI as shown in Eq. 2: Similarly, from O iJ (f ), the site effect of the J-th site G J (f ) can be calculated using the source spectra of the i-th event S i (f ) and the Q factor as shown in Eq. 3: Then combining Eqs. 2 and 3 with Eq. 1, we can get the spectra of the target record O I J (f ) as shown below: It is sufficient that the j-th site and i-th event are involved in it through records O kj (f ) and O il (f ), (k, l) = (I, J ). For example, we can reproduce ground motion of a major event that occurred in the old days O I J (f ) using the earthquake source spectrum obtained from a old dataset at different location O Ij (f ), a recent record O iJ (f ) at the J-th site, and the result of SI, namely, S i (f ), G j (f ), and the Q factor based on the strong motion records of recent years.

Seismic networks
The National Seismological Network (RSN) is a collaborative program between the Seismic and Volcanic Auscultation Area from the Costa Rican Electricity Institute (ICE) and the Central American Geology School from the University of Costa Rica (UCR). For this study, we used strong motion observation sites, 19 REFTEK accelerometers with a sampling rate of 200 Hz and a 24bit A/D converter (Trimble 2014) from the ICE network, and five stations accelerometerseismometer Sixaola 3 V2 (Raspberryshake 2019), installed in 2015 into the UCR network referred to as Sixaola sites (Fig. 2, Tables 2 and 3). The ICE stations are classified according to the Costa Rica's seismic code (Colegio Federado de Ingenieros y Arquitectos de Costa Rica (CFIA) 2010). Each Sixaola station has a sampling rate of 100 Hz, employs 14-bit A/D converter, and has a sensitivity as fine as 2.28E+03 counts/(m/ 2 ), reason why these stations can only record the events where the PGA is larger than a few tens cm/s 2 . Therefore, these five stations were used only for the reproduction of the records during a major event as described in Section 4.4.1.
Also, to constrain the determination of the source spectra of the reference events, we included one broadband station from the Incorporated Research Institutions for Seismology/International Deployment of Accelerometers (IRIS/IDA) (2019). We expect that the records of a broadband station have lower self-noise level in the low frequency range in comparison with strong motion records. For this study, we select the station JTS000, administrated by the Costa Rican Vulcanology and Seismology Observatory (OVSICORI) located in Las Juntas, Guanacaste, which is located at a rock site with low amplification (Záček et al. 2012). The station has a sampling rate of 40 Hz (http://ds.iris.edu/wilber3/find stations/ 10485770). Hereafter, this site will be referred to as JTS000. The final inversion was conducted with 13 stations (Table 3).

Selection of events
We selected events from the RSN seismic catalog, that surrounded Costa Rica, that had magnitudes between 5.1 and 7.6 for a proper wave path coverage. Using this criteria, we initially selected a database of nineteen earthquakes from 2009 to 2019. The information on focal mechanism and the moment magnitude were obtained from the catalogs of United States Geological Survey (USGS) (2019) catalog and the Global Centroid Moment Tensor (CMT) (2019). Due to their low magnitude size, for the events smaller than 5.2M w this information was not available in either catalog. Therefore, we calculated it according to the equation (Kanamori 1977), using the value of magnitude from the RSN catalog. The correction factor of the radiation pattern for SV and SH waves from the focal mechanism were calculated for all the events (Table 2).
From the original data set, we selected 14 events with a magnitude between 5.1 and 7.6, expecting that their corner frequency would be lower than the target frequency range, and that their flat acceleration Location and focal mechanism of the selected major events spectra be within this range. Seven events were excluded during the analysis as described in Section 4.1 (Fig. 3, Tables 2, and 3)

Data processing
From the 19 events, we collected a total of 184 strong motion records from ICE, 213 records from UCR, and 10 broadband records at JTS000. We resampled them at a sampling rate of 100 Hz, and then rotated the waveform data of the horizontal components to the orientation north-south/east-west orientation if needed (NIED 2020).
To extract the main part of the S waves, we applied the Husid plot to the two horizontal components, with a time window between 5% and 95% of the energy release (Husid 1969). We carefully reviewed the output from this automatic process to assure the onset and end of the S waves.
Next, we calculated the Fourier spectral density of acceleration using Fast Fourier Transform. From the analysis, we rejected some records of ICE sites, where the signal to noise ratio (SNR) were not sufficient, log 10 (SNR) < 0.4. The SNR is defined here as the larger PGA of the two horizontal components divided by the noise level (RMS amplitude of the two horizontal components in the time domain at the pre-P waves part of the records).
In this study, we used the Sixaola sites shown in Table 1 only to reproduce the acceleration spectra at sites related to the mentioned intensity anomaly or located above the source of the 2012 Sámara earthquake.

Reference events and preliminary analysis
As mentioned before, we incorporated the broadband station JTS000 because it has a lower noise level than strong motion observation stations and is expected to   Figure 4a shows the observed acceleration spectra of the ten events from the original data set recorded at JTS000. The four events shown by broken curves have a lower limit frequency of the acceleration flat level higher than the target frequency of 1 Hz. Their seismic moments and magnitudes given in Table 2 seem too big to be consistent with these corner frequencies. This implies that they do not follow the ω −2 model and may have an intermediate frequency range of source spectra between the corner frequency of displacement spectra and the acceleration spectra corner frequency. Therefore, we did not use these four events to prevent unnecessary complications because the model spectral shape in this intermediate frequency range, between the corner frequency of displacement and that of the acceleration spectra, not yet well established. The lower limit frequency of the other six events was sufficiently lower than the target range, the reason why the flat level of acceleration may cover the target frequency range. On the other hand, it becomes challenging to determine the corner frequency accurately, due to the relatively low S/N ratio outside the frequency range of analysis. Therefore, we directly estimated the acceleration flat level of these six events, select them as reference events, and reduce the data set into six reference events and the 2012 Sámara earthquake. This is different from the usual application of SI, where both seismic moment and the corner frequency are determined based on an accurate spectral estimate over a broad frequency band both below and above the corner frequency.
In this study, the low S/N ratio of accelerograms, except at station JTS000, does not allow us to take this usual way. Then, we focus mainly on the acceleration flat level.
We conducted the preliminary analysis of SI using event #10 2017 Jacó, 6,5 Mw (Fig. 4b) as the unique reference event and the records of the events shown in Fig. 4a at JTS000 plus eighteen ICE sites, ignoring the stations with anomalous high site effect, to obtain a tentative attenuation model Q(f ) = 243f 0.4206 . The source spectra of the reference event #10 were parameterized by the flat-level A 0 (#10) 7.0E+18 (NM/s 2 ) calculated by A 0 = (2ρf c) 2 Mo using its seismic moment M o =7.87E+18 (NM) given in Table 2 and the corner frequency fc (#10) of 0.15 Hz determined by fitting the w −2 model of source spectral curve to the observed one. Figure 4b compares the observed spectra (solid curve) with the reproduced ones (broken curve) calculated using the given source spectra and the tentative attenuation model in the frequency range lower than 2 Hz for the corner frequency. We considered the above values of A 0 (#10) and f c (#10) just tentative ones and revised them during the analysis using the lowest site effect station described in Section 4.4, because it was not believable that JTS000 would have a surface layer of density r 2,700 kg/m3 and Vs 3,600 m/s, which are values used for the determination of seismic moment by USGS. These values are used only as a reference for the constrain of the reproduced curve. The real values can only be obtained by a specific site study. In the case of the effect of source rupture directivity on the corner frequency and acceleration flat level,  showed that the directivity effect above the source corner frequency diminishes with increasing frequency for small earthquakes. The dataset of this study consists of moderate-sized events with a working frequency range from 0.5 Hz. However, the catalog of the used events does not contain information on the rupture direction of the seismic source except for the 2012 Samara earthquake. Thus, any correction that requires this information could not be conducted.
To constrain SI more strongly, we set the six reference events shown by the solid curves in Fig. 4a and regulated their acceleration at flat levels of source spectra as follows. First, using the tentative attenuation model, we converted the observed spectra at JTS000 to those at the unit distance. Then, we estimate the correction factors of the five events based on the converted spectra of event #10 (2017 Jacó, 6.5 Mw) except for the frequency range under the influence of κ-decay, and their acceleration flat level. We used these correction factors in further analysis, confining an ambiguity only in the flat level of acceleration spectra of event #10 that hereafter is denoted A 0 (#10).
Despite using the six reference events during the analysis, a considerable instability of outputs was observed: unacceptably large site effects at some sites, source spectra quite different from the w −2 model, and so on. This suggested that the input data may have been inaccurate or corrupted. Then we established the following way to check the correctness of the input and the SI output quality.

Consistency check using the individual site effects
We introduced two consistency checks and applied them to all twenty sites. In SI, it is assumed that each site must have its unique site effect, although various earthquakes were observed there. For each of the six reference events, it is possible to derive individual site effect by calculating the ratio of the observed spectra to the input at the seismic bedrock (i.e., the product of the given source effect and the path effect of the tentative attenuation model). Hereafter this is termed the "individual site effect". In our SI, the individual site effect of the used sites at each site was checked, and the outlier curves were detected as shown in Fig. 5.  Figure 5 shows an example of the deviation of the individual site effect that was considered acceptable, roughly within the factor of 4 or 5 within the frequency range from 1 Hz to 5 Hz, judged visually. On the other hand, Fig. 5b shows an example of the outlier (black broken curve), that looks biased constantly over the frequency. As the correct values are unknown, we eliminated the records that correspond to the outlier curves from the input data of SI. Eliminating a few records from SI sometimes resulted in a different type of instability of the output. Then, we had to repeat this sequence of checking, elimination, and SI until the output results were stabilized.
However, this consistency check cannot detect the following problem shared by all reference events. For example, in Fig. 5b, after the elimination of the black broken curve, the remaining three curves show a site effect as big as 60 at the peak around 3 Hz. This large value could not be accepted without a careful discussion, so we conducted the second consistency check described below.

Consistency check using synthetic spectra
Using Eq. 4, synthetic spectra of the I-th event at the J-th site (target record) can be calculated from source and site effects resulting from SI and two-element records, namely, the i-th event recorded at the J-th site and the I-th event recorded at the j-th site. It is possible to compare the synthetic spectra with the observed one if an existing record is used as the target. It is also possible to synthesize the ground motion at a site  We compare the synthetic spectra using various combinations of two-element records and look for outlier curves. For example, event #6 (2016 Capellades, 5.3 Mw) recorded at JTS000 was synthetically reproduced by two sets of element records: one for #06 at PPA800 and #10 (2017 Jacó, 6.5 Mw) at JTS000, the other for #06 at LPPG00 and #10 at JTS000 (Fig. 6). The synthetic one by the former pair (gray curve) shows an acceptable similarity with the observed one (black solid curve), however, the latter pair shows a drastic discrepancy (broken curve in black). This seems to be caused by the record of #6 at LPPG00 because the other element record was shared. Thus, this record was eliminated from SI. This second consistency check was conducted for all possible target records where observed records exist.

Consistency check working flow and final inversion with the 2012 Sámara earthquake
After several iterations of the two consistency checks, the outputs of SI using the six reference events were stabilized. Then, we conducted the final trial, but added the records of event #1, the mainshock of the 2012 Sámara earthquake, to the input and X TOES00 X X J T S 0 0 0 X X X X X X Fig. 8 A) Site effects as the result of the spectral inversion, and B) normalized ones by the PPA800 curve ran the two consistency checks again (Fig. 7). For the final interaction, we only used the reviewed data set, excluding the records with suspicious data (Table 3).
In parallel, we revised A 0 (#10) (2017 Jacó, 6.5 Mw) using site information for station PPA800 because this showed the lowest site effect curve and was located on rock exposed for the construction of a dam in a mountainous region. An in situ survey provided a Vs of 1,200 m/s on exposed lutite and sandstone (Piedra and Climent 2017). We estimated the amplification of vertically propagating SH waves applying the wave propagation theory based on the ratio of the impedance for SH waves, i.e., the square root of the product of density and Vs. The total amplification is given by the consecutive product of this impedance ratio of the boundaries that SH waves pass through. Thus, the impedance of the intermediate boundaries is canceled out, and the ratio of the crust to the top layer remains. This means the reverberation effect is neglected as usual for rock sites. We use a V s = 3, 600m/s and a density of 2,700 kg/m3 for the crust following USGS, whereas for the site with V s = 1, 200m/s and an assumed density of 2,000 kg/m3. According to Wohlenberg (1982) this is almost sandstone's lower limit. If there were actual density measurements, the site effect could be estimated more accurately. The role of density of the surface layer of rock sites has not been sufficiently examined so far, but it should. Then, we multiplied it by the surface effect 2.
For a better understanding in the earthquake engineering community, we regulated the value of A 0 (#10) so that the site effect of PPA800 in the target frequency range becomes this estimate through SI (Fig. 6). Figure 8a shows the site effects of the results. These thirteen sites successfully passed the two consistency checks. Site amplification at PPA800, which has the lowest site effect curve, was almost flat between 2 Hz and 5 Hz. The dotted black horizontal line shows the estimate mentioned above. In the frequency range higher than 5 Hz, a decay was observed commonly over sites. This seemed an effect of κ-decay. We discuss it in Section 5.3. Figure 8b shows the site effect normalized by the curve of PPA800.
The attenuation model of these results is Q = 179f 0.5598 . We checked the performance of SI by the comparison of the observed and the synthetic spectra. For JTS000 for event #10 (2017 Jacó, 6.5 Mw), for example, in Fig. 9 top-left, the thin solid curve shows the acceleration spectra observed, the thick solid curve is the synthetic one at the bedrock using the results of SI, and the thin broken curve is the reproduced ones using Eq. 1, as shown in Fig. 8a. Figure 9 top-right, shows another example (i.e., the case of SGBA00 for event #01 2012 Sámara earthquake 7.6 Mw). The acceptable similarity between the observed spectra and the reproduced ones confirms the success of our SI. Figure 9 bottom-left shows the source effect of event #1 determined by SI about A 0 =1.0E+19 (nm/s 2 ) and those of the other six reference events, for which acceleration flat levels were listed in Table 3. Figure 9 bottom-right shows the obtained values of the quality factor and its model compared with that of Moya (2009).

Reproduction of Fourier amplitude spectra of the 2012 Sámara earthquake
Using the results of SI, we applied Eq. 4 for synthetic spectra to reproduce the 2012 Sámara earthquake (Mw 7.6) records in the frequency domain at the five Sixaola sites (Table 2) where observed records did not exist during the mainshock. Namely, we used ICE stations as the core to produce a set of source, path, and site effects, and acceleration spectra at Sixaola sites as element records to prevent unnecessary complications due to their 14 bits A/D conversion. We could then post-dict the acceleration spectra using Eq. 4, as shown below. At station JTS000 located in Las Juntas (Fig. 10), the synthetic spectra of event #1 (2012 Sámara, 7.6 Mw) (Fig. 11 top-left) are almost ten times as large as those observed during event #10 (2017 Jacó, 6.5 Mw) (Fig. 4a) with a standard deviation of the order of 2. The observed spectra at SGBA00 (black broken curve) are shown as the reference of MMI VI. The similarity of the amplitude level of the two spectral curves shows a good consistency of the synthesis with observed intensity distribution, although a discrepancy is shown in the frequency range of the site effect of JTS000 at around 4 Hz (Fig. 8a).
Then, we reproduced the records of event #1 at two Sixaola sites just above the seismic source (RSN 2012): NICO00 at the colonial city of Nicoya, and PAQE00 in Paquera to the southeast of the Nicoya peninsula ( Fig. 10 for their location, and Fig. 11 topright and middle-left for their reproduced acceleration spectra, respectively). Both sites were in the zone of MMI VII, but the reproduced spectral level differed. The synthetic spectra at NICO00 have an amplitude of about 1.0 m/s 2 s (100 gal s). These are slightly bigger than the observed spectra at SGBA00 (black broken curve) and lower than those of PAQE00, which shows a broad peak of the amplitude 4.0 m/s 2 s (400 gal s) around 2 Hz with a side peak at 5 Hz, that implies the existence of thick soft sediments. This difference can be consistent with the surface geology that Nicoya city is on basalts from the Nicoya Complex, whereas PAQE00 is characterized by alluvial deposits (Denyer et al. 2013a, b), although more detailed physical information of the sites is not available at the present.
For the anomalous intensity MMI VII area observed around Naranjo town during event #1, we reproduced the acceleration spectra at RAMO00 in San Ramón and POAS00 in San Pedro de Poás to compare them with this anomalous intensity distribution ( Fig. 10 for their location, and Fig. 11 middleright and bottom-left for their reproduced spectra, respectively). The reported intensity was MMI VII at POAS00 in San Pedro, whereas at surroundings stations were MMI VI or lower (RSN 2012). The synthetic spectra of RAMO00 showed a broad peak from 1 Hz to 2.5 Hz, and its amplitude 1.0 m/s 2 s (100 gal s). The amplitude level was similar to the observed spectra at SGBA00 and comparable with or a little smaller than that of NICO00. POAS00 showed a broad peak from 1.3 to 2 Hz and its amplitude 3.0 m/s 2 s (300 gal s). It was bigger than NICO00 and a little smaller than PAQE00, consistently with MMI VII reported. POAS00 is the closest site to the Poás volcano. According to RSN (2012), this zone is composed of weathered pyroclastic deposits. It is thought that the seismic waves were amplified by this soft deposit, which may be thicker than surrounding areas.
Finally, we reproduced the acceleration spectra at another Sixaola site MEXI00, located at the central part of the capital city San José in a basin composed of volcanic sediments typical of the central part of Costa Rica and surrounded by major active faults (Montero 2001 ; Fig. 10). The synthetic spectra in Fig. 11 (bottom-right) showed amplitude similar or slightly smaller than the observed spectra at SGBA00. These sites presented a peak of 1.3 Hz, followed by a strong decay at the higher frequencies. This suggests the same or lower intensity in MEXI00, consistently to the location of San José at the border of MMI VI and V. Our SI results could explain the anomalous intensity distribution.

On spectral inversion
In the final SI, we used seven events, six for reference, and one major earthquake with an anomalous intensity (2012 Sámara earthquake, 7,6 Mw). Their spatial distribution covers Costa Rica almost evenly, except for the Caribbean coast: two at the border with Panamá, two at the central part, two at the Nicoya peninsula, and one at the border with Nicaragua (Fig. 3). Except for event #1, they were recorded at JTS000. However, for the final inversion, four events were excluded from the analysis due to their spectral shape (Fig. 4a  Fig. 11 Examples of reproduction of Fourier Acceleration Spectra of event #1 the 2012 Sámara earthquake at sites where any strong motion record was not obtained at that time: JTS000 (top-left), Sixaola stations NICO00 (top-right), PAQE00 (middle-left), RAMO00 (middle-right), POAS00 (bottom-left) in the intensity anomaly, and MEXI00 (bottom-right) in the national capital San José. Solid black curves and gray shadows show the average and the standard deviation, respectively. Observed spectra of SGBA00 (broken curve) are shown together as a reference for MMI VI broken curves) because of one of two possible reasons. The first case has a splitting of the corner frequencies by the appearance of an intermediate segment in the spectral curve, which implies a complex source process for earthquakes even as small as the magnitude 5 class. The other case is the group of events with an exceptionally high corner frequency, which means a higher stress drop than usual. Considering the complex tectonic setting (e.g., around a triple junction and the Panamá Fracture Zone), the stress drop, among other factors around this area. These high stresses drop events are also crucial for earthquake engineering because of their capacity to emit a higher acceleration than usual events in the same magnitude class. Further observations are necessary to reveal the background source process for them.
The screening using the two consistency checks introduced in this study showed a good performance. We detected several records that cause "individual site effect" or synthetic spectra significantly different from others (Figs. 5 and 6,respectively). For the former, we took the majority, whereas, for the latter, we could rely on the observed ones, and this choice stabilizes SI well. Being in the majority is, however, not always a criterion for correctness. It may be better to know the reason for their run-off from the majority. Theoretically, an erroneous scaling factor of the data file can cause it, though this is unlikely. The tentative attenuation model is a candidate for the cause, considering the discrepancy over a wide frequency range. The sparse and uneven distribution of records over the used events and sites used, can cause a trade-off between source and site effects. However, in this study, we could not identify the cause and had to leave it as a future issue.
Thirteen ICE sites passed successfully the two-step screening. The excluded six ICE sites were not concentrated in any region: one in the north, three in the central part, and two in the south (Fig. 2) With the almost even distribution of the used events, we could use a suitable spatial configuration for SI, except in the southern part where only one site LPPG00 was involved in the final SI. The records of the broadband site JTS000 played an important role, their existence was a fortune, and further studies for Costa Rica will bridge the seismological observations and strong ground motion observation needed for earthquake engineering.
The resulting attenuation model Q = 179f 0.5598 was considered representative for the northern and central part of Costa Rica. This is consistent with the value of Moya (2009) in the range from 1 Hz and 2 Hz, in which the author estimated a value of Q = 131.6f 1.1 for the central part of the country (Fig. 9b). As its exponent of the frequency is greater than the unity, the attenuation factor exp(−ρf R/Qβ) for the Q model of Moya (2009) is a monotonically increasing function for the frequency within the target frequency range, whereas that of this study shows decreasing feature and seems more likely. Therefore, in this study, we could extend the applicable region for the attenuation model. The differences between the two studies can be attributed to the difference in the area coverage, since Moya (2009) focused only in the Central Valley, which is mainly a sedimentary basin, while we cover a bigger area with more geological variations. Attenuation in the crust may be related to heterogeneities caused by the volcanism and the crustal deformation. In this context, the southern part can have an attenuation different from the rest of the country, as Moya (2009) suggested. The extension of coverage to the southern part of Costa Rica and the northern part of Panamá, and regional difference of attenuation models may be interesting research topics. However, they are beyond the scope of this study and are thus left for the future.
We obtained the site effect of JTS000 and thirteen ICE sites. As expected, the site effect of PPA800, classified as S1 of the site classification of the seismic code of Costa Rica (CFIA 2010), presented the lowest curve, whereas the highest site effect with values up to 10 times as high (Fig. 8 left and right). On the other hand, JTS000 shows a non-negligible site effect, that is not biased over the target frequency range, but a peak at around 4.0 Hz. According to the URL of Project IDA (https://ida.ucsd.edu/?q=station/ station-jts-las-juntas-de-abangares-costa-rica/ fdsn-station-page-jts) for this site, its geological materials are coral and limestone. As the very low ambient noise level of this site shows the absence of perturbation by human activities, we imagine that this site effect may be due to sedimentary rock layers such like limestone, perhaps thicker than several hundred meters. We wait for further research on it in the future.

On reproduction of acceleration spectra
By the nature of SI, the outputs shown in Fig. 11 are obtained assuming the elastic soil behavior. In general, soft soil behaves non-linearly under strong ground motion. An additional procedure based on the soil dynamics may be necessary to convert these outputs to the ground motion at the engineering bedrock using calculation of elastic response and again to convert this bedrock motion to obtain the strong motion at the ground surface computing nonlinear response. The information on underground velocity structure and the dynamic property of each site's soil are indispensable. The intensity at these five Sixaola sites, however, equal to or below MMI VII. A strong nonlinear behavior may not take place, but just a weak one. Therefore, we believe that the synthetic spectra shown above may provide the upper bound and be used for the safe side consideration in the engineering sense.
We did not involve the Sixaola sites directly in SI. Their usability, however, has not yet been tested and may be worth trying after sufficient accumulation of significant strong motion records. From the viewpoint of seismic disaster prevention, the existing Sixaola sites show their importance of their location at fire stations in busy cities. This study showed one of the ways how we could integrate the low-cost accelerographs for disaster prevention (IT accelerographs, that is Sixaola) with those for scientific or technological purposes as ICE sites.
It may be worth mentioning that the Fourier acceleration amplitude spectra can be converted to the input energy spectra with a mathematical exactness (Ordaz et al. 2003). This shows one way in which the results of this study can provide useful information for earthquake engineering.

On K-decay
As shown in Figs. 4a and 9 (top), κ-decay exists in observed acceleration spectra and the determined site effect curves. In this study, we conducted SI using reference events, assuming that their source spectra follow the w −2 model, namely, that the source spectra do not have high frequency decay. Therefore, κ-factor of site effect curves may include the site specific κ 0 and the source contribution κ s . Oth et al. (2011) demonstrated that κ-decay should be related to attenuation within the near-surface weathered layers using the relation log 10 S(f > f E ) = S 0 − log 10 exp(−πκf ), where S and S 0 denote acceleration spectra, and a flat level f E the lower frequency bound, considered 10 Hz or 15 Hz without significant differences in their results. They compared the source spectra from the results of SI using reference site and strong motion records of surface sensors of KIK-net (NIED 2020) with those of borehole sensors. The decay represented by κ (κ−decay) appeared stronger in the SI results. The average κ-factor were 0.029s for surface sensors and 0.015s for borehole sensors. They proposed a correction for the source spectra in terms of the ω −2 model. Studies of κ in the context of the site effect have concentrated on the distance-independent term κ 0 , where κ = κ 0 + κ s +κ R (R) and κ R (R) denotes the distance-dependent component (Anderson 1991;Ktenidou et al. 2014;Ktenidou et al. 2015;Parolai 2018). κ 0 means κ value at zero distance, and represents the attenuation just beneath the site that the shear waves suffer at the final stage of the propagation (Ktenidou et al. 2015). The rest of the effect represented by κ s may be due to the seismic source, and can be resolved for dataset of adequate sources and stations (Ktenidou et al. 2014). If the separation of the three effects by SI had been perfect, κ R (R) would have been allocated to the path effect. Figure 12 (top-left) shows the decay part of the site effect curves (Fig. 9a) but normalized by their own values at 5.0 Hz. The broadband site JTS000 shows smooth and the sharpest decay. Among a few curves of weak decay, the rock site PPA800 decays smoothly. In determining the κ-factors for these two sites, we extracted their curves in the frequency range from 5 Hz to 7 Hz and conducted regression analysis using an exponential function. The results are shown in Fig. 12 (top-right) for PPA800 with κ=0.043 s. This implies that its cause may lie in sedimentary surface layers below the exposed rock layer or at a seismic source. Figure 12 (middle-left) shows those for JTS000 with κ=0.168 s. The peak of the site effect at around 4 Hz shown in Fig. 8 (left) is suspected as the cause of this large value, and a significant occupancy of the site specific κ 0 in it.
We checked their dependence on distance and magnitude using the observed spectra at JTS000 (Fig. 4a), using the same procedure as above for determination of κ-factor of site effects but with the correction using the path effect determined in this study. Only one exception is event #16, of which frequency range for determination started from 5.5 Hz due to a small bump just above 5 Hz. A negative and small trend of distance dependence (Fig. 12 bottom) may show that the assumption we made may be appropriate, namely κ-factor's dependence on distance was almost allocated into the path effect. However, due to a small trend shown and significant deviation for its magnitude dependence, it may be difficult to conclude with certainty about the source contribution to κfactor ( Fig. 12 bottom-right). More precise estimate of κ-factor using more events and stations distributed appropriately may enable to reveal the source contribution κ s in the future.

Final considerations and conclusion
We conducted SI to separate source, path, and site effects using mainly the strong motion records obtained at the sites of the ICE-UCR network covering all of Costa Rica at essential electricity supply facilities, such as power distribution centers and hydroelectric dams. We also used the broadband records at the site JTS000 (JTS of IRIS/IDA).
For SI to suppress the instability of outputs, we introduced two-step screening. One step is the Fig. 12 Analysis of the decay at high frequency. Top-left: Decaying part of the site effect curves normalized by their values at 5 Hz. Top-right: Determination for the site effect of PPA800 κ = 0.043; middle-left: κ=0.168 for the site effect of JTS000. Middle-right: Decaying part of the observed spectra at JTS000 that are corrected by the path effect and normalized by their values at 5 Hz. Bottom-left: κ-factor's dependence on distance for the observed spectra of JTS000 corrected by the path effect shown with that of the site effects at PPA800 (thin broken line) and JTS000 (thin solid line); bottom-right: that on magnitude. Double arrows in the panels top-left and middle-right show the frequency range used to determine the value of κ. The broken curves in the panels top-right and middle-left show the calculated ones using determined values of κ consistency check of the inputs using the "individual site effect", which is derived from the observed spectra tentative source and path effects based on the result of preliminary inversion. The other step is the consistency check of the outputs, comparing observed spectra with synthetic ones based on the site and path effects of outputs from SI. This screening performed well in being used iterative to check, eliminate records, and do an inversion. We also applied the regulation of source spectra of reference events considering the estimate of site effect for one-dimensional SH wave propagation at the site of the lowest site effect (PPA800).
Then, in the frequency domain, we applied the results to several Sixaola sites to reproduce the strong ground motion during the 2012 Sámara earthquake at specific sites that were not monitored during the event. We obtained the synthetic spectra consistent with the reported intensity distribution in and around the anomaly surrounding the POAS00 site near Naranjo township.
Besides application to a specific data set of Costa Rica, the two methods we introduced in this study: consistency checks, and spectral synthesis, have a generality that can be applied to other data sets in other regions to separate source, path, and site effect from strong motion records, and to post-dict the acceleration spectra for past major earthquakes at sites that are not monitored then. In this study, we find them a useful tool to find and discard records that have been corrupted and it was not detectable with any other method.