###### Abstract  
The paper is devoted to the problem of numerical modeling of earthquake response of porous saturated soil deposits to seismic waves propagation. Site-specific earthquake response analysis is a necessary and important component of seismic hazard assessment. Accounting for the complex structure of porous saturated soils, i.e., the content in them, in addition to the solid matrix, pore water, gas mixture and ice, is especially important for the water areas in the zones of continuous or sparse permafrost, as well as the massive release of bubble gas from bottom sediments. The purpose of this study is to introduce an algorithm and its Matlab implementation for numerical modeling of the nonlinear response of porous saturated soil deposits to vertical P- and SH-waves propagation. The presented MatNERApor package consists of a set of Matlab scripts and functions. The package was tested and verified using the records of vertical seismic arrays of the Kik-net network. In addition, the records of local earthquakes obtained by ocean bottom seismographs in the Laptev Sea in 2019-2020 were used to demonstrate the effect of the water layer above the seabed sites on the reduction of vertical motions spectra. The results of the calculations showed good agreement with the data obtained from real seismic records, which justifies the correctness of the theoretical basis of the presented algorithm and its software implementation. title = [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [  
Seismic hazard assessment is a necessary component of research work in the design and operation of structures in seismically hazardous areas. Site-specific earthquake response analysis is an important component of seismic hazard assessment. According to national and international building code regulations, one of the main results of seismic hazard assessment are sets of accelerograms (time series of ground accelerations), consisting of two horizontal and vertical components. The influence of soil layers on the propagation of seismic waves (seismic microzonation) should be taken into account. Experimental seismic 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 since long-term instrumental seismic observations in sparsely populated areas or in water areas are technically difficult and expensive. This problem can be solved by numerical simulation of the seismic wave propagation through soil deposits with specified parameters. Many algorithms for soil response modeling have been developed over several last decades, that can be divided into two general groups: equivalent linear [2; 3; 4] and nonlinear models [5; 6; 7]. Verification studies have shown that equivalent linear models are suitable when shear strains in the soil remain low while nonlinear models are more accurate for simulating response for soft soil sites that experience higher strains and stronger ground motions [8; 9]. It is commonly believed that the main contribution to the destruction of structures during earthquakes is made by shear waves since usually their amplitude is maximum during strong ground motions and they have a characteristic spectral composition. However, records from recent earthquakes show that many structures have experienced significant damage because of the high intensity of vertical earthquake motions [10; 11; 12]. Vertical ground motions also significantly affect slope stability [12; 13]. The most common algorithms assume homogeneous solid soil structure while the real soils are porous and saturated by multicomponent fluids. In addition to water, the pore space of soils can contain a significant amount of gas mixtures and ice.