Localized stationary seismic waves predicted using a nonlinear gradient elasticity model

This paper aims at investigating the existence of localized stationary waves in the shallow subsurface whose constitutive behaviour is governed by the hyperbolic model, implying non-polynomial nonlinearity and strain-dependent shear modulus. To this end, we derive a novel equation of motion for a nonlinear gradient elasticity model, where the higher-order gradient terms capture the effect of small-scale soil heterogeneity/micro-structure. We also present a novel finite-difference scheme to solve the nonlinear equation of motion in space and time. Simulations of the propagation of arbitrary initial pulses clearly reveal the influence of the nonlinearity: strain-dependent speed in general and, as a result, sharpening of the pulses. Stationary solutions of the equation of motion are obtained by introducing the moving reference frame together with the stationarity assumption. Periodic (with and without a descending trend), as well as localized stationary waves, are found by analyzing the obtained ordinary differential equation in the phase portrait, and integrating it along the different trajectories. The localized stationary wave is in fact a kink wave and is obtained by integration along a homoclinic orbit. In general, the closer the trajectory lies to a homoclinic orbit, the sharper the edges of the corresponding periodic stationary wave and the larger its period. Finally, we find that the kink wave is in fact not a true soliton as the original shapes of two colliding kink waves are not recovered after interaction. However, it may have high amplitude and reach the surface depending on the damping mechanisms (which have not been considered). Therefore, seismic site response analyses should not a priori exclude the presence of such localized stationary waves.


Introduction
For the prediction of the so-called seismic site response -the response of the top soil layers of the earth -induced by seismic waves, typically 1-D models are employed [1]. The so-called equivalent linear scheme is used very often, in which the soil stiffness and damping are modelled assuming constant shear modulus and material damping ratio, respectively [2]. The actual values of the shear modulus and damping ratio of a specific layer in the soil profile are based on the maximum level of strain observed inside that layer; this requires iteration as it is not possible to determine the maximum strain level a priori. For high maximum strain levels in the soil layers, such equivalent shear modulus and damping ratio cannot accurately represent the behavior over the entire duration of a seismic event, as the strains vary significantly. In such cases, a nonlinear time domain solution is typically used to account for the variation of the shear modulus and damping ratio during shaking (e.g., [1]).
Contemporary research into nonlinear time-domain models for seismic site response analysis is mostly focused on the development of advanced constitutive models so as to capture important features of the soil behavior such as anisotropy, pore water pressure generation and dilation [3]. However, limited research has been devoted to the fundamental nonlinear-dynamics aspects of the seismic site response. Employing the commonly used hyperbolic constitutive model [4,5], softening behaviour and super-harmonic resonances were recently demonstrated for a superficial soil layer under uniform harmonic excitation at the lower boundary [6]. However, the possibility for localized stationary waves such as solitons to propagate through the soil column and reach the surface has not been widely recognized in the seismological literature and neither in the geo-technical literature, although a number of publications hint at the possibility of Love-type surface solitary waves, solitary waves along faults, and the importance of nonlinearity in general [7,8,9,10,11,12,13,14].
In this paper, we therefore investigate the existence of localized stationary waves in the shallow subsurface with the constitutive behaviour governed by the hyperbolic model, implying that the shear modulus is strain dependent (i.e., non-polynomial nonlinearity). As the classical wave equation with this particular nonlinearity has non-physical discontinuous solutions, we employ a nonlinear gradient elasticity model, which sometimes is also called a higher-order gradient continuum or a micro-structured solid; the equation of motion is of the Boussinesq-type. Compared to the classical continuum, the stress-strain relation has higher-order gradient terms to capture the effect of small-scale soil heterogeneity/ micro-structure (yet keeping the description of the material homogeneous), which introduces dispersive effects particularly for the shorter waves [15,16,17,18,19,20,21,22,23]. Such higher-order gradient terms are naturally obtained using asymptotic homogenization techniques for periodically inhomogeneous media [15,24,25,26,27]. Dispersive effects significantly influence the behaviour of localized stationary waves, as such waves exist exactly because of the balance between dispersive and nonlinear effects, allowing their propagation without distortion. The dispersion prohibits the formation of jumps, which leads to physically realisable solutions. So-called higher-order dispersion correction and higher-order dispersive nonlinearity has been discussed in the context of the well-known Korteweg-de Vries equation [28]. For higher-order gradient elasticity continua, the existence and properties of localized stationary waves have been studied for polynomial nonlinearity in both the macro-scale and the microscale terms (i.e., in the higher-order derivative terms) [13,29,30,31,32,33,34,35]. However, such waves have never been studied in the context of higher-order gradient elasticity continua that are dictated by hyperbolic nonlinearity, which is done in this paper.
Sticking to the 1-D assumption, we derive the equation of motion for the nonlinear gradient elasticity model based on Newton's second law and Eringen's general stress-strain relation [36], which allows introducing the higher-order derivatives in an unambiguous manner (Section 2). An ordinary differential equation from which stationary solutions for the nonlinear gradient elasticity model can be obtained, is derived in Section 3. Periodic (with and without a descending trend) as well as localized stationary waves are found (Section 5.2); the latter is in fact a kink wave and is obtained by integrating along a homoclinic orbit in the phase portrait [18], like it can be done for the well-known sine-Gordon equation. Section 4 presents a novel numerical scheme to solve the nonlinear equation of motion in space and time. It exploits the structure of the partial differential equation in order to simplify the computation of the spatial finite-difference approximations. The pseudospectral scheme presented in [29] cannot be directly applied for the hyperbolic nonlinearity, since the roots of a polynomial with non-integer exponents would have to be determined; this cannot be done analytically and is therefore numerically expensive. For that reason, the introduced numerical scheme directly considers the hyperbolic nonlinearity. Another advantage of the numerical scheme is that, because it employs the finite-difference method, we are not limited to periodic boundary conditions. The scheme is used to study the propagation of arbitrary initial pulses (Section 5.1) as well as to check the stationarity of the stationary waves identified (Section 5.2). In addition, it allows to demonstrate that the kink wave is in fact not a true soliton. This is done by means of a numerical collision experiment in which two kink waves propagate in opposite direction and pass each other (Section 5.3); after interaction, their original shapes are not recovered.
Even though not being a true soliton, the kink wave, which may have high amplitude, can propagate through the soil column and potentially reach the surface depending on the strength of the material and geometrical damping mechanisms (which have not been considered). Therefore, seismic site response analyses should not a priori exclude the presence of such localized stationary waves.

