Exact seismic response of smooth rigid retaining walls resting on stiff soil

The assessment of forces exerted on walls by the backfill is a recurrent problem in geotechnical engineering, owing to its relevance for both retaining systems and underground structures. In particular, the work by Arias and colleagues, and later also the one by Veletsos and Younan, among others, becomes pertinent when considering pressure increments on underground structures triggered by seismic events. As a first step, they studied the response of a rigid retaining wall resting on rigid bedrock subjected to SV waves, introducing some simplifying assumptions. This paper presents the exact solution to this reference problem. The solution is given in horizontal wavenumber domain; hence, it comes in terms of inverse Fourier transforms, which can be approximated numerically in Mathematica, which in turn are verified against finite‐element simulations. Specific features of this exact solution that were not captured by prior engineering approximations are highlighted and discussed.


INTRODUCTION
The study of the forces that soil exerts over retaining structures dates back centuries. A first clear distinction can be made on the basis of the methods that have been employed for these studies: • Methods based on plastic behavior of soils, foremost, limit-state theory, 1 for so-called yielding walls, that is, walls that move enough so as to elicit plastic response in the soil. • Methods formulated on the basis of the theory of linear elasticity to describe both soil and structure, 2 which presupposes small deformation in both soil and wall. The wall in these cases is termed unyielding.
The proper approach for the analysis does depend on the typology of the wall that is to be assessed or designed. Static load scenarios can be designed following classic methods 3 ; on the contrary, when it comes to considerations in the dynamic earthquake setting, the right choice remains unclear to this day.
On one hand, yielding walls require consideration of the soil plastic behavior in the time domain, and thus, researchers have proposed pseudo-static methods, which, in spite of intrinsic limitations, deliver satisfactory, even conservative, outcomes for free-standing gravity or cantilever walls (except if resting on stiff rock). On the other hand, approaches based on linear elasticity are suitable to ascertain seismic pressures on unyielding walls (as, e.g., restrained basement walls or rigid-enough underground structures), as long as the structure whose walls are being analyzed does not experience other soil-structure interaction (SSI) effects (rocking). The National Earthquake Hazard Reduction Program (NEHRP) recommendations 4 sanction this classification, while also suggesting that a wall whose top deflects around 0.2% of its height already triggers plastic behavior in its soil vicinity, ergo the threshold to consider a wall as yielding is in practice fairly narrow.
The influence of soil-structure interaction, the effect of water in the backfill, the role of soil cohesion on the seismic response, and other topics are still subjects of ongoing study.

