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 behavior 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.

out 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.
Keywords Nonlinear gradient elasticity model · Stationary waves · Localized kink wave · Homoclinic orbit · Wave interaction

Introduction
For the prediction of the so-called seismic site response -the response of the top soil layers of the earthinduced 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 modeled assuming constant shear modulus and material damping ratio, respectively [2]. The actual val-ues 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 behavior 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 Lovetype 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 behavior 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 gra-dient terms are naturally obtained using asymptotic homogenization techniques for periodically inhomogeneous media [15,[24][25][26][27]. Dispersive effects significantly influence the behavior 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 realizable solutions. So-called higher-order dispersion correction and higher-order dispersive nonlinearity has been discussed in the context of the wellknown Korteweg-de Vries equation [28]. For higherorder gradient elasticity continua, the existence and properties of localized stationary waves have been studied for polynomial nonlinearity in both the macro-scale and the micro-scale 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 (Sect. 2). An ordinary differential equation from which stationary solutions for the nonlinear gradient elasticity model can be obtained is derived in Sect. 3. Periodic (with and without a descending trend) as well as localized stationary waves are found (Sect. 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 for example, 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 (Sect. 5.1) as well as to check the stationarity of the stationary waves identified (Sect. 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 (Sect. 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] ρ 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 stressstrain 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 medium, and 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 [6,38]

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 integrodifferential 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 (secant) 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., [39]). 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 Eqs. (3) and (6) into Eq. (1) yields the corresponding equation of motion: For simplicity, we now relate the strain-dependent additional elastic moduli to the conventional straindependent 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. Throughout the paper, the subscript "0" indicates that the quantity relates to small-strain/linear behavior. Using Eqs. (8) and (9), Eq. (7) can be written as which is the final form of the equation of motion of the nonlinear gradient elasticity model. In this work, we use the hyperbolic soil model typically employed for seismic site response analyses, with the following expression for the strain-dependent 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. (8)) Using the same definition for T (Eq. (9)), the stressstrain relation Eq. (6) reduces to (the z, t dependence is omitted for brevity) The corresponding equation of motion reads as follows: This linear equation was derived before by Metrikine and Askes [40] from a discrete model using a continualization procedure. A similar equation was used by Georgiadis et al. [41] 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 : 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 Sect. 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 Eqs. (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 ): 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 third-order 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 third-order 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 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 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 Table 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 an 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 Sect. 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 behavior. 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 Gaussian-cosine 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 inter- act 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 Sect. 4 needs also an initial 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, 9, 10, and 11. In each figure, the temporal evolution of the solution is shown and the initial condition is compared to the solution at the end of the simulation. In each computation, the velocity of the waves is chosen as c = 100 m/s. For the computation of each solution, periodic boundary conditions have been used. 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 [42,43], see "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 behavior, 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 [44, e.g.,]. 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 shortcoming 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 shortcoming too. For the linear problem, it can be overcome by adding a term with fourth-order 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 nonlinear 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 behavior 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 Boussinesqtype. The higher-order gradient terms (compared to those in the classical wave equation) capture the effect of small-scale soil heterogeneity/microstructure, which introduces dispersive effects particularly for the shorter waves; the dispersion prohibits the formation of jumps, which leads to physically realizable 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: strain-dependent 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 [45,46].
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Com- mons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

A Numerical scheme for the linear equation of motion waves
In order to solve Eq. (14), which is expressed as a finite difference scheme is introduced. Applying explicit finite difference approximations, Eq. (34) can be written as  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.

B Stationary solution of the sine-Gordon equation
A well-known equation which exhibits the so-called kink solitons is the sine-Gordon equation [42,43]. In this appendix, we demonstrate that such a solution can be obtained by introducing a moving reference frame, assuming stationarity, and integrating the thus obtained equation along a homoclinic orbit in the phase portrait. In this paper, we applied essentially the same procedure to get localized solutions of the equation describing sta-tionary waves in the nonlinear gradient elasticity model (Eq. (18)).
We show in the following relations between kink solutions of the equation of motion of the nonlinear gradient elasticity model (Eq. (10)) and kink solutions of the sine-Gordon equation.
Since it was realized that the sine-Gordon equation led to kink solitons, the importance of this equation greatly increased. The sine-Gordon equation has several important physical applications [47][48][49][50]. It can be stated as u ,tt − u ,zz + sin u = 0, where z ∈ R denotes the space variable, t ∈ R denotes time, and u := u(z, t).
Using the same procedure as in Sec. 3, we apply the transformation ξ = z − ct in order to determine stationary solutions of the sine-Gordon Eq. (36). With this Fig. 13 Homoclinic orbit for the ordinary differential equation (37) with saddle point S and fixed point P transformation, we obtain from Eq. (36) the following ordinary differential equation For Eq. (37), the homoclinic orbit on a cylindrical phase space, which is 2π -periodic with respect to u, is given by (cf., [51]) This homoclinic orbit is shown in Fig. 13. Starting at the upper or lower part of this homoclinic orbit close to the saddle point S and solving Eq. (37) numerically, we obtain the so-called kink or anti-kink solution, respectively. By solving Eq. (37) one can also show that the analytical expression for the kink and anti-kink solutions is given by where ξ 0 denotes the center position. A numerically obtained anti-kink solution of the sine-Gordon Eq. (37) is shown in Fig. 14 for the case c = 0.8, and it clearly coincides with the corresponding analytical solution from Eq. (39).
These kink or anti-kink solutions are localized stationary solutions of the sine-Gordon Eq. (36). 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.