## Abstract

Coupled-mode equations describing the propagation and scattering of internal waves over large-amplitude arbitrary topography in a two-dimensional stratified fluid are derived. They consist of a simple set of ordinary differential equations describing the evolution of modal amplitudes, based on an orthogonality condition that allows one to distinguish leftward- and rightward-propagating modes. The coupling terms expressing exchange of energy between modes are given in an analytical form using perturbation theory. This allows the derivation of a weak-topography approximate solution, generalizing previous linear solutions for a barotropic forcing that were described in 2002 by Llewellyn Smith and Young . In addition, the orthogonality condition derived is valid for a different set of eigenmodes defined on a sloping bottom, which shows a better convergence rate when compared with the standard set of modes. The work presented here provides a useful and simple framework for the investigation of internal wave propagation in an inhomogeneous ocean, along with theoretical insight.

## 1. Introduction

Ocean internal waves and internal tides have been identified as important constituents of the oceanic circulation (e.g., Wunsch and Ferrari 2004). They play an important role in the dissipation of energy and the mixing of water masses, with a major impact, as yet not fully understood, on the meridional overturning circulation (MOC; Kunze 2017). For instance, internal tides generated by barotropic tidal currents over topographic features are thought to account for 1 TW of energy dissipation in the deep ocean (Egbert and Ray 2001) of the 2.1 TW required to support the MOC (Munk and Wunsch 1998; Wunsch and Ferrari 2004).

One of the difficulties in understanding tidal-driven mixing and dissipation lies in their propagation, thus transporting and redistributing energy over large distances, with large variations in the distance traveled and hence large uncertainties in the localization of the ultimate breaking events, if any. Low-mode (modes 1 and 2) internal waves can propagate over very long distances, while higher modes generally dissipate and break locally. In addition, internal tides interact with their environment (topography, varying stratification and mesoscale currents), resulting in scattering and reflection, thus adding variability to their propagation (Garrett and Kunze 2007).

The conversion of the barotropic tide has been investigated by many authors (Cox and Sandstrom 1962; Baines 1973; Bell 1975a,b), and theoretical estimates, based on a small-topographic-amplitude approximation among others, have been formulated by Llewellyn Smith and Young (2002) and St. Laurent and Garrett (2002). These have been used to elaborate global maps of the estimated internal tide generation (Nycander 2005; Melet et al. 2013; Falahat et al. 2014), with reasonable accuracy at least in the deep ocean (Garrett and Kunze 2007). The subsequent propagation and scattering of the internal tide are less well documented, especially with regard to the theoretical aspects. The propagation and diffraction of internal wave beams were investigated by Rainville and Pinkel (2006), and the scattering was studied—again in the limit of weak topography—by Müller and Xu (1992), but questions remain, particularly in regions of large-amplitude topography such as coastal shelves (Kelly et al. 2012; Klymak et al. 2016).

Theoretical models addressing propagation for finite-amplitude topography can be classified into two groups, often formulated in a two-dimensional domain. The first uses the method of characteristics and conformal mapping to build solutions for isolated topography (Müller and Liu 2000; Bühler and Holmes-Cerfon 2011; Maas 2011). It is applicable to horizontally homogeneous stratification and a limited class of vertical profiles, although these can be general enough to apply to the ocean. Another approach is based on a vertical mode decomposition of the solution, which may be convenient because these modes propagate independently in a horizontally homogeneous medium over a flat bottom, with known horizontal wavenumbers and group velocities. In addition, such a decomposition is particularly relevant in the context of surface (barotropic)-internal-tide decomposition, and the generation of internal tides, where modal decomposition is a natural choice and has been used either for theoretical predictions (Llewellyn Smith and Young 2002) or for diagnosing observations or numerical simulations (e.g., Zhao et al. 2019). It was used by Griffiths and Grimshaw (2007) to investigate the generation of internal tide at continental shelves, a formulation reused later on by Griffiths (2011) to propose a global model of the internal tides. Maugé and Gerkema (2008) used a similar formulation in a more complex (and three-dimensional) domain, including background flow and weakly nonlinear effects, yet under the approximation of small topographic height. Kelly et al. (2012) derived an energy equation to diagnose the internal tide scattering on a continental slope, and Shimizu (2011) proposed a slightly different formulation based on the linearized shallow water equations with multiple layers. An intermediate method, introduced by Pétrélis et al. (2006) in the context of the barotropic tide conversion and applied more recently to the scattering of internal tide (Mathur et al. 2014), is based on representing the solution using a Green’s function and solving the resulting integral equation along the bottom boundary.

Using a vertical-mode decomposition in a horizontally inhomogeneous medium results in nonconstant modal amplitudes as well as nonzero horizontal derivatives of the eigenfunctions. The latter have to be evaluated numerically to be able to solve the equations that determine the former. Kelly et al. (2013) introduced a semianalytical two-dimensional model that confines the coupling at the interface of horizontal steps with homogeneous medium properties (equivalent depth), by matching the solutions between adjacent steps. This approach is a generalization of earlier work by Sjöberg and Stigebrandt (1992) where each step’s contribution was added independently (thus discarding backward reflection and interferences). Within this formulation, horizontal derivatives of the modes are replaced by integrals representing the product of mode functions between adjacent steps.