Approaches based on elasticity
When straining of both the soil and wall remains relatively small, linear elasticity theory provides a proper framework to describe its evolution. The only completely rigorous solution (from the mathematical standpoint) using this approach was obtained by Wood. 5 Wood considered a 2D finite stratum of soil on rigid bedrock confined between walls deforming under linear-elastic plane-strain conditions. Although this configuration may resemble a different system at first, logic dictates that if the distance between walls increases, the effects of one wall over the other one must fade, leaving the system virtually equivalent to the one-wall configuration insofar the displacements field around the wall are concerned. Hence, the problem we are to discuss may well be dubbed the singular Wood problem.
The solution to this twin-wall problem has also been extended recently to the case of inhomogeneous soil. 6 Wood's solution was also extended by Papazafeiropoulos and Psarropoulos,7 who introduced the concept of "distress spectrum." Approximate models for one-wall systems have been developed throughout the years in order to circumvent excessive mathematical complexity, but every effort of tackling this very problem relied on introducing some form of simplification.
The first attempt to consider a one-wall system directly traces back to the work by Matsuo and Ohara. 8 They assumed a confined solution 9 in which no vertical displacement develops anywhere in the soil domain, not even at regions close to the wall. This solution was shown to be unbounded in the incompressible-material limit (undrained soil conditions). Arias et al 10 developed an approximate solution wherein the stresses in the soil were modeled using a spring-like behavior assumption.
Two classic solutions for the one-wall case were derived by Veletsos and Younan. 11,12 In the first solution, a variation of the confined solution was utilized: instead of neglecting a component of the displacement field, these researchers assumed that no normal vertical stress develops anywhere in the soil ( yy = 0 everywhere) in order to simplify and solve the equations of motion, thus retrieving the displacement field over the whole domain. Regrettably, satisfying vertical equilibrium and the boundary condition xy | w = 0 became impossible due to the introduction of this, among others, severe simplification. In the second solution, 11 a substantial improvement on Scott's model 13 was presented. It was assumed that the soil stratum behaves as an elastically supported, semi-infinite horizontal viscoelastic bar with distributed mass and that the horizontal gradient of vertical displacements in the calculation of shear stresses is negligible. A solution was thus derived for the medium impedance, which related the stress field on the wall to the simpler displacement field in the far-field (referred as "1D soil column" or "shear beam"). Later, 14 their analysis was refined to consider different wall typologies, yet the foregoing simplifications were kept. Nevertheless, this work became the first one in attempting to include the effect on wall flexibility on the magnitude and distribution of pressures, and others followed the path it opened. 15 Other models have relied on different approaches, 16 but current one-walled state-of-the-art models, as those proposed by Kloukinas et al, 17 later used in the already-mentioned "kinematic framework" proposed by Brandenberg et al, 18 have been directly influenced by the Veletsos-Younan landmark work 11 and have inherited the same pivotal assumption to simplify the equations of motion: normal vertical stresses are forced to be zero everywhere, yy = 0.
A more detailed summary of these and other previous studies can be found in a recent paper. 19 Despite simplifications, the aforementioned models have proved their reliability against both numerical simulations 20 and analytical solutions of the two-walled system with large spacing between them. 5 In this paper, we provide the exact solution of the problem, computed with no either physical or mathematical simplifications. The solution procedure is inspired by those recurrently furnished on the basis of the Elastic Waveguide Theory 21 and does not resort to any unconventional mathematics. Results are given in the form of algebraic expressions in the horizontal wavenumber domain and compared against results of finite-element simulations in the spatial domain. We also discuss features of the exact solution that we have identified were absent from simplified one-wall solutions and ponder in what cases these features become relevant.
Finally, we must acknowledge the inherent limitations of this simplified system: most retaining walls do not rest directly on rigid bedrock but on compliant soil. However, removing such simplification would prevent reaching an analytical solution for the whole displacement field. This simplified model remains useful because it allows to obtain insight on how earthquakes influence pressures on retaining walls and how different contact conditions between wall and soil can affect these forces.

DERIVATION OF THE EXACT SOLUTION IN HORIZONTAL WAVENUMBER DOMAIN
The system shown in Figure 1 comprises a rigid wall connected to a rigid base (bedrock) and restraining a layer of soil. No specific boundary conditions, other than rigidity, were used for the wall in Veletsos and Younan's original work, as the simplifying constraint they introduced, yy = 0, precluded those considerations anyway. Following Wood, 5 we shall consider that the wall is smooth, that is, frictionless, and thus, no shear stress develops at the soil-wall interface.

Presentation: Total displacements
Consider an semi-infinite soil layer, resting on rigid bedrock and bounded to the left by a rigid (thus, it moves with the bedrock base) and smooth (thus, no shear stress develops at the interface) wall.
The excitation of the system is an imposed time-varying displacement at the base g (t). Assuming that the soil behaves as a homogeneous, density , isotropic linear-elastic solid (characterized by Lamé parameters , shear modulus, and ), undergoing small deformations, the equations that govern the dynamic response of the material are the so-called Cauchy-Navier equations 9 : subjected to the following boundary conditions: whereas at y = H and at x = 0 Including hysteretic damping in the material response can be achieved (once we move to the frequency domain) by simple substitution of the real shear modulus by (1 + i d ), where d can be referred as "damping ratio." We also assume the Poisson ratio to be independent of damping; hence, damping also affects the second elastic constant:

Introduction of relative displacements
Both u t and v t represent total displacements (horizontal and vertical, respectively), measured with respect to some fixed reference. This particular problem lends itself to easier treatment in terms of relative displacements. Thus, the change of variable u t (x, , t) = u(x, , t)+ g (t), v t (x, y, t) = v(x, y, t) was introduced in Equation (1) to obtain the equations of motion in terms of relative displacements: As for Equation (2a), whereas Equation (2b) does not change insofar the stresses do not depend on the displacements themselves but on their spatial gradients: and similarly, Equation (2b) becomes Note that the original problem can be classified as "unforced wave propagation with inhomogeneous mixed boundary conditions," whereas this modified one corresponds to "forced wave propagation (because now there is an external body force d 2 g ∕dt 2 in Equation 3a) with homogeneous mixed boundary conditions." Purposefully, nothing has been said about initial conditions because we shall focus on steady-state response to a harmonic load. Considering a transient regime governed by the initial conditions can be achieved by appealing to Laplace's transform in time. This could be easily done, but it would entail extra muddling of the problem, as one would have to invert yet another integral transform. Moreover, such endeavor is unnecessary if we choose to focus on obtaining transfer functions linking the external load to the response.