General framework
The starting point to derive the equation of motion of the gradient elasticity model is Newton's second law. For transverse waves propagating in the vertical direction z and considering the one-dimensional situation, it reads [37] ρ ∂ 2 u ∂t 2 = ∂σ zx ∂z . (1) Here, ρ denotes the material density, u(z, t) the horizontal displacement in the x direction, t is time, and σ zx is the shear stress. In a nonlinear system, the stress-strain relation can generally be written as follows [36]: This relation expresses that the stress σ zx at location z and time t generally depends on the strain ε zx at all points of the. Moreover, on the entire strain history (note that ζ and τ denote auxiliary space and time variables, respectively); the specific nonlocality and history dependence is contained in the kernel function g(z, t, γ). The strain depends on the displacement in the following way [37]: The nonlinearity comes into play through the γ dependence of the kernel function, which expresses that the elasticity operators, i.e. the elastic coefficients, are strain dependent. In this work, we relate γ to the deviatoric component of the so-called octahedral strain, which for the one-dimensional case leads to [38,39]

Equation of motion for nonlinear gradient elasticity model
If the distribution of the kernel function g(z, t, γ) in z and t is arbitrary, the equation of motion, as obtained by substituting Eq. 2) into Eq. (1), will be of the integro-differential type. However, as our aim is to derive the equation of motion for a gradient elasticity model, which is a partial differential equation, we restrict ourselves to Dirac delta functions. Apart from the conventional term, we include terms with double space and double time derivatives: where δ(...) denotes Dirac's delta function, and (...) ,ζζ and (...) ,τ τ denote double partial differentiation with respect to ζ and τ , respectively (this notation is used throughout the paper, when useful, to denote partial derivatives); G(γ) is the strain-dependent shear modulus, G (L) (γ) and G (T ) (γ) are additional strain-dependent elastic moduli (related to higher-order derivative terms, as shown below), and γ = γ(ζ, τ ). The semi-local operators with derivatives of the Dirac function come with the time (T ) and length (L) scales that characterize the specific history dependence and nonlocality of the medium, respectively; it is well-known that particular nonlocal and history effects can indeed be captured by using higher-order derivative terms in the equation of motion (e.g., [40]). The kernel function with the two additional terms (with double space and double time derivatives) and corresponding signs is chosen in accordance with the linear model (see Eq. (14) below), for which it ensures unconditional stability as well as realistic lower and upper bounds for the speed of energy transfer of the propagating wave [15]. Now, inserting Eq. (5) into the Eq. (2), we obtain the stress-strain relation of the nonlinear gradient elasticity model: Substituting the Eqs. (3) and (6) into Eq. (1) yields the corresponding equation of motion: ∂u ∂z .
(7) For simplicity, we now relate the strain-dependent additional elastic moduli to the conventional strain-dependent shear modulus G(γ) using dimensionless constants B 1 and B 2 : Without loss of generality, we interrelate the characteristic length and time scales: Here, c 2 0 = G 0 /ρ, with G 0 being the well-known small-strain shear modulus from linear elasticity; c 0 is the corresponding shear-wave speed. Tthroughout the paper, the subscript "0" indicates that the quantity relates to smallstrain/linear behavior. Using Eqs. (8) and (9), Eq. (7) can be written as In this work, we use the hyperbolic soil model typically employed for seismic site response analyses, with the following expression for the straindependent shear modulus [4]: Here, γ ref denotes a reference shear strain, and β is a dimensionless constant (0 < β < 1).

