## Abstract

A layer of homogeneous, isotropic, elastic material overlays a substrate of similar material. The shear wavespeed within the layer is less than that of the substrate causing waves to be trapped within the layer. At the interface a long inclusion, that grows gradually until it reaches a constant thickness, is introduced. The inclusion is composed of a material whose shear wavespeed is less than that in the layer; it is described as slow. It is imagined that the lowest surface-wave mode of the structure is incident to the growing inclusion. Numerical calculations show that the growth of the slow inclusion brings the wavenumber of this lowest mode into an interval where it is close to that of the second mode, thus exciting it. This process is repeated when the wavenumber of the second mode is brought close to that of the third. Within these intervals, energy is exchanged among the coupling modes. Outside of these localized intervals, the modes propagate independently of one another and their amplitudes vary such that the flux of energy in each mode is conserved; they are said to propagate adiabatically. Reflections are also excited, but are shown to be very small in magnitude.

## 1. Introduction

The goals of the present paper are to derive a system of coupled local-mode equations to describe coupling among elastic surface waves perturbed by a slow, interfacial inclusion, such as that shown in figure 1; and to use asymptotic and numerical calculations to approximate their solution so that the coupling among the modes can be explicitly studied. The lowest surface-wave mode is imagined to be incident to the inclusion from the uniform region, from the left, in figure 1; it is subsequently coupled to various higher modes, in the region to the right, by the presence of an inclusion with a very slow shear wavespeed (it is described as slow) and a variable thickness. The curves and *x*_{2}=0 bound the inclusion, and the parameter *ϵ*, where *ϵ*≪1, controls its rate of growth. Interest in problems of this kind arose from an effort to model a potential method to detect ultrasonically a defect in the adhesion of a thin film to its substrate. For this reason the shear wavespeed of the inclusion was taken as much smaller than those of the surrounding materials: it was imagined to be composed of weak, perhaps unconsolidated, material.

The local modes are determined by choosing a reference waveguide. In this paper it is a set of uniform layers whose thicknesses match those of the varying structure, in figure 1, at each cross-section; thus, for an inclusion at the interface, a two layered waveguide on a substrate with a rigid base is used as the reference for each cross-section A–A. The modes for this reference are described as *local*. These modes are calculated using the Thomson–Haskell, or propagator matrix, technique (Aki & Richards 2002, pp. 261–282).

The propagation of guided waves in a slowly changing, horizontally inhomogeneous environment is frequently studied using coupled modes: It is used to study acoustic propagation underwater (Brekhovskikh & Godin 1999, pp. 243–282), the propagation and scattering of electromagnetic surface waves (Shevchenko 1971, pp. 22–142), and the propagation of elastic waves trapped in slow surface layers (Malischewsky 1987, pp. 54–69; Kennett 1984). More specifically, coupled local modes are used for elastic-wave problems by Odom (1986), Maupin (1988) and Tromp (1994). Odom uses local modes to describe the propagation of antiplane shear waves in a surface waveguide of variable thickness, and includes the continuous spectrum in his coupled-mode theory. Maupin (1988) gives a formal development of the coupled local-mode approximation for inplane elastic waves, and Tromp (1994) extends her work to waves in three dimensions; neither author works out a particular case in detail. Maupin (1988) introduces the idea of using an equivalent body force to take account of the failure of the local modes to satisfy fully the continuity conditions at sloping boundaries. This difficulty was noted previously by Rutherford & Hawker (1981) in the context of underwater acoustics.

The formulation of the coupled local-mode equations assumes that wavefields are continuous throughout the region of interest. That implies that no singularities in stress can be present. The function is selected to fulfil this condition (Bogy & Wang 1971). If the horizontal inhomogeneity introduces sudden changes, as it might if the inclusion began with a step, then a different approach is needed. Two such approaches are developed by Malischewsky (1974) and by Pagneaux & Maurel (2001).

The theoretical development in this paper is close to that of Maupin (1988); nevertheless, parts of this development are re-examined in this paper. In doing so, the possibility that local modes with complex horizontal wavenumbers contribute to the coupled local-mode equations is noted, though not explored. The presence of a continuous spectrum is avoided by imagining that a rigid base is present. Explicit calculations are confined to the surface-wave modes so that it is simply a formal mechanism used to make the complete set of eigenmodes discrete. Its depth *d* is made very large when the surface-wave modes are numerically calculated. The present work differs from the earlier work of Maupin (1988) in that, by seeking to work through a representative problem, explicit expressions describing the coupling coefficients are calculated. By doing so it becomes apparent (i) that the coupling coefficients may need to be approximated to *O*(*ϵ*^{2}) and (ii) that the coupling coefficients can become very large if the related wavenumbers approximate one another. If the wavenumbers coalesce, the coupled local-mode equations must themselves be re-formulated.