In this paper, we present a method based on modal decomposition that relies on perturbation theory to formulate a set of “coupled-mode equations” describing the evolution of modal amplitudes. In this approach, quantities *F*(*x*, *y*, *z*) are represented in the form

where the modes *F*_{m} are defined by one-dimensional eigenproblems with “frozen” properties; that is, horizontal variables enter only parametrically so that the eigenproblems are local with respect to *x* and *y*. The quantities *A*_{m}(*x*) are range-dependent modal amplitudes the horizontal evolution of which is obtained. This method is known as the “reference waveguide method” (Brekhovskikh and Godin 1992) or “coupled-mode theory” (Desaubies and Dysthe 1995, and references therein). The result is a set of coupled-mode equations which contains coupling coefficients accounting for exchanges between modes. Obtaining the coupled-mode equations requires extensive integration by parts to account for the possibly different boundary conditions satisfied by *F* and *F*_{m}, as will be shown.

The general approach works in three dimensions with variable bathymetry, rotation, compressibility, stratification, and background velocity fields. In this paper we limit ourselves to the two-dimensional case for an incompressible ocean with rotation and vertical stratification *N*(*z*) so that the background state is at rest, so as to be able to focus on the details of the method and compare with previous results. We derive a set of eigenmodes and corresponding orthogonality condition that make explicit the distinction between counterpropagating modes, allowing us to project the initial linearized equation onto the basis of vertical modes in a concise way. In addition, the orthogonality condition enables us to use a different set of vertical modes that match the bottom boundary condition with finite-amplitude topography, thereby exhibiting more rapid convergence with *n* compared to the standard modes.

The use of perturbation theory allows one to recast the horizontal derivatives of modes as a function of the values of the modes themselves and the gradient of the topography. This manipulation is exact provided the topography is a differentiable function of *x*. Using a WKBJ approximation for the vertical modes, we can further derive analytical expression for the coupling between modes for both amplitude and energy flux. Injecting these coupling terms in the coupled-mode equations, using a weak-topography approximation, leads to analytical solution for the scattering of an incoming wave with arbitrary modal constituents.

The paper is organized as follows. Section 2 presents the linearized equations and the vertical modes used for projection. Some of the properties of these modes that are useful for the subsequent calculations, in particular results from perturbation theory, are described. Then, in section 3, the coupled-mode equations are derived and discussed in the context of energy exchange between modes. Section 4 examines two examples of internal wave scattering over topography, including solutions based on the coupled-mode equations. In section 5, we introduce a different set of vertical modes that have the boundary slope embedded in the bottom boundary condition, and discuss changes to the results that follow from this new set. Conclusions are given in section 6, and the appendix contains some calculations whose results are needed in the main text.

## 2. Problem formulation and normal mode decomposition

### a. Governing equations and boundary conditions

We consider the linearized equations for a stratified fluid in rotation under the Boussinesq approximation. We restrict ourselves to perturbations that are independent of one of the horizontal coordinates, thus considering a two-dimensional domain bounded by a rigid lid at *z* = 0 and bottom at *z* = −*H*(*x*). The equations are formulated in the frequency domain:

where *N*(*z*) is the buoyancy frequency and standard notation is used for the variables (e.g., Vallis 2006). The following boundary conditions hold:

The horizontal boundary conditions depend on the details of the configuration investigated, and will be specified later on.

The problem can be recast in terms of the velocities *u* and *w*, the variables present in the vertical boundary conditions. The equations become

with *α*(*z*) = (*N*^{2} − *ω*^{2})/(*ω*^{2} − *f*^{2}) being the square of the inverse characteristic slope, which has no horizontal variations in the problem considered.

We introduce the “state vector” **Ψ** = (*u*, *w*)^{T}, which allows us to rewrite the problem in the following matrix form:

with

### b. Vertical normal modes

To apply a modal decomposition to the linear problem (8), we first need to define a complete set of modes. We consider the range-independent linear equations [i.e., with frozen properties with respect to *x* so that *H*(*x*) = *H* is treated as a constant] and apply a Fourier transform along the *x* direction. We obtain the following eigenvalue problem, with *ϕ*_{n} and *φ*_{n} being the vertical functions for *u* and *w*, respectively:

with eigenvalue *k*_{n}. Here and in the rest of the paper, prime denotes a vertical derivative. This set of equations is usually formulated for one variable or the other; that is,

As above we introduce the state vector for the normal modes **Φ**_{n} = (*ϕ*_{n}, *φ*_{n})^{T} that will be used below.

The boundary conditions corresponding to the rigid lid and flat bottom are *φ*(0) = *φ*(−*H*) = 0, implying *ϕ*′(0) = *ϕ*′(−*H*) = 0. With this choice of boundary conditions, both equations in (12) lead to an eigenproblem with real eigenvalue $kn2$. Then the eigenfunctions *ϕ*_{n} and *φ*_{n} are associated with the following orthogonality conditions:

where *j*_{n} is a real arbitrary normalization constant and *δ*_{n,m} is the usual Kronecker symbol. For simplicity, and without loss of generality, we will set *j*_{n} = 1 in the rest of the paper.

The right-hand sides in the orthogonality conditions (13) and (14) reflect the fact that counterpropagating modes with identical vertical structure are degenerate. However, one can take a combination of the two conditions to distinguish between the modes. We thus define modes corresponding to the positive and negative roots by

With this notation, positive and negative roots are associated with modes propagating toward increasing and decreasing *x*, respectively, and a change of sign of *φ*. Unlike the usual convention, “mode” in this paper will thus refer to the state vector **Φ**_{n}; that is, counter-propagating modes have different structures, although they are composed of similar functions (*ϕ*_{±n}, *φ*_{±n}). A “generalized” orthogonality condition that makes the distinction between right-going and left-going modes can be formed by taking $\u222b\u2212H0\u2061[\u2061(\varphi n\phi m)\u2032\u2212\u2061(\varphi m\phi n)\u2032]\u2009dz$. Using (10) and (11), this gives

By virtue of the boundary conditions, the LHS vanishes. Anticipating further calculation, we introduce the matrix

and the notation $\u2329\cdots \u232a\u2261\u222b\u2212H0\cdots \u2009dz$ for conciseness. The orthogonality condition is thus

where $\Phi n\u2020=\u2061(\varphi n*,\phi n*)=\u2061(\varphi n,\u2212\phi n)$ is the conjugate transpose of **Φ**_{n} = (*ϕ*_{n}, *φ*_{n})^{T}, since choosing *j*_{n} > 0 implies that the vertical modes *ϕ*_{n} are real while *φ*_{n} are imaginary.

### c. Perturbation analysis and normal mode properties

When using modal analysis in an inhomogeneous domain, one uses locally defined normal modes to project the equations. This procedure results in terms that contain horizontal derivatives of the eigenmodes (see the next section and the appendix), for which no analytical expression exists in general, thus requiring numerical evaluation. Here, we use perturbation theory (Courant and Hilbert 1966) to express these derivatives as a function of the eigenmodes themselves and the horizontal derivative of inhomogeneities of the medium. Interestingly, this method has been routinely used in the framework of ocean and atmospheric acoustics (e.g., Brekhovskikh and Godin 1992), but has not been used for internal wave propagation previously, to our knowledge.

The perturbation theory is applied as follows: consider the eigenproblem (10) and (11), in two media with slightly different properties: *H*_{n} = *H*, *H*_{m} = *H* + *δH* and *α*_{n} = *α*, *α*_{m} = *α* + *δα*. At this stage, the physical origin of the perturbations (e.g., horizontal inhomogeneities) does not matter. We have introduced here a perturbation in *α* for generality and further computation, although we will not use it to include variations in *N*^{2} because these are not addressed in this paper and would require a different set of linear equations. Perturbations of *α* may instead be caused by variations in *ω*, and this is used below to obtain an expression for the modal group velocity. As a result of these perturbations in *α*, we have *k*_{m} → *k*_{m} + *δk*_{m}, *ϕ*_{m} → *ϕ*_{m} + *δϕ*_{m}, and *φ*_{m} → *φ*_{m} + *δφ*_{m}. Taking the first variation of (15) gives

The LHS vanishes at the upper boundary but not at the bottom. Expanding the lower boundary condition *φ*_{m}(−*H*_{m}) = 0 leads to $\delta \phi m\u2061(\u2212H)=\phi m\u2032\u2061(\u2212H)\delta H=\u2212ikm\varphi m\u2061(\u2212H)\delta H$. Hence

Since *N*^{2} does not vary with *x* in this paper, we first take *δα* = 0 and obtain the following equation:

We can introduce the physical origin of the perturbation [here, *δH* = (*dH*/*dx*)*δx*] and take the limit of infinitely small perturbation to construct the horizontal derivative of the eigenmodes. This calculation is exact provided that *H*(*x*) is continuous and differentiable. Setting *n* = *m* in (20) gives the following expression for the horizontal derivative of the eigenvalues:

Furthermore, combining (17) and (20) for *n* ≠ *m* leads to an expression for the horizontal derivative of the eigenvectors:

Note that, while we restrict ourselves to variations of the topography along a single dimension here, the use of perturbation theory in 3D configurations (although using a different orthogonality condition, since right-/left-propagating mode separation would not hold) and for variations of other physical quantities, such as properties of the medium that would enter the eigenproblem as parameters, is straightforward.

This method also allows one to compute the modal group velocity: for this purpose, we repeat the above procedure with variable *ω* (and, thus, variable *α*) and *δH* = 0, to obtain

where the result

## 3. Coupled-mode equations

### a. Derivation of the coupled-mode equations

We now derive the coupled-mode equations that describe the evolution of the modal amplitude. To do so, we assume that normal modes actually depend on the horizontal coordinates in a parametric way (e.g., Desaubies and Dysthe 1995; Griffiths and Grimshaw 2007) by expanding the variables as

which is the appropriate version of (1). At each location, **Φ**_{m} are defined by the eigenproblem (10) and (11) with frozen (i.e., the local value of) *H*. The initial system of equations in (8) is projected onto the basis of vertical modes, with the help of the equalities for the horizontal variations of these modes [(21) and (22)] obtained as a result of perturbation theory developed above. This method is known as the “reference waveguide method” (Brekhovskikh and Godin 1992).

Projection of the equations onto mode *n* is carried out by multiplying (8) by $\Phi n\u2020P$ and integrating over [−*H*(*x*), 0]. Since Ψ and the **Φ**_{m}s do not obey the same condition at *z* = −*H*, one must integrate by parts and insert the modal expansion at an intermediate stage of the procedure (e.g., Llewellyn Smith and Young 2002; Kelly 2016). The procedure is detailed in the appendix. One finally obtains a set of first-order ordinary differential equations (ODEs), the so-called coupled-mode equations, in the form

This system of equations describes the evolution of the modal amplitudes over range. The right hand side terms account for energy exchange between modes expressed by the coupling coefficients *g*_{nm}. These are given by

These expressions require the use of the results in the appendix along with (22). With the normalization used in the paper, *ϕ*_{n}*ϕ*_{m} scales like *H*^{−1} and *g*_{nm} scales like *H*^{−1}*dH*/*dx*.

The barotropic mode, with *n* = 0, is unaffected by the internal modes, since *k*_{0} = 0 because of the rigid lid condition. Thus, it can act only as a forcing generating internal modes, while *g*_{00} ensures mass conservation, that is, that ⟨*u*_{0}⟩ is constant with respect to *x* (noting that $\varphi 0=1/H$ and *φ*_{0} = 0).

The set of equations in (26) combined with the normal mode definition (10) and (11) and the coupling coefficients (27) are a formally *exact* reformulation of the initial problem (8) using modal decomposition, provided that the depth profile *H*(*x*) is continuous and differentiable. It has the advantage of reducing the horizontal dependence of the solution to a simple system of ODEs where the distinction between leftward and rightward propagating mode appears explicitly. A similar approach was used in Griffiths and Grimshaw (2007) and in Kelly et al. (2016) in a more general setting (time-domain formulation, and varying stratification in a three-dimensional domain and background currents), in the form of coupled shallow-water equations. A closer formulation in a two-dimensional domain is the coupling equations for linear tides (CELT) model (Kelly et al. 2013), in which the horizontal domain is divided into sections with horizontally constant properties (*N*, *H*, and *f*) and the solution is constructed by applying matching conditions at the interfaces.

The coupled-mode equations formulation derived in this paper is different in several ways: first, it is continuous with respect to *x*. Therefore, one cannot derive the coupling terms for discontinuous topography (since perturbation theory does not apply exactly). However, the ODEs in (26) can be solved separately on each side of the discontinuity and the full solution can be obtained by applying dynamical conditions (on the pressure and the velocity) at the interface, following, for example, Robinson (1969). When doing so, the coupled-mode equations yield a linear matrix problem, akin to stepwise formulations as in the CELT model.

In addition, the horizontal derivatives of the eigenfunctions and eigenvalues are recast in terms of the variations of *H*, thus providing an analytical expression for the mode coupling. Note that this was partially done in the derivation by Griffiths and Grimshaw (2007), in a formulation that was valid for variations of *H* only (through the chain rule), whereas perturbation theory is applicable to variations of other properties [e.g., *N*(*z*, *x*), *f*] as well. We mention, however, that this could imply a different set of initial linear equations and additional complexity: for instance, nonhomogeneous stratification is associated with geostrophic currents and additional interaction (advective) terms.

### b. Analytical expression for mode coupling and energy exchange

We now give expressions for the horizontal energy flux and the energy exchange between modes based on the coupled-mode equations derived previously. The horizontal flux of energy *I* is given by

To obtain an expression for the modal energy flux, we now derive an expression for the modal projection of the pressure. From the vertical component of (2) and (3), we have

The procedure is now similar to other derivations given in this paper. To summarize, we multiply both sides of the above equality by *φ*_{n}, expand $w=\u2211mAm\phi m$, use the second orthogonality condition (14), and integrate the left-hand side by parts, expanding the pressure as $p=\u2211mpm\varphi m$, with boundary terms vanishing. Then, using (10) and the first orthogonality condition (13), we obtain

Note that leftward and rightward modes are degenerate in the orthogonality conditions (13) and (14) and are thus combined in (30). However, inserting the modal decomposition for pressure and velocity from (28) into (30), we obtain (after reordering the terms in the sum and using *k*_{−n} = −*k*_{n})

in which one recognizes the contribution of fluxes associated with left- (*n* < 0) and right-going (*n* > 0) waves. From this expression we can define the modal energy flux for mode *n* as

Note that the barotropic energy flux is infinite according to this expression, which is a well-known limitation of the rigid-lid approximation. As mentioned previously, however, the barotropic mode only satisfies mass conservation and is an adiabatic reservoir of energy flux for the internal tides (i.e., forcing) in this limit.

From this expression, one sees that the exchange of energy between modes depends on the variations of the depth of fluid and on the values of the normal modes at the bottom (because *g*_{nm} does). Moreover, the (*k*_{n} − *k*_{m})^{−1} dependence shows that adjacent and copropagating modes are more prone to exchange energy with each other than counterpropagating modes.

We can obtain more useful information about the strength of the coupling (and, thereby, mode scattering) by using analytical expressions for the eigenmodes. We thus use the WKBJ approximation in the eigenproblem (10) and (11), which gives (e.g., Gill 1982)

where $mn=kn\alpha $ and $\beta \u2061(z)=H\alpha /\u2329\alpha \u232a$. The eigenmode is associated with the eigenvalue

We recall, for clarity, that *k*_{n}, *m*_{n}, and *β* are also functions of *x* (whereas *α* is not).

Injecting this form into the expression for the coupling coefficient (27), we obtain

This expression can be substituted into (33) to provide an analytical expression for the modal energy exchange:

## 4. Solutions to the coupled-mode equations

The set of coupled equations in (26), even when analytical expressions such as (36) are available from exact (or approximate WKBJ) solutions for the vertical modes, has no general analytical solution. In this section, we first derive approximate solutions valid in the limit of weak topographic coupling. Then numerical solutions of the coupled-mode equations are presented. These are used to illustrate the method presented in this paper and also to provide a reference for evaluating the accuracy of the approximate solution. We should emphasize that none of the solutions (approximate or numerical) here relies on using WKBJ approximation for the vertical modes. The horizontal (i.e., coupled-mode equations) and vertical (normal mode) solutions are distinct problems. The WKBJ approximation is only necessary if a fully analytical solution is desired (this will be discussed explicitly in what follows).

We first recast the coupled-mode equations by using a “phase compensated” amplitude:

where we have introduced the phase function $\theta n\u2061(x)=\u222bxkn\u2061(x\u2032)\u2009dx\u2032$. Inserting this into (26) gives

where $\theta nm=\u222bx\u2061[kn\u2061(x\u2032)\u2212km\u2061(x\u2032)]\u2009dx\u2032$. Hence *B*_{n} remains constant if *g*_{nm} = 0; that is, changes in the variables *B* are associated with scattering. In the following, we will solve this new set of ODEs instead of (26).

Specifying the horizontal boundary conditions leads to a boundary-value problem. One can specify the value of all modes, that is, incoming and outgoing at a physical boundary, or implement a radiation condition by specifying *B*_{n} = 0 at the entrance boundary for each mode: the left and right boundary for *n* > 0 and *n* < 0, respectively.

### a. Approximate solutions for weak coupling

As shown in section 3b, the strength of the coupling scales like the relative topographic variations and the inverse mode number difference. This can be used to derive approximate solutions when the coupling is expected to remain weak.

#### 1) Adiabatic approximation

The so-called adiabatic approximation describes the evolution of modes when coupling terms, and thereby mode scattering, are neglected (Milder 1969; Brekhovskikh and Godin 1992; Desaubies and Dysthe 1995). In this case, energy flux conservation is ensured by retaining the diagonal terms *g*_{nn}. The solution to (39) is then

where the exponential term accounts for the normalization of the vertical modes to ensure energy flux conservation. This solution is equivalent to the WKBJ solution of the depth-separated linear equations for baroclinic motions (in two dimensions) of Rainville and Pinkel (2006). Using in addition the WKBJ approximation for the vertical modes introduced above, (34), and assuming that *β*(*H*) is nearly constant (this is exactly true for stratification profiles with *α* = [*a* exp(*Hz*) + *b*]^{2}, with *a* and *b* constant), the adiabatic solution can be further simplified to

Under the adiabatic approximation, the modal amplitudes adjust to variations of the depth of fluid so that the energy flux is conserved mode-wise. This solution is not of great interest in the context of oceanic internal wave scattering, since scattering is explicitly neglected.

#### 2) Weak topography approximation

Let us introduce the following asymptotic expansion: $Bn=\u2211j\epsilon jBn\u2061(j)$, where the small parameter *ε* is a nondimensional measure of the topographic amplitude *δH*/*H*. Since *g*_{nm} scales with *H*^{−1}*dH*/*dx*, the coupled-mode equations in (39) then become

where the zeroth-order amplitude $Bn\u2061(0)$ is constant and set by the boundary conditions. Note that the previous adiabatic solution can be seen as a zeroth-order approximation, in which the diagonal terms in *g*_{nn} are explicitly retained.

Integrating (42) to first order in the expansion gives

To go further, we can also expand the eigenvalues and eigenfunctions in powers of *ε*, and retain only their leading-order, *x*-independent part [e.g., if *H* ≈ *H*_{0} + *εδH*(*x*), then *ϕ*, *φ*, and *k* are function of *H*_{0} only at leading order]. The first-order solution then becomes

where the 0 subscript refers to quantities evaluated at *H* = *H*_{0}. This expression provides an approximate semianalytical solution for the mode amplitude evolution over arbitrary topography. The first term corresponds to the adiabatic correction to the modal amplitude, while the second accounts for energy transfer between modes. The solution obtained is a generalization (in a one-dimensional horizontal domain) of the linear, weak topography solution of Llewellyn Smith and Young (2002) for the conversion of the barotropic tide. Here, any incident barotropic or baroclinic mode with amplitude $Bn\u2061(0)$ can act as a “forcing” tide, while the amplitude of the scattered waves $Bn\u2061(1)$ are first order. It is worth mentioning that this first-order solution is accurate as long as the amplitude of $Bn\u2061(1)$ remains small, so that the small-parameter expansion remains valid. Since the small parameter is *ε* ∝ *δH*/*H*, this depends on the horizontal extension of the topographic irregularities (as the amplitude of the first-order solution can grow and accumulate over large distances) in addition to their magnitude. Using the WKBJ approximation for the vertical modes, as in section 1, it is possible to derive analytical solutions for a variety of profiles *H*(*x*) and boundary conditions.

### b. Numerical solution: Method and example

For large-amplitude topography, the approximate solution derived above is likely to be inaccurate, and one is left with the full boundary value problem (39) to be solved numerically. To illustrate the method described in this paper, we consider an idealized configuration in which a first baroclinic mode, coming from the left, encounters a ridge. The latter is modeled by a cosine bump with

for |*x* − *x*_{0}| < *L*_{h} and *H*(*x*) = *H*_{0} otherwise, where *H*_{0} is the depth away from the ridge and 2*L*_{H} the ridge width. We use two different values of the ridge amplitude while keeping the far-field depth *H*_{0} = 3 and slope criticality max(*dH*/*dx*)(*α*)^{1/2} = 0.8 constant. The width of the topography changes accordingly. We take *α* = 1 (the results can be rescaled to realistic values by changing the aspect ratio of the domain and the amplitude of the velocities accordingly). A numerical solution is obtained and compared with the approximate solution derived above. It will be further compared with an exact 2D numerical solution in the next section.

Instead of solving the BVP (39) directly, we reduce the numerical complexity by turning the problem into a set of iterative initial value problems, using the fact that coupling terms between counterpropagating modes are weaker than for copropagating modes. The following set of equations is obtained:

At each iteration, the above set of coupled ODEs is solved in both directions separately, with the contribution from counterpropagating modes taken from the previous iteration.

Boundary conditions are specified on the zeroth-order amplitude $Bn\u2061(0)$ while $Bn\u2061(i>0)=0$ is imposed at the “entrance” boundary (i.e., left and right boundary for positive and negative *n*, respectively). Note that the barotropic amplitude is given by the adiabatic solution and thus $B0\u2061(i>0)=0$, because it is not altered by the baroclinic modes as shown in section 3. The numerical solution is computed using a finite number of modes *N* in each direction and a finite number of iterations. When reconstructing the physical solution, we only use the first 3*N*/4 modes because the highest modes resolved are inaccurate: indeed, they should exchange energy with higher modes, but these latter modes are not present due to truncation, so that the corresponding coupling terms are missing. This is similar to the aliasing effect due to nonlinearities in pseudospectral methods. Since coupling between modes is most important for closest modes, the lower modes are only weakly affected by this truncation effect.

The solution to (46) is computed using 50 vertical modes (in each direction) for a “moderate amplitude” ridge Δ*H*/*H* = 1/6 with $B1\u2061(0)=1$ on the left boundary, thus forcing the incident first baroclinic mode. Only a few iterations (two here) are necessary to obtain an accurate approximation to the solution, which is shown in Fig. 1. One sees the formation of beams emanating from the topography in both directions: 14% of the incoming flux (contained in the first baroclinic mode) is forward-scattered into higher copropagating modes, while 5% is reflected (i.e., backscattered). As can be seen in the right panel of Fig. 1, the weak-topography solution (44) provides a reasonable approximation to the numerical solution, predicting overall 15% of energy flux forward scattered and 4% reflected. This agreement is lost for higher values of the amplitude of the topography, as expected. The solution for Δ*H*/*H* = 1/2 is shown in Fig. 2. In this case, much more of the incoming energy flux is forward scattered into higher copropagating modes (80%) which thus carry most of the energy at the right boundary (nearly 4 times as much as the remaining first baroclinic mode). In contrast the amplitude of the reflected beam is much weaker and accounts for only 0.4% of the incoming energy flux. By comparison, the weak-topography approximate solution predicts 125% and 0.01% forward scattered and backscattered, respectively. It greatly overestimates the generation of copropagating low-order modes (mode 2) while failing to predict large-enough amplitudes for the higher counterpropagating modes.

## 5. Normal modes over a sloping bottom

### a. Definition of the sloping-bottom modes

The set of normal modes *ϕ*_{n} and *φ*_{n}, defined in section 2b and used in the subsequent derivation and solution of the coupled-mode equations, forms a complete basis that is uniformly convergent over [−*H*, 0] (Jackson 1914, his theorem 1). Because the boundary conditions are different from those of the solution to (8), the convergence rate of the basis scales as 1/*N*, where *N* is the total number of modes. In this section, we introduce a different set of vertical modes with modified boundary conditions that incorporate the bottom slope. This decomposition can be used in place of the previous one, although the coupling coefficients are different. We discuss the implication and potential utility of this different decomposition, its range of validity and the corresponding modal coupling coefficients.

Here, one can recognize the “natural” boundary condition (5) by identifying *γ* with the local value of *dH*/*dx*. The corresponding eigenproblem is solved either numerically using a collocation method with Chebyshev polynomials or analytically for particular stratifications (e.g., constant *α* in the examples treated in this paper). The corresponding solution is shown in Fig. 3, compared with the standard, flat-bottom modes defined in section 2b. Note that in the limit of vanishing slope (*γ* → 0), both solutions are identical. Despite the different boundary condition, the bracket in the LHS of (15) still vanishes, so the solutions still obey the orthogonality condition (17) and define a complete basis suitable for modal expansion. However, the modified eigenproblem is not in standard Sturm–Liouville form. In particular, the eigenmodes as well as the eigenvalues now take complex values.

### b. Properties of the modal expansion

This new basis may be useful for several reasons, which are discussed here, among with some noteworthy properties of the “sloping bottom” normal modes. First, the convergence of the corresponding series is expected to be more rapid with *n* than for the “flat bottom” set, because the boundary condition at the bottom is the same for the basis function and the solution. Figure 4 shows the rate of convergence computed on the basis of an exact solution in a wedge with *α* = 1 and two different slopes. The solution (shown in Fig. 4, upper-left panel) corresponds to a first baroclinic “wedge” mode propagating upslope in the apex (Wunsch 1969):

with *q* = 2*nπ*/log[(*α*^{−1} + *γ*)/(*α*^{−1} − *γ*)] and *γ* being the bottom slope. The convergence of the series is measured as the normalized absolute value of the error for *u* and *w* given a number of modes *N* used in the projection (Fig. 4, upper-right panel). (The latter is performed in a standard way, numerically evaluating $An=\u2329\Phi nP\Psi wedge\u232a$ using Clenshaw–Curtis quadrature.) For the flat-bottom mode decomposition, *w* exhibits the highest error and lowest convergence rates: slightly slower that 1/*N*. On the contrary, the sloping-bottom mode decomposition exhibits higher convergence rate, close to 1/*N*^{2}, especially for highest bottom slope, and a significant decrease of the error associated with *w*. This is further illustrated in Fig. 4 (lower-left panel) by the magnitude of the modal amplitude resulting from the projection on both type of eigenmodes, which shows that higher modes have much weaker amplitude for the sloping-bottom modes than for the flat-bottom modes.

While the orthogonality condition (17) is common to the two sets of basis function, it is worth noting that (13) and (14) are not verified by modes on a sloping bottom. Instead we have

This means that the vertically integrated modal energy density (or horizontal flux), given by similar integrals, has no meaning within this decomposition. In other words, sloping-bottom modes are intrinsically coupled with each other, in the sense that one cannot define an energy per mode unless *dH*/*dx* = 0 (i.e., the bottom is flat).

To gain further insight in the structure of these modes, we use the WKBJ approximation. Then, incorporating the boundary condition at the top, the solution takes the form

with $m\u2261k\alpha $, and *ϕ*_{n} follows from (11). Note that this solution is exact for some stratification profiles (Baines 1973). The quantities *m* and *k* take discrete value as a result of the boundary condition at the bottom, which becomes

This equality shows that the eigenvalue with the sloping bottom condition is equal to the flat-bottom eigenvalue plus an additional imaginary part. This imaginary part is associated with enhanced amplitude of the modal structure near the bottom (cf. Figs. 3 and 4, lower-right panels). This behavior is consistent with the physics of internal wave beam reflection over a slope, which is known to exhibit amplification of the energy density near the slope as the slope get higher, in particular when the angle becomes close to criticality.

A major caveat of the sloping-bottom boundary condition is that it can give rise to singularities. While it is beyond the scope of this manuscript to derive a complete spectral theory for the eigenproblem with sloping-bottom boundary condition, this singular behavior can be discussed using the above WKBJ solution. The relation (50) has an inverse hyperbolic tangent singularity [with ℑ(m) → ∞] for $\alpha \u2061(\u2212H)\gamma =1$, which coincides with the well-known criticality condition of internal waves on a slope. This is associated with increasing confinement of the eigenmode structure, with very large amplitude in the vicinity of the bottom (due to the increasing value of the imaginary part of *m*) for $\alpha \u2061(\u2212H)\gamma \u21921$. Away from criticality, the WKBJ solution is well defined. This criticality issue is absent when using the standard modes, as the bottom slope does not enter their definition.

### c. Coupling coefficients and coupled-mode equations for sloping-bottom modes

We now list the counterparts to the formulas established in the text in the presence of a bottom slope. The perturbation analysis now includes a change in the slope, *γ*_{n} = *γ* and *γ*_{m} = *γ* + *δγ*. We deduce from (47)

and

The projection of the linear system of (8) onto the new vertical modes is carried out in the same fashion as over a flat bottom, considering a frozen value of *γ* = *dH*/*dx* in addition to *H* itself. In this case, however, the boundary conditions for expansion functions and independent variables are the same, so there is no need to use integration by parts. The counterpart to (A3) is

After some calculation, one obtains the same coupled-mode equations that are in (26) as before, but with the following coupling coefficients:

### d. Comparison of the solutions using flat-bottom and sloping-bottom modes

A comparison of the solution to the coupled-mode equations and the direct two-dimensional reconstructed fields is shown in Fig. 5. It is further compared to a numerical solution of the 2D initial problem (8) using a collocation method with Chebyshev polynomials in both directions (with 80 collocations points in the vertical and 60 in the horizontal), after recasting the problem in a rectangular domain using a stretched vertical coordinate. The solution is computed with the same topography (with Δ*H*/*H* = 1/6) and in the same domain as in Fig. 1 but, for simplicity, we consider instead an initial value problem in which a first baroclinic mode is coming from the left and no waves exit the domain through the left boundary (i.e., *A*_{n} = 0, *n* ≠ 1). This implies that weak-amplitude (left-propagating) modes are entering the domain at the right boundary and vanish above the topography through interactions with each other and right-propagating modes. Both methods yield solutions very similar to the numerical solution. As previously shown in Fig. 4, the solution above the topography is better reproduced with sloping-bottom mode decomposition, especially near the bottom. Yet, the modal amplitudes at the right boundary (where *dH*/*dx* = 0) are very similar with both types of eigenmodes and compare very well with the amplitudes computed from the projection of the numerical solution onto the set of eigenmodes, up to *n* ≈ ±20. Thus, while using sloping-bottom modes has the advantage of better reconstructing the dynamical fields above nonflat topography, especially near the bottom, they do not provide significant improvement on the modal amplitudes.

## 6. Conclusions

We have derived a representation of linear internal waves propagating in a domain with nonuniform, arbitrary bathymetry, based on a normal mode expansion and coupled-mode equations. It constitutes a simpler, yet exact, reformulation of the initial set of equations, provided variations in range of the environment are smooth (i.e., differentiable). While we have restricted the analysis to variations of bottom depth, the method can be generalized to horizontal variations of any background fields, provided an orthogonality condition can be found and the background state is stationary.

This formulation relies on a complete basis of vertical normal modes. This basis is defined locally, and the propagation and interaction between modes are expressed using variational calculus. This representation has several advantages. The coupled-mode equations express the modal amplitude evolution in a straightforward manner, making explicit the distinction between modes propagating in opposite directions and displaying the coupling term between modes. Using perturbation theory, we derived an analytical expression for the latter, from which an equation for the energy exchange between modes, that is, wave scattering, follows. Our results show that the strength of the energy exchanges scales like the relative topography variations (*H*^{−1}*dH*/*dx*) and the inverse mode number difference. Besides providing a useful framework for the analytical investigation of internal wave scattering, it also saves us from evaluating the horizontal gradients of modal structures numerically, and hence can reduce the numerical cost of the method.

We obtained an approximate solution for the scattering of incoming internal tide by isolated topography by solving the coupled-mode equations in the limit of weak topography. This solution is an extension to baroclinic incoming tide to the solution derived by Llewellyn Smith and Young (2002) in the context of barotropic tide conversion.

Last, we introduced a different vertical mode decomposition that includes the bottom slope in the boundary condition. The basis functions in this representation have the advantage of satisfying the same boundary conditions as the physical fields. While these vertical modes can exhibit a singularity that seems to coincide with slope criticality, the corresponding series converge faster than the usual vertical modes, and they provide a more accurate reconstruction of the dynamical fields especially near the bottom. A potentially interesting property of this decomposition is the intensification of the mode amplitude near the bottom, which matches with physical observations at near-critical slopes and can be related to processes like boundary streaming. In combination with an appropriate parameterization, it could also provide a way to describe bottom-intensified mixing with no need to capture the smallest vertical wavenumbers generated. However, the energy flux can no longer be decomposed into separate modes.

The coupled-mode formulation derived in this paper can be useful for analytical and numerical investigation of internal tide dynamics, as well as for diagnostics of observations and outputs from numerical models. Future work will address the generalization of the present formulation to nonlinearities, three-dimensional configurations and cases with a background geostrophic flow.

## Acknowledgments

This work was supported by ONR Award N00014-13-1-0347. We are grateful to an anonymous reviewer and to Jonas Nycander for suggestions that helped to improve the paper.

### APPENDIX

#### Series Expansion and Galerkin Integration

Let us consider two functions *f*(*x*, *z*) and *g*(*x*, *z*). We are concerned with series expansions of the form $g=\u2211mAmgm$ in the vertical integrals ⟨*f*∂_{x}*g*⟩ and ⟨*f*∂_{z}*g*⟩. If *g* and *g*_{m} do not satisfy the same boundary conditions, one must perform the integrations with care in a Galerkin fashion to avoid exchanging differentiation and summation. This is done using Leibniz’s rule or integrations by parts back and forth, with the series expansion performed at an intermediate stage. We find

For the terms explicitly involved in the derivation of the coupled-mode equations, this gives for the projection of the LHS of (8):

The second term on the RHS is given by (22) for *n* ≠ *m*. For *n* = *m*, we make use of Leibniz’s rule again to obtain

Last, for the projection of the RHS of (8), one has

## REFERENCES

_{2}internal tide field

## Footnotes

© 2020 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).