Limit case
For completeness, we here consider the limit case of γ/γ ref → 1, which reduces the model to a linear gradient elasticity model. In that case (cf. Eq. Using the same definition for T (Eq. (9)), the stress-strain relation Eq. (6) reduces to whereby . The corresponding equation of motion reads as follows: This linear equation was derived before by Metrikine and Askes [41] from a discrete model using a continualization procedure. A similar equation was used by Georgiadis et al. [42] to investigate the existence of horizontally polarized surface waves.

Stationary wave solutions
It is possible to determine stationary solutions of the equation of motion Eq. (10) that do not change their shape while propagating in the nonlinear medium. This means that, starting at an initial condition, the stationary solution does not change with respect to the coordinate ξ = z − ct, which moves with the velocity c ∈ R. Thereby, c can be arbitrarily chosen and leads to a specific solution. In order to find possible stationary solutions, we apply the transformation ξ = z − ct and assume stationarity, which yields Here, we recall that these particular subscripts (time or space variable(s) preceded by a comma) denote partial derivative operators. Substituting expressions of Eq. (15) into Eq. (10), we get an ordinary differential equation for the determination of stationary solutions to Eq. (10): In order to evaluate ∂ 2 ∂ξ 2 (G(γ)u ,ξ ) in Eq. (16) for the hyperbolic soil model from Eq. (11), it has to be noted that the absolute value function g(x) = |x| is not differentiable for x = 0. However, it is weakly differentiable with sgn(x), which denotes the sign function, as weak derivative. Using this, we specifically obtain Integration of Eq. (17) with respect to ξ, grouping terms related to u ,ξξξ , and setting y := u ,ξ ,() := ∂ ∂ξ yields This is a nonlinear second-order ordinary differential equation, which can be solved for stationary wave solutions of Eq. (10). Specific stationary solutions are determined in Section 5.

Numerical scheme
In order to solve Eq. (10) in space and time, we use a newly developed fully implicit scheme for the numerical solution of partial differential equations of the form where G ∂u ∂z denotes that G can depend on ∂u/∂z in an arbitrary manner. It is assumed that the solution u(z, t) of Eq. (10) exists in time t ∈ [0, T] and space z ∈ [Z , Z h ]. Therefore, a grid in time and a grid in space are introduced. Using Eq. (19) can be written as By using Eq. (22) and (23), the structure of the considered partial differential equation is exploited. This will simplify the computation of the spatial finite difference approximations, as can be seen in the further course of this section. Assuming that the solution is known at the timepoints t n−1 and t n and replacing the time derivative by a finite-difference approximation, Eq. (23) results in the nonlinear equation with and where f i is the i-th component of f , u n i is a grid function approximating the solution at time t n and space z i , i.e. u n i ≈ u(z i , t n ), and u n i,z approximates u ,z (z i , t n ). Equations (24) and (26) approximate Eq. (23) (evaluated at t = t n ) up to an accuracy of O(∆t 2 ), whereby O(·) represents the big O notation. Next, the space derivatives are discretized. For this, standard finite-difference approximations are used again. In order to simplify the notation, the time index n is omitted in the following, such that u n is written as u or u n i is written as u i . This leads to the approximations, which have all an accuracy of O(∆z 2 ): (29) Now, the advantage of the exploitation of the structure of Eq. (19) and the introduction of the function h(u ,z ) can be seen. In Eq. (19) the thirdorder derivative of G(u ,z ) with respect to z appears. However if the hyperbolic soil model G(γ) from Eq. (11) is used for G(u ,z ) in Eq. (19), then G(u ,z ) = G(γ) contains the absolute value function, which is only one time weakly differentiable. A direct finite difference approximation of the thirdorder derivative of G(u ,z ) must be able to deal with this problem. By using Eq. (28) and (29), the missing differentiability in the case G(u ,z ) = G(γ) is circumvented.
Since Eq. (24) is nonlinear, a numerical scheme has to be used in order to calculate u n+1 i iteratively. For this, Newton's method is used, for which the computation of the Jacobian matrix is necessary.
Using |x| ,x = sgn(x), sgn(x)x = |x| and defining k(x, y) := it follows that This concludes our new numerical scheme for the computation of solutions to Eq. (19). In order to simulate the solution on an open domain, absorbing boundary conditions are used. The numerical solution of the linear equation of motion (Eq. (14)) can be computed in a similar manner (using a finite difference scheme). For details, see appendix Appendix A.

Numerical results
In this section we show and discuss numerical results for the nonlinear Eq. (10) with the choice of parameter values indicated in Tab. 1. The values of G 0 , ρ, β and γ ref have been chosen to represent soil, and the values of the parameters related to the higher-order gradient terms are similar to the ones used in [15].

Solutions for a Gaussian pulse
First of all, the temporal evolution of a specific solution is studied, where as initial condition a Gaussian pulse is used, i.e.
Here, the amplitude u 0 = 0.016 m and standard deviation σ = 30 m have been chosen to obtain a relatively high strain level, in accordance with [1]. In order to compute the numerical solution using the scheme described in Section 4, an initial condition at time point t −1 (i.e. u −1 ) has to be chosen. In this study, u −1 = u 0 is used, resulting in a solution with zero initial velocity. The resulting numerical solution can be seen in Fig. 1. The initial pulse divides into two parts, which propagate in opposite directions. We can also observe that the nonlinear solution does not have sharp edges, which would be the case if the higher-order derivative terms were omitted in Eq. (10); hence the gradient elasticity model yields physically admissible behaviour. In order to see the effect of the nonlinear terms in Eq. (10), the temporal evolution of the corresponding solution for the linear case (Eq. (14)) is shown in Fig. 2. There is a clear difference between the linear and nonlinear solutions. It can be seen that, although the general behavior of the solution stays the same, the nonlinear solution is much sharper compared to the linear solution (but still smooth). In order to make this more clear, Fig. 3 shows both solutions at the end of the simulation time.
Another effect of the nonlinearity can be observed in Fig. 4, which shows the solution of the nonlinear equation of motion Eq. (10) with a Gaussiancosine pulse as initial condition, which is given by  Thereby, an amplitude of u 0 = 0.016 m and standard deviation of σ = 30 m have been chosen. It can be well observed that, due to the nonlinearity, different parts of the pulses having different slopes (i.e., strain levels), travel with different speed. In order to study the interaction of two solutions, two Gaussian pulses starting at z 1 = −500 m, z 2 = 500 m with amplitudes u 1 = 0.008 m, u 2 = 0.016 m, respectively, and the same standard deviation σ = 30 m are considered. The corresponding initial condition has been calculated by adding both pulses. The obtained temporal evolution of the corresponding solution is displayed in Fig. 5. It is shown that both pulses divide directly in two parts, propagating in different directions. Thereby, the interaction of the resulting waves does not destroy the structure of each wave but does lead to a very small change in their amplitude, which can be verified by comparing the waves that did interact with the ones that did not. It can thus be concluded that such waves stay localized after collision with each other. At the same time, these waves are not stationary and their dispersion can be clearly seen in Fig. 5, so the waves are in principle not even expected to regain their shape after interaction. We therefore study results for stationary solutions in the next sections.

Results for stationary wave solutions
In order to obtain the stationary wave solutions, we solve Eq. (18) numerically for y = u ,ξ . Thereby, the traveling velocity c can be arbitrarily chosen. An integration with respect to ξ yields finally the solution u, which is used as initial condition for Eq. (10). Since the scheme described in Section 4 needs also an initial condition at the time point t −1 , the corresponding values at this time have to be computed as well. Because the stationary solution propagates with the velocity c, the solution u −1 at time point t −1 can be computed by shifting the stationary solution in space by c t −1 .
From Eq. (18), we obtain the corresponding phase portrait for u ,ξ and u ,ξξ as shown in Fig. 6 for c = 100 m/s, where P 1 and P 2 denote the two existing fixed points and S denotes the saddle point. Starting from this saddle point, we can find two homoclinic orbits, which we refer to as left homoclinic orbit and right homoclinic orbit.
In the following investigation of stationary solutions, we look at four cases: • Stationary solution, which is located in the phase portrait "outside homoclinic orbit, close" (in the following denoted as case 1 ).
• Stationary solution, which is located in the phase portrait "outside homoclinic orbit, far" (in the following denoted as case 2 ).
• Stationary solution, which is located on the left homoclinic orbit (in the following denoted as case 3 ).
• Stationary solution, which is located in the phase portrait "inside homoclinic orbit" (in the following denoted as case 4 ).
The four corresponding trajectories in the phase portrait are shown in Fig. 7a. Thereby, the blue colored trajectory corresponds to case 1, the black colored trajectory to case 2, the dashed red colored trajectory to case 3 and the yellow colored trajectory case 4. Fig. 7b shows the behavior of the trajectories close to the origin of the phase portrait. The trajectories of case 1 and case 3 do not lie on top of each other. In case 1 the trajectory encircles both homoclinic orbits, whereas in case 3 it is located on the left homoclinic orbit. Next, we show the temporal evolution of the four solutions related to cases 1-4 in detail. The corresponding results are shown in Figs. 8-11. In In these results, we can identify two periodic solutions, which correspond to case 1 and case 2 . Clearly, the periodic pattern converges to a more edged pattern as the periodic solution gets closer to the homoclinic orbits. Starting at one of the homoclinic orbits, we obtain the so-called kink solution or kink wave, which also appears on a homoclinic orbit for the sine-Gordon equation [43,44], see appendix Appendix B. This solution is not periodic anymore and stays constant before and after the kink in solution. Another solution is obtained when starting inside a homoclinic orbit, like in case 4 . Then, we obtain a descending periodic pattern in the stationary solution if starting inside the left homoclinic orbit and an ascending periodic pattern, if starting inside the right homoclinic orbit. For case 4, modified periodic boundary conditions have been used in order to periodically continue the pattern at the boundaries. For this case, the derivatives of the solution at the left and right end of the spatial computation domain are equal and the difference between the corresponding function values of the solution are constant in time. Based on this, the difference at the right end of the spatial domain is added such that there is a periodic transition from the right to the left end value of the solution.
The amplitude of the stationary waves u(ξ) can be controlled using the parameter c. In addition, the width of the stationary waves depends on the distance to a homoclinic orbit. The closer the phase-portrait trajectories (u ,ξ , u ,ξξ ) are to a homoclinic orbit, as shown in Fig. 6, the wider the stationary waves u(ξ) are.
As can be seen in all figures, the form of the used initial condition is maintained in the corresponding temporal evolution. In other words, different parts of the solution do not interact with each other and stay stationary, which is to be expected. This implies that such waves in reality may be able to propagate over large distances and, when excited by an earthquake, even reach the surface. This is particularly relevant for the kink wave, which may be excited due to a sudden permanent local displacement associated with sliding along a fault.

Interaction of kink solutions
The result of interaction of two kink wave solutions is shown in Fig. 12, for c = 100 m/s. We observe an interesting nonlinear behaviour, which is different from classical soliton collision experiments where localized solutions regain their complete shape after collision, as for example is the case for the Korteweg-de Vries equation [e.g., 45]. Compared to the initial condition, the displacement amplitude at the end of the simulation has become slightly larger at the left and right sides of the considered spatial domain; the displacement has the value of 3.98 mm compared to 3.84 mm for the initial condition. Moreover, as shown in Fig. 12c, the plateau at the end of the simulation is slightly bent downwards compared to the initial condition, which is flat (i.e., zero) in the spatial interval from z = −9 m to z = −6 m. This difference compared to the corresponding values of the initial condition (i.e., stationary waves) has been analyzed for different spatial and temporal discretizations; for each of the considered discretizations, the difference did not vanish. Therefore, we can conclude that the kink solutions change their shape due to interaction. Hence, the kink wave cannot be referred to as a true soliton.

Discussion
In this section, we discuss two aspects regarding the model adopted in this paper: a potential improvement to overcome the presence of evanescent waves and a potential way to obtain closed-form expressions for the additional elastic moduli in the equation of motion. Regarding the first point, it can be verified that the linear version of our model is stable for all possible initial conditions (i.e., for all wavenumbers) and that the group velocity associated with the propagating wave is finite for all wavenumbers. However, the model does allow the presence of an evanescent wave through which information can travel at infinite speed. This is a short-coming of the model as it does not comply with the so-called strict/Einstein's causality condition [15]. It is expected that the nonlinear gradient elasticity model has this short-coming too. For the linear problem, it can be overcome by adding a term with fourthorder time derivative in the equation of motion so as to include the effect of micro-scale inertia (thus balancing the highest order of spatial and temporal

derivatives)
; then, all waves appear to have finite propagation speed [15]. Such a term is often included in the equation of motion of nonlinear gradient elasticity models too [32].
Secondly, we recall that the equation of motion of the nonlinear gradient elasticity model has been obtained using an assumed strain dependence of the additional elastic moduli (Eq. (8)). Generally speaking, the equation of a higher-order gradient elasticity model can be derived starting from a periodically layered classical continuum and applying an asymptotic homogenization procedure [24,25]. By applying such a procedure for a non-linear periodically layered classical continuum, closed-form expressions may be found for the additional elastic moduli (of the homogenized medium) in terms of the strain-dependent shear moduli of the original layers and their densities. It may also reveal how the time and length scales should be chosen, and with that the coefficients B 1 and B 2 . Finally, it may naturally lead to the inclusion of the fourth-order time-derivative term representing the effect of micro-scale inertia referred to above.