Two works addressing this latter issue are a study of mode coupling between two modes in a slowly tapered earth–ionosphere waveguide (Budden 1975), and a study of linearly coupled oscillators with slowly varying frequencies (Grimshaw & Allen 1979). The work of Budden (1975) is the more relevant: When two wavenumbers coalesce the coupled local-mode equations are re-formulated using combinations of the coupling modes; from this a two-by-two system with a single turning point is isolated; this is reduced to Airy's equation. If there are several intervals where different wavenumbers coalesce, then a new calculation is needed for each interval. In the representative case worked through in this paper, there are two coupling intervals. However, the wavenumbers do not coalesce (for real *x*_{1}), though they come close; as a consequence, in these intervals, though the coupling coefficients become large, they remain well defined. In these intervals the coupled local-mode equations are numerically integrated.

Elsewhere the modes propagate adiabatically; that is, they propagate independently and their amplitudes vary so as to conserve the flux of energy in each mode. In these intervals the coupled local-mode equations are solved using a WKBJ approximation. WKBJ methods are widely used for many problems in wave propagation. Three references, that may be useful to the reader are: Berry & Mount (1972), Steele (1976) and Fedoryuk (1993). Fedoryuk gives an example of a first-order system with turning points (the Stueckelberg system, pp. 346–351).

*Notes on scaling*: the equations of inplane, linear, isotropic elasticity arewhere *α*, *β*, *γ*=1,2 and summation over repeated indices is used. The *T*_{αβ} are the components of stress, *U*_{α} are the components of particle displacement and *G*_{α} are the components of a body force per unit mass. The time-dependence will be taken as exp(−i*ωT*). The material properties are the density *ρ*, and the elastic or Lamé parameters *λ* and *μ*. Within each layer the material properties are constant; however, a dependence on (*x*_{1}, *x*_{2}) may be exhibited to indicate that they change from layer to layer. The material properties *ρ*_{0}, *λ*_{0} and *μ*_{0} characterize the substrate and are constant. The shear modulus *μ*_{0} is used to scale the stress components; thus, *τ*_{αβ}=*T*_{αβ}/*μ*_{0}.

A reference angular frequency *ω*_{0} is used to scale the time. The shear wavenumber *k*_{T0}=*ω*_{0}/*c*_{T0}, where *c*_{T0}=(*μ*_{0}/*ρ*_{0})^{1/2}, is used to scale lengths. The scaled coordinates are *x*_{α}=*k*_{T0}*X*_{α} and the scaled time is *t*=*ω*_{0}*T*. The scaled angular frequency becomes *Ω*=*ω*/*ω*_{0}, and the time-dependence thus becomes exp(−i*Ω**t*), though hereafter it is suppressed. The scaled particle displacement components are *u*_{α}=*k*_{T0}*U*_{α}. The scaled thickness of the inclusion is *h*_{l}=*k*_{T0}*H*_{l}, that of the layer is *h*_{0}=*k*_{T0}*H*_{0} and that of the substrate is *d*=*k*_{T0}D.

The following combinations of material parameters prove helpful in writing the equations(1.1)Note that in the substrate *p*=1 and *m*=1.

Also of use in later work is the relation(1.2)with the notation ∂_{α} meaning ∂/∂*x*_{α}.

## 2. Horizontal equations of motion

### (a) Equations of motion

Expressing the equations of motion in terms of the components of particle displacement results in a system of equations that contains the terms i*β*_{p}∂_{2}, in addition to . This is awkward because *β*_{p} is the eigenvalue. The equations can, instead, be written as a system that isolates the eigenvalue by eliminating the stress component *τ*_{22}, given by equation (1.2), from both the equations of motion and the constitutive equations.