Dealing with the time domain: Assume harmonic loading and steady-state response
The disclosed objective of this work is, chiefly, finding the transfer function for the earth thrust on the wall. Therefore, a more convenient, equivalent way of proceeding is to assume a harmonic decomposition of both the load and the response and, invoking superposition, focus onto one sole harmonic. In summary, this means assuming u = ûe i t , v =ve i t as well as d 2 ∕dt 2 =Ẍ g e i t , where can be any excitation frequency and û,v must be understood as complex numbers, whose modulus represents the magnitude of each displacement and its argument represents the phase lag between load and response. Introducing these changes into Equation (3) yields with boundary conditions (from Equations 4a to 4c):

Nondimensionalization and key symmetry argument
Working with dimensionless equations facilitates keeping track of the parameters of the problem, and also eases the burden of algebraic manipulations. Regardless of its convenience, this step is optional and does not bear any special significance within the solving procedure. For the purpose of writing the last equations in dimensionless form, the following dimensionless variables and parameters are introduced: where represents the dimensionless horizontal coordinate, the dimensionless vertical coordinate,ũ andṽ the dimensionless horizontal and vertical displacements, respectively, r the dimensionless excitation frequency, and c = c( ) = c s ∕c p = √ 2(1 − )(1 − 2 ) the ratio between P-wave and S-wave propagation velocities, which, in the case of isotropic linear-elastic solid, is a function of the Poisson ratio, , solely. Once again, using this change of variables in Equation (5), the equations of motion become At this moment, one must consider how to deal with the unboundedness of the horizontal coordinate. Because x ∈ [0, +∞), an approach based on applying Laplace's transform on displacements may seem on point at first; however, opting for this approach requires knowledge of both the value of the displacement and its first derivative at the wall, what is not the case: the value ofũ is known yetũ∕ is not, and conversely,ṽ∕ is known whereasṽ is not. Nevertheless, the conditions seem ideal to apply Fourier sine transform and Fourier cosine transform, yet this option appears confusing as it would require using different transforms for each displacement field component. The fitting approach for this specific problem is appealing to symmetry and then applying the standard Fourier transform. The symmetry argument reads: this type of one-wall system formulated in terms of relative displacements, amounts to applying a symmetry condition on an infinite layer wherein the body force suddenly changes sign, but not magnitude, at = 0, see Figure 2.
Thus, in lieu of the semi-infinite soil strip domain, let us consider an infinite domain ∈ (−∞, +∞). At = 0, there is a discontinuity in the external body force: originally, it was set to −1 and defined over x ≥ 0; now, it follows a step function, 1 − 2 ( ), defined over ∈ R ( ( ) is the Heaviside function centered at = 0). Hence, the system that will finally be solved is governed by the following equations and boundary conditions (equivalent to Equation 8 with boundary conditions of Equation 9): subject to the boundary conditions at = 0: We resort to the Fourier transform 22 at this point: where k is the dimensionless horizontal wavenumber (its relation to the physical wavenumber k being kH = k); hence, Equation (10) becomes likewise, Equation (11) becomes the following: Equation (14) represents the problem in the original vertical (nondimensional) variable and in the wavenumber (k) domain, assuming harmonic loading and steady-state conditions.

