## 1 Introduction  
Site-specific response analysis is a necessary and important component of seismic hazard assessment made for the purposes of residential and industrial facility construction in any tectonically active region [1, 2]. An accuracy of geohazard assessment for such installations as nuclear power plants, mineral production and manufacturing facilities, bridges, etc., is crucial due to high environmental risks. Especially, it concerns the structures that have direct contact with open water areas, such as offshore oil platforms, underwater pipelines, marine terminals, as the sea environment has a high vulnerability and ability to quickly transport pollution. The need for site response analysis is declared by national and international building code regulations, but they usually do not indicate a specific algorithm that must be used for such assessment. Experimental recordings are more preferred; however, the surveyors very often face the lack of records collected at the observation site, associated with the earthquakes with appropriate parameters. Long-term instrumental seismic observations at the seafloor are technically more difficult and expensive than on land; this exacerbates the problem of insufficient experimental data of strong motion records. This challenge is usually overcome by numerical simulation of seismic wave propagation through soil deposits with specified parameters. Many algorithms for soil response modeling have been developed over several last decades. Most of the site response models can be divided into two general groups: equivalent linear and nonlinear models. Equivalent linear models (EQL) approximate stress-strain hysteresis during loading cycles with an secant shear modulus, while the nonlinear models (NL) simulate stress-strain curves using different constitutive soil models. Such models were developed for the case of SH-wave propagation. EQL methods are frequency-domain and firstly implemented in the SHAKE code [3]. The following popular computer programs, such as EERA [4] and STRATA [5], utilized the same concept with some improvement in computational efficiency, graphical user interface, the possibility to account for multiple input motions and uncertainties of soil parameters. NL methods are usually time-domain and differ mainly by constitutive soil models. For example, NERA code [6] is a time-domain algorithm, which uses IM-model introduced independently by [PERSON] and [PERSON] [7,8]. The IM model is based on simulating nonlinear hysteretic stress-strain curves using a series of mechanical elements that have different stiffness and sliding resistance. One of the most popular codes, DEEPSOIL [9], performs both frequency- and time-domain analysis with hyperbolic constitutive models. It uses the modified Kondner-Zelasko model (MKZ) [10] and general quadratic/hyperbolic (GQ/H) [9]. Hysteretic stress-strain behavior is taken into account by defining model parameters using an automatic curve fitting procedure. Verification studies have shown that equivalent linear models are suitable when shear strains in the soil remain low. However, they are not applicable for the pore pressure modeling or permanent displacement calculation and also can lead to spurious resonances, while nonlinear models are more accurate for simulating response for soft soil sites that experience higher strains and stronger ground motions. At the same time, such models usually require much more input parameters and are computationally less efficient [1,11]. The model types mentioned above can be used for one-, two-, and three-dimensional simulation [12]. Two-dimensional analyses are commonly used for geotechnical structures having lateral dimensions that are large compared to the wavelength of significant earthquake ground motion; the out-of-phase effects result in ground motions spatial variations [12,13]. Three-dimensional simulation is even more complicated and used in cases when the motions or boundary conditions vary significantly in three dimensions.