Using the vector(2.1)where the superscript *T* indicates the transpose, the equations of motion are then written as(2.2)where the operator is defined by(2.3)A symbol such as ∂_{2}[*a*(*x*_{1}, *x*_{2}), when operating on *u*_{2}, means ∂_{2}[*a*(*x*_{1}, *x*_{2})*u*_{2}]. * F* is the dimensionless body force and is described bywhere

*f*

_{α}=

*pG*

_{α}/(ω

_{0}

*c*

_{T0}).

Figure 1 indicates the direction of the unit normal at the sloping interface. The components of traction acting on it are(2.4)where . Note the presence of the *τ*_{11} term. The boundary and continuity conditions for the problem are thus(2.5)

The inner product used in this paper is defined as(2.6)The vectors * V* and

*both have the form shown in equation (2.1). The interval of integration for equation (2.6) needs to be broken up if discontinuities arise in the integrand, as they do in calculating the coupling coefficients (see (A 1)).*

**U**### (b) The eigenvalue problem

Set *ϵx*_{1}=const. so that a specific cross-section, such as *A*−*A* in figure 1, is being considered. Next consider a layered structure where the interfacial layer has a uniform thickness *h*_{l}(const.), so that the material properties depend only on *x*_{2}. Thus * U* can be written as

*=*

**U**

**u**_{p}(

*x*

_{2})exp(i

*β*

_{p}

*x*

_{1}). Substitution into equation (2.2), with

*set to zero, gives(2.7)The accompanying boundary and continuity conditions are identical to those stated in equation (2.5) except that the first simplifies to(2.8)Equations (2.7) and (2.8), together with the remaining conditions given by equation (2.5), are the eigenvalue problem giving the eigenmodes*

**F**

**u**_{p}and eigenvalues

*β*

_{p}corresponding to this uniform layering.

Setting * V*=

**u**_{q}exp(i

*β*

_{p}

*x*

_{1}) and

*=*

**U**

**u**_{q}exp(i

*β*

_{q}

*x*

_{1}), substituting into equation (2.6), integrating by parts, and invoking the conditions (2.5) as modified by equation (2.8), gives(2.9)where the superscript * indicates complex conjugation. This is the orthogonality condition. The modes are labelled as described in Auld (1990, pp. 155–160): Modes occur in pairs and are labelled so that one member of each pair carries energy or decays exponentially in the +

*x*

_{1}direction, while the other member does so in the −

*x*

_{1}direction. Subscripts

*p*are integers and take plus and minus values, but not 0; a plus value is taken for propagation in the +

*x*

_{1}direction and a minus is taken for propagation in the −

*x*

_{1}direction. Thus

*β*

_{−p}=−

*β*

_{p}.

Note that the eigenvalues can be complex. This is a characteristic feature of elastic-wave problems: it is known that the dispersion relation both for propagation in an elastic plate with traction-free surfaces (Auld 1990, fig. 10.14), and for propagation in a layer on a rigid base (Kirrmann 1995) each have complex branches. Though the dispersion equation for the modes used in this paper is not explored in sufficient detail to uncover the complex branches, they are likely to be present. The complex eigenvalues are included in the early development of the coupled local-mode equations so that it may be assumed that the eigenmodes form a complete set. Kirrmann (1995) and Besserer & Malischewsky (2004) discuss this issue in greater detail.

By taking the complex conjugate of equation (2.7), and giving close attention to the mode labeling just noted, it can be shown thatMoreover, by asking that the modes of equations (2.7) and (2.8) be *symmetric upon reflection in ϵx*_{1}=const., it can be shown that, if , then(2.10)The modes are normalized so that the conditions of (2.10) are satisfied. One could also ask that the modes be antisymmetric, or a combination of the symmetric and antisymmetric modes. Moreover, *u*_{2} must be assigned a value at *x*_{2}=−*h*_{0}; this is discussed in §4*c*.

If *β*_{p,q} are real and *β*_{q}=*β*_{p}, then either *p*≡*q*, or the eigenmodes **u**_{p} and **u**_{q} are degenerate. In the former case, which is the usual one, in each of the previous statements, beginning with equation (2.9) *q* is replaced by *p*, so that *P*_{pp}≠0, and so on. In the latter case, *P*_{pp} has a very particular behaviour, as was pointed out by Budden & Eve (1975) and Budden (1975) in a study of degenerate modes in the earth-ionosphere waveguide. Assume that *β*_{p}≈*β*_{q}. Recall that the eigenmode **u**_{p}(*x*_{2}) is a function of its eigenvalue: **u**_{p}(*x*_{2}):=* u*(

*x*

_{2},

*β*

_{p}), where

*is the general solution to equation (2.7). For*

**u***β*

_{p}close to

*β*

_{q}, and

*an analytic function of*

**u***β*It then follows that(2.11)where

*P*

_{qp}=0 because

*(*

**u***x*

_{2},

*β*

_{q}) is orthogonal to

*(*

**u***x*

_{2},

*β*

_{p}), provided

*β*

_{q}≠

*β*

_{p}. If

*β*

_{p}→

*β*

_{q}, then

*(*

**u***x*

_{2},

*β*

_{p})→

*(*

**u***x*

_{2},

*β*

_{q}) and the two modes become degenerate. In this limit,

*P*

_{pp}=

*P*

_{qq}=

*P*

_{qp}=0; this is the self-orthogonality noted by Budden and Eve, and Budden, that arises when eigenmodes are degenerate.

### (c) The local eigenmodes and eigenvalues

The thickness of the layering does change as one moves from cross-section to cross-section; *ϵx*_{1} is a parameter identifying a particular uniform layering at each cross-section. Therefore, and *β*_{p}(*ϵx*_{1}) are the *local* eigenmode and corresponding eigenvalue for each , and they are the solution to the eigenvalue problem posed by (2.7) and (2.8). Moreover, *P*_{pq}(*ϵx*_{1}), defined by (2.9), changes as one moves from cross-section to cross-section. The symmetries indicated by (2.10) do not change. In the problem solved in this paper, as *ϵx*_{1} varies, a real is brought to an interval where it approximates a real (*p*≠*q*), and there (2.11) describes the behaviour of *P*_{pp}; however, the two modes *never* become degenerate.

## 3. The coupled local-mode equations

### (a) The sloping interface

Equation (2.5) indicates that the solution to the problem of interest demands that *t*_{1} and *t*_{2} both be continuous across . However, equation (2.8) indicates that the local modes only enforce the continuity of *τ*_{12} and *τ*_{22}; *τ*_{11}, and therefore *t*_{1}, is discontinuous at the sloping interface. Because , where , and . The discontinuity in *t*_{1} is then equal towhere and the + defines the side to which the normal points. However, this discontinuity has an effect upon the solution that is equivalent to that of a body force (−* F*) (Aki & Richards 2002, pp. 38–42). Therefore, to enforce the requirement that

*t*

_{1}be continuous, the equal and opposite, body force

*is introduced so as to cancel the effect of the discontinuity. The force*

**F***is given by(3.1)*

**F**### (b) The evolution equations

At this point it is useful to introduce the slow variable ; it is this variable, rather than *x*_{1}, that governs the horizontal variation. A solution to equation (2.2), with the body force given by equation (3.1), is sought in the form(3.2)where(3.3)and *c*_{0}≡0. The local modes **u**_{r}(*y*_{1}, *x*_{2}) satisfy equation (2.7), at each *y*_{1}, and all but the first condition of (2.5); the first condition is enforced by introducing equation (3.1).

Substituting the sum (3.2) into (2.2), and using (2.7), giveswithUsing equation (2.9), one notes that **u**_{q} is orthogonal to **u**_{r} ∀*r*≠*p*, and that *P*_{qp}≠0 for (no modes are allowed to be degenerate). This gives the coupled local-mode equations(3.4)The interactions among the modes are represented, to , by the coupling coefficients *K*_{qr}; these are given by(3.5)Noting the symmetries given by equation (2.10), it can be shown that(3.6)and, from the definition of *P*_{qr}, that(3.7)The calculation of equation (3.5) follows the same pattern as that found in Odom (1986) for his scalar problem. Appendix A outlines the intermediate steps. The outcome is(3.8)The term [*p*] means (*p*^{+}−*p*^{−}), and the + defines the side to which the normal points; the [*c*], [*a*], [1/*m*] and [*b*/*m*] have a similar meaning. Note that the coupling coefficients depend critically on *f* and , and that the term in square brackets is evaluated along *x*_{2}=*h*_{l}(*y*_{1}).

## 4. A representative inclusion

### (a) The inclusion

Consider the profile of an inclusion that grows smoothly, beginning from the origin, to a thickness Δ*h* over a distance *l*, and thereafter remains a uniform thickness (see figure 1). In §3*a* it was asked that . Set , with *ϵ*≪1. To ensure that no singularities in the stress components arise, select *h*_{l}(*y*_{1}) such that *h*_{l}(0)=0, *f*(0)=0, and *f*(Δ*h*)=0. A profile that satisfies these criteria isThen(4.1)Other profiles could have been picked without affecting the theoretical formulation provided that they do not introduce stress singularities. Lastly recall that the inclusion is composed of a material having a small shear wavespeed.

### (b) Simplifying the coupled, local-mode equations

Mode 1 of the uniform structure on the left, in figure 1, is imagined to be incident to the inclusion after having propagated from a distant source. Then *β*_{1} must be real if mode 1 is not to have a negligible amplitude. It will be assumed that *β*_{1} does not approaches a complex *β*_{p}, *p*≠1 and that it does not become complex, as *y*_{1} varies. Accordingly, modes with complex *β*_{p} are not likely to be excited by the incident mode; therefore, from this point forward only modes with real *β*_{p} are considered. Moreover, only transmitted or reflected modes with amplitudes larger than *O*(*ϵ*^{2}) are considered, because the coupling coefficients are only calculated to *O*(*ϵ*^{2}).

The coupled, local-mode equations can now be put into a simpler form. Equation (3.7) is used to eliminate the diagonal terms from the right-side of equation (3.4). Further, *P*_{pp} can be combined with the *c*_{p} by defining(4.2)and(4.3)Recall from equation (2.11) that, if *β*_{p}≈*β*_{r}, the denominator of equation (4.3) is proportional to (*β*_{p}−*β*_{r}). The function *f*, given by equation (4.1), is embedded in *Q*_{pr}. The coupled, local-mode equation (3.4) now takes the form(4.4)Recall that *q*_{0}≡0. The term *δ*_{rp} is the Kronecker delta. The combination of symbols *r*∈{*β*_{r} real} indicates that the sum is taken only over eigenmodes with real eigenvalues. The term coupling coefficient will continue to label the whole symbol *Q*_{pr}/(*β*_{p}−*β*_{r})^{2}.

Equation (4.4) is not an optimal set of equation to solve numerically because they form a two-point boundary-value problem, not an initial-value one. For the transmitted modes, where *p* is a positive integer, all the *q*_{p} are specified at *y*_{1}=0, whereas for the reflected modes, where *p* is a negative integer, the *q*_{p} are specified at *y*_{1}=Δ*h*.

### (c) A numerical investigation

To understand how the solutions to equation (4.4) might behave, a numerical investigation of the coupling coefficients is undertaken. Numerical values for the various parameters are set, assuming that these are representative. In figure 1 the material of the substrate is labelled 0, that of the inclusion 1 and that of the layer 2. The parameters are given the following valuesTo calculate *a*, *b*, *m* and *c* use the parameterand *c*_{T}=(*μ*/*ρ*)^{1/2}. In all cases *κ*=(3)^{1/2}. Note that . Setting *h*_{0}=1 determines the reference angular frequency *ω*_{0}; for example, if the physical thickness of the layer is *H*_{0}=1.0 mm, then *ω*_{0}=2.5×10^{6} s^{−1}.

Figure 2 shows the normalized slowness *s*_{p}=*c*_{T0}/*c*_{p} graphed against *Ω*, at different points in the interval *y*_{1}∈[0,0.5], for the first four modes. Recall that *β*_{p}=*Ω**s*_{p}. Note that the dispersion curves for modes 1 and 2 approach one another, and subsequently those for modes 2 and 3 do so as well. Figure 2*d* shows the approach of modes 1 and 2 in detail. These graphs suggest that there are intervals in which modes couple strongly because their wavenumbers approximate one another.

To show how the wavenumbers vary with *y*_{1}, for fixed *Ω*, figure 3 shows and for *y*_{1}∈[0,0.5] and *Ω*=10. If the condition defining the coupling intervals is set as , then the coupling interval for modes 1 and 2 lies between *a*_{1}≈0.17 and *a*_{2}≈0.23, and that for modes 2 and 3 between *b*_{1}≈0.26 and *b*_{2}≈0.34. The choice of these intervals is somewhat arbitrary and other reasonable choices are possible, including those that cause the intervals to overlap; however, such choices do not change the underlying arguments.

Next consider the behaviour of the coupling coefficients *Q*_{12}/(*β*_{1}−*β*_{2})^{2} and *Q*_{23}/(*β*_{2}−*β*_{3})^{2}, and their numerators *Q*_{12} and *Q*_{23}, for *y*_{1}∈[0,0.5] and *Ω*=10. *Q*_{pr} is defined by (4.3). The increase of the coupling coefficient *Q*_{12}/(*β*_{1}−*β*_{2})^{2} in the coupling interval [*a*_{1},*a*_{2}] is very marked, whereas *Q*_{12} hardly varies in the same interval. The changes in the second coupling interval [*b*_{1},*b*_{2}] are similar. It is clear then that the variation in the difference of the wavenumbers is the cause of the sudden increase in the coupling coefficients.

*Notes about the numerical work*: The eigenvalues *β*_{r} and the eigenmodes **u**_{r} are calculated at 100 points of a grid, from *y*_{1}=0 to *y*_{1}=Δ*h*, using the propagator matrix technique. This grid remained fixed throughout the calculations. The system of equations that were solved is given in appendix B by equation (B 1); the propagator matrix elements are listed in eqn (7.45) of Aki & Richards (2002, p.273), though they need to be scaled to agree with the non-dimensional variables used here. There are only two layers and a substrate in the *x*_{2} direction, so that only three propagator matrices are needed. The propagator-matrix technique was implemented in Matlab. The eigenvalues *β*_{p} are calculated by the finding zeros of the dispersion equation as a function of *β*, for various values of *Ω* at *y*_{1}=0, 0.2, 0.3, and at other values not shown in figure 2. Fixing *Ω*=10, the eigenvalues, the eigenfunctions and the coupling coefficients were calculated at each point of the *y*_{1} grid.

In doing so, *u*_{p2} is assigned a value 10^{−3} at *x*_{2}=−*h*_{0} for all *p*. Once this is done *u*_{p1} is calculated at *x*_{2}=−*h*_{0}; the stresses *τ*_{p12} and *τ*_{p22} are zero there. Setting *u*_{p2} to be real causes the symmetries indicated in equation (2.10) to be enforced. The eigenfunctions are then calculated for arbitrary *x*_{2} using the propagator matrices, that propagate the values at the surface into the interior.

Recall from equation (4.4) that is the term of interest rather than **u**_{p}; the normalization, at *x*_{2}=−*h*_{0}, taken for the numerical calculations is an intermediate step. Equation (2.11) suggests that *P*_{pp} becomes small when *β*_{p}≈*β*_{q}; however, note that 〈**u**_{p}, **u**_{p}〉=*P*_{pp} and *P*_{pp}≠0.

In the course of the formal development of the coupled, local-mode equations it was convenient to imagine that there is a rigid base at *x*_{2}=*d*. To calculate numerically the surface-wave modes the *d* is set to 0.02*k*_{T0} (*D*=0.02 m). For *ω*_{0}=2.5×10^{6} s^{−1} and *c*_{T0}=2.5×10^{3} m s^{−1}, *d*=20. For *s*_{p}>1 and *Ω*>1, this was sufficient to ensure that any up-going disturbances excited at *x*_{2}=*d* had effectively zero amplitudes, for the surface-wave modes used in these calculations. A grid of 500 points from *x*_{2}=−*h*_{0} to *x*_{2}=*d* was used to calculate *P*_{pp}.

## 5. Asymptotic and numerical approximations

Figures 2–4 indicate that mode 1, incident from the left, couples strongly to mode 2 in [*a*_{1},*a*_{2}], and mode 2 couples strongly to mode 3 in [*b*_{1},*b*_{2}]. Outside of these coupling intervals the modes propagate adiabatically. All the remaining modes, including the reflected ones, are of negligible magnitude. This interpretation is here supported by asymptotically approximating the *q*_{p} in the intervals of adiabatic propagation, and by truncating the coupled local-mode equations in the coupling intervals to a system that includes only the two coupling modes, and then integrating them numerically.

Recall that the *q*_{p} with positive subscripts are the amplitudes of the normalized modes propagating to the right, while those with negative subscripts are those propagating to the left. Let identify the set of all integers, the set of positive integers, and the set of negative integers. A lower case subscript *p* will indicate that .

### (a) The transmitted modes in [0, *a*_{1}), (*a*_{2}, *b*_{1}) and (*b*_{2}, Δ*h*]

In this subsection, a capital subscript *P* indicates that . Consider first the interval [0,*a*_{1}). The *q*_{r} are expanded aswhere *q*_{10}=*q*_{1}(0). Mode 1 is incident from the left so that *q*_{1}(0)≠0, while *q*_{R}(0)=0, *R*≠1. These expansions are substituted into equation (4.4) and the equations for integrated from 0 to *y*_{1}, where *y*_{1}∈[0,*a*_{1}). This gives(5.1)The exponential term, within the integrand, has an argument multiplied by the large parameter indicating that the integral can be asymptotically approximated (Copson 1971, pp. 21, 29–34). Introducing the new variable of integrationthe integral is approximated as(5.2)Recall that *f*(0)=0, so that *Q*_{Pr}(0)=0. This indicates that the phase term in the coupled local-mode equation (4.4), introduces an order *ϵ* term, when the first-order system is integrated, provided |*β*_{P}−*β*_{r}|=*O*(1) and the functions multiplying the exponential term vary smoothly.

Examining equations (5.1) and (5.2), it is determined that(5.3)In deriving this expression the reflected terms have been assumed to be of *O*(*ϵ*). This must be the case because, as *ϵ*→0, there can be no reflections. Section 5*e* gives estimates of the reflected terms. Returning to the original expansion (3.2), it is now possible to write(5.4)This is equivalent to the first term of the WKBJ approximation.

Coupling of modes 1 and 2 occurs in [*a*_{1},*a*_{2}] so that *q*_{2} emerges from the coupling interval being of the same order as *q*_{1}; numerical integration indicates that this is the case. An analysis very similar to that just undertaken shows that(5.5)The remaining .

Coupling of modes 2 and 3 occurs in [*b*_{1},*b*_{2}] so that *q*_{3} emerges from the coupling region being of the same order as *q*_{2}. Again repeating the previous analysis gives(5.6)The remaining . Thus, for *y*_{1}∈(*b*_{2},Δ*h*],(5.7)This, again, is equivalent to a WKBJ approximation. The wavefield (5.4) enters the sloping region and the wavefield (5.7) emerges.

Note that throughout the sloping region, from (5.3), (5.5) and (5.6), truncating the system of equations (4.4) to a three-by-three system for the first three modes will give a system whose solution is accurate to order *ϵ*. Moreover, because *f*(Δ*h*)=0 and hence *Q*_{Pr}(Δ*h*)=0, all , where *P*≠1, 2, 3.

### (b) Coupling of two modes in [*a*_{1},*a*_{2}] and [*b*_{1},*b*_{2}]

In [*a*_{1},*a*_{2}] (*β*_{1}−*β*_{2}) varies rapidly and becomes of order at its minimum, while the remaining *β*_{r} remain separated from *β*_{1,2}: the integral in (5.2) cannot be approximated as indicated.

Definingthe coupled local-mode equations in this first coupling region, truncated to *O*(*ϵ*), may be written as(5.8)whereUsing the symmetries indicated by equation (2.10), note that *Q*_{12}=−*Q*_{21} for *β*_{1,2} real. The matrix * B* is thus Hermitian. The initial condition is determined from equation (5.3) and is given asThis system of equations is numerically integrated.

Equation (5.8) could also be written by setting(5.9)and solving for *q*_{1} and *p*_{2}. The lower limit of the integrations in the exponential terms of * B* would begin at

*a*

_{1}rather than 0.

A very similar two-by-two system of equations describes the coupling of modes 2 and 3 in [*b*_{1},*b*_{2}]. It is given by equation (5.8) if subscripts 1 and 2 are replaced by 2 and 3. The initial conditions for the second coupling region are

### (c) Numerical integration

Figure 5 illustrates the outcomes of the arguments just advanced. The magnitudes |*q*_{p}|, Re(*q*_{p}) and Im(*q*_{p}), for *p*=1, 2, 3 are graphed against *y*_{1} in figure 5*a–c*, respectively. Figure 5*c* also indicates the coupling intervals. The solid curves show a solution composed of equations (5.3), (5.5) and (5.6) in the intervals of adiabatic propagation, and of a numerical solution to equation (5.8) in [*a*_{1},*a*_{2}] and a numerical solution to a similar two-by-two system in [*b*_{1},*b*_{2}]. The dashed curves show a numerical solution to (4.4), taken throughout [0,Δ*h*], truncated to a three-by-three system. The previous analysis indicates that this is also accurate to *O*(*ϵ*). Figure 5*b*,*c* indicate how the phases of the leading order terms change. Writing *q*_{2} as indicated by equation (5.9) shows that propagation terms are part of this term so that the rapid changes in both Re(*q*_{2}) and Im(*q*_{2}) in figure 5*b*,*c* are to be expected. Similar comments apply to *q*_{3}.

The solid curves demonstrate that the coupling in each of the coupling intervals brings a lower order term to the same order as the incident one. Moreover, they indicate that, in similar problems, identifying the coupling intervals with a calculation such as that shown in figure 3, and subsequently numerically integrating the truncated two-by-two system gives a satisfactory solution. The dashed curves indicate that a numerical integration throughout *y*_{1}∈[0,Δ*h*] of the three leading order terms gives a solution that is consistent with the asymptotic analysis and the choice of coupling intervals.

*Notes about the numerical work*: A Matlab routine for a system of stiff differential equations was used in all the numerical integrations. The numerical information needed for this routine was available from the work of §4.

### (d) Energy flux

Returning to (4.4), and using (5.3), it is straightforward to show that(5.10)where *q*_{1}(0)=1. This states that the energy contained in the first three right-propagating eigenmodes is conserved to *O*(*ϵ*^{2}). Looking at figure 5*a* it is seen that this relation is followed to very good accuracy.

The dimensionless, time-averaged flux of energy through a vertical cross-section, such as A−A of figure 1, carried by the *r* th mode iswhere (4.2) defining *q*_{r} has been used. The transmission coefficients for the energy flux, with respect to *q*_{1}(0)=1, are thus defined asThese terms can be calculated from figure 5*a* for the specific case considered here; for example, .

### (e) An estimate of the reflections

A capital subscript *P* indicates that , in contrast to the convention used in §5*a*. The boundary conditions at *y*_{1}=Δ*h* are that *q*_{p}(Δ*h*)=0. Moreover, it is known from §5*a* that , for *p*>3, so that *q*_{p}(Δ*h*), *p*=1, 2, 3, are the leading order terms in the initial conditions for the reflected modes.

Consider first the interval (*b*_{2},Δ*h*]. The *q*_{r} are expanded aswhere *q*_{10}=*q*_{1}(*a*_{2}), *q*_{20}=*q*_{2}(*b*_{2}) and *q*_{30}=*q*_{3}(*b*_{2}). These expansions are substituted into equation (4.4) and the equations for integrated from Δ*h* to *y*_{1}, where *y*_{1}∈(*b*_{2},Δ*h*]. This givesAgain, the feature to note is the presence of the *ϵ*^{−1} multiplying the argument of the exponential term in the integrand. At this point the analysis is identical to that used to arrive at equation (5.2). The integral contributes an *ϵ* multiplying other terms; thus the *q*_{P}(*b*_{2}) are *O*(*ϵ*). The coupling does not increase the order of the reflected terms, though it certainly mixes them, so that the *q*_{P}(*y*_{1}) are *O*(*ϵ*) in *y*_{1}∈[0,*a*_{1}).

## 6. Concluding remarks

Coupled local-mode equations describing propagation in an environment similar to that shown in figure 1 have been investigated. Two limitations to their use have become apparent: (i) To explicitly calculate the coupling coefficients the slope must be *O*(*ϵ*), as indicated by equation (3.1); the equation (3.4) is thus known only to *O*(*ϵ*^{2}). (ii) The equations fail if two eigenmodes become degenerate because the modes become self-orthogonal. Nevertheless, the coupled local-mode equations have proven to be a useful starting point to examine both adiabatic and non-adiabatic behaviour in a slowly changing environment. When two wavenumbers approximate one another, but do not coalesce, strong mode coupling occurs; this can be exhibited by numerically integrating a truncated two-by-two system solely through the coupling interval. Elsewhere the propagation is adiabatic. Lastly, launching a single lowest mode in a thin film, but receiving it along with several higher modes with strong amplitudes may be a useful indication of the de-adhesion of the film from its substrate, and from the travel times of the higher modes the location of the de-adhesion could be estimated.

## Acknowledgments

We thank the NSF (CMS-9812820) for supporting initial parts of the work, John Dempsey for the reference describing the elastic singularities of an embedded wedge, M. V. Perel for the reference describing coupling in a tapered earth-ionosphere waveguide, and a referee for several other references on mode coupling.

## Footnotes

- Received January 16, 2005.
- Accepted June 23, 2005.

- © 2005 The Royal Society