Conclusions
The aim of this paper was to investigate the existence of localized stationary waves in the shallow subsurface with constitutive behaviour governed by the hyperbolic model, implying that the shear modulus is strain dependent (i.e., non-polynomial nonlinearity). To this end, we derived a novel equation of motion for a nonlinear gradient elasticity model, which is of the Boussinesq-type. The higher-order gradient terms (compared to those in the classical wave equation) capture the effect of small-scale soil heterogeneity/micro-structure, which introduces dispersive effects particularly for the shorter waves; the dispersion prohibits the formation of jumps, which leads to physically realisable solutions. In order to obtain stationary solutions of the equation of motion, we introduced the moving reference frame together with the stationarity assumption, which yields an ordinary differential equation. We also presented a novel numerical scheme to solve the nonlinear equation of motion in space and time up to an order of accuracy of O(∆t 2 + ∆z 2 ), which exploits the structure of the partial differential equation in order to simplify the computation of the spatial finite-difference approximations.
Periodic (with and without a descending trend) as well as localized stationary waves were found by analyzing the above mentioned ordinary differential equation in the phase portrait, and integrating it along the different trajectories. The localized stationary wave is in fact a kink wave and was obtained by integration along a homoclinic orbit. Generally speaking, the closer the trajectory lies to a homoclinic orbit, the sharper the edges of the corresponding periodic stationary wave and the larger its width (i.e., period). The numerical scheme was used to verify the stationary character of the waves. In addition, it was used to study the propagation of arbitrary initial pulses, which clearly reveals the influence of the nonlinearity: straindependent speed in general and, as a result, sharpening of the pulses. Finally, we simulated using the numerical scheme a collision experiment in which two kink waves propagate in opposite direction and pass each other. We found that, after interaction, their original shapes are not recovered. Therefore, we conclude that the kink wave identified in this work is in fact not a true soliton.
Even though not being a true soliton, the kink wave, which may have high amplitude, can propagate through the soil column and potentially reach the surface depending on the strength of the material and geometrical damping mechanisms (which have not been considered). Therefore, seismic site response analyses should not a priori exclude the presence of such localized stationary waves. Follow-up research should be devoted to quantifying the influence of the damping mechanisms on the decay of the waves for specific soil profiles; here, the Masing model for hysteresis could serve to include the material damping [46,47]. (A.2) Collecting all unknown values of u at the timepoint t n+1 at the left-hand side and the rest at the right-hand side leads to a system of linear algebraic equations, which has to be solved for each time step. In order to simulate the solution in an open domain, absorbing boundary conditions have been used.  Similarly as for the equation of motion (Eq. (10)) of the nonlinear gradient elasticity model, they are located at a homoclinic orbit of the corresponding ordinary differential equation for the transformed coordinate ξ = z − ct.