Exact solution in the horizontal wavenumber space
The system of Equations (13) can be condensed even more by introducing additional auxiliary variables  = V∕ , subject to the boundary conditions equivalent to Equation (14): Introducing a vector of unknowns, X = , Equation (15) can be written simply as where The eigenvalues of this matrix are where each solution corresponds to either S waves (first two solutions) and P waves (last two solutions) propagating within the bulk of the stratum. Because we do not pursue the inversion of the solution in this paper, we shall not delve in the intricacies of properly choosing branches of these multivalued functions. 23 The eigenvectors (normalized to 1 in the fourth entry) are thus, the solution must be or introducing the following shorthands = where A, B, C, D are "constants" (which nevertheless depend on the parameter k) to be determined from the boundary conditions, starting from the easiest ones, V(0) = 0 and U(0) = 0 followed by the most convoluted conditions at = 1, xy = 0 and yy = 0 Recast Equations (23), (24), (25c), and (26b) in matrix form the value of the constants is found after solving (inverting) this linear system. The actual expressions are computed by recourse to Mathematica. 24 The expressions are frankly convoluted, so denominator and numerator are shown separately. For instance, the coefficient A is expressed as A = A N ∕A D . Therefore, and the corresponding denominators There is no tabulated formula to invert the integrals where the coefficients appear. Analytical inversions must be performed by recourse to contour integration and residue calculus, 23 which turns out to be an extensive and challenging task that the authors have not succeeded in yet.

EVALUATING THE SOLUTION IN MATHEMATICA
Instead of taking that path, one can approximate the integral numerically, as the coefficients can be evaluated in Mathematica with relative ease.
We use the built-in function "NIntegrate" 25 with default parameters, which, based on the actual form of the integrand, automatically selects the most efficient discretization of the integration interval and the most advantageous quadrature rule.
When approximating these integrals numerically, the integration has to be stopped at some point; it is known (from previous study carried out by the first author in the context of his doctoral thesis 26 ) that the horizontal wavelengths ∼ H are those that control the response close to the wall and are also the shortest wavelengths that are present. Therefore, one can limit the integration to |k| ∈ [0, 10]. See Figure 3 for an example of the aforementioned coefficients, while Figure 4 displays en example of the numerical integration. Finally, Figure 5 shows how the implementacion produces one of the graphics in the paper.
The notebook used by the authors to generate the data in the following figures is available in GitHub (see Supplementary Material section), so that any reader can download it, reproduce the results or generate new ones if so desires.
Additionally, a Matlab script that generates the figures in the text can be found in the same GitHub repository.

VERIFYING THE SOLUTION
Once these coefficients have been computed, the solution in the wavenumber domain has been found. Results derived from the exact solution are compared to FEM analysis performed on Abaqus/Standard. 27 This package relies on "infinite elements" to model semi-infinite domains, yet this type of elements are known to cause issues in plane strain problems. 28 For this reason, as the unboundedness in the horizontal direction cannot be accounted formally, a very slender model of length L (L∕H = 275∕2) was used. This approach mirrors the one advocated recently by Durante et al, 29 which delivered satisfactory results in the quasi-static regime, yet it had not been tested in frequency domain dynamic analysis previously.  The parameter values are encapsulated in the following table (neither Poisson's ratio nor damping factor are shown as they will be varied across simulations, their values will be pointed out explicitly when displaying results). Regarding units, these are not a concern because we are dealing with a linear problem and all the results will be expressed in dimensionless form.
Relevant information concerning the model: • Mesh: the mesh is refined around the wall, the target element size being H∕100 whereas the mesh coarsens in the far-field to a element size H∕4. The refinement is maintained up to a distance H∕2. The transition zone between wall mesh and far-field mesh is stretched horizontally for a distance 4.5H. Element type used: "CPEG8R" (quadratic-interpolation plane-strain elements 8-node quadrilateral, reduced integration). Quadratic elements are deemed more efficient as they capture the far-field response properly with just four elements across the layer width. • Loading: a load type "Gravity" is used to model the body forces elicited by the shaking. The intensity is 0.73, and the acceleration is directed towards the wall. • Damping: the damping type is "Structural," which works exactly as in the analytical derivation when we substitute by (1 + i d ); that is, it introduces a damping force phased ∕2 radians with respect to the elastic forces, and whose magnitude is scaled by the factor d .
The details on the analyses carried out in Abaqus follow: • First step: "Frequency" perturbation analysis uses the Lanczos method to extract the vibration modes in the interval [0, 350] cycles per time unit (this range is chosen so as to capture the first modes of the corresponding shear beam given the parameter values in Table 1). The frequency domain is uniformly subdivided using 500 points. • Second step: "Steady-state Dynamics -Modal" analysis. Using the newly obtained information on the natural modes of the system, the transfer function, in that frequency range, for the any response is obtained.
The static analysis used to generate the FEM data in Figure 11 is carried out with the same mesh and parameters, but the analysis type is "Static, General" instead.

Vertical displacement at upper-left corner
One can then move to calculate, for instance, the vertical displacement at the top of the wall,ṽ( = 1, = 0) =ṽ top (see Figure 6), by inverting the transform: where A, B, C, D, , and are functions of the dimensionless wavenumber k and the integration is limited to k ∈ [−10, 10]. The comparison of Equation (30b), as evaluated in Mathematica, to the results obtained in Abaqus, once properly expressed in nondimensional form, is displayed in Figures 7 and 8. The characteristic displacement U =Ẍ g H 2 ∕ is the same that was used to write displacements in dimensionless form in Equation (7)

Earth thrust over the wall
The earth thrust represents the integral of the stresses that develop at the soil-wall interface. In the case of the smooth wall, harmonic loading and steady-state conditions, it can be written as (the subindex w indicates evaluation at the wall) ergo the amplitude of the thrust can be expressed in terms of the solution coefficients in the horizontal wavenumber domain as the inverse Fourier transform below: (32) Likewise, the pressures themselves can also be expressed in dimensionless fashion in terms of the coefficients: Similarly to Equation (30b), Equation (32) can be evaluated numerically to produce Figures 9 and 10. For = 1∕3 and d = 0.01, see Figure 9. The plot includes numerical results from Equation (32) in Mathematica and two previous research results: the thrust provided by Younan and Veletos 11 (evaluated using the five first modes), based on the assumption yy = 0, and its direct state-of-the-art heir by Kloukinas et al, 17 which considers only the fundamental mode. There seems to be some numerical artifact the authors have not yet identified, which makes the numerical result oscillate mildly between the first two spikes. Other than that, the agreement is practically global (there is a meager mismatch of amplitudes at the first and third resonant spikes). All the features predicted by the exact solution are present in the FEM analysis, including a substantive "dip" (deamplification) right before the second peak. As reported in Figure 7 and figures thereafter, this resonance happens at ≈ 2 s = c s = p , because c p ∕c s = 2 as = 1∕3. This secondary resonance would have been overlooked had we assumed the thrust resonates at the same natural frequencies as the 1D soil column in the far-field. 11,30 The effect of damping on the thrust is addressed in Figure 10. The figures correspond to = 0.1 (so one should expect, again, to find resonance at ∕ s = 1, 1.  Again, there is an excellent compliance of numerical to analytical results, in the three cases. It is instructive to note that in this case, = 0.1, the dent at the "intermediate resonance" seems not to point up but down and that increasing damping smooths out the two second resonant peaks and the aforementioned dip that happens before the resonance associated with p .
The physical meaning of this "intermediate" resonance has not been spelled out yet: it remains to demonstrate if it is an artifact of mode conversion or of a surface wave propagating along the wall-soil interface, and how pernicious its effects can be over the wall integrity. Something is sure: it had been unaccounted for in the previous engineering solutions, which only displayed resonance at the wall at the same natural frequencies as the far-field.

In-depth analysis of intermediate resonance
This section is aimed to explain the causes of the extra set of resonance frequencies. Let us begin by noting that the normal horizontal stress developing at the soil-wall interface, and the thrust itself by extension, can be decomposed in two parts. From Equation (31), introducing the linear-elastic constitutive law, and writing it in dimensionless form as in Equation (32): The first addend corresponds to the stress induced by gradients of horizontal displacement, viz. the straining of the soil bordering the wall in the horizontal direction. The second corresponds to the gradients of vertical displacement over the wall, namely, the straining along the vertical direction. See how in previous engineering solutions 11,17 both components were linked through the simplifying assumption yy = 0, which forces them to have the same phase in all the cases, regardless of excitation frequency or Poisson's ratio of the soil, and also linked their amplitude through the so-called compressibility factor. The exact solution renders a more complex picture, wherein the two components are independent and behave in significantly different fashion. Let us showcase this.
We already haveṽ top expressed in terms of an inverse Fourier transform of the coefficients, Equation (30b); similarly, the integral over the wall of the horizontal gradient of the horizontal displacement can be expressed as Recall that an inertial force oriented towards the wall (e.g., towards the negative abscissa) is just a proxy for the imposed movement along positive abscissa that appeared when making the change of variable from total to relative displacements. Hence, the compressive force induces always a relative displacement towards the wall; thus, the baseline displacement is negative and its phase corresponds to zero at low frequency because this negative displacement is in phase with a negative force. Nevertheless, the direction of the vertical displacement does depend on the compressibility of the soil: in the case = 0.1, the soil slides down compressing the soil at the interface vertically, whereas in the case = 1∕3, the soil slides up, stretching it. The threshold value from one behavior to the other (the one corresponding to no vertical displacement whatsoever) seems to be = 0.15. This is easy to visualize attending to the low-frequency response (quasi-static regime, corresponding to ∕ s ≪ 1) as depicted in Figure 11: We posit that each component plays a different role in furnishing the total thrust and that the phase lag between the two is a critical factor. In order to check this assertion, the phase of the two addends in Equation (34b), Equations (30b) and (35), is compared next using the analytical solution (see Figure 12).
We acknowledge that the vertical displacement at the top corner of the layer is in phase with the horizontal gradient of horizontal displacement when = 0.1; hence, once the second resonance frequency is excited, the phase lag between the two is so big that the contribution of v top , which clearly increases more due to this resonance, diminishes the total amplitude of the two combined addends. Conversely, in the case = 1∕3, the vertical displacement is positive when the horizontal displacement gradient is negative; hence, the former is delayed radians, half a period, with respect to the latter; nevertheless, when the second resonance frequency is excited, both components contribute positively to the response as both their phases have the same sign (see Figure 13 for the final effect on the thrust itself).
In conclusion, the appearance of the new resonant spikes and their different amplitude, which depends on the compressibility of the soil (the Poisson ratio), is elicited by the different behavior of the two components of the thrust on the

Simultaneous resonance
It has already been discussed that there are two contributions to the thrust, each with different traits: • One component that depends on the integral over the wall of the horizontal gradient of horizontal displacement, Equation (35), which resonates at the same frequencies as the far-field, odd multiples of s = c s ∕2∕H. • A second component proportional to the vertical displacement at the top of the layer, v top , which resonates at the same frequencies as the far-field (just like the prior term) and at a second set of frequencies corresponding to odd multiples of p = c p ∕2∕H = c( )c s ∕2∕H = c( ) s . At each frequency, the lag between the second and first components depend on the compressibility of the soil.
Previous engineering results 11,17,30 only captured the first set of resonances. Assume that we are to consider a set of different soils having the same density, , and shear stiffness, , but different compressibility, , the location of the first of set of natural frequencies, those that coincide with the "shear beam" resonance modes that are proportional to s , is the same in every instance, but the location of the second set of frequencies will depend on the value of Poisson's ration through the factor c( ) = c p ∕c s = √ 2(1 − )∕(1 − 2 ). See that if = 7∕16 = 0.4375, then c = 3, whence p = 3 s : the first resonance mode of the second set coincides with the second mode of the first set. This happens above the aforesaid threshold ≈ 0.15; thus, based on the discussion in the previous section, the two contributions of the thrust reinforce each other.  Figure 14 displays the unconservative limitations of prior approximate solutions: even though the amplitude of the first peak is even slightly overestimated, the second resonance, which in this case includes also the first peak corresponding to the second set of modes, is overlooked.
This simultaneous resonances are not restricted to a narrow band of values of around 0.4375, but it happens for other values. Figure 15 displays a comparison between 0.42 , 0.4375 , 0.45. In all three cases, remarkable increments with respect to the values delivered by the engineering approximations are found in the region surrounding ≈ 3 s .

Thrust eccentricity
The resulting moment at the base of the wall produced by the seismic pressures,M, is another important design variable. See that, given the following definition of this variable, one can obtain the second form through integration by parts: where y ′ represents a "dummy" integration variable; in dimensionless form: the auxiliary variable q( ) represents the dimensionless partial integral of the stresses (from zero up to height ), where ′ = y ′ ∕H and obviouslyQ∕Ẍ g H 2 = q(1). Having both the expressions for the thrust and the moment, the thrust eccentricity is defined as and so it is also calculated in Mathematica (see Supplementary Material). Clearly, the last expression will be a function of the dimensionless frequency, r, of the damping factor, d , and also of the Poisson ratio, the latter being a feature that was not present in either Younan and Veletsos' results 11 (insofar we assume and to be the two independent elastic constants), which does depend on r and d , or in Kloukinas et al, 17 which remains constant independently of any of the aforesaid parameters. The following figures explore the influence of these dimensionless parameters on the eccentricity, and include comparisons to previous results. We acknowledge from Figure 16, on one hand, that the eccentricity proposed by Younan and Veletsos resembles the exact one for both below-the-threshold value of Poisson's ratio ( = 0.1) and above-the-threshold values ( = 1∕3, = 7∕16), yet their model misses significant oscillations in the intermediate region of the plot in the above-threshold case. Moreover, we also note that for the "simultaneous resonance" discussed in Section 4.2.2, = 7∕16, the exact eccentricity matches well the one of the simplified model at the frequency where the phenomenon occurs ( ≈ 3 s ). Finally, let us highlight how the value 2∕ ≈ 0.64 can be used as an approximate eccentricity for practical purposes. 17

Stress distribution on the wall
Not only the stress resultant and the eccentricity are of interest, but also the actual distribution of stresses across the height of the wall deserves attention. This distribution has been extensively studied in previous research, 6,7,11,12,16 because it represents a parameter that can usually is compared to experimental centrifuge data. 31 Figure 17 shows the distribution of stress magnitude for three values of Poisson's ratio, = 1∕4, 1∕3 , 7∕16 at six frequencies, from s ∕2 to 3 s spaced uniformly by s ∕2 increments. Damping is d = 0.10 in all plots.
Features that worth highlighting in Figure 17: • In tile (A), the distribution resembles the one obtained using the quasi-static approximation, and the expected trend (larger stresses as Poisson's ratio increases) can be acknowledged. • Tile (B) displays the stress distributions at the first resonance peak, which features substantially larger stresses (see that the values at the top are about five times the characteristic valueẌ g H whereas the previous quasi-static results barely exceeded 1.5 times that value). • Tile (F) depicts the stress distribution that a load exciting the second resonance frequency of the far-field shear-beam would elicit. See that for both = 1∕4 and = 1∕3 the magnitude remains relatively small but for = 7∕16 (the value that makes the first frequency of the second set coincide with the second one of the first set) the stresses around the base increase substantially. The latter trait agrees with the observations in Figure 14, which shows a substantial thrust increment around that frequency, and to Figure 16B, where the eccentricity for = 7∕16 corresponding to this frequency is smaller than the one for lower frequencies, what is consistent with large stresses located closer to the base of the wall.

CONCLUSIONS AND FUTURE WORK
In this text, we derive the steady-state harmonic solution of the problem consisting of a long (so 2D plane-strain assumptions apply) rigid smooth wall, which rests directly on rigid bedrock, subjected to SV waves. The solution can be expressed in terms of elementary functions in the horizontal wavenunber domain. By evaluating the solution numerically, it has been shown that there exist resonance phenomena that current simplified models do not capture, as the latter presuppose that the natural frequencies of the displacement field close to the wall resemble those of the displacement at the far-field. It has been also proven that this new resonance can lead to large forces that were previously unaccounted for in simplified models, in particular when ≈ 7∕16, and that these forces can provoke large moments at the base because the eccentricity does not decrease dramatically. Comparing the exact solution to the one-wall state-of-the-art model 17 also yields as conclusion the latter being safe as long as ≤ 0.15.
To formally close the problem, the inverse Fourier transform involving the coefficients in Equation (15) should be computed to obtain the exact displacement field in every point of the domain.
This exact solution can potentially be generalized to the case of rough wall (the soil being bonded to the wall), by mildly adapting the boundary conditions at the wall, and applying the same procedure outlined in this text; this particular case possess great interest as removing the capacity of the soil bordering the wall to deform vertically should make the second set of natural frequencies vanish. An intermediate case between the two in which the wall is neither smooth nor rough-but is governed by a certain friction law-may be also amenable utilizing a slightly different approach. Another interesting extension is the flexible wall case, which, our intuition suggests, could be considered in conjunction with the symmetry argument if the wall deflection was averaged on a per-cycle basis (hence, this averaged relative displacement would be zero all over the wall). This may provide new insight on the relation between wall stiffness and seismic earth pressures so as to confirm the insight obtained by Veletsos and Younan through simplified models and simulations. 14