The technique referred as Geometrical Optics entails considering the wave propagation in a heteroge- neous medium as if it happened with inﬁnitely small wavelength. This classic simpliﬁcation allows to obtain useful approximate analytical results in cases where complete description of the waveform behavior is virtually unattainable, hence its wide use in Physics.This approximation is also commonly termed Ray Theory, and it has already been thoroughly applied in Seismology. This text presents an application of Geometrical Optics to 1D Site Response (1DSR): it is used herein to, ﬁrst, explain and elucidate the generality of some previous observations and results; second, to partially settle an open question in 1DSR, namely “what are the equivalent homogeneous properties that yield the same response, in terms of natural frequencies and resonance amplitude, for a certain inhomogeneous site?”, provided few assumptions. Keywords Geotechnical Engineering Site Response Analysis Geometrical Optics Ray Theory 1 Introduction Site conditions and soil layering have paramount inﬂuence on seismic ground motions [ Benuska, 1990 ]. The evolution of mechanical properties in the soil can induce a concentration of seismic energy in the softer layers [ Towhata, 1996 ], which usually are those closer to the surface, and in turn translates into higher risk for man-made above-ground and underground structures. A number of scholars have been occupied with studying the modiﬁcation of ground seismic response due to the presence of inhomogeneous stratiﬁcation in the soil by means of simpliﬁed one dimensional (1D) wave propagation models under the assumption of linear-elastic medium [ Dobry et al., 1971 ,Gazetas, 1982 ,Vrettos, 2013 ,Vrettos, 1990 ,Zhao, 1997 , Afra and Pecker, 2002, Travasarou and Gazetas, 2004, Semblat and Pecker, 2009], which happens to be suitable if 1)the material response is linear-elastic (there is an absence of either non-linear material response or large deformations), 2)the soil is horizontally layered in depth (thus the problem reduces to one-dimensional, vertical, wave propaga- tion), 3) SH waves moving vertically conform the excitation, 4) material damping is frequency-independent. The customary approach focuses on ﬁnding steady-state solution in frequency domain, which happens to be the dominant component of the response, except under some speciﬁc circumstances [ Sarma, 1994 ], and it is also known to provide a good approximation to Ratios of Response Spectra (RRS) [ Dobry et al., 1992 ,Joyner et al., 1992 ]. This APREPRINT - FEBRUARY 17, 2020 simple approach has been favored in seismic building codes [ Urzúa et al., 2017 ], even for sites displaying strong non-linear response [Ordaz and Meli, 2004], as a ﬁrst step towards a more realistic assessment. Although the problem can be solved analytically for some types of cs(y), a preferred approach is to deﬁne an equivalent uniform shear-wave (also termed “effective velocity") but how to deﬁne this parameter has remained an open debate. In particular, the usage of the harmonic average of the velocity distribution has been an option pursued since the 70s [ Madera, 1970 ]. Despite the fact that this average overestimates (underestimates) the fundamental frequency of the site in cases having surface-to-bedrock decreasing (increasing) stiffness [ Urzúa et al., 2017 ], it has been known [Zhao, 1997 ] to yield excellent estimates of the higher mode natural frequencies in sites having continuous proﬁle given by the so-called “generalized parabola" [ Rovithis et al., 2011 ]. A mathematical explanation of this fact was provided through recourse to the asymptotic behavior (for large frequencies) of the exact solution by the authors of this manuscript: it was found that in the high-frequency limit (corresponding to the resonance mode as high as desired) the exact equivalent velocity was indeed the harmonic mean of the proﬁle. This text seeks to ﬁt the aforementioned result within the general theory of ray approximation for wave propagation inhomogeneous media [ Kravtsov and Orlov, 1990 ] and to generalize the previous result, delivering a clear answer to the question “when should one use the harmonic mean of the shear wave velocity distribution as equivalent homogeneous velocity to calculate resonance modes in 1DSR?”, provided the proﬁle fulﬁlls two straightforward conditions. Answering the prior question can have important consequences, since equivalent velocity in terms of harmonic average is currently a widely-used parameter when it comes to account for soil heterogeneity, both in Site Response Analisys proper (see for instance the deﬁnition of Vs30in ASCE 7/16 [ ASCE, 2016 ] or the discussion in [ Urzua, 1974 ]) and Soil-Structure Interaction [NIST, 2012]. 1.1 Governing equation cs,top cs,basecsy u = 0𝑦τxz = 0 HzX ሷg cs,top cs,basecsy 𝑦τxz = 0 Hz ut=Xgeiϖt Figure 1: Scheme of the system considered in this study. The harmonic wave-propagation problem to be considered, in terms of the total displacement amplitude ^ut, is [Rovithis et al., 2011] d dz (z)d^ut dz +$2^ut= 0 s.t. ^ut(z=H) =Xg ,d^ut dz z=0= 0; (1) wherezis the coordinate stretching from the free surface to the rigid base, represents the density of the soil (which is assumed to be constant), the seismic S-wave excitation is modeled as a uniform displacement at the soil base Xgei$t, and(z)represents the variation of shear modulus with depth. Damping is not shown explicitly in eq. (1) as it is not of the viscous type but hysteretic, hence the shear modulus is interpreted as an storage modulus (real number) plus a loss modulus (imaginary number) as (z)!(z)(1 + id), wheredis referred as coefﬁcient of material damping. As per usual arrangement, the density is taken as constant throughout the deposit. For the sake of clarity, let us restate that the time variation comes given by ei$t, hence both Xgand^utare the amplitudes that accompany the phase: Xgis a positive real number, whereas ^utrepresents a complex number whose modulus is 2 APREPRINT - FEBRUARY 17, 2020 the total displacement magnitude and whose argument represents the time shift between load (input) and displacement (response, output). 1.2 Previous results In a previous work [ Garcia-Suarez et al., 2020 ], the authors focused on a speciﬁc, commonly used, distribution of S-wave velocity, the so-called “generalized parabola” (density is assumed to be constant): cs(z) cs;base=s (y) base= b+ (1b)z Hn = 1 + 1=n1y Hn ; (2) wherey=Hzis the complementary depth coordinate, b= 1=n= (cs;top=cs;base )1=nandnis termed “inho- mogeneity factor” (also referred to as “inhomogeneity parameter” [ Dakoulas and Gazetas, 1985 ]). The parameter b controls the maximum variation of cs, which in this case always corresponds to top-to-bottom or vice versa, since the function only allows for monotonic evolution. The dimensionless expression of the problem, in terms of relative displacements, was used to re-expresses (1) in the following manner: d d f2()d~u d +r2~u= 1;s.t. ~u(= 0) = 0 ,d~u d =1= 0; (3) wheref() =f(y=H) =cs(y=H)=cs;base ,=y=H ,r=$H=cs;base and~u=c2 s;base ^u=($2XgH2),^u=Xg^ut being the amplitude of the relative displacement (relative to the base movement). f()was assumed to be given by eq. (2) for all the developments carried out in the aforementioned work. By analyzing the asymptotic behavior of the exact solution to (3)[Rovithis et al., 2011 ], the authors found that the base-to-top dynamic ampliﬁcation could be expressed as A( ;n;r ) =p cs;base=cs;top cos $ cs;eq=H+Ocs;eq=H $ ; (4) withcs;eqgiven by cs;eq=H ZH 0dy cs(y)!1 ; (5) that is, the harmonic mean of the S-wave velocity distribution over the deposit. Since the “error” term in eq. (4), O(:::), decays as the inverse of the frequency, its magnitude is guaranteed to be negligible if $>>cs;eq=H, hence retaining the ﬁrst addend would be enough as to approximate the overall response in this regime. Likewise, in this regime the resonant peaks behaved as it the system was homogeneous having an equivalent material damping given by d;eq=dq cs;top=cs;base: (6) Bear in mind, throughout this whole document, that any result concerning “high excitation frequency” could also be expressed in terms of “long wavelengths” of wave propagating within the soil. The qualitative behavior put forward by eqs. (4) and (5) had been established previously [ Zhao, 1997 ], but it had not been justiﬁed mathematically in terms of asymptotic behavior of the exact solution. The previous result, on one hand, provided a simple estimate of the frequency response, in terms of amplitude and natural frequencies, and, on the other one, entailed the possibility of unequivocally deﬁning a “homogeneous-equivalent” properties so as to obtain a response similar to the one of the exact solution of (3)in the high-frequency regime, provided the S-wave velocity distribution being furnished by eq. (2). The obvious next move was to ascertain the scope of these results, verifying under what conditions they could be generalized to other proﬁles f(). This text is occupied with exactly that task. 2 Application of Geometrical Optics For the purpose of generalizing the aforementioned results, only two requirements will be demanded of the function f(): 3 APREPRINT - FEBRUARY 17, 2020 1)The proﬁle is continuous and bounded with continuous and bounded ﬁrst derivative, mathematically: f(y)2 C1([0;H]). This precludes considering layered proﬁles with sharp transitions (inter-layer impedance). 2)Only positive values of velocity are allowed, that is, f(y)>08y2[0;H]. This assumption prevents considering the case of zero stiffness at the free-surface [ Rovithis et al., 2011 ] (this case, it will be shown, can be considered as a generalization of the results to be obtained though). It is to be proved that an expression equivalent to eq. (4) holds for any distribution f()verifying just these two conditions. To do so, it turns out to be more convenient re-adopting the formulation in terms of total displacements (1), which, in dimensionless form reads d d f2()d~ut d +r2~ut= 0; (7a) or f2()d2~ut d2+ 2f()f0()d~ut d+r2~ut= 0; (7b) where in this case the non-dimensional total displacement is ~ut= ^ut=Xg, and all the other dimensionless parameters and variables have been presented earlier. The boundary conditions of (1) are, now in dimensionless form, ~ut(= 0) = 1 andd~ut d =1= 0: (8) Recall the high-frequency regime represents the focus of this work; such a regime corresponds to r=$H=cs;base 1, therefore= 1=r 1may be considered as a small parameter [ Simmonds and Mann, 2013 ]. Thus, rewrite eq. (7b) multiplying by 2, and also dividing by f2() = 0 : 2d2~ut d2+ 22f0() f()d~ut d+1 f2()~ut= 0; (9a) and, for the sake of brevity, call q() = 1=f2()(see that this function is always positive and bounded based on the restrictions of positive and non-zero wave velocity) and m() =f0()=f()(ﬁnite and positive too) 2d2~ut d2+ 22m()d~ut d+q()~ut= 0: (9b) This is a singular perturbation problem with spatially-changing coefﬁcients. The classic WKB method [Bender and Orszag, 2013 ] is the tool-of-choice to ﬁnd an approximate solution. See that the condition f()6= 0 8y2[0;H]is critical as it precludes the existence of turning points [Bender and Orszag, 2013]. In a nutshell, the WKB method advocates the use of the “ansatz” ~ut()exp 1 1X n=0nSn()! (10) in eq. (9b). A painstaking derivation can be found in the appendix at the end of this text. Let us proceed to state the main results that utilizing the method yields: S0() =iZ 1p q(0)d0=iZ 1d0 f(0); (11) recall thatq() = 1=f2() = (cs;base=cs(y))2; moreover S1() =C+ ln(f()1=2); (12) whereCrepresents a constant. Combining this two contributions of the phase (the only ones that do not vanish in the high-frequency limit) one obtains the physical optics approximation, which for the base-to-top ampliﬁcation yields A($) =f(= 1)1=2 cos $H cs;eq=p cs;base=cs;top cos $H cs;eq; (13) wherecs;eqis the harmonic mean deﬁned in eq. (5). This is the same expression as eq. (4), neglecting the “error” terms, therefore this result implies that the expression is valid for any proﬁle abiding by the two constraints enunciated at the beginning of the section, regardless of other ancillary features. It is instructive to note that utilizing just S0yields the denominator and the harmonic mean homogeneous-equivalent velocity whereas including S1adds the amplitude correction associated to the heterogeneity. The formal validity of the approximation is discussed in the appendix. If one wants to display the damping factor, d, explicitly in eq. (13), just substitute cs;eqbycs;eqp1 + id. 4 APREPRINT - FEBRUARY 17, 2020 2.1 Generalization to evolving density A similar result would have been obtained had we allowed the density to be a function of depth. Had we done so, considering=(y) =baseg(),=(y) =baseh()(hencef()2=g()=h()), satisfying the same conditions as discussed for f(), the auxiliary coefﬁcients m()andq()would have become m() =1 2g0() g(); (14a) q() =h() g(): (14b) These expressions yield a generalization of eq. (13): A($) =(h(= 1)g(= 1))1=4 cos rR1 0d f(); (15) reversing the non-dimensionalization, A($) =basebase toptop1=4 sec!H cs;eq ; (16) wherecs;eqis still being given by the harmonic mean of f() =p g()=h(), eq. (5), and now the change of amplitude depends on the base-to-top relative evolution both of stiffness and the density and density (observe that eq. (13) simply corresponds to base=top). 3 Physical signiﬁcance 3.1 Interpretation The physical interpretation of these results is straightforward once they are framed in terms of wavelength: high frequencies entail relatively short wavelengths, short with respect to H, propagating inside the material bulk. Short wavelengths tend to “smooth out” those spatial features of the medium larger than themselves, and thus, one may say, the wave does not “see” any shift in stiffness and it behaves as if there was only uniform soil all around, while long wavelengths “see” every stiffness alteration, and hence it tends to displace “chunks” of truly heterogeneous soil. The reason for obtaining the harmonic mean instead of the arithmetic mean derive from the fact that in this problem the travel distance is kept constant, not the travel time. Let us illustrate this point carefully: the stratum in Figure 1 can be conceptualized as a series of inﬁnitely many slices each with the same inﬁnitesimal thickness, h=H=N where the number of slices, N, tends to inﬁnity, hence we calculate the average velocity of the propagating wave as the ratio between the total travel distance over the total time that the journey takes cs;mean =PN i=1hiPN i=1Ti=HPN i=1(H=Ncs;i)!1RH 0d(y=H ) cs(y)=cs;eq (17) asNtends to inﬁnity. The subscript idenotes properties concerning the i-th discrete slice. This means that the proper average to consider in this kind of problem is the harmonic mean. This interpretation is by no means original, it is simply an adaptation of well-known results in Classical Optics (i.e., propagation of light as a electromagnetic wave) which date back to the works of Fermat and Huygens. 3.2 This result in the broad context of Ray Theory The previous results can be considered as the outcome of application of so-called ray approximation methods to the problem (1). These employ the short-wavelength (high-frequency) limit to analyze wave propagation phenomena (electromagnetic, acoustic or elastic). Ray methods are a powerful tool that has already been used in Elastodynamics; see for example the monograph by Prof. Achenbach and coworkers [ Achenbach et al., 1982 ], wherein these methods are applied to the study of the effect of cracks over waves propagating within homogeneous elastic bodies. The physical intuition behind these results has been common currency in the physics community, including seismology, see Aki and Richards’ treaty on Quantitative Seismology [ Aki and Richards, 2002 ], in particular, Chapter 4, section 5 APREPRINT - FEBRUARY 17, 2020 4, concerning ray methods as a means to consider seismic wave propagation in the heterogeneous earth. However, these appraisals had not been passed on to the geotechnical engineering community at large nor, it appears, explic- itly spelled out in the mathematical terms of the WKB approximation applied to considerations of unidimensional propagation in a ﬁnite domain, although similar expansions are usually discussed in the geometrical optics literature [Achenbach et al., 1982 ,Kravtsov and Orlov, 1990 ], including the renowned Debye expansion , which does not rely on the deﬁnition of small parameters directly. This intuition also allows to reply to the following question: “what does it happen when f(y)=2C1([0;H])but f(y)2C0([0;H])(continuous function, non-continuous ﬁrst derivative) or even f(y)2C1([0;H])(non-continuous function)?”. The reasoning that was followed in section 2 is forbidden since f0()cannot be deﬁned. The getaway is found in Fermat’s Principle: basically, Fermat’s Principle states that the problem we are solving possesses a variational structure, and hence the equations we are solving are nothing but its corresponding Euler’s equations [Gelfand et al., 2000 ]. If the functions describing the parameter evolution are not continuous enough, this may prevent us from deﬁning the Euler’s equations properly, but the underlying stationarity condition [ Kravtsov and Orlov, 1990 ] remains in effect. The phenomenon at the intersection of different mediums can be described by an analogue of the Snell’s law in Optics. The only “exotic” feature that this speciﬁc instantiation introduces concerns where the inhomogeneity is localized. Normally, the wave propagates in a medium with varying refraction index, hence the elliptic operator that appears in Euler’s equations is simply a Laplacian, which in the 1D-case boils down to a neat second derivative in space, @2()=@y2, where the dot represents the dependent variable; in this case, at least part of the homogeneity is conﬁned “within the operator” due to the change in stiffness of the medium, hence the operator one has to deal with is @((y)@()=@y)=@y, hence this problem can be considered as an instantiation of Sturm-Liouville Theory [Weinberger, 2012]. In order to perform a one-to-one mapping between the classic theory and prior results, let us state that the ﬁrst component to the phase, S0is equivalent to the eikonal function (is actually equal to it times an extra factor i), the refraction index corresponds to the inverse of the velocity distribution, and, just like the general theory foresees, the eikonal function is proportional to the integral of the refraction index over the arc-length of the trajectory followed by the ray [ Kravtsov and Orlov, 1990 ], which in this case happen to be straight vertical lines. The second component of the phase,S1, yields the leading-order approximation of the amplitude, which is provided by solving a transport equation. Moreover, the ampliﬁcation factor in both eq. (16) and eq. (13) can be retraced to the ray divergence parameter, which is but a consequence of the energy ﬂux conservation in the system. See [Kravtsov and Orlov, 1990] for details. Before concluding this digression, let us mention, on one hand, that the case corresponding to d=d(y) 1 can be assimilated into a classic scenario dubbed as medium with weak absorption , and, secondly, that the case of vanishing stiffness at the top of the layer can be understood as the presence of a caustic surface at the top of the stratum [Stavroudis, 2012]. 4 Comparison to numerical results Let us illustrate how and when the approximate formula eq. (13) agrees well with an exact solution that could still be obtained numerically for cases of complicated proﬁles. In order to illustrate this, two “toy” proﬁles shall be utilized. 4.1 Increasing contrast First, let us consider the inﬂuence of sudden mechanical property changes on the ﬁtness of the approximation. For such purpose, consider f() = 13 8 1 + tanh p 1 3 : (18) This proﬁle gives a contrast between base and top of 1=4, where all the variation is concentrated at = 1=3. The parameterpcontrols the gradient in that zone, the higher pthe quicker the evolution and the more sudden the change. In the limitp!1 the function becomes an step function and the gradient is inﬁnite, so the condition f()2C1([0;1]), upon which eq. (13) was obtained, stops holding. Three different cases have been considered using Mathematica [Wolfram, 2000], by numerically solving eq. (7b) for three increasing values of p(5, 10, and 100). We observe in Figure 2 how the approximate formula eq. (13) successfully captures the behavior in most of the frequency range, but the agreement worsens as the smooth change becomes a localized snap at p= 100 . This result illustrates how it behooves the gradient to remain bounded everywhere for the approximation to match the exact solution. 6 APREPRINT - FEBRUARY 17, 2020 4.2 Proﬁle reversal This time, let us consider a proﬁle with reversals (softer regions in between stiffer ones). For this task, use f() = 1 +  sin K 2 ; (19) where controls the amplitude of the reversals whereas Kdoes the same for the number of reversals that happen within the layer. This example is intended to conﬁrm the lack of inﬂuence of the monotonicity of f()over the result eq. (13). On can acknowledge in Figure 3 that the number of oscillations, K, and their amplitude , in and of itself, do not inﬂuence the convergence of the approximation to the exact value, but the gradient, which in this case is linearly proportional to both Kand, is the parameter to control, and obviously the more turns within the ﬁnite layer and the acuter these are have a direct translation on the gradient magnitude. This explains the deterioration from the ﬁst case to the last one. This display embodies the irrelevance of monotonicity requirements on the proﬁle evolution for eq. (13) to hold. 7 APREPRINT - FEBRUARY 17, 2020 0 1 2 3 4 5 60153045 a) p=5 A()Asymptotics Exact 0 0.5 101 y/H 0 1 2 3 4 5 60153045 b) p=10 A() 0 0.5 101 y/H 0 1 2 3 4 5 6 /s,base0153045 c) p=100 A() 0 0.5 1 cs(y) / cs,base01 y/H Figure 2: Base-to-top dynamic ampliﬁcation A($)(with corresponding vertical S-wave proﬁle cs(y)), comparison between numerical solution of eq. (1) ( Exact ) and eq. (13) ( Asymptotics ) 8 APREPRINT - FEBRUARY 17, 2020 0 2 4 6 8 100153045 a) =0.5 ; K=3 A()Asymptotics Exact 0 0.5 1 1.5 201 y/H 0 2 4 6 80153045 b) =0.9 ; K=5 A() 0 0.5 1 1.5 201 y/H 0 2 4 6 8 /s,base0153045 c) =0.7 ; K=8 A() 0 0.5 1 1.5 2 cs(y) / cs,base01 y/H Figure 3: Base-to-top dynamic ampliﬁcation A($)(with corresponding vertical S-wave proﬁle cs(y)), comparison between numerical solution of eq. (1) ( Exact ) and eq. (13) ( Asymptotics ) 9 APREPRINT - FEBRUARY 17, 2020 5 Implications for NIST and ASCE recommendations Current recommendations [ NIST, 2012 ] by the National Institute of Standards and Technology (NIST) regarding Soil-structure Interaction (SSI) instruct the use of the harmonic mean of shear-wave velocity distribution, for any value of excitation frequency , when considering SSI (rocking, swaying and other phenomena) of foundations resting over heterogeneous soil. This recommendation may be re-evaluated in light of the results displayed previously, as it appears to be unﬁtting to use the harmonic mean of the velocity distribution as the equivalent one for the whole range of excitation frequencies; such a value should be limited to speciﬁc combinations of frequency, soil mechanical properties csand depth to bedrock H so that$H=cs 1. When it comes to deﬁne equivalent homogeneous property to consider the amplitude of the fundamental mode in a 1D site response scenario that includes heterogeneity, our results indicate that using harmonic average to be unﬁtting except in cases having very minor heterogeneity. Other results published by the authors suggest, contrari- wise, properties of the soil layer in contact with bedrock should be used as these control the fundamental response [Garcia-Suarez and Asimaki, 2019]. 6 Conclusions and future work 6.1 Conclusions This work has been occupied, chieﬂy, with a generalization of previous results obtained under the assumption of S-wave velocity proﬁle abiding by eq. (2) [ Garcia-Suarez et al., 2020 ] to any proﬁle meeting two simple requirements: 1) continuity and boundedness of ﬁrst derivative of f(), 2) positive (non-zero) values of csthroughout the deposit, including at the free-surface. It has been shown that this result stems from a classical result in Physics, the so-called Geometrical Optics approximation. These results also settle partially the debate over the equivalent-homogeneous properties that should be used to describe the response of inhomogeneous soil deposits overlying rigid bedrock: provided $H=cs;base 1and the two constraints in Section 2 are satisﬁed, eq. (5) described the equivalent-homogeneous S-wave velocity, whereas the resonance peaks behave as if the material damping came given by eq. (6). In spite of the application , if the amplitude as small as for the linear-elastic response assumption to remain in effect, the approximation eq. (13) may still be applied to a portion of the earthquake spectrum. 6.2 Future work Since the approximation is not limited to ﬁnite layers with zero impedance at the bottom, the previous results could be extended to waves propagating vertically in a halfspace, hence generalizing the results to those scenarios wherein radiational damping [Kramer, 1996] has to be accounted for. 7 Acknowledgments The authors are thankful to Prof. Michael Ortiz for his input regarding connections to Geometrical Optics Theory at large. The format used in this preprint has been dispensed by Mr. George Kour through Overleaf under Creative Commons CC BY 4.0.