###### Abstract
It is important to consider such a soil structure if the site under study is located in the water area within the zones of continuous or sparse permafrost, as well as the massive release of bubble gas from bottom sediments. Such zones are particularly extensive in the Arctic seas of Russia, especially on the Arctic shelf of Eastern Siberia, where a large number of discovered methane seeps are located [14; 15; 16]. The water content of the soils does not affect horizontal ground motion significantly since water has no resistance to shear stress [17]. Conversely, porous soil structure and position of the water table considerably affect the vertical response [12; 18; 19]. Some investigations based on seafloor seismic records analysis show that the vertical component of seafloor motions are significantly lower than those of onshore motions near the P-wave resonant frequencies, caused by the water layer above the station [20; 21; 22]. Some works justify the applicability of the well-known codes, which are usually used for SH-wave, to the P-wave propagation case [12; 19]. When solving the wave equation, the velocity of the P-wave is used instead of the S-waves. Additionally, constrained moduli are used instead of shear moduli. However, for the P-wave case, there are no indications of a specific well-known code in the literature [1]. The purpose of the present study is to introduce an algorithm and its Matlab implementation for numerical modeling of the nonlinear response of porous saturated soil deposits to P- and SH-waves propagation, based on the existing literature on this topic. We consider one-dimensional vertical propagation of P- and SH-waves in porous saturated soil layered system. The basics of the NERA algorithm [5] were used to model the nonlinear behavior of soils. The same applies to spatial and temporal discretization for the numerical solution of analytical wave equations. The [PERSON] [23; 24] and [PERSON] [25] theories for seismological low-frequency range were used to substantiate the applicability of the derived equations of motion for the case of porous saturated soil. The procedure for passing from the solution for SH-waves to the solution for P-waves is based on solving the wave equations using the P-wave velocity instead of S-wave velocity and degradation curves of constrained moduli for saturated soils instead of shear moduli according to [12]. Simulation of the reduction of P-waves spectra due to the passage of part of their energy into the water column (water layer above seabed site) and resonance in it, as formulated in [21], is carried out through the use of the series of band-stop filters around the resonant frequencies. The algorithm presented in this paper was implemented in the Matlab programming environment. The 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 [26] since the vertical arrays in water areas are still exotic. The records for two sites were used: the TCGH15 site is characterized by sandy soils lying on the bedrock, and the KSRH10 site is characterized by clayey soils. 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 column on the reduction of vertical motion 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.