A conservation law for testing methods of prediction of the seismic wave response of a protuberance emerging from flat ground
aa r X i v : . [ phy s i c s . g e o - ph ] J a n A conservation law for testing methods of prediction ofthe seismic wave response of a protuberance emergingfrom flat ground
Armand Wirgin ∗ January 22, 2020
Abstract
We establish the equations which translate a conservation law for the problem of the seismicresponse of an above-ground structure (e.g., building, hill or mountain) of arbitrary shape andinquire whether both the implicit (formal) and explicit (numerical) solutions for the responseobey this law for the case of a cylindrical, rectangular protuberance. Both the low-order (poorapproximations of the response) as well as higher-order (supposedly better approximations) turnout to satisfy the conservation of flux relation, which means that the satisfaction of this relationis a necessary, but not sufficient, means for determining whether a solution to the scatteringproblem is valid.
Keywords: seismic response, conservation law, energy, flux, above-ground feature.Abbreviated title: Conservation of flux relative to the seismic response of a protuberanceCorresponding author: Armand Wirgin,e-mail: [email protected] ∗ LMA, CNRS, UPR 7051, Aix-Marseille Univ, Centrale Marseille, F-13453 Marseille Cedex 13, France,( [email protected] ) ontents M for fixed frequency and shear modulus in Ω . . . 239.3 Effect of the variation of the frequency f for various M and µ [1] . . . . . . . . . . . . 24
10 Conclusions 26 Introduction
Seismic waves are known to damage buildings or industrial facilities located either on flat orhilly ground and thought to provoke land and rock slides on hills and mountains [5, 6, 7, 9, 10,11, 13, 14, 15, 18, 20, 23, 25, 27, 29, 31, 34, 45]. These man-made or natural features can begrouped into the term: above-ground structure (AGS) or protuberance for short. The seismicresponse of protuberances is an ongoing research theme in theoretical and applied seismology,mainly because of the variety of AGS’s, the social and economic impact of the damage issue anda relatively-poor understanding of the empirical evidence that the seismic wave field is amplifiedon the stress-free boundary of the protuberance relative to this field on flat ground. For thesereasons, a great variety of (mainly-numerical) solutions to specific AGS seismic response problemshave been proposed [2, 3, 4, 5, 6, 8, 10, 12, 13, 15, 17, 19, 21, 22, 23, 24, 25, 27, 28, 30, 31, 32, 33,35, 36, 37, 38, 39, 40, 42, 43, 44], but the question that is often ignored or eluded is: how good arethese solutions?If the underlying boundary-value problem (BVP)is correctly formulated then the ultimate testof whether a solution is valid or not is whether it satisfies the equations inherent in the BVP, butto carry out this test is often very painstaking because it requires generating the displacement fieldand its gradient at all points of the ground and on the upper, curved, boundary of the protuberance,not to speak of the field everywhere within and below the protuberance. Another manner to testthe solution is by finding out if it satisfies a conservation law (such as that of energy). Since weshall be concerned with frequency-domain formulations of the BVP, the conservation of energy lawtakes the form of what we term the conservation of flux law which states that the input flux equalsthe scattered flux plus the absorbed flux [41].We show that, to employ this law in its first, abstract form, requires computing the scatteringamplitude in the far-field zone (whether the protuberance filler medium is lossy or non-lossy) aswell as the displacement field at all points within the protuberance (only when the filler is lossy).The latter task (for lossy fillers) is likewise painstaking, so that it is opportune to show (as is doneherein) that the said task can be replaced by computing the field and its gradient only on the lower,flat (ground-level) boundary of the protuberance. To show that the so-contrived conservation (offlux) law makes sense, we apply it to the case of a cylindrical protuberance of rectangular shapesubmitted to a shear-horizontal plane body wave.The result of this operation is that the said solution (deriving from a domain decompositionseparation of variables (SOV) technique) indeed satisfies the conservation law for all orders ofapproximation of the solution which suggests that the SOV solution is, in a sense, exact. However,:i) if this law is not satisfied the solution cannot be reasonably qualified as exact (i.e., might even befallacious if the difference between the input and output (scattered plus absorbed) fluxes is large,and ii) even if this law is satisfied, the solution may be incorrect (this is demonstrated numericallyherein). This means that the conservation law is a necessary, but not sufficient condition for testingthe validity of a solution to the seismic response problem.
In the first approximation, the earth’s surface is considered to be (horizontally-) flat (termed”ground” for short) and to separate the vacuum (above) from a linear, isotropic, homogeneous(LIH) solid (below), so as to be stress-free. In the second approximation the flat ground is locally3eformed so as to penetrate into what was formerly the vacuum half space. We now define theprotuberance as the region between the locally-deformed stress-free surface and what was formerlya portion of the flat ground. This protuberance is underlain by the same LIH solid as previously, butthe solid material within the protuberance is now assumed to be only linear and isotropic (i.e., nothomogeneous). In fact, we consider the specific case in which the material within the protuberanceis in the form of a horizontal bilayer so as to be able to account for various empirically-observedeffects that are thought to be due to inhomogeneity of the protuberance material. Furthermore, weassume that: the protuberance is of infinite extent along one of the cartesian coordinates, and itsstress-free boundary to be of arbitrary shape (in its cross-section plane). The underlying problemof much of what follows is the prediction of the seismic wave response of this earth model.The earthquake sources are assumed to be located in the lower half-space and to be infinitely-distant from the ground so that the seismic (pulse-like) solicitation takes the form of a body (plane)wave in the neighborhood of the protuberance. This plane wavefield is assumed to be of the shear-horizontal ( SH ) variety, which means that: only one (i.e., the cartesian coordinate z ) componentof the incident displacement field is non-nil and this field does not depend on z .We shall assume that the protuberance boundary does not depend on z and that the (oftenrelatively-soft) medium filling the protuberance as well as the (usuallly relatively-hard) mediumbelow the protuberance are both linear and isotropic. Furthermore the medium of the below-ground half space is assumed to be homogeneous, whereas that of the protuberance to be piecewisehomogeneous (however, this heterogeneity is such as to not depend on z ). It ensues that thescattered and total displacement fields within and outside the protuberance do not depend on z .Thus, the problem we are faced with is 2D ( z being the ignorable coordinate), and it is sufficientto search for the z -component of the scattered displacement field, designated by u sz ( x ; ω ) in thesagittal (i.e., x − y ) plane, when u iz ( x ; ω ) designates the incident displacement field, with x = ( x, y )and ω = 2 πf the angular frequency, f the frequency. The temporal version of the displacementfield is u z ( x ; t ) = 2 ℜ R ∞ u iz ( x ; ω ) exp( − iωt ) dω wherein t is the temporal variable.Fig. 1 describes the scattering configuration in the sagittal plane. In this figure, k i = k i ( θ i , ω )is the incident wavevector oriented so that its z component is nil, and θ i is the angle of incidence.The portion of the ground outside the protuberance is stress-free but since the protuberanceis assumed to be in welded contact with the surrounding below-ground medium, its lower, flat,boundary is the locus of continuous displacement and stress, as is requisite for the incident field tobe able to penetrate into the protuberance and then be scattered outside the protuberance in theremaining lower half space.Most of what is offered in this study does not imply any restrictions either on the shape of thestress-free portion of the boundary of the protuberance or on the (lossy or non-lossy) nature of themedium filling the protuberance.The three media (other than the one of the portion of the space above the protuberance, beingoccupied by the vacuum, is of no interest since the field cannot penetrate therein) are M [ l ] ; l =0 , , µ [ l ] ; l = 0 , , β [ l ] ; l = 0 , , β [ l ] = β ′ [ l ] + iβ ′′ [ l ] , with β ′ [ l ] ≥ β ′′ [ l ] ≤ β [ l ] = q µ [ l ] ρ [ l ] , and ρ [ l ] the (generally-complex) mass density. The shear-wave velocity β [0] is assumed to be real, i.e., β ′′ [0] = 0. 4igure 1: Sagittal plane view of the 2D scattering configuration. The protuberance occupies theshaded areas and the medium within it is a horizontal bilayer. The protuberance occupies (in the sagittal plane (SP)) the finite-sized region Ω S Ω . Thebelow-ground half-space occupies the region Ω . Ω is entirely filled with M [0] whereas Ω is filledwith M [1] and Ω with M [2] .Always in the sagittal plane, the flat ground is described by Γ G , with x, y the cartesian coordi-nates in the SP) and is composed of three segments; Γ l , Γ m , and Γ r , which designate the left-hand,middle, and right-hand portions respectively of Γ G . The protuberance is an above-ground structurewhose upper and lower boundaries (in the SP) are Γ p and Γ m , the latter being of width w .The analysis takes place in the space-frequency framework, so that all constitutive and fieldvariables depend on the frequency f . This dependence will henceforth be implicit (e.g., u z ( x ; f ),with x = ( x, y ), will be denoted by u ( x )).The seismic solicitation is an incident shear-horizontal (SH) plane wave field of the form u i ( x ) = a i exp( i k i · x ) = a i exp[ i ( k ix x + k iz y )] , (1)wherein a i = a i ( ω ) is the spectral amplitude of the seismic pulse, k i = ( k ix , k iy ), k ix = k [0] sin θ i , k iy = k [0] cos θ i , k [ l ] = ω/β [ l ] ; l = 0 , , as: u ( x ) = u [ l ] ( x ) ; ∀ x ∈ Ω l , l = 0 , , , (2)with the understanding that these fields satisfy the 2D SH frequency domain elastic wave equation(i.e., Helmholtz equation) (cid:16) △ + (cid:0) k [ l ] (cid:1) (cid:17) u [ l ] ( x ) = 0 ; ∀ x ∈ Ω l , l = 0 , , , (3)with the notations △ = ∂ ∂x + ∂ ∂y in the cartesian coordinate system of the sagittal plane.In addition, the field u [0] satisfies the radiation condition u [0] ( x ) − u i ( x ) − ∼ outgoing wave ; k x k → ∞ . (4)due to the fact that Ω is unbounded (i.e., a semi-infinite domain).The stress-free nature of the boundaries Γ l , Γ p , Γ r , entail the boundary conditions: µ [0] u [0] ,y ( x ) = 0 ; ∀ x ∈ Γ l + Γ r , (5) µ [2] u [2] ,y ( x ) = 0 ; ∀ x ∈ Γ p , (6)wherein u ,ζ denotes the first partial derivative of u with respect to ζ .The fact, that the horizontal segment Γ between the two media filling the protuberance isassumed to be an interface across which two media are in welded contact, entails the continuityconditions: u [2] ( x ) − u [1] ( x ) = 0 ; ∀ x ∈ Γ , (7) µ [2] u [2] ,y ( x ) − µ [1] u [1] ,y ( x ) = 0 ; ∀ x ∈ Γ . (8)Finally, the fact, that Γ m was assumed to be an interface across which two media are in weldedcontact, entails the continuity conditions: u [0] ( x ) − u [1] ( x ) = 0 ; ∀ x ∈ Γ m , (9) µ [0] u [0] ,y ( x ) − µ [1] u [1] ,y ( x ) = 0 ; ∀ x ∈ Γ m , (10)The purpose of addressing such a boundary-value (direct) problem is to determine u [ l ] ( x ); l =0 , , l ; l = 0 , ,
2. The principal ambition of the analysis which follows is rather to establish aconservation law governing u [ l ] ( x ) ; l = 0 ,
1, this being done, in the first part of our study, withoutactually solving for u [ l ] ( x ) ; l = 0 ,
1. However, to show that this conservation law makes senseand is useful, we shall, in the second part of this study, appeal to a separation-of-variables (SOV)solution of the problem in which the protuberance is of rectangular shape.6
Basic ingredients of the conservation of flux law
Eq. (3) yields (cid:16) △ + (cid:2)(cid:0) k [ l ] (cid:1) (cid:3) ∗ (cid:17) u [ l ] ∗ ( x ) = 0 ; ∀ x ∈ Ω l , l = 0 , , , (11)wherein ( X + iY ) ∗ = X − iY . It follows (with d Ω the differential surface element) that Z Ω l u [ l ] ∗ x ) n(cid:16) △ + (cid:0) k [ l ] (cid:1) (cid:17) u [ l ] ( x ) − u [ l ] ∗ ( x ) (cid:16) △ + (cid:2)(cid:0) k [ l ] (cid:1) (cid:3) ∗ (cid:17) u [ l ] ( x ) o d Ω = 0 ; ∀ x ∈ Ω l , l = 0 , , , (12)or Z Ω l n u [ l ] ∗ ( x ) △ u [ l ] ( x ) − u [ l ] ( x ) △ u ∗ [ l ] ( x ) o d Ω+ Z Ω l n(cid:0) k [ l ] (cid:1) − (cid:2)(cid:0) k [ l ] (cid:1) (cid:3) ∗ o k u [ l ] ( x ) k d Ω = 0 ; ∀ x ∈ Ω l , l = 0 , , . (13)We want to apply Green’s second identity to the first integral, and to do this we must define the(closed) boundaries of Ω l . We already know that the boundaries of Ω ; l = 1 , was not closed. To close it, we imagine a semicircle Γ R , of large radius R (taken tobe infinitely large in the limit), centered at the origin O , to be drawn so as to intersect the groundat x = ±R and to intersect the y axis at y = −R . Thus, designating by Γ ± p the upper(lower)portions of Γ p , the closed boundaries of Ω l ; l = 0 , , ∂ Ω = Γ + p ∪ Γ , ∂ Ω = Γ m ∪ Γ − p ∪ Γ , ∂ Ω = Γ l ∪ Γ b ∪ Γ r ∪ Γ R→∞ , (14)and we shall designate by ν l the unit vector normal to ∂ Ω l that points towards the exterior of ∂ Ω l .We now apply Green’s second identity to obtain (with d Γ the differential arc element) Z ∂ Ω l n u [ l ] ∗ ( x ) ν l · ∇ u [ l ] ( x ) − u [ l ] ( x ) ν l · ∇ u ∗ [ l ] ( x ) o d Ω+ Z Ω l n(cid:0) k [ l ] (cid:1) − (cid:2)(cid:0) k [ l ] (cid:1) (cid:3) ∗ o k u [ l ] ( x ) k d Ω = 0 ; ∀ x ∈ Ω l , l = 0 , , . (15)or, in condensed form (on account of the non-lossy nature of the solid filling Ω and the homogeneousnature of the solids filling Ω l ; l = 1 , ℑ Z ∂ Ω u [0] ∗ ( x ) ν · ∇ u [0] ( x ) d Γ = 0 , (16) ℑ Z ∂ Ω l u [ l ] ∗ ( x ) ν l · ∇ u [ l ] ( x ) d Γ − ℑ (cid:2)(cid:0) k [ l ] (cid:1) (cid:3) Z Ω k u [ l ] ( x ) k d Ω = 0 ; l = 1 , . (17)More explicitly, and on account of (14), ℑ Z Γ l +Γ r u [0] ∗ ( x ) ν · ∇ u [0] ( x ) d Γ + ℑ Z Γ m u [0] ∗ ( x ) ν · ∇ u [0] ( x ) d Γ+ ℑ Z Γ ∞ u [0] ∗ ( x ) ν · ∇ u [0] ( x ) d Γ = 0 , (18)7 Z Γ m u [1] ∗ ( x ) ν · ∇ u [1] ( x ) d Γ + ℑ Z Γ u [1] ∗ ( x ) ν · ∇ + ℑ Z Γ − p u [1] ∗ ( x ) ν · ∇ u [1] ( x ) d Γ u [1] ( x ) d Γ+ ℑ (cid:2)(cid:0) k [1] (cid:1) (cid:3) Z Ω k u [1] ( x ) k d Ω = 0 , (19) ℑ Z Γ u [2] ∗ ( x ) ν ·∇ u [1] ( x ) d Γ+ ℑ Z Γ p u [2] ∗ ( x ) ν ·∇ u [2] ( x ) d Γ+ ℑ (cid:2)(cid:0) k [2] (cid:1) (cid:3) Z Ω k u [2] ( x ) k d Ω = 0 . (20)Due to ν ( x ) = − ν ( x ) ; ∀ x ∈ Γ b , and the boundary and continuity conditions, we have: u [2] ∗ ( x ) ν · ∇ u [2] ( x ) = − µ [1] µ [2] u [1] ∗ ( x ) ν · ∇ u [1] ( x ) ; ∀ x ∈ Γ , (21)so that, since µ [1] µ [2] was assumed to be real, (20) becomes: − µ [1] µ [2] ℑ Z Γ u [1] ∗ ( x ) ν · ∇ u [1] ( x ) d Γ + ℑ (cid:2)(cid:0) k [2] (cid:1) (cid:3) Z Ω k u [2] ( x ) k d Ω = 0 . (22)Due to the stress-free boundary condition, (19) reduces to ℑ Z Γ m u [1] ∗ ( x ) ν · ∇ u [1] ( x ) d Γ + ℑ Z Γ u [1] ∗ ( x ) ν · ∇ + ℑ (cid:2)(cid:0) k [1] (cid:1) (cid:3) Z Ω k u [1] ( x ) k d Ω = 0 , (23)so that the linear combination of these last two equations gives rise to: ℑ Z Γ m u [1] ∗ ( x ) ν · ∇ u [1] ( x ) d Γ + µ [2] µ [1] ℑ (cid:2)(cid:0) k [2] (cid:1) (cid:3) Z Ω k u [2] ( x ) k d Ω + ℑ (cid:2)(cid:0) k [1] (cid:1) (cid:3) Z Ω k u [1] ( x ) k d Ω = 0 . (24)Due to ν ( x ) = − ν ( x ) ; ∀ x ∈ Γ b , and the boundary and continuity conditions, we have: u [0] ∗ ( x ) ν · ∇ u [0] ( x ) = − µ [1] µ [0] u [1] ∗ ( x ) ν · ∇ u [1] ( x ) ; ∀ x ∈ Γ m , (25) u [0] ∗ ( x ) ν · ∇ u [0] ( x ) = 0 ; ∀ x ∈ Γ l + Γ r , (26)so that (18) becomes − µ [1] µ [0] ℑ Z Γ m u [1] ∗ ( x ) ν · ∇ u [1] ( x ) + ℑ Z Γ ∞ u [0] ∗ ( x ) ν · ∇ u [0] ( x ) d Γ = 0 , (27)whence the linear combination of (18) and (24) gives rise to: ℑ Z Γ ∞ u [0] ∗ ( x ) ν ·∇ u [0] ( x ) d Γ+ µ [2] µ [0] ℑ (cid:2)(cid:0) k [2] (cid:1) (cid:3) Z Ω k u [2] ( x ) k d Ω+ µ [1] µ [0] ℑ (cid:2)(cid:0) k [1] (cid:1) (cid:3) Z Ω k u [1] ( x ) k d Ω = 0 . (28)Eqs. (27) and (28) are alternate expressions of the same sought-for conservation law.8 Separation-of-variables representation of the field in the halfspace underneath the protuberance
The fact that the conservation law involves the field on Γ ∞ = Γ R→∞ means that we mustdispose of an expression for the field u [0] ( x ) in the half space (and notably in the far-field zonethereof) underneath the protuberance. To do this, we are not obliged to solve the forward scat-tering problem, but only obtain a representation of u [0] ( x ) that obeys the Helmholtz equation, theradiation condition and the stress-free boundary condition on Γ l + Γ r . To do this, we employ theseparation-of-variables (SOV) technique in terms of the cartesian coordinates to obtain u s = u [0] ( x ) − u i ( x ) − u r ( x ) = Z ∞−∞ B ( k x ) exp[ i ( k x x − k [0] y y )] dk x k [0] y ; ∀ y ≤ , ∀ x ∈ R , (29)wherein k [0] y = q(cid:0) k [0] (cid:1) − (cid:0) k x (cid:1) , ℜ k [0] y ≥ ∀ ω ≥ , ℑ k [0] y ≥ ∀ ω ≥ . (30)and u r ( x, y ) = u i ( x, − y ) . (31)It follows that: u i ( x,
0) + u r ( x,
0) = 2 u i ( x,
0) = 2 a i exp[ ik x x ] , (32) u i,y ( x,
0) + u r,y ( x,
0) = 0 , (33)whence u [0] ( x,
0) = 2 u i ( x,
0) + u s ( x, , (34) u [0] ,y ( x,
0) = u s,y ( x,
0) = − i Z ∞−∞ B ( k x ) exp[ i ( k x x ] . (35)Fourier inversion then yields B ( k x ) = i π Z ∞−∞ u [0] ,y ′ ( x ′ ,
0) exp[ − i ( k x x ′ ] dx ′ . (36)We assume that Γ m , whose width is w , extends from x ′ = − w/ x ′ = w/
2, and make use of thestress-free boundary condition on Γ l + Γ r to obtain, with the change of variables k x = k [0] cos φ : B ( φ ) = π B ( k x ) = i π Z w/ − w/ u [0] ,y ′ ( x ′ ,
0) exp[ − ik [0] x ′ cos φ ] dx . (37)The introduction of (36) into (29) gives, after the interchange of the orders of integration u s ( x ) = i π Z w/ − w/ u [0] ,y ′ ( x ′ , H (1)0 ( k [0] R ) dx ′ ; ∀ y ≤ , ∀ x ∈ R . (38)wherein R = p ( x − x ′ ) + y and H (1)0 ( k [0] R ) = 1 π Z ∞−∞ exp[ ik x ( x − x ′ ) − k y y )] dk x k [0] y ; ∀ y ≤ , ∀ x ∈ R (39)9s the zeroth-order Hankel function of the first kind [1].Since we are particularly interested in the field in the far-field zone, we appeal to the well-known[1] asymptotic form of the Hankel function H (1)0 (cid:0) k [0] k x − ( x ′ , k (cid:1) = H (1)0 (cid:0) k [0] R (cid:1) ∼ (cid:16) πk [0] R (cid:17) / exp[ i ( k [0] R − π/ k [0] R → ∞ . (40)The change of variables x = r cos φ , y = r cos φ , x ′ = r cos φ ′ (with the understanding that φ ′ = 0for x > φ ′ = π for x <
0) gives rise to R = p (( r cos φ − r ′ cos φ ′ ) + ( r sin φ ) ∼ r − x ′ cos φ ; r ′ r << , (41)the condition r ′ r << | x ′ | = r ′ ≤ w/ r → ∞ . Consequently, u s ( x ) ∼ B ( φ ) (cid:16) πk [0] r (cid:17) / exp[ i ( k [0] r − π/ | x ′ | = r ′ ≤ w/ , k [0] r → ∞ . (42)wherein B ( φ ) is as in (37).What this all means is that the scattered field behaves asymptotically (i.e., for k [0] r → ∞ ) likea cylindrical wave whose complex amplitude is B ( φ ). We now dispose of the means for evaluatingthe integral on Γ ∞ in the conservation of flux relations. Now let us return to (27)-(28) which can be written either as I − J = 0 . (43)or I + K = 0 . (44)with I = ℑ Z Γ ∞ u [0] ∗ ( x ) ν · ∇ u [0] ( x ) d Γ , (45)and let us examine I more closely. On account of (34) we have I = ℑ Z Γ ∞ [ u ∗ i + u ∗ r + u ∗ s ] ν · ∇ [ u i + u r + u s ] d Γ = ( I ii + I rr ) + ( I ir + I ri ) + ( I is + I si ) + ( I rs + I sr ) + I ss , (46)or, due to the facts that: d Γ (cid:12)(cid:12) Γ R = R dφ and ν · ∇ u [0] (cid:12)(cid:12) Γ R = u [0] ,r ( R , φ ), we obtain I = ℑ lim R→∞ Z ππ [ u i ∗ ( R , φ ) + u r ∗ ( R , φ ) + u s ∗ ( R , φ )][ u i,r ( R , φ ) + u r,r ( R , φ ) + u s,r ( R , φ )] R dφ , (47)in which we must adopt the following cylindrical coordinate representations of u i and u r , wherein x = r cos φ and y = r sin φ : u i = a i exp[ ir ( k ix cos φ + k iy sin φ )] , u r = a i exp[ ir ( k ix cos φ − k iy sin φ )] . (48)10t is easy to show [41] that: I ii + I rr = I ir + I ri = I is + I si = 0 . (49)In [41] we showed that the employment of the stationary phase technique enables to find: I rs + I sr = 4 ℜ (cid:2) a i ∗ B (3 π/ θ i ) (cid:3) . (50)The employment of (42) yields (see [41]) I ss = 2 π Z ππ k B ( φ ) k dφ . (51)The last step is to carry out the sum in (46 so as to obtain I = 2 π Z ππ k B ( φ ) k dφ + 4 ℜ (cid:2) a i ∗ B (3 π/ θ i ) (cid:3) , (52)which is the detailed expression of the term I in the conservation law which we wrote either as I − J = 0 or I + K . Let us now examine in detail the term J of this law. The definition of J is J = µ [1] µ [0] ℑ Z Γ b u [1] ∗ ( x ) ν · ∇ u [1] ( x ) d Γ , (53)but, as underlined earlier, our ambition was not to solve the boundary-value problem, namely forthe field u [1] ( x ) within the basin, so that we cannot go beyond expressing the conservation law as I − J = 2 π Z ππ k B ( φ k dφ + 4 ℜ (cid:2) a i ∗ B (3 π/ θ i ) (cid:3) − µ [1] µ [0] ℑ Z Γ b u [1] ∗ ( x ) ν · ∇ u [1] ( x ) d Γ = 0 , (54)whose meaning, is unfortunately not clear for the moment.To cope with this problem, we make use of the alternative expression (27) of the conservationlaw I + K = ℑ Z Γ ∞ u [0] ∗ ( x ) ν · ∇ u [0] ( x ) d Γ + X j =1 µ [ j ] µ [0] ℑ (cid:2)(cid:0) k [ j ] (cid:1) (cid:3) Z Ω j k u [ j ] ( x ) k d Ω = 0 , (55)which authorizes us to write (28) as I + K = 2 π Z ππ k B ( φ ) k dφ + 4 ℜ (cid:2) a i ∗ B (3 π/ θ i ) (cid:3) + X j =1 µ [ j ] µ [0] ℑ (cid:2)(cid:0) k [ j ] (cid:1) (cid:3) Z Ω j k u [ j ] ( x ) k d Ω = 0 . (56)This expression informs us that if the media within the protuberance domain Ω ∪ Ω are non-lossy, then ℑ k [1] = ℑ k [2] = 0 so that K = P j =1 µ [ j ] µ [0] ℑ (cid:2)(cid:0) k [ j ] (cid:1) (cid:3) R Ω j k u [ j ] ( x ) k d Ω = 0, which entails11 = µ [1] µ [0] ℑ R Γ b u [1] ∗ ( x ) ν ·∇ u [1] ( x ) d Γ = 0. This means that J accounts for damping (i.e., absorption)in Ω since J vanishes in the absence of a loss mechanism, as expressed by c [ j ] ; j = 1 , k [ j ] ; j = 1 ,
2) being complex. A valid question is then: damping/absorption of what? Itis not energy because the units of J = − K are not those of energy, but rather something relatedto energy, which we shall name ’flux’. If we reason in terms of energy, the only means by whichenergy can be lost in a scattering problem such as ours is by radiation damping, which is themechanism by which (scattered) energy escapes to the outer reaches of Ω , i.e., escapes to r → ∞ in the half-space y < y . As statedearlier, radiation damping is exclusively related to the scattered wave portion of the total field,and the function that expresses this relation to the scattered field is B ( φ ), so that the the termexpressing radiation damping must be π R ππ k B ( φ ) k dφ in I . We are therefore authorized to callthis term the ’scattered flux’. Now to continue to reason in terms of energy, it is pertinent to ask:if ’lost energy’, is the sum of absorbed and radiation damping energies, then what is the ’providedenergy’, if we accept the fact that energy must be conserved, i.e., ’lost energy’=’provided energy’ ?In our problem, the means by which energy is provided is obviously via the incident plane wave,or in other terms: the means by which flux is provided is via the incident plane wave and the onlyterm in our ’conservation of flux’ relation that can account for this is 4 ℜ (cid:2) a i ∗ B (3 π/ θ i ) (cid:3) in I .Since this term also contains a quantity related to the scattered field (via B (3 π/ θ i ) ) it is morecoherent to normalize the expression of the conservation of flux law I − J = 0 by dividing it by − ℜ (cid:2) a i ∗ B (3 π/ θ i ) (cid:3) so as to obtain − π R ππ k B ( φ ) k dφ ℜ (cid:2) a i ∗ B (3 π/ θ i ) (cid:3) + µ [1] µ [0] ℑ R Γ b u [1] ∗ ( x ) ν · ∇ u [1] ( x ) d Γ4 ℜ (cid:2) a i ∗ B (3 π/ θ i ) (cid:3) = 1 , (57)which can be written as S + A = I , (58)wherein S = − R ππ k B ( φ ) k dφπ ℜ (cid:2) a i ∗ B (3 π/ θ i ) (cid:3) = −JK , (59)is the so-called ’normalized scattered flux’, A = µ [1] µ [0] ℑ R Γ b u [1] ∗ ( x ) ν · ∇ u [1] ( x ) d Γ2 ℜ (cid:2) a i ∗ B (3 π/ θ i ) (cid:3) = LK , (60)is the so-called ’normalized absorbed flux’, and I = 1 . (61)is the so-called ’normalized incident flux’. In black-box language, we can say the the conservationlaw expresses the fact that the ’input flux’ I equals the ’output flux’ S + A , the latter being thesum of the scattered flux S and the ’absorbed flux’ A , wherein, for convenience we have droppedthe term ’normalized’ which is henceforth implicit.Note that (58)-(61) are in agreement with the conservation of flux relation previously obtainedin [41] for the case of a basin filled with a lossy or non-lossy medium.Note also that in the case the protuberance media are non-lossy, the conservation law does notdepend explicitly on any of the constitutive parameters of the filler, i.e., c [ j ] = c ′ [ j ] ; j = 1 , µ [ j ] ; j = 1 ,
2. However it does depend implicitly on these parameters via B ( φ ).12 Demonstration that the conservation of flux relation is satisfiedby the formal solution of the scattering problem when the pro-tuberance is of rectangular shape and is composed of two lossyor lossless media
From now on, the option is to completely solve the forward scattering problem. Fig.2 describesthe scattering configuration in which the bilayer protuberance outer boundary (in the sagittalplane) is rectangular. As previously, the width of the protuberance is w , and its other characteristicFigure 2: Sagittal plane view of the 2D rectangular protuberance scattering configuration. Notethat now the boundary Γ p of the above-ground feature is composed of three connected portions,Γ g , Γ s and Γ d .dimensions are the bottom ( h ) and top ( h ) layer thicknesses, with h = h + h being the heightof the protuberance. What was formerly Γ p is now Γ g ∪ Γ s ∪ Γ d , wherein Γ g is the leftmost verticalsegment of height h , Γ s is the top segment of width w and Γ d is the rightmost vertical segment ofheight h . Everything else is as in fig.1. Owing to the fact that the configuration comprises three distinct regions, each in which theelastic parameters are constants as a function of the space variables, it is opportune to employ13omain decomposition and separation of variables (DD-SOV). Thus, as previously, we decomposethe total field u as: u ( x ) = u [ l ] ( x ) ; ∀ x ∈ Ω l , l = 0 , , , (62)with the understanding that these fields satisfy the 2D SH frequency domain elastic wave equation(i.e., Helmholtz equation) and u [0] satisfies the radiation conditionThe stress-free nature of the boundaries Γ g , Γ d , Γ s , Γ l and Γ r entail the boundary conditions: µ [0] u [0] ,y ( x ) = 0 ; ∀ x ∈ Γ l + Γ r , (63) µ [ l ] u [ l ] ,x ( x ) = 0 ; ∀ x ∈ Γ g + Γ d , l = 1 , , (64) µ [1] u [1] ,y ( x ) = 0 ; ∀ x ∈ Γ s . (65)Finally, the fact that Γ and Γ m were assumed to be interfaces across which two media are inwelded contact, entails the continuity conditions: u [0] ( x ) − u [1] ( x ) = 0 ; ∀ x ∈ Γ m , (66) µ [0] u [0] ,y ( x ) − µ [1] u [1] ,y ( x ) = 0 ; ∀ x ∈ Γ m , (67) u [1] ( x ) − u [2] ( x ) = 0 ; ∀ x ∈ Γ , (68) µ [1] u [1] ,y ( x ) − µ [2] u [2] ,y ( x ) = 0 ; ∀ x ∈ Γ . (69) As in the case of arbitrarily-shaped protuberances, the SOV technique gives rise to the fieldrepresentation u [0] ( x ) = u i ( x ) + u r ( x ) + u s ( x ) (70)wherein (in a somewhat simpler notation) u s ( x ) = Z ∞−∞ B ( k x ) f ( k x , x ) exp( − ik y y ) dk x k [0] y , (71)with: f ( k x , x ) = exp( ik x x ) . (72)Note that the scattered field u s is expressed as a sum of plane waves, some of which arepropagative (for real k [0] y ) and the others evanescent (for imaginary k [0] y ).Within the rectangular protuberance, the same SOV technique, together with the boundaryconditions (64)-65), give rise to the field representations: u [1] ( x ) = ∞ X m =0 h a m exp (cid:0) ik [1] ym y (cid:1) + b m exp (cid:0) − ik [1] ym y (cid:1)i f m ( x ) , (73) u [2] ( x ) = ∞ X m =0 d m cos (cid:2) k [2] ym ( y − h ) (cid:3) g m ( x ) , (74)14herein: f m ( x ) = exp( ik xm x ) , g m ( x ) = cos[ k xm ( x + w/ , (75) k xm = mπw , k [ l ] ym = q(cid:0) k [ l ] (cid:1) − (cid:0) k xm (cid:1) , ℜ k [ l ] ym ≥ , ℑ k [ l ] ym ≥ ω ≥ , l = 1 , . (76)The functions f ( k x , x ) and g m ( x ) satisfy the orthogonality relations:12 π Z ∞−∞ f ( k x , x ) f ( − K x , x ) dx = δ ( k x − K x ) ; K x ∈ R , , (77)1 w Z w/ − w/ g m ( x ) g l ( x ) dx = δ lm ǫ l ; l = 0 , , , ... , , (78)in which δ ( k x − K x ) is the Dirac delta distribution, δ lm the Kronecker delta symbol, and ǫ l theNeumann symbol. The boundary condition (63) and the continuity condition (67) entail µ [0] π Z ∞−∞ u [0] ,y ( x, f ( − K x , x ) dx = µ [1] π Z w/ − w/ u [1] ,y ( x, f ( − K x , x ) dx ; ∀ K x ∈ R , (79)whereas the the continuity conditions (66), (68) and (69) give rise to:1 w Z w/ − w/ u [0] ( x, g l ( x ) dx = 1 w Z w/ − w/ u [1] ( x, g l (( x ) dx ; l = 0 , , , .... , (80)1 w Z w/ − w/ u [1] ( x, h ) g l ( x ) dx = 1 w Z w/ − w/ u [2] ( x, h ) g l (( x ) dx ; l = 0 , , , .... , (81) µ [1] w Z w/ − w/ u [1] ,y ( x, h ) g l ( x ) dx = µ [2] w Z w/ − w/ u [2] ,y ( x, h ) g l (( x ) dx ; l = 0 , , , .... . (82)The introduction of the field SOV representations into these four expressions lead, after use is madeof the orthogonality relations (77)-(78), to: a l + b l = 2 a i ǫ l I + l ( k ix ) + ǫ l Z ∞−∞ B ( k x ) I + l ( k x ) dk x k [0] ; l = 0 , , , .... , (83) B ( K x ) = − w π µ [1] µ [0] ∞ X m =0 ( a m − b m ) k [1] ym I − m ( K x ) ; ∀ K x ∈ R , (84) a l exp( ik [1] yl h ) + b l exp( − ik [1] yl h ) = d l cos( k [2 yl h ) ; l = 0 , , , .... , (85) a l exp( ik [1] yl h ) − b l exp( − ik [1] yl h ) = d l i µ [2] µ [1] k [2] yl k [1] l sin( k [2 yl h ) ; l = 0 , , , .... , (86)15n which K [1] y = q(cid:0) k [1] (cid:1) − (cid:0) K x (cid:1) and: I ± m ( k x ) = Z w/ − w/ exp( ± ik x x ) cos[ k xm ( x + w/ dxw , (87)and it is easy to show that I ± m ( k x ) = i m (cid:2) ( ± k x + k xm ) w (cid:3) + ( − i ) m (cid:2) ( ± k x − k xm ) w (cid:3) . (88)We thus have at our disposal four coupled expressions (i.e., (94)-(86) which should enable us todetermine the four sets of unknowns {B ( k x ) } , { a m } , { b m } , { d m } . Eqs. (85)-(86) readily yield: a l = d l exp( − ik [1] yl h )2 iµ [1] k [1] yl h iµ [1] k [1] yl cos( k [2 yl h ) + µ [2] k [2] yl sin( k [2 yl h ) i ; l = 0 , , , .... , (89) b l = d l exp( ik [1] yl h )2 iµ [1] k [1] yl h iµ [1] k [1] yl cos( k [2 yl h ) − µ [2] k [2] yl sin( k [2 yl h ) i ; l = 0 , , , .... , (90)whence a l + b l = d l κ l , a l − b l = − id l σ l ; l = 0 , , , .... , (91)with κ l = cos( k [1] yl h ) cos( k [2] yl h ) − µ [2] k [2] yl µ [1] k [1] yl sin( k [1] yl h ) sin( k [2] yl h ) , (92) σ l = sin( k [1] yl h ) cos( k [2] yl h ) + µ [2] k [2] yl µ [1] k [1] yl cos( k [1] yl h ) sin( k [2] yl h ) . (93)Finally,the introduction of (91) into (94)-84)results in: d l κ l ǫ l = 2 a i I + l ( k ix ) + Z ∞−∞ B ( k x ) I + l ( k x ) dk x k [0] ; l = 0 , , , .... , (94) B ( k x ) = iw π µ [1] µ [0] ∞ X m =0 d ( M ) m σ m k [1] ym I − m ( k x ) ; ∀ k x ∈ R . (95)16 .6 Approximations of the sets of equations Until now everything has been rigorous provided the equations in the statement of the boundary-value problem are accepted as the true expression of what is involved in the seismic response ofthe protuberance. In order to actually solve for the sets { d m } , and then for { a m } , { b m } , { d m } , {B ( k x ) } (each of whose populations were considered to be infinite until now) we must now resortto approximations.The approach is basically to replace (79)-(86) by the finite system of linear equations a ( M ) l + b ( M ) l = 2 a i ǫ l I + l ( k ix ) + ǫ l Z ∞−∞ B ( M ) ( k x ) I + l ( k x ) dk x k [0] ; l = 0 , , , ...M. , (96) B ( M ) ( K x ) = − w π µ [1] µ [0] M X m =0 ( a ( M ) m − b ( M ) m ) k [1] ym I − m ( K x ) ; ∀ K x ∈ R , (97) a ( M ) l exp( ik [1] yl h ) + b ( M ) l exp( − ik [1] yl h ) = d ( M ) l cos( k [2 yl h ) ; l = 0 , , , ....M , (98) a ( M ) l exp( ik [1] yl h ) − b ( M ) l exp( − ik [1] yl h ) = d ( M ) l i µ [2] µ [1] k [2] yl k [1] l sin( k [2 yl h ) ; l = 0 , , , ....M , (99)from which we obtain, as by the previous steps, a ( M ) l + b ( M ) l = d ( M ) l κ l , a ( M ) l − b ( M ) l = − id ( M ) l σ l ; l = 0 , , , ....M , (100)and B ( M ) ( k x ) = iw π µ [1] µ [0] M X m =0 d ( M ) m σ m k [1] ym I − m ( k x ) ; ∀ k x ∈ R . (101) d ( M ) l κ l ǫ l = 2 a i I + l ( k ix ) + Z ∞−∞ B ( M ) ( k x ) I + l ( k x ) dk x k [0] ; l = 0 , , , ....M , (102)in which the superscript ( M ) signifies the M -th order approximation of the indicated quantity, theprocedure being to increase M so as to generate the sequence of solutions for M = 0, M = 1, etc.until the values of the first few members of these sets stabilize and the remaining members becomevery small.It is important to underline the fact that the approximate solutions { a ( M ) m } , { b ( M ) m } , { d ( M ) m } , {B ( M ) ( k x ) } , together with the associated approximate field representations u s ( M ) ( x ) = Z ∞−∞ B ( M ) ( k x ) f ( k x , x ) exp( − ik y y ) dk x k [0] y , (103) u [1]( M ) ( x ) = M X m =0 h a ( M ) m exp (cid:0) ik [1] ym y (cid:1) + b ( M ) m exp (cid:0) − ik [1] ym y (cid:1)i f m ( x ) , (104) u [2]( M ) ( x ) = M X m =0 d ( M ) m cos (cid:2) k [2] ym ( y − h ) (cid:3) g m ( x ) , (105)17atisfy all the conditions of the boundary-value problem (i.e., Helmholtz equations, radiationscondition, boundary conditions and continuity conditions) for every M ≥
0, but the associatedfield representations are mathematically not complete which is the reason why these solutionsare qualified as approximate. We show hereafter that, in spite of this incomplete nature of ourapproximate solutions, the latter satisfy the conservation of flux relation for all M ≥ On account of the last statement, the demonstration of the conservation of flux law for all M ≥ J ( M ) + K ( M ) = L ( M ) ; ∀ M = 0 , , , ... , (106)is satisfied.Let us begin with L ( M ) , which, due to the continuity relations across Γ m , is L ( M ) = µ [1] µ [0] ℑ Z w/ − w/ u [1]( M ) ∗ ( x, u [1]( M ) ,y ( x, dx = − ℑ Z w/ − w/ u [0]( M ) ∗ ( x, u [0]( M ) ,y ( x, dx . (107)The key point is to invoke the boundary conditions on Γ l and Γ r so as to obtain L ( M ) = − − ℑ Z ∞−∞ u [0]( M ) ∗ ( x, u [0]( M ) ,y ( x, dx . (108)Eq. (103) informs us that: u [0]( M ) ∗ ( x,
0) = 2 u i ∗ ( x, u s ( M ) ∗ ( x,
0) = Z ∞−∞ h a i ∗ δ ( k x − k ix )+ B ( M ) ∗ ( k x ) k [0] ∗ y i exp( − ik x x ) dk x , (109)and u [0]( M ) ,y ( x,
0) = u s ( M ) ,y ( x,
0) = − i Z ∞−∞ B ( M ) ( k x ) exp( ik x x ) dk x , (110)whence, after the interchange of integrals, L ( M ) = − ℑ Z ∞−∞ dk x h a i ∗ δ ( k x − k ix ) + B ( M ) ∗ ( k x ) k [0] ∗ y i Z ∞−∞ dk ′ x B ( M ) ( k x ) ik [0] y Z ∞−∞ dx exp[ i ( k x − k ′ ( x )) x ] . (111)However, Z ∞−∞ exp[ i ( k x − k ′ ( x )) x ] = 2 πδ ( k x − k ′ x ) dx . (112)whence, by use of the sifting property of the Dirac delta distribution, L ( M ) = ℜ " πa i ∗ B ( M ) ( k ix ) + π Z ∞−∞ kB ( M ) ( k x ) k dk x k [0] ∗ y . (113)With the change of variables k x = k [0] cos φ and k ix = k [0] cos φ i , and making use of previously-evoked definitions B ( k [0] cos φ ) = B ( φ ) π , and B ( k [0] cos φ i ) = B ( φ i ) π = B ( θ i + π ) π , we get L ( M ) = 2 ℜ a i ∗ B ( M ) (cid:18) θ i + 3 π (cid:19) + 1 π Z ππ k B ( M ) ( φ ) k dφ , (114)18nd since, by definition, 2 ℜ a i ∗ B ( M ) (cid:0) θ i + π (cid:1) = K ( M ) and π R ππ k B ( M ) ( φ ) k dφ = J ( M ) , we con-clude that L ( M ) = J ( M ) + K ( M ) ; ∀ M = 0 , , , ... , (115)which means that our formulation, which involves only the continuity relations across Γ m and theplane wave representation of u [0]( M ) , is such as to obey the conservation of flux law for all M ≥ L ( M ) = − µ [1] µ [0] ℑ Z w/ − w/ u [1]( M ) ∗ ( x, u [1]( M ) ,y ( x, dx . (116)Eq. (105) informs us that u [1] ∗ ( M ) ( x,
0) = M X m =0 h a ( M ) ∗ m + b ( M ) ∗ m i f m ( x ) , (117)and u [1]( M ) ,y ( x,
0) = M X m =0 ik [1] ym h a ( M ) ∗ m − b ( M ) ∗ m i f m ( x ) , (118)which, on account of (100), give rise to u [1] ∗ ( M ) ( x,
0) = M X m =0 d ( M ) ∗ m κ ∗ m f m ( x ) , (119)and u [1]( M ) ,y ( x,
0) = M X m =0 d ( M ) m k [1] ym σ m f m ( x ) . (120)It follows, after sum and integral changes, that L ( M ) = − µ [1] µ [0] ℑ M X m =0 d ( M ) ∗ m κ ∗ m M X n =0 d ( M ) n k [1] yn σ n Z w/ − w/ dxf m ( x ) f n , (121)or, on account of the orthogonality relation relative to f m , L ( M ) = − µ [1] µ [0] ℑ M X m =0 d ( M ) ∗ m κ ∗ m M X n =0 d ( M ) n k [1] yn σ n w δ mn ǫ m , (122)which, because of the sifting property of the Kronecker delta, becomes L ( M ) = − w µ [1] µ [0] ℑ M X m =0 k d ( M ) ∗ m k κ ∗ m k [1] ym σ m ǫ m . (123)Now, by virtue of the continuity of displacement across Γ m and the boundary condition on Γ l andΓ r , we get Z ∞−∞ h u i ( x,
0) + u s ( M ) ( x, i exp( − ik x x ) dx = Z w/ − w/ u [1]( M ) ( x,
0) exp( − ik x x ) dx ; ∀ k x ∈ R , (124)19hich, by virtue of (109), (117), (119), (112), and the sifting property of the Dirac delta distribution,takes the form 4 πa i δ ( k ix − k x ) + 2 πk [0] y B ( M ) ( k x ) = w M X m =0 d ( M ) m κ m I − m ( k x ) ; ∀ k x ∈ R , (125)from which it follows that Z ∞−∞ kB ( M ) ( k x ) k dk x k [0] ∗ y = − a i ∗ Z ∞−∞ B ( M ) ( k x ) δ ( k ix − k x ) dk x + iw π Z ∞−∞ dk x B ( M ) ( k x ) M X m =0 w π d ( M ) ∗ m κ ∗ m I + m ( k x ) , (126)or on account of (101) Z ∞−∞ kB ( M ) ( k x ) k dk x k [0] ∗ y = − a i ∗ B ( M ) ( k ix ) + iw (2 π ) µ [1] µ [0] M X m =0 d ( M ) ∗ m κ ∗ m M X n =0 d ( M ) ∗ n σ n k [1] n Z ∞−∞ dk x I + m ( k x ) I − n ( k x ) . (127)Owing to the orthogonality relations satisfied by f ( k x , x ) and f m ( x ) we find Z ∞−∞ I + m ( k x ) I − n ( k x ) dk x = 2 πw δ mn ǫ m , (128)so that Z ∞−∞ kB ( M ) ( k x ) k dk x k [0] ∗ y = − a i ∗ B ( M ) ( k ix ) + iw π µ [1] µ [0] M X m =0 k d ( M ) m k κ ∗ m σ n k [1] n ǫ m . (129)Taking the real part of this expression finally yields π Z k [0] − k [0] kB ( M ) ( k x ) k dk x k [0] y + π ℜ h a i ∗ B ( M ) ( k ix ) i = − w µ [1] µ [0] ℑ M X m =0 k d ( M ) m k κ ∗ m σ n k [1] n ǫ m . (130)which, on account of (123) and the change of variables k x = k [0] cos φ again yields1 π Z ππ k B ( M ) ( φ ) k dφ + 2 ℜ (cid:20) a i ∗ B ( M ) ( θ i + 3 π (cid:21) = L ( M ) , (131)wherein we recognize that 2 ℜ (cid:2) a i ∗ B ( M ) ( θ i + π ) (cid:3) = K ( M ) and π R ππ k B ( M ) ( φ ) k dφ = J ( M ) , so thatonce again we obtain J ( M ) + K ( M ) = L ( M ) ; ∀ M = 0 , , , ... , (132)thus showing that the M -th order approximate solutions indeed satisfy the conservation of flux lawfor all M ≥ M ≥
0, and (ii) is valid both inthe absence of losses (in which case L ( M ) = 0) and the presence of losses (in which case L ( M ) = 0).20 Verifications of the conservation of flux law by means of the nu-merical solutions for the case of a lossy or lossless protuberanceof rectangular shape
Until now, we treated the the forward-scattering problem in a formal manner, i.e., based on theformal (i.e., not numerical) solution of the four equations (96)-(99). We now address the problem ofhow to actually (i.e., numerically) solve this system of four equations and then employ the numericalsolutions to see if they are such as to (numerically) satisfy the conservation of flux relation.
We already showed that this system can be reduced to two coupled systems of equations (101)-(103), which we re-write here for convenience: B ( M ) ( k x ) = iw π µ [1] µ [0] M X m =0 d ( M ) m σ m k [1] ym I − m ( k x ) ; ∀ k x ∈ R . (133) d ( M ) l κ l ǫ l = 2 a i I + l ( k ix ) + Z ∞−∞ B ( M ) ( k x ) I + l ( k x ) dk x k [0] ; l = 0 , , , ....M , (134)The final trick is to introduce(133) into (134) so as to give rise to the system of linear equationsfor the { d l } M X m =0 E ( M ) lm d ( M ) m = c l ; l = 0 , , , ...M , (135)in which: E ( M ) lm = δ lm κ l ǫ l − iw π µ [1] µ [0] k [1] ym σ m J lm , c l = 2 a i I + l ( k ix ) , J lm = Z ∞−∞ I + l ( k x ) I − m ( k x ) dk x k [0] y . (136)As it stands, the matrix equation E ( M ) d ( M ) = c is not particularly-appropriate for the determi-nation of the diffraction coefficient vector d ( M ) . The reason for this is that certain elements of thematrix E ( M ) become very large for large M so as to make the intversion of E ( M>> problematic.The way to resolve this problem is actually quite simple: in (135), divide E lm by σ m andmultiply d m by σ m so as to obtain ∞ X m =0 E ( M ) lm F ( M ) m = G l ; l = 0 , , , ... , (137)in which: E lm = E lm ǫ l σ m = δ lm κ l σ l − ǫ l iw π µ [1] µ [0] k [1] ym J lm , G l = c l ǫ l = 2 a i ǫ l I + l ( k ix ) , F m = d ( M ) m σ m . (138)From the numerical point of view, there is now no difficulty in solving for F m via (137).The fulll numerical procedure is then to increase M so as to generate the sequence of numericalsolutions { F (0)0 } , { F (1)0 , F (1)1 } ,....until the values of the first few members of these sets stabilize and21he remaining members become very small. This is usually obtained for reasonably-small values of M , especially in the low frequency regime of interest in our seismic response problem.The problem that was ignored until now is that of providing a suitable means for numericallyevaluating J lm . This problem be treated either in the manners explained in ([39, 40]) or as follows.From the definition (87) of I ± lm and that (39) of the Hankel function, it follows, after changes inthe order of integration, that J lm = πw Z w/ − w/ dx ′ cos[ k xl ( x ′ + w/ Z w/ − w/ dx cos[ k xm ( x + w/ H (1)0 ( k [0] | x ′ − x | ) , (139)this double integral (over finite limits) being accessible to standard (e.g., Simpson) quadraturetechniques. All the following numerical examples apply to the case of a small (in height) double-layer hillsubmitted to an obliquely-incident seismic plane body wave: a i = 1 ( a.u. ), θ i = 70 ◦ , h = 75 m , h = 75 m , w = 750 m , µ [0] = 6 . GP a , β [0] = 1629 . ms − , β [1] = 1300 − i ms − , µ [2] = 2 GP a , β [2] = 1000 − i ms − . Thus, the only changes from one example to another concern M , µ [1] or f .22 .2.1 Effect of the variation of M for fixed frequency and shear modulus in Ω Fig. 3 is relative to a variation of M . S A O ( b ) ,I (r) Figure 3: The three panels depict fluxes as a function of M . The upper, middle and lower panelsare relative to the normalized scattered flux S ( M ) , the normalized absorbed flux A ( M ) and thenormalized output flux O ( M ) = S ( M ) + A ( M ) respectively. The blue curves are all the result ofnumerical solutions and the red line in the lower panel depicts the normalized input flux I = 1which is the goal for O if the conservation law O = I is to be satisfied. Case f = 2 Hz and µ [1] = 4 GP a S ( M ) and A ( M ) are seen (at least graphically) to stabilize starting with M = 2, O ( M ) is seen to depart progressively (although slightly) from I = 1, the reason for this beingthe numerical errors in the computation of J lm (because when the latter quantities are computedwith more accuracy it was found that the difference of O ( M ) from I = 1 decreases). In spite of this,the conservation law is seen to be satisfied with an error of less than a half percent.An important feature of fig. 3 is that the conservation law is numerically satisfied for all M even though the lower-order (i.e., M <
2) solutions are in obvious error (this meaning that theyare not stabilized, assuming that the correct values of S [ M ] and A [ M ] are their stabilized (large M )values). f for various M and µ [1] Figs. 4-6 are relative to a variations of f . S A O ( b ) ,I (r) mu2=(2000000000,2000000000,1) c2r=(1000,1000,1) c2i=(−10,−10,1) M=(0,0,1) N=(400,400,1)J=(201,201)K=(101,101) ||Bi||=0.19446 Figure 4: The three panels depict fluxes as a function of f . The upper, middle and lower panelsare relative to the normalized scattered flux S ( M ) , the normalized absorbed flux A ( M ) , and thenormalized output flux O ( M ) = S ( M ) + A ( M ) respectively. The blue curves are all the result ofnumerical solutions and the red line in the lower panel depicts the normalized input flux I = 1which is the goal for O if the conservation law O = I is to be satisfied. Case M = 0, µ [1] = 4 GP a .24 .5 2 2.5 30.70.80.91 f (Hz) S A O ( b ) ,I (r) Figure 5: Same as fig. 4 except that now M = 5 and µ [1] = 4 GP a . S A O ( b ) ,I (r) Figure 6: Same as fig. 4 except that now M = 5 and µ [1] = 2 GP a .25n these three examples, the conservation law is again seen to be satisfied with an error of lessthan a half percent for all M . The comparison of fig. 4 (relative to the M = 0 solution) and fig. 5(relative to the M = 5 stabilized, presumably correct, solution) shows that the M = 0 solution isacceptable only at very low frequencies f < . Hz . This latter finding is bad news for those whooffer approximate solutions similar to our M = 0 solution as explanations of how an above-groundstructure responds, over a rather large range of frequencies, to a seismic wave.
10 Conclusions