## 4 Approach of the Algorithm with the Kik-Net Vertical Arrays Data  
To verify the modeling results, the records of vertical seismic arrays of the Kik-net network [26] were used. Figure 7 shows the \(V_{P}\) and \(V_{S}\) velocity profiles at TCGH15 and KSRH10 sites. The TCGH15 site is characterized by sandy soils lying on the bedrock at a depth of 22 m. The KSRH10 site is characterized by clayey soils lying on the bedrock at a depth of 36 m [37]. Figure 6: Block-scheme of MatNERApor Matlab package for numerical modeling of nonlinear response of porous saturated soil deposits to seismic body waves propagation. Numbers in blue circles—subsets of the package (see description in text). For most cases, only velocity profiles are available in the database [26]. Other geotechnical parameters were calculated according to various equations (Tables 1 and 2). The bulk density can be calculated as a function of the P-wave velocity [38] or as a function of the S-wave velocity [39]. Figure 7 shows that at a certain depth there can be a sharp increase in \(V_{P}\) with a relatively small increase in \(V_{S}\). Such a horizon can be interpreted as the level of groundwater. The bedrock is characterized by a significant increase in both \(V_{P}\) and \(V_{S}\). In accordance with the accepted model of the elementary volume of soils with different moisture content (Figure 2), the dependences of density on \(V_{P}\) are better suited for finding the bulk density of non-cohesive soils, since an increase in moisture of such soils should lead to an increase in bulk density. Conversely, to find the bulk density of cohesive soils, the dependences of density on \(V_{S}\) are better suited since an increase in moisture of such soils should lead to a smaller change in density or to its decrease. Two earthquakes were chosen to model the site response as they resulted in approximately the same PGA at both sites (\(\sim\)60 Gal):  
(1) \(\mathrm{Mw}=5.7\), origin time 6 October 2004 14:40:38 UTC (for TCGH15 site) [40];  
(2) \(\mathrm{Mw}=5.2\), origin time 10 May 2008 18:24:02 UTC (for KSRH10 site) [40]. Figure 8 shows degradation curves of shear modulus (from literature), constrained moduli above the water table and under the water table (calculated) and corresponding damping ratio curves used for modeling the responses for TCGH15 and KSRH10 Kik-net sites. We used the shear modulus degradation curves from [41] for sandy sites, from [42] for clayey sites and from [2] for stiff bedrock. Figures 9 and 10 show a comparison of the response spectra modeled by MatNERApor and recorded accelerograms on the surface of soil layers at TCGH15 and KSRH10 sites, respectively. Calculations were carried out for all three components of records. It can be seen that the simulation results basically correspond to real records. Figure 7: Velocity profiles at sites TCGH15 (**a**) and KSRH10 (**b**) of the Kik-net network. Modeling accuracy can be assessed by using the residuals as follows [12]:  
\[ \begin{array}{l}R_{SA}(T)=\ln[SA(T)]_{record}-ln[SA(T)]_{model}\\ R_{SV}(T)=\ln[SV(T)]_{record}-ln[SV(T)]_{model}\\ R_{SD}(T)=\ln[SD(T)]_{record}-ln[SD(T)]_{model}\end{array}  \tag{40}\]  
where \([SA(T)]_{record}\), \([SV(T)]_{record}\), \([SD(T)]_{record}\) are spectral acceleration, velocity and displacement of the record for period T, \([SA(T)]_{model}\), \([SV(T)]_{model}\), \([SD(T)]_{model}\) are spectral acceleration, velocity and displacement of the modeled accelerograms for period T.