Unintended and unwanted high frequency oscillation motion is commonly observed in small-scale experiments and in numerical simulations when soil is subjected to harmonic input motions. This high frequency motion has been often attributed to the drawbacks of actuating systems in experimental setups and to numerical noise in computational analyses. This work presents introductory consideration supporting the hypothetical idea that the recorded and the computed high frequency oscillation motion can possibly be the consequence of an unrecognized before physical phenomenon of soil elastic waves released upon unloading due to soil inherent hysteretic stress-strain behaviour and affecting the steady state response to harmonic excitation. To this aim, simplified numerical studies representative of the most basic soil mechanical properties are carried out. The results reveal the presence and the potential importance of soil-released unloading elastic waves when understanding soil free field response and structural response in numerical simulations of small-scale experimental setups. Chosen numerical cases are compared with available examples of experimental works on dry sand from the literature. The discussion provides further comments, including the potential importance of the soil-released elastic waves in real earthquakes of large magnitudes. Keywords: soil dynamics, wave propagation, soil nonlinearity, finite element modelling, soil constitutive modelling, soil elastic waves. 1. Introduction The recent advances in seismic geotechnical engineering are provided by high quality small-scale experimental works in flexible soil containers and advanced numerical studies. The most important aims of the physical modelling are twofold: improved understanding of soil-structure interaction during seismic events and creating reference baseline for validation of the numerical modelling studies, thus, allowing their use in the subsequent predictive studies of full-scale problems. Although, advanced physical and numerical modelling are considered as reliable ways of investigating seismic behaviour of soil, both these methodologies are often affected by the presence of unexpected high frequency oscillations in recorded or computed motions. Physical modelling of soil seismic behaviour is typically achieved in flexible soil containers placed on shaking tables in 1g (e.g. Durante et al. [1]) or centrifuge experimental setups (e.g. Lanzano et al. [2]). Such experimental studies often record unexpected and undesired high frequency content even though they used simplified sinusoidal input motions of a single frequency applied at base. For example, high frequency content in a form of a regular pattern of higher harmonics, e.g. 3ω, 5ω, 7ω etc. (where ω is the driving frequency), was observed by many research centres participating in the LEAP- GWU-2015 Project (Kutter et al. [3]) and LEAP-UCD-2017 Project (Kutter et al. [4]) on saturated soil slopes. Experimental work on dry soil also revealed the presence of undesired high harmonics for an intended perfect sinusoidal input motion (e.g. Lanzano et al. [2]; Conti & Viggiani [5]; Conti et al. [6]; Abate & Massimino [7]). Sometimes the high harmonics of a regular pattern 3ω, 5ω, 7ω etc. were also recorded at base level. For example, Madabhushi [8] showed example input motion recorded at base in a centrifuge test to comprise harmonics in the regular form of ω, 3ω, 5ω, 7ω etc. As a result of this observation, the presence of high harmonics is often attributed also to the actuating systems (Brennan et al. [9]). In more detail, Kutter et al. [3] and Manandhar et al. [10] indicated difficulties in controlling high frequency content in 2 actuators, Lanzano et al. [2] pointed out at potential resonance with high natural modes of the shakers, whereas Yao et al. [11] showed that high harmonics are recorded on an empty shaking table for a very high amplitude input motions and they attributed these high harmonics to nonlinearities in the actuating system. The presence of high frequency content recorded in the experimental work has also been explained as a result of physical phenomena. In studies regarding saturated soil, high frequency content could be associated with ‘de-liquefaction shock waves’ (Kutter & Wilson [12]) or in other words with dilation, sudden reduction in excess pore pressure and consequent sudden increase in soil stiffness (Bonilla et al. [13]; Roten et al. [14]; McAllister et al. [15]; Wang Gang et al. [16]). On the other hand, high frequencies occurring in the experimental setups containing dry soil were attempted to be explained by ‘soil fluidisation’ (Dar [17]), shear band development at shallow soil depths (Gajo & Muir Wood [18]) or attributed to be a proof of ‘pounding’ between soil and piles (Chau et al. [19]). In general, one is reminded that any nonlinearity in a physical system is expected to generate high harmonics for a single driving frequency of ω (e.g. Nekorkin [20]). For instance, high harmonics for perfect sine input motion were observed in frictional base isolation systems (Kelly [21]; Fan et al. [22]; Wiebe & Christopoulos [23]) or in a harmonic oscillator with sliding friction (Vittorino et al. [24]), where abrupt changes in stiffness and changes in the direction of the friction force induced nonlinearity and led to the generation of high frequency motion of the oscillator. Some numerical studies in the past showed that soil inherent nonlinearity can also generate high harmonics for single sine input motions (Pavlenko [25]; Pavlenko & Irikura [26]; Mercado et al. [27]; Veeraraghan et al. [28]). In detail, Pavlenko [25] and Pavlenko & Irikura [26] used a simple hysteretic soil model and monochromatic seismic input motions to show regular high frequency patterns for an example soil profile and for the 1995 Kobe earthquake site, respectively. A recent study by Mercado et al. [27] showed similar conclusions with the amount of harmonics being dependent on the level of nonlinearity based on comparisons of experimental data from centrifuge on saturated soil with single-phase numerical simulations obtained with a simplified hyperbolic, backbone soil model. The presence of high harmonics was shown to follow the pattern of exponential decay in the amount of the consecutive harmonics as a result of the distortion of a ‘sinusoidal wave’ towards a ‘square wave’. Finally, Veeraraaghan et al. [28] approximated a soil column as a single nonlinear element without considering wave propagation and showed that the shape of the stress-strain hysteresis determines the presence of high frequency components. To sum up, some physical explanations to the observed high harmonics have been attempted by the previous studies, however the theoretical analyses proposed so far are not fully consistent with each other, are based only on simplified soil models, do not include suitable comparisons with experimental evidence. Moreover, these studies explain only the presence of high frequency components occurring in the evaluated spectral responses; however, they do not explain the origin of the potential presence of additional waves generated in soil. The occurrence of high frequency content appears also in numerical studies. For instance, the performance of the numerical predictions in the LEAP-UCD-2017 project was presented by Manzari et al. [29]. This study highlighted the fact that some constitutive models were able to represent the high frequency content in a similar way as obtained in the centrifuge tests. Some other constitutive models were able to damp out high frequencies. Therefore, the presence or absence of high harmonics was attributed more to the implemented numerical or viscous damping and less to the constitutive models themselves. To the best of the authors’ knowledge, majority of the predicting teams used, as input motion at the base of the soil container, the signal recorded in centrifuge tests, thus, including high harmonics. Similarly, Bilotta et al. [30] showed that different predicting teams chose different methods to remove considered as ‘spurious’ high frequencies, e.g. by means of result filtering or application of viscous damping, and were able to obtain comparable response spectra at soil surface with experimental data. A more specific study, dealing strictly with the occurrence of high frequency content in numerical studies, was presented by Tsiapas & Bouckovalas [31]. They observed that the source of the high frequency oscillations can be attributed to ‘abrupt stiffness changes due to dilation and/or unloading-reloading’, attributed such frequencies to numerical noise and proposed a filtering method to remove such noise from the computations of dry and saturated soil. The problem of high frequency content, in recorded or computed motions, appears in both, physical and numerical modelling studies. This paper shows a potential explanation to the sources of such high frequencies based on the hypothetical idea of the release of soil elastic waves upon unloading in nonlinear hysteretic soil. To recall, the presence of elastic waves in the transient solution of the dynamic response of a system subjected to harmonic input motion is well recognized (e.g. Kramer [32]). Such waves are needed in the response of the system to satisfy the initial conditions, appear solely in the transient part of the response during which are progressively damped out prior to the steady state response is reached where the solution contains only the frequency of the input motion. The occurrence of elastic waves in the steady state solution has not been recognized in any dynamic system, to the best of the Authors’ knowledge. On the other hand, this work points out towards the hypothetical possibility of the release of soil unloading elastic waves in the steady state solution as a result of abrupt changes of stiffness from loading to unloading in nonlinear hysteretic soil. In general, the 3 importance of unloading waves is well recognized in solid mechanics (e.g. Nowacki [33], Wang [34]). On the other hand, the importance of unloading waves in wave propagation in soil has not been studied explicitly. Some very initial numerical studies on the effects of unloading waves were limited only to a soil column analysed by Fellin [35] and a brief study by Song et al. [36]. The former considered a compressive wave propagation under a single cycle of loading and unloading and showed that a shock wave forms as a result of soil nonlinearity in oedometer-like conditions. The latter study (Song et al. [36]) considered a semi-infinite column with a rectangular shock pulse with loading and unloading fronts, which interact with each other due to the unloading wave being faster and chasing the loading wave. However, the semi-infinite column was modelled with linearly hardening plasticity on loading and pure elasticity on unloading, thus not being representative of real soil behaviour. To sum up, these studies did not: address in detail the propagation of unloading waves in soil, account for inherent depth-dependence of soil stiffness or compare the simulations with experimental data. This paper suggests a potential explanation of the presence of high frequency oscillations in experimental and numerical studies as possibly related to the hypothetical, unknown before, physical phenomena of soil elastic waves released upon unloading (i.e. soil unloading elastic waves) and ‘trapped’ in a soil column modelled with a nonlinear hysteretic material. The presented study comprises primarily numerical studies compared in a few cases with experimental examples from the literature (Durante [37], Dar [17]). In more detail, Section 2 of the paper is dedicated to the methodology, including introduction to the conducted numerical studies and a short description of the experimental setups available to the authors. The results are presented in Section 3. Firstly Section 3.1 shows a numerical study of a semi- infinite column modelled with a simplified depth independent constitutive law and subjected to sinusoidal input motion at base. This section is a general introduction to the propagation of unloading waves in a material of a hysteretic type of stress- strain behaviour. It allows to build a link between previous works and the findings of the presented paper, mainly regarding which high frequency components in the spectral response are due to the numerical issues as a result of the occurrence of strain discontinuity, which due to wave distortion towards a ‘square wave’ as a result of soil inherent nonlinearity (i.e. as per previous works [25], [27]) and which can possibly be representative of the actual additional waves generated in soil due to the physical response of a nonlinear hysteretic soil. Section 3.2 presents the numerical studies of a soil column of the height representative of the soil container used in the experimental works. Herein, mainly numerical consideration, supported where available by experimental data from literature, is given to the hypothetical idea of the existence of soil unloading elastic waves, primarily in the steady state solution to harmonic excitation, but also possibly in scaled earthquake input motions. Finally, Section 3.3 shows a numerical study on a boundary value problem with some comparisons with the available experimental data of a similar testing setup to reveal how structural response may be affected due to the release of soil unloading elastic waves. The discussion in Section 4 follows with further comments on the consequences of the potential existence of soil unloading elastic waves in the dynamic response of soil, including their potential importance when understanding some observations in large magnitude earthquakes from the past. Finally, the most important findings are listed in Section 5 with the conclusions of this paper. 2. Methodology 2.1. Numerical studies The numerical studies in this paper are based primarily on simplified soil constitutive models (Section 3.1 and 3.2). An advanced soil constitutive model, implemented in the multiaxial stress space, is used when comparing numerical results of a boundary value problem with a pile (Section 3.3). The simplified 1D soil constitutive models, incorporating the basic ingredients of the advanced model, have been chosen to limit constitutive features only to those sufficient to discuss the nature of, firstly, the propagation of unloading waves in a semi-infinite column (Section 3.1), secondly, the release and ‘entrapment’ of soil unloading elastic waves in a soil column (Section 3.2) The numerical models in Section 3.2 include a simple geometry of a 0.8m high soil column representative of the soil specimen in small scale experimental works carried out in the past. The soil density has been set to 1332kg/m3 unless specified otherwise. The boundary conditions allow full wave reflections at the base since the soil container is placed on a rigid shaking table. The side nodes of the soil column are linked by tie connectors ensuring the same lateral and vertical displacements of the nodes at the same height. The soil column has been discretized into 32 eight-node, quadratic finite elements of an equal size of 0.025m. This mesh refinement ensures that the element size is less than the standard ‘rule of thumb’ and accounts for the refined mesh size criteria proposed by Watanabe et al. [38]. In more detail, taking the highest frequency of interest to be 80Hz as per the one of the reference experimental works (Durante [37]), the slowest elastic shear wave velocity at shallow depth being approximately 40m/s, one obtains the maximum distance in grid spacing of 0.05m (i.e. the quadratic element size of 0.1m) for the standard discretization rules. Watanabe et al. [38] showed that such grid 4 spacing/element size may need to be reduced to account for soil nonlinearity and reduction in the shear wave velocity due to developing plasticity. As such, the quadratic element size of 0.025m adopted in this study can be considered much smaller than the standard discretization practise and within the updated discretization rules proposed by Watanabe et al. [38]. The time increment has been chosen to ensure the fastest elastic shear wave in the soil column (i.e. around 120m/s at soil base) does not reach two consecutive nodes within the same time increment. As such, the time increment of 2·10-5 seconds has been chosen for the advanced soil constitutive model implemented with the implicit integration scheme, and 1·10-5 seconds for the simplified soil constitutive model, since the constitutive law is integrated with an explicit integration scheme. The numerical analyses have assumed a homogeneous soil layer within the soil container. This is admissible because the difference between the mechanical properties of the two layers in the experimental setup by Durante [37] is relatively small and the experimental work by Dar [17] used a single homogenous soil layer. The soil column in the numerical studies is subjected to horizontal acceleration time histories applied at base. The analysed input motions include sinusoidal input motions and scaled earthquake input motion. Perfectly sinusoidal input motions (i.e. of a single harmonic only) are introduced in a smooth manner using cycles with steadily increasing amplitude in order to limit the introduction of spurious waves due to the transient response prior to the steady state response is reached. 2.1.1. Simplified constitutive models Two simplified constitutive models are used in this work. First constitutive model is used in a semi-infinite soil column thus disregarding depth effects. Second constitutive model is used in a soil column of 0.8m height, thus representative of the height of soil in the small-scale experimental works [17, 37]. Firstly, Section 3.1 of this paper shows a simplified numerical study on the propagation of a sinusoidal wave in a semi-infinite soil column. For this purpose, taking inspiration from the typical stress-strain response of soils, a simple, 1D constitutive model based on a hyperbolic response and leading to nonlinear hysteretic behaviour was implemented as a user- defined subroutine (UMAT) of the FEM code Abaqus (Dassault Systèmes [39]). The column is modelled with a material of the following constitutive law: 𝜏̇=ቐ(ିఛౢ౟ౣିத)మ ଶ·஻·ఛౢ౟ౣ 𝛾̇,𝑓𝑜𝑟 𝛾̇<0 (ఛౢ౟ౣିத)మ ଶ·஻·ఛౢ౟ౣ 𝛾̇,𝑓𝑜𝑟 𝛾̇>0, (eq. 1) where τ lim is the limiting (positive) value of the stress, which is assumed constant within the soil column and B is the constitutive parameter defining the stiffness. The model response can be imagined to be representative of the shear stress versus shear strain response in a simple shear test, in which the volumetric behaviour is not considered. For the chosen example values of τ lim=4000Pa and B=0.001, typical stress-strain response is plotted on Figure 1. The semi-infinite column (Section 3.1) was modelled with the same input parameters and assuming density of 2000kg/m3. The second simplified 1D law was defined with the following incremental expression, i.e. amended expression (1) to include depth dependence of stiffness and strength with a square root of depth: 𝜏̇=ቐ(ିఛ೗೔೘(ௗ)ିఛ)మ ଶ·஻·ఛ೗೔೘(ௗ) 𝛾̇,𝑓𝑜𝑟 𝛾̇< 0 (ఛ೗೔೘(ௗ)ିఛ)మ ଶ·஻·ఛ೗೔೘(ௗ) 𝛾̇,𝑓𝑜𝑟 𝛾̇> 0, where 𝜏௟௜௠(𝑑)=√ௗ √ௗ೘ೌೣ𝜏௟௜௠ (eq. 2) where d is the current depth in the soil column, d max is the maximum depth, i.e. in this work 0.8m for the simulated soil column in Section 3.2, τ lim and B are the model parameters, chosen in this work to be τ lim=6000Pa and B=0.00016. The simplified model can be deemed representative of any family of constitutive models in the finite element method (e.g. originating from elastoplasticity or hypoplasticity) which predict a hysteretic type of stress-strain behaviour. The simplified model with depth dependence (eq. 2) has been calibrated in a way to represent approximately the soil response in the small-scale experiments. Typical stress-strain response of the model at two chosen depths is shown on Figure 2. The initial G 0 shear stiffness profile with the representative shear wave velocity profile predicted by the constitutive law (eq. 2) is shown on Figure 3. It can be observed that the G 0 profile fits closely the well-known empirical evaluation from Hardin & Drnevich [40]. The fastest time of a wave travelling from the bottom to the top of the soil column 5 for such definition of G 0 profile is approximately 0.0085sec. The first natural frequency of the soil column for an average shear wave velocity of 95m/s is slightly below 30Hz as evaluated from the well-known analytical expression: f=v s/4H (Kramer [32]), where v s is the average shear wave velocity, H is the height of a soil column. This assessment has been verified numerically as shown in Appendix B, where a slightly higher value of 33Hz has been indicated as the natural frequency of the soil column when subjected to an input motion of a flat frequency spectra (i.e. a ‘white noise’ type of motion). In addition, the second soil natural frequency was evaluated to be around 88.5Hz as shown in Appendix B, which is close again to the value of 89Hz computed using the simple empirical approach. This evaluation of soil natural frequencies will be shown important when presenting the idea of the release of soil unloading elastic waves. Note that the authors in this work preferred simplified constitutive models for the sake of identifying the basic ingredients needed to observe the phenomenon of soil unloading elastic waves, however the general patterns of the shown results can easily be replicated when using other constitutive approaches provided those approaches account for soil inherent irreversibility in a form of a material hysteresis (e.g. by using the Severn-Trent model by Gajo [41], a hypoplastic sand model by Von Wolffersdorff [42] or one of the SANISAND models by Dafalias & Manzari [43], as shown initially by Kowalczyk [44]). Figure 1. Typical response of the simplified constitutive model (eq. 1) under ramped up sinusoidal input motion. Figure 2. Typical response of the simplified 1D constitutive law (eq. 2) under ramped up sinusoidal input motion: a) at mid-depth, b) at base of the 0.8m high soil column. 6 Figure 3. Stiffness characteristics of the modelled soil column of 0.8m height: a) the initial G 0 shear stiffness profile (compared with the empirical evaluation by Hardin & Drnevich [40]) and b) shear wave velocity profile. 2.1.2. Advanced soil constitutive model An advanced soil constitutive model, namely the Severn-Trent sand model (Gajo [41]) has been used in the numerical simulations of the small-scale experiments with piles of Durante [37]. Such constitutive model incorporates all the features of the simplified model, but in addition offers much more accurate representation of the real soil behaviour under cyclic loading. Short description of the Severn-Trent model used herein is given below. Note that other advanced soil constitutive models, such as a hypoplastic sand model (Von Wolffersdorff [42]) or one of the SANISAND models (i.e. Dafalias & Manzari [43]), could equally be used in predicting similar trends of the dynamic response of soil, as shown in (Kowalczyk [44]). The Severn-Trent sand model in the version proposed by Gajo [41] is a classical kinematic hardening elastoplastic model. Its current model formulation is based on the original formulation by Gajo & Muir Wood [45][46]. The model focuses on deviatoric soil response, thus neglecting grain crushing under high pressures, and is based on the concepts of classical kinematic hardening (arranged in such a way to mimic the bounding surface plasticity approach) to keep some memory of the past loading history. Importantly, the Severn-Trent model is based on the well-known concepts such as: the critical state, the Mohr-Coulomb failure, the state parameter Ψ (Been & Jefferies [47]) to determine the dependence of soil state on density and pressure, and plastic hyperbolic stiffness dependence on the distance from the failure surface. Two key features of the Severn-Trent model are the use of a hyperelastic formulation for representing the soil small-strain stiffness response and the numerical implementation based on an implicit, back-Euler integration method. The original formulation of the model [41] includes an evolving fabric tensor accounting for changes of elastic stiffness due to plastic strain. This feature has been disabled, i.e. fabric tensor is kept constant, in this study for the sake of computational efficiency. Further details of the model formulation are not presented herein for the sake of brevity and the interested Reader is addressed to the cited works. The calibration of the constitutive model is based on the general guidelines given by Dietz & Muir Wood [48] for G/G0 curves in small scale experiments carried out in dry sand. These guidelines were supported by Seed & Idris limits [49] and laboratory experimental data obtained by Kokusho [50], both carried out on different sands, therefore can be deemed appropriate for sands in general. The calibration of the model parameters (Table A1) is shown on Fig. A1 in the Appendix A together with the validation of the calibration with a single element cyclic test compared with relevant laboratory data. 2.2 Benchmark experimental work used in comparisons Chosen numerical cases in this paper are compared with available experimental example tests of dry sand placed in a flexible soil container. The brief presentation of the reference experimental works is provided below. The work by Durante [37] was aimed primarily at investigating seismic soil-structure interaction for a group of five piles, however; free field measurements were also obtained. The soil container of 0.8m height was filled by using a pluviation method with two layers of Leighton Buzzard (LB) sand, the top ‘softer’ layer of LB sand fraction E, the bottom ‘stiffer’ layer of LB sand fractions B+E. The initial density of the top layer is 1332kg/m3 and the void ratio of around 0.9, whereas for the bottom layer it is 1800kg/m3 and 0.53, respectively. Nevertheless, note that the difference in terms of the 7 relative density between the two soil layers is rather small (approximately 25% for the top layer, 40% for the bottom layer). Therefore, the numerical studies assume modelling this experimental setup as a single homogeneous soil layer. The experimental measurements were filtered with a lowpass filter, 80Hz, 5th order. The soil natural frequency was estimated to lie within a range of 25-30Hz depending on the amplitude of the input motion varying from 0.01g to 0.1g. The applied input motions considered various input frequencies and moderate amplitudes of up to around 0.15g. More details on the experimental work can be found in the cited works (Durante [37], Durante et al. [1]) and are not repeated here. The experimental work by Dar [17] was carried out in a flexible soil container of 0.8m height as the experimental work of Durante [37]. The main difference with respect to the experiments by Durante [37] is that the soil container was filled with a single nonhomogeneous soil layer of LB sand, fraction B, deposited by using the pluviation technique. The achieved void ratio was evaluated to be around 0.6 resulting in the relative density representative of dense sand. The density of the deposited sand was around 1700kg/m3. The experimental measurements included free field response as well as investigation of soil-structure interaction. The natural frequency of the soil was evaluated by Dar [17] to be around 23Hz, for the amplitude of motion of 0.1g. The amplitudes of the applied input motions included strong motions as large as 1.0g. 3. Results 3.1 Semi-infinite soil column The aim of the numerical study in this section is to show the general effects of the propagation of fast unloading waves on the dynamic response of a semi-infinite column modelled with a material exhibiting hysteresis in the stress-strain behaviour (eq. 1). In case of a semi-infinite column the response is not affected by wave reflections as presented in the subsequent Section 3.2. This will be shown crucial for the actual observation of soil unloading elastic waves ‘trapped’ in the soil column representing the soil specimen in a soil container. Herein, the semi-infinite column is subjected to the input motion of the driving frequency of 25Hz and the maximum amplitude of 0.08g. The computed results are shown on Figure 4. First of all, distortion in the sine wave towards a ‘square wave’ can be seen in the computed accelerations with the increasing propagation length. The distortion of the sine wave is due to the discontinuity occurring in the computed strains. Such discontinuity is a result of fast unloading waves which erase the slowly propagating loading waves (i.e. the amplitude of the computed strains decreases with the increasing propagation length). The occurrence of strain discontinuity can be noticed at 1m and 2m heights and leads to a ‘jump’ in the computed accelerations. In addition, very high frequency oscillations (Fig. 4d) can be observed as typically expected in the finite element computations due to the occurrence of strain discontinuity (thus the very high frequency oscillations can be considered as numerical noise). The spectral response computed for the motion at the height of 2m and shown on Figure 5 reveals the presence of high harmonics of two origins. Firstly, high harmonics, i.e. high components of Fourier series, are needed to represent the wave distortion towards a ‘square wave’. In this case, the amplitude of the consecutive harmonics decreases with exponential decay (similarly to the observations by Mercado et al. [27]). The increased amount of very high frequency harmonics between 1500-2500Hz is representative of the numerical oscillations due to the occurrence of strain discontinuity. To sum up, this section has shown how the propagation of fast unloading waves results in the occurrence of strain discontinuity and the distortion of a sinusoidal wave towards a ‘square wave’ in the computed accelerations. Moreover, the strain discontinuity causes numerical problems and results in very high frequency numerical oscillations in the computed accelerations. Importantly for what will be discussed in Section 3.2, the high harmonics in the presented spectral response on Fig. 5 are not due to physically present additional waves but represent either the distortion of the sinusoidal wave (harmonics of the pattern ω, 3ω, 5ω etc, from 25Hz onwards with an exponential decay in their amounts) or computational problems with representing strain discontinuity (high frequencies from around 1500 to 2500Hz). 8 Figure 4. Horizontal acceleration and shear strain records: a) input motion, b) height 1.0m, c) height 2.0m, d) zoom in on the numerical oscillations computed for the height of 2.0m 9 Figure 5. Spectral horizontal acceleration evaluated at the height of 2m : a) up to 3000Hz (with very high frequency components due to numerical oscillations), b) up to 500Hz (to represent solely the wave distortion towards a ‘square wave’). 3.2 Soil column of 0.8m height This section presents results of the numerical studies obtained mainly with a simplified constitutive law defined by equation (2) when representing a soil column of 0.8m height. However firstly, the employed numerical algorithms are checked while modelling the soil column with linear elastic properties equivalent to the shear stiffness profile shown on Figure 3a (i.e. depth dependent stiffness) without any addition of damping in the system. Such a linear elastic soil column, of the first natural frequency of 33Hz is subjected to 10Hz input motion. For the smoothly introduced input motion (Fig. 6a), the response at the top of the column is slightly affected by the elastic waves induced due to the transient response. The elastic waves affecting the motion at the top of the column do not produce repetitive sine cycles (i.e. the consecutive sine cycles differ from each other). Since no damping is present in the model, the induced elastic waves remain ‘trapped’ in the column during the constant amplitude part of the input motion and after the driving force is ceased. On the other hand, when the motion is introduced in an abrupt manner (Fig. 6b), the computed response at the top of the column is heavily affected by the elastic waves of the transient response. In fact, after the driving force is ramped down to nul, the remaining elastic waves are characterised by large amplitude and more than a single natural frequency component (i.e. higher modes of free vibrations are present in the solution). This is confirmed when investigating the evaluated spectral response for both cases in the steady state and the ‘coda’ part of the motion (Fig. 6c and Fig. 6d). In general, the computed responses presented on Fig. 6 are expected, however here are useful for comparisons with responses given by a hysteretic material model as shown below. Nextly, the soil column is analysed under the same input motions but with a hysteretic material (eq. 2) with the G 0 initial shear stiffness as shown on Fig. 3a. Note that the amplitude of the input motion has been chosen in such a manner that the hysteresis loop is relatively narrow (i.e. little nonlinearity develops), thus remains fairly comparable with the linear elastic material used for the computations shown on Fig. 6. In contrast to the linear elastic calculations, such change in the material definition results in apparently much less expected response computed at the top of the soil column. Fig. 7a presents smoothly introduced input motion. The response at the top of the soil column is not strongly affected similarly to the case shown on Fig. 6a. The computed amplitudes are very similar for both, although somehow unexpectedly, the actual amplification of motion is slightly higher when using the nonlinear model when comparing with the linear model (however the very first cycles of response, for which the nonlinear model remains practically almost linear, the amplification appears to be the same for both). The very beginning of the motion is in fact similar to the one computed on Fig. 6a. However, later the material hysteretic damping appears somehow ‘inefficient’ in removing the high frequency oscillations. This high frequency motion appears now in the consecutive sine cycles with the same repetitive pattern, i.e. although containing high frequency oscillation waves, a ‘sort of’ steady state is apparently reached under the constant amplitude of the input motion. When the driving force is ceased, the high frequency oscillation waves are still present, however, in this situation the material damping appears ‘sufficient’ to slowly reduce the amplitude of the remaining waves. Fig. 7b shows the case of the input motion introduced in an abrupt manner. Again, the very first cycle of motion is comparable with this predicted by the elastic material (Fig. 6b), with strong presence of high frequency motion, thus possibly elastic waves due to transient response. Nevertheless, this time this high frequency motion diminishes slightly in the consecutive sine cycle, such as this time the material hysteretic damping in the soil column was ‘sufficient’. After a couple of cycles, again an apparent steady state solution emerges with no further damping out of the high frequency oscillation waves. The computed steady state 10 solutions for both cases (Fig. 7a and b) appear to be the same thus independent on whether the input motion is introduced smoothly (Fig. 7a) or abruptly (Fig. 7b). Consequently, so is the coda part of the signal after the driving force is damped out. The spectral evaluation of the computed responses is shown on Fig. 7c and 7d and confirms that the responses in the steady state and in the coda part of the computations, in fact, perfectly match for both, smoothly and abruptly, introduced input motions. Moreover, the coda part of the motion reveals that the remaining waves in the system appear to be the soil elastic waves, i.e. soil natural frequencies of 33Hz and 88Hz have been identified, with frequencies of 144Hz and 199Hz being probably roughly representative of higher modes of free vibrations (see Fig. 6c and 6d for comparisons). These elastic waves appear as ‘trapped’ in the soil due to reflections from the top free end and the base fixed end, and slowly being damped out. The peculiarities observed in the response of the soil column modelled with a hysteretic model (eq. 2) will be subject to further investigation when looking at loading cases representative of those experienced in a flexible soil container. Figure 6. Soil column modelled with elastic depth dependent material under harmonic (10Hz) input motion: a) computed response for smoothly introduced input motion, b) computed response for abruptly introduced input motion, c) evaluated spectral response when harmonic force is applied, d) evaluated spectral response in the coda part of the motion with harmonic force ceased. 11 Figure 7. Soil column modelled with hysteretic depth dependent material under harmonic (10Hz) input motion: a) computed response for smoothly introduced input motion, b) computed response for abruptly introduced input motion, c) evaluated spectral response in the steady state, d) evaluated spectral response in the coda part of the motion. The next numerical analyses are focused in more detail on what may be possibly experienced in small-scale experiments. To this aim, horizontal acceleration input motions of various driving frequencies and amplitudes are used herein in order to present how soil elastic waves can hypothetically be released in hysteretic nonlinear soil upon unloading during the steady state response of soil placed in a flexible soil container. Firstly, results for 10Hz input motion of the maximum amplitude 0.137g are shown. Fig. 8 reveals distortion of the sinusoidal wave, both in terms of the computed accelerations (Fig. 8a) as well as the computed shear strains (Fig. 8b). High frequency components of motion can be seen in the spectral response in Fig. 8c. Importantly, the harmonic of 30Hz (i.e. close to the soil first natural frequency) has a large amplitude, moreover, the harmonic of 90Hz (i.e. close to the soil second natural frequency of soil) is also significant, thus can possibly be representative of soil elastic waves as indicated for the previous example shown on Fig. 7. Note that the remining harmonics would be representative of the distortion of a sinusoidal wave towards a ‘square wave’ as shown by others (e.g. Pavlenko [25], Mercado et al. [27]) and in Section 3.1 of this work. On the other hand, the high harmonics in Fig. 8c are not of the ‘exponential decay’ pattern (e.g. herein the harmonic 90Hz is larger than the one 70Hz), thus the evaluated pattern of high harmonics cannot be representative solely of the distortion towards a ‘square wave’ but apparently could possibly be associated to the presence of soil elastic waves. The pattern of high harmonics in the evaluated response spectra is always in a form of ω, 3ω, 5ω etc since represents the steady state (i.e. repetitive sine cycles), thus the exact values of soil natural frequencies (i.e. 33Hz and 88Hz) cannot be shown here. Figure 8d shows that the presence of high frequency waves is also visible in the computed relative displacements at 12 the top of the soil. The hysteretic nonlinearity in soil and the abrupt changes from loading to unloading stiffnesses are shown on Figure 9. Finally, Figure 10 compares the computed accelerations at depth of around 40mm (filtered with 80Hz Butterworth filter as per the experimental measurements) and compares those with the available experimental data (Durante [37]) for the same input motion of 10Hz driving frequency and the amplitude of 0.137g. It can be observed that the amplitude of the motion has been captured correctly together with only little phase shift between the motion at base and top (as expected for the case of the driving frequency below the first soil natural frequency). High frequency motion is present in both; however, is apparently much more visible in the computations, possibly due to different damping characteristics in the system, for example damping induced by the physical presence of the soil container, not modelled numerically, may play a role. Moreover, the experimental records experience nonsymmetric response within a single sine cycle with stronger distortion observed on the downward branch of the measurements, possibly due to the accumulation of plastic shear strain from the previous loading scenarios, not accounted for in the simplified numerical studies. Figure 11 shows the results for 10Hz input motion with the amplitude of motion increased up to 0.2g. Such amplitude results in even more clearly visible high frequency oscillations, i.e. distortions in the computed accelerations (Fig. 11a) and shear strains (Fig. 11b), and the presence of high harmonics again not representative of the exponential decay pattern (see increased amounts of frequencies 90Hz, 130Hz and 150Hz on Fig. 11c). Herein, further observation can be made, i.e. the occurrence of strain discontinuity somehow similar to the one shown in a simpler case of the propagation of unloading waves in a semi-infinite column analysed in Section 3.1. The occurrence of such strain discontinuity (herein strong discontinuity) is shown in more detail on the enlarged view on Figure 12. Although very high frequency numerical oscillations (due to the problems of standard numerical algorithms to represent strain discontinuity as shown in Section 3.1) appear in the computed accelerations (Fig. 12a) and shear strains (Fig. 12b), regular ‘jumps’ in these quantities can also be observed. It can be noticed that the time distance between these ‘jumps’ at the soil column top, is approximately equal to the time of 0.017sec needed for the fastest wave (i.e. elastic wave) to travel from the top to the base, and reflecting back to the top. Note that following such reflection at base the shear strain moving upwards changes sign due to the phase change upon the reflection from the fixed end. The level of the induced nonlinearity for 10Hz 0.2g input motion shown on Figure 13 is higher than shown for the case of 10Hz 0.137g. Note that the amplitude of the reflected elastic wave is such that it causes unloading-reloading loops in the stress-strain response of soil at the top of the column (highlighted with circles in Figure 13a and shown in more detail in Figure 13c). Figure 14 shows the computed accelerations when the soil column is subjected to even higher amplitude of input motion of 0.63g. The computations are compared with the experimental data from the past (Dar [17]) since such high amplitude input motions were not introduced in more recent experimental works (e.g. Durante [37]) but where apparently the release of soil unloading elastic waves is expected to be more pronounced. The comparison between the computations and measurements indicates generally good agreement in terms of the amplitude and the presence of high frequency motion in the steady state response, with both, the simulations and experiments, being heavily affected by repetitive high frequency motion. The patterns of high frequencies are not fully consistent with each other; however, the observed responses are obtained for input motions of a single driving frequency (i.e. perfectly sinusoidal in the numerical study, and, from visual inspection, apparently perfectly sinusoidal also in the experiment by Dar [17]). Finally for the case of the 10Hz input motion, it can be observed that low amplitude input motion of 0.03g (Fig. 15) results in the response much less affected by the high frequency motion (due to the almost linear behaviour). In addition, the spectral response evaluated from the response in the time interval 0.3-0.8sec (Fig. 15c) reveals slight presence of the exact harmonic of 33Hz likely due to the initial transient response and little damping, thus this time not fully representative of the steady state response. The almost linear response in the soil column and little damping are evident when investigating the stress-strain behaviour shown on Fig. 16. The next analysed case deals with an input motion of the driving frequency of 25Hz and the maximum amplitude of 0.077g. In the computed accelerations and shear strains on Figures 17a and 17b, respectively, a ‘jump’ in the accelerations occurs due to the discontinuity in strains (thus also the occurrence of very high frequency numerical oscillations). Apparently for higher frequency of the input motion (25Hz) even a lower amplitude than for the 10Hz input motion (i.e. Fig. 11) is enough to obtain a strain discontinuity. On the other hand, this time the elastic wave reflected at the top towards the base of the column does not return back to the top of the column in the steady state cycles. This can be seen in the computed accelerations and shear strains but also in the relative displacements (the latter one this time in a form of a single frequency). The spectral response of the accelerations computed at the top shows high harmonics however these are due to the representation of the ‘jump’ in the accelerations and not due to physically present high frequency waves, for this reason the consecutive harmonics occur with exponential decay in their amounts. The presence of the soil elastic waves is revealed only in the coda part of the motion when the driving force is ceased. For example, the period of oscillations in the 13 lateral displacement changes from 0.04sec to 0.03sec (i.e. representative of an elastic wave of around 33Hz frequency). The fact that the steady state response in this case is not affected by the wave reflection is due to the fact that the time to travel down and up again is longer than the subsequent unloading happening in the column, which is the consequence of applying motion with higher frequency of 25Hz rather than 10Hz. Finally, Figure 18 compares the filtered computed accelerations with the recorded ones (Durante [37]). Generally, the amplitude of motion is similar as is the phase shift with some difference observed in the form of distortion of the sinusoidal wave. In this case, the distortion in the sine wave seems to appear in different places of the acceleration branches and may in fact be representative of two different aspects, i.e. wave distortion towards a ‘square wave’ (i.e. in the simulations), elastic wave arrival or some other experimental effect (in the experimental measurements). To complement the use of simplified sinusoidal input motions where the identification of soil unloading elastic waves might be easier, the authors present also a case of earthquake input motion. The chosen transient loading history is a scaled Tolmezzo earthquake as analysed by Durante [37]. The analysed duration includes the first 1.5sec of the time history. The computed and measured time histories are shown on Figure 19a and 19b. The experimental data is available only in a filtered format thus it is not possible to clarify if the generated high frequencies (Fig. 19f) are representative of the wave distortion towards a ‘square wave’ or due to the release of soil elastic waves. On the other hand, the computed accelerations can represent in a certain part of the motion (i.e. when the incident motion is of a lower frequency and amplitude) the presence of soil elastic waves ‘trapped’ in a soil column as indicated on Fig. 19c for the motion between 0.5-0.7sec. For completeness, the detail of the time history between 0.5-0.7sec is also shown for the experimental measurements (Fig. 19d) when one may also find patterns of wave distortion, resembling those potentially representative of soil elastic waves. Finally, the evaluated spectral responses at the top of the soil column confirm the presence of some high frequency components of motion even though these were absent at the base level. Note that the analysed case of a scaled earthquake does not present any longer the steady state response as per the previously analysed cases with harmonic excitation, however it reveals how soil elastic waves may be released (possibly also due to soil nonlinear hysteretic behaviour and not only due to the transient response to non-harmonic input motion) in ‘real earthquakes’ and remain ‘trapped’ in a stratum if a rigid or much stiffer formation constrains this stratum at base. To sum up, the results presented in this section reveal how high frequency motion in the presented computations can possibly be representative of an unrecognized before hypothetical physical phenomenon of soil unloading elastic waves released in a soil column due to soil inherent nonlinearity and irreversibility in a form of material hysteresis. The identified phenomenon is more prominent for input frequencies below the first soil natural frequency when the ‘entrapment’ of soil elastic waves in the soil column is observed. In fact, one could analyse a case with input motion of even lower driving frequency, which would allow for an increased number of elastic wave reflections in the soil column (omitted here for the sake of consistency with typical input motions applied in experimental works). In general, the release of soil unloading elastic waves can be deemed somehow similar to the presence of elastic waves in the transient response in any dynamic system under harmonic excitation (e.g. Kramer [32]) where the solution contains two types of frequencies of motion (i.e. the driving frequency and the natural frequencies). Apparently, unloading in nonlinear hysteretic soil and abrupt change from loading to unloading stiffness may possibly lead to a similar effect induced by the transient response and resulting in the release of elastic waves also in the steady state response. 14 Figure 8. Results for 10Hz input motion of the maximum amplitude of 0.137g: a) accelerations, b) shear strains, c) spectral response for the steady state response, d) relative displacements. 15 Figure 9. Stress strain curves for 10Hz input motion of the maximum amplitude of 0.137g up to 0.5s of the computed time history: a) top of the soil column, b) base of the soil column. Figure 10. Filtered horizontal accelerations at depth of 40mm for 10Hz input motion of the maximum amplitude of 0.137g: a) numerical computations, b) experimental measurements, c) numerical computations (zoom in), d) experimental measurements (zoom in). 16 Figure 11. Results for 10Hz input motion 0.2g: a) accelerations, b) shear strains, c) spectral response in the steady state, d) relative displacements. 17 Figure 12. Results for 10Hz input motion of 0.2g amplitude in zoom in: a) accelerations, b) shear strains. Figure 13. Stress strain curves for 10Hz input motion 0.2g up to 0.5s of computed time history: a) top of the soil column, b) base of the soil column, c) top of the soil column (zoom in on unloading-reloading strain cycle). 18 Figure 14. Horizontal accelerations at depth of 170mm for 10Hz input motion of the maximum amplitude of 0.63g: a) computations, b) experimental measurement data digitalized from (Dar [17]). 19 Figure 15. Results for 10Hz input motion 0.03g: a) accelerations, b) shear strains, c) spectral response, d) relative displacements. Figure 16. Stress strain curves for 10Hz input motion of varying amplitudes up to 0.5s of computed time history: a) top of soil column, b) base of soil column. 33Hz 20 Figure 17. Results for 25Hz input motion 0.077g: a) accelerations, b) shear strains, c) spectral response, d) relative displacement. 21 Figure 18. Horizontal accelerations at depth of 40mm for 25Hz input motion of the maximum amplitude of 0.077g comparison with experimental data Durante [37]. 22 Figure 19. The horizontal acceleration at the top and the bottom of the soil column for the scaled earthquake input motion: a) numerical computations (filtered), b) experimental measurements by Durante [37] (filtered), c) numerical computations between 0.5-0.7s (unfiltered), d) experimental measurements between 0.5-0.7s (filtered), e) evaluated spectral response for computations, f) evaluated spectral response for experimental measurements. 3.3 Boundary value problem of soil-structure interaction on piles Finally, this section presents a comparison of the numerical and experimental results for a boundary value problem on an example case of soil-structure interaction which reveals how soil generated high frequency motion, possibly representative of soil elastic waves, can impact on structural response. In this case, slender piles embedded in soil with a small mass attached to their tops (as a result of the attached fixing and measuring devices) are subject to a simple analysis. The soil and pile discretization is shown on Figure 20. Only a half of the flexible soil container has been modelled in order to reduce the computational times. The soil has been modelled assuming a homogenous soil layer, similarly to the studies shown in Section 3.2. The soil element size for the full 3D analysis has been optimised for the sake of the computational time. Instead of 32 equal size elements of 0.025m along the 0.8m height (as in the free field studies in Section 3.2), herein 22 quadratic brick elements have been used with the element size varying from 0.02m at the top to 0.06m at the bottom, therefore, accounting for slower waves propagating in soil at the top of the soil container. The chosen range of element size remains well within the standard meshing procedures as explained in Section 2.1. The boundary conditions have been specified at soil side and base nodes. The base nodes were constrained in the vertical direction with horizontal acceleration time history being applied. The long soil side nodes have been constrained in the horizontal direction perpendicular to the long side plane. Finally, the nodes on both short sides have been tied together 23 to ensure the same horizontal displacements on both sides to mimic the presence of the lateral constraint provided by the flexible soil container. Note that the 3D numerical study accounts for the presence of a single pile only, whereas the experimental setup of Durante [37], to which the numerical results are compared, contained five piles. The reason of modelling a single pile only is to show to the Reader that the observed oscillations in the accelerations of a pile can be related only to the high frequency motion generated in soil and not to the interaction between the piles. Nevertheless, previous numerical studies (Kowalczyk [44]) showed that the same effects can also be observed if the whole group of five piles is modelled or if a different seismic input motion is analysed. Figure 21 shows the comparison of the computed and measured horizontal accelerations on piles when soil specimen is subjected to 10Hz input motion of 0.137g. It can be observed that the high frequency oscillation motion is present in both, the computations and the experiments (Durante [37]), in the form of a regular repetitive pattern in, what appears to be, the steady state response to the single harmonic input motion. Note that experimental work contained only fractional amounts of high harmonics recorded at base (the discussion provides further comments on this aspect), whereas the numerically applied input motion contained solely 10Hz frequency. The pattern of high frequency motion registered experimentally on piles is not fully consistent with respect to those predicted by the simulations. Potential reasons for such discrepancy could lie in the approximately evaluated mass of measuring devices placed on the top of the piles or plastic shear strain accumulation in the experimental work due to the numerous dynamic tests performed on the same soil specimen thus possibly affecting the response of piles embedded in soil. Nevertheless, the fact of observing high frequency motion on the piles in the numerical and experimental studies, even if of different patterns, can be considered as further potential evidence supporting the idea that high frequency motion oscillations can be generated in the dynamic response of nonlinear hysteretic soil and affect the structural response. In other words, the piles with a mass on their top can be thought to act as ‘elastic inclusions’ placed in soil which act as ‘additional measuring instrumentation’, pick up and amplify the high frequency motion generated within the soil mass. Figure 20. Mesh discretization of the example of the soil-structure interaction. 24 Figure 21. The horizontal acceleration records in the large boundary value problem: a) numerical simulation, b) experiment. 4. Discussion The results of Section 3 have presented introductory consideration supporting the hypothetical idea of the phenomenon of the release of soil unloading elastic waves due to soil nonlinear hysteretic behaviour which results in the observation of high frequency oscillation motion in the soil dynamic response. These results, if confirmed in more detailed dedicated research in the future, may have significant importance on various aspects of earthquake geotechnical engineering, including: pointing out gaps in the theory regarding stress wave propagation in materials similar to soils, recognizing the potential origin of high frequency motion measured in small-scale experiments, understanding the source of very high frequency numerical noise in the finite element numerical analyses, showing new insights on seismic soil- structure interaction and, finally, possibly changing the way we think of stress wave propagation in soil in real earthquakes. First of all, this study shows the general importance of the propagation of unloading waves in soils, which is typically not acknowledged in soil dynamics. The importance of unloading waves is well recognized within solid mechanics (e.g. Nowacki [32], Wang [33]). In the light of unknown analytical solution for the wave propagation problems in the nonlinear hysteretic materials, the propagation of unloading waves have been investigated in this work using simplified soil constitutive models in the framework of the finite element method. To the best of the Authors’ knowledge, the propagation of unloading waves have not been studied before in pressure dependent materials such as soils. Moreover, the release and the ‘entrapment’ of soil unloading elastic waves in a soil column modelled with a hysteresis material has also not been shown before in stress wave propagation problems, thus indicates potential gaps in the theory of wave propagation. Note that the numerical studies presented in Section 3.2 have chosen the depth dependent material definition (eq. 2) to be representative of real soil behaviour. Nevertheless, the release of unloading elastic waves can be shown also in a hysteretic depth independent material in a column of a fixed height, thus the presented phenomenon would be applicable to general mechanics, in this work omitted for the sake of compliance with the real soil behaviour and the direct applicability of the obtained results to general experimental observations. The presented hypothetical idea of the release of soil unloading elastic waves may cast new light on the potential source of high frequency motion oscillations registered often in experimental works when a soil specimen is subjected to simplified sinusoidal input motions. Some explanation to the occurrence of the high harmonics in the evaluated spectral responses was pointed out in the past towards soil nonlinearity which results in distortion of a sine wave towards a ‘square wave’ (e.g. Mercado et al. [27]). Nevertheless, apparently experimental works do not register high frequency components represented by the high harmonics of the pattern of an exponential decay (i.e. representative of the distortion towards a ‘square wave’). In the opinion of the Authors, such wave distortion towards a ‘square wave’ would also rather not lead 25 experimentalists to filter experimental measurements. On the other hand, this work has initially indicated a new potential explanation of irregular patterns of high harmonics in the spectral response of experimental measurements. It has been shown how an increased amount of high harmonics close to soil natural frequencies can be observed in the dynamic steady state response of soil, due to the release of soil unloading elastic waves (Section 3.2). Moreover, the numerical studies in this work suggest that high harmonics may not be limited to those below 80Hz as shown in the experiments (Durante [37]). Thus, it would be interesting to filter experimental data with a higher cut-off frequency in the future works, in order to validate the speculative findings presented in this paper. A part of the numerical and experimental comparisons in Section 3.2 is based on the experiments conducted by Dar [17]. These experiments, in general, were often characterized by ‘double-peak’ acceleration records in the top of the soil container when subjected to high amplitude input motions (i.e higher than 0.15g not used in the work by Durante [37]). Originally, the explanation of such ‘double peaks’ was attributed to ‘soil fluidisation’ potentially taking place in soil (Dar [17]), or, possibly, to strain localisation due to shear band development (Gajo & Muir Wood [18]) even though the maximum lateral displacements were reported to be around 2-3 mm (Gajo & Muir Wood [18]), thus likely too low to represent the shear band development. In addition to the peculiarities observed in the horizontal accelerations, some further unrecognized phenomena on soil surface was observed in the vertical displacements measured by Dar [17] and Durante [37] as mentioned by Gajo & Muir Wood [18] and Kowalczyk [44], respectively. These irregularities may possibly be related to the shear localisation (Gajo & Muir Wood [18]) or be representative of the soil-released unloading elastic waves ‘trapped’ in the soil specimen. Clear distinction on whether shear localisation, soil unloading elastic waves, or combination of these phenomena leads to ‘double peaks’ in the horizontal accelerations at soil surface under high amplitude input motions, could be confirmed in future dedicated experimental works. Importantly, regarding experiments in flexible soil containers, some consideration is worth to be given to the recorded input motions in experimental works. Namely, a difference between small-scale tests performed at 1g (e.g. Durante [37]) and in centrifuge (e.g. Madabhushi et al. [51]) is sometimes observed in terms of the recorded input motions. In particular, high harmonics are apparently more prominent at the soil base in centrifuge tests, notwithstanding the intended input motion should consist only of a single harmonic. For this reason, the high frequencies are sometimes considered to be generated by the actuator itself (Brennan et al. [9]; Yao et al. [11]). The possibility that the high harmonics at base are due to the actuators themselves cannot be excluded. This, indeed, might be the case in centrifuges at certain geotechnical research centres (e.g. Kutter et al. [3]). On the other hand, it is worth observing that, the soil specimen in a soil container placed on a shaking table are one complex dynamic system, in which energy induced by the actuators is reflected from the soil specimen back to the shaking table and its actuators. Moreover, the differences between the dynamic impedances of the soil specimen in soil containers on shaking tables are expected to vary in a wide range for various experimental setups, especially in 1g apparatuses as compared to the apparatuses in centrifuge tests (for obvious reasons of payload). Thus, the high frequency motion generated within the soil mass is expected to have different effects on ‘heavy’ 1g experimental setups (as the one used herein for comparisons [17], [37]) with respect to the ‘light’ centrifuge setups, potentially inducing unwanted high harmonics even at the level of the shaking table and controlling systems. This work has revealed how the hypothetical release of soil unloading elastic waves can possibly be a source of high frequency motion observed in simplified experimental setups subjected to sinusoidal input motions. Note that there still might be other possible sources of high frequency motion in soil response to dynamic excitation. Separation of the elastic part of the wave (elastic precursor) from the plastic part of the wave can be such an example (Kowalczyk [44]). Such phenomenon can take place when there is no smooth, continuous change in the soil stiffness between the elastic and the plastic part of the stress-strain behaviour. This might be a case of soil, especially when we recall that the elastic behaviour is associated often with an undisturbed soil fabric, whereas plastic behaviour is associated with changes in the soil fabric. Therefore, intuitively it may appear possible that such initiation of change in soil fabric may result in a non-perfectly- smooth transition between the elastic and plastic stiffnesses, thus allowing the elastic precursor separation and resulting in the ‘entrapment’ of separated soil elastic waves. Although this phenomenon was indicated as potentially possible to occur in soil placed in flexible soil containers (Kowalczyk [44]), this work shows that the release of soil unloading elastic waves would rather be dominant in small-scale experimental studies. On the other hand, in real earthquakes, the separation of the elastic precursor cannot be excluded. Nevertheless, this has not been a subject of further consideration in this paper. The results in this paper are mainly related to simplified input motions used in the numerical and experimental studies. On the other hand, one may want to ask if soil unloading elastic waves can also be released in real earthquakes. Some indication to this point has been shown in Section 3.2 of this paper where scaled earthquake input motions were analysed and revealed generation of high frequency components in the spectral response. Indeed, there can be further consequences of the release of soil unloading elastic waves in strong motion earthquakes where nonlinear soil behaviour is 26 expected to be experienced. For example, the results shown in this paper can contribute to the discussion of double-peak spectra earthquakes (e.g. Gallegos & Saragoni [52]) where the two peaks are attributed to represent the source (i.e. the driving frequency) and the site (i.e. the soil natural frequency) response suggesting that the response of a site bounded at base by a stratum of large stiffness contrast may be heavily affected by soil released elastic waves ‘trapped’ in this stratum. In fact, free vibrations (i.e. soil elastic waves) have been explicitly observed in some earthquakes (Ruiz & Saragoni [53]) including the site of the Mexico City, where a soft clay layer is underlined by much stiffer rock. Regarding the seismic response of the Mexico City, there is a number of unexplained observations in relation to the seismic response recorded in the Mexico City, including ‘anomalous’ records on the rock sites contaminated by unexpected long-period vibrations or unusually long duration of motion in the soft clay layer characterized sometimes apparently by the presence of soil elastic waves and recorded at the stations located on the soft clay (e.g. station Aux as shown recently in [54]). According to some previous research works (e.g. Singh et al. [55]) the ‘anomalous’ records on the rock sites could be representative of the underlying stratigraphy of low-velocity layers, however based on the findings in this work, an alternative explanation of the ‘anomalous’ records could be speculated as those ‘anomalies’ possibly being representative of the overlying stratigraphy, i.e. the soft clay soil layer in the lake basin, where soil elastic wave may get ‘trapped’ as in the case of the simplified numerical studies of the soil container. The presented work has also important implications on the numerical modelling of soil response under dynamic excitation. It has been shown that the standard formulation of the finite element method encounters difficulties in representing the propagation of soil unloading elastic waves for relatively large amplitude motions resulting in the occurrence of strain discontinuity and generation of spurious very high frequency numerical oscillations. The presence of soil unloading elastic waves can explain the numerical difficulties that are typically encountered in computational modelling of seismic geotechnical problems. In any case, the occurrence of strain discontinuity might demand the use of more sophisticated numerical techniques (such as the discontinuous Galerkin method or higher order finite elements, e.g. Semblat & Brioist [56]; or Godunov’s method, e.g. Fellin [34]) than those currently used in practice and employed in this work. It would be presumed that the use of such techniques would allow for further developments of the numerical findings shown in this paper. In addition, this work may point out on the importance of two other aspects in soil numerical modelling. The first one is the appropriate definition of soil elastic stiffness through an appropriate constitutive model. Only constitutive models which define very small strain initial stiffness of soil in a reliable way, can be expected to compute soil response comprising high frequency motion of accurate representation. The second numerical aspect to consider carefully is the amount of viscous damping often used in numerical studies. It has not been shown in this work, however excessively large viscous damping may remove any high frequency motion from the computations including the one representing soil elastic waves. This study has presented initial numerical and some potential experimental evidence supporting the hypothetical idea of the release of soil unloading elastic waves in nonlinear hysteretic soil. Herein, the results have been based on s-wave propagation in dry soil subjected primarily to harmonic input motions. In the light of a substantial number of recent experimental works dealing with saturated soil and liquefaction cases, the authors of this work were somehow ‘lucky’ to have access to examples of the experimental research carried out in dry sand, for which the numerical studies were developed, and where the steady state solution under harmonic motion can apparently be reached (i.e. the plastic deformation and the change in soil stiffness characteristics after each loading cycle is negligible). Nevertheless, the presence of soil elastic waves would equally be expected in the response of saturated soil. However, in that case, the steady state response under harmonic excitation is typically not reached due to the changes in the stiffness of saturated soil as a result of pore pressure generation, thus the presence of soil unloading elastic waves is expected to be much more difficult to be recognized. Finally, it has to be recalled that, although the numerical studies have rather convincingly identified the presence of soil unloading elastic waves, the experimental examples used here for the purpose of comparisons with the numerical studies are more speculative. In fact, these experimental works were focused on other aspects of soil dynamic response (i.e. the response of piles in Durante [37] or the development of a soil container in Dar [17]), and, moreover, could also be influenced by other phenomena resulting in observation of high frequency motion, e.g. wave reflections due to bi- layered soil profile or wave scattering from piles in the work of Durante [37], or shear band development in the work of Dar [17]. Therefore, it appears that there is need for further investigation of the presented idea, including experimental studies dedicated towards clear confirmation of the potential existence of the phenomenon of soil unloading elastic waves in the dynamic response of soil. 27 5. Conclusions This paper presented initial consideration supporting the hypothetical idea of the release of soil elastic waves upon unloading in nonlinear hysteretic soil (i.e. soil unloading elastic waves), mainly in the steady state response when subjected to harmonic excitation, but also possibly in the response to earthquake-like input motions. The potential existence of soil unloading elastic waves in the steady state response was shown, primarily, in the finite element numerical simulations, and secondarily, by comparisons with past experimental works on shaking tables. The paper showed that soil unloading elastic waves can be a potential explanation to high frequency motion registered in experimental works in the steady state solution and can potentially be important when analysing real earthquake records and the structural response. Acknowledgments The research presented herein has been initiated as a part of the XP-Resilience project. This project received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 721816. The authors would like to thank Professor Colin Taylor (University of Bristol) for his insightful remarks regarding the complexity of physical modelling of soil behaviour in small-scale experiments. The authors would like to thank to Dr Maria Giovanna Durante (University of Texas) and Professor Luigi Di Sarno (University of Liverpool) for providing examples of raw experimental data of the small-scale experimental work on piles. Appendix A Severn-Trent constitutive model- calibration and performance This appendix presents the input parameters (Table A1) and the chosen calibration (Fig. A1) for the Severn-Trent constitutive model together with an example of its response when simulating a cyclic simple shear test on Toyoura Sand compared with experimental data (Shahnazari & Towhata [57]) shown on Fig. A2. Figure A1 Calibration of Severn-Trent model compared with guidelines by Dietz & Muir Wood [48], Seed & Idriss limits [49] and laboratory data by Kokusho [50]. 28 Table A1 Input parameters for the Severn-Trent sand model Parameter Description νΔ Intercept for critical-state line in ν-ln p plane at p=1Pa 2.194 Δ Slope of critical-state line in ν-ln p plane 0.0267 φcv Critical-state angle of friction 33° m Parameter controlling deviatoric section of yield surface 0.8 k Link between changes in state parameter and current size of yield surface 3.5 A Multiplier in the flow rule 0.75 kd State parameter contribution in flow rule 1.3 Bmin Parameter controlling hyperbolic stiffness relationship 0.0005 Bmax Parameter controlling hyperbolic stiffness relationship 0.002 α Exponent controlling hyperbolic stiffness relationship 1.6 RR Size of the yield surface with respect to the strength surface 0.02 ER Fraction of G 0 used in the computations 1.0 Figure A2. Typical model responses for the simulation of a drained cyclic simple shear test on loose sand, corresponding to the calibration of the constitutive parameters in Table A.1 compared with laboratory results: a) numerical simulation, b) experiments by Shahnazari & Towhata [57]. 29 Appendix B Ormsby wavelet identification This appendix presents the numerical dynamic identification of an 0.8m high soil column modelled with a simplified constitutive model (eq. 2) used in the numerical simulations shown in Section 3.2 of this paper. The chosen input motion is an Ormsby wavelet (Fig. B1a) of a flat spectral response between 10Hz and 110Hz (Fig. B1b) and the maximum amplitude of 0.0001g. The soil column subjected to such input motion responses at its top as shown on Fig. B1c. The evaluated natural frequencies of the soil column are around 33.0Hz and 88.5Hz (Fig. B1d). Fig. B1 Dynamic identification of an 0.8m soil column be means of an Ormsby wavelet. 30