Spectral methods for coastal-trapped waves

Continental Shelf Research 31 (2011) 1481–1489

Contents lists available at ScienceDirect

Continental Shelf Research journal homepage: www.elsevier.com/locate/csr

Research papers

Spectral methods for coastal-trapped waves E.R. Johnson , J.T. Rodney Department of Mathematics, UCL, Gower Street, London WC1E 6BT, UK

a r t i c l e i n f o

abstract

Article history: Received 7 August 2010 Received in revised form 3 June 2011 Accepted 15 June 2011 Available online 23 June 2011

This note describes efficient and accurate spectral numerical schemes to compute both propagating and evanescent free baroclinic coastal-trapped waves over general depth profiles for arbitrary vertical density profiles in horizontally semi-infinite domains. The general problem is recast into a linear eigenvalue problem for the along-shore wavenumber k, which can be solved directly, without initial guesses or searching, using any standard linear eigenvalue package to find real and complex eigenmodes simultaneously. An equivalent recasting gives the linearised eigenvalue problem for the frequency o. A novel, nonlinear, boundary condition is derived that is particularly effective for modes whose offshore decay is weak, as in the long-wave limit. The resulting nonlinear eigenvalue problem is solved by a highly efficient Newton–Kantorovich algorithm. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Coastal trapped waves Evanescent modes Semi-infinite domains Spectral accuracy

1. Introduction Discussions of low-frequency coastal motions often require knowledge of the structure of the local coastal-trapped wavefield. In their discussion of the response of the Gulf of California to remote forcing (Martinez and Allen, 2004) apply the widely used and time-tested FORTRAN code of Brink and Chapman (1987) (BC herein). Webster and Holland (1987), WH herein, describe a method for which they note two advantages over BC: first, the calculation of complex wavenumbers directly, as would result from the determination of evanescent modes or in solutions including bottom friction, and second, determination of forced waves that are nonlinear in the eigenvalue. It is the purpose of this note to use subsequent advances (Trefethen, 2000; Shen, 2001; Weideman and Reddy, 2000) in the exposition and availability of spectral routines in MATLAB to give highly accurate, easy-to-implement schemes for both propagating and evanescent CTWs in nonlinear eigenvalue problems. Wang and Mooers (1976), Huthnance (1978) and BC find real modes individually by iteration of the nonlinear eigenvalue problem. WH transform the problem to a linear eigenvalue problem and describe iterative methods of solution. All the methods approximate the solution on a finite-difference grid truncated at some fixed offshore position with an assumption on the vanishing of the normal derivative of the offshore or tangential component of velocity, or of the vertical structure of the solution there. Boyd (2008) (Chapter 17) points out that an alternative to domain truncation is to expand the

solution in basis functions intrinsic to the infinite interval. Thus here the basis functions in the semi-infinite horizontal offshore direction are chosen to be the Laguerre functions which automatically decay exponentially at infinity, i.e. the solution is approximated by a Chebyshev expansion in the vertical and a Laguerre expansion in the horizontal. The Laguerre expansion requires no specific offshore boundary condition and can be optimised (Shen, 2001) to provide resolution only where the pressure field differs significantly from zero. The spectral expansions concentrate grid points close to boundaries leading to highly economical representations of the solutions and significantly smaller matrices. The linearised eigenvalue problem can thus be solved immediately, without initial guesses or searching for eigenvalues, using any standard package to yield all real and evanescent modes simultaneously. These advantages become more marked for higher modes, where the accuracy of initial guesses in BC and WH becomes critical, and at stronger stratifications where direct implementation of the methods in BC and WH fails. In the limiting case of extremely long waves the offshore decay of the pressure field is slow. The direct spectral method described above remains more accurate than earlier methods but a novel, nonlinear, boundary condition can be formulated that removes the slowest decaying mode and allows very accurate eigenvalues to be obtained straightforwardly in these limiting cases. The algorithm for this new boundary condition is discussed in greater detail in Section 4.3.

2. Equations of motion  Corresponding author.

E-mail addresses: [email protected] (E.R. Johnson), [email protected] (J.T. Rodney). 0278-4343/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.csr.2011.06.009

Consider a semi-infinite rotating incompressible ocean with uniform Coriolis frequency f and depth H bounded by a monotonic

1482

E.R. Johnson, J.T. Rodney / Continental Shelf Research 31 (2011) 1481–1489

z

0

y

z=0

x

z = −h(y)

oð1o2 Þb2 B2 pz þ ohy py ¼ hy kp, ðz ¼ hðyÞ, yr 1Þ,

ð11bÞ

opy ¼ kp, ðy ¼ 0Þ,

ð11cÞ

pz ¼ 0,

ð11dÞ

ðz ¼ 1, y 4 1Þ,

a2 pz þ B2 b2 ð0Þp ¼ 0,

ð11eÞ

System (11) is closed by a far-field (y-1) boundary condition on the pressure field. In the flat outer ocean region the pressure field can be expanded as

z = −H

L

ðz ¼ 0Þ:

Fig. 1. The flow domain and pressure field for mode 2 over depth profile (29) for uniform stratification, B¼0.8 and o ¼ 0:3. The isobars are equally spaced with negative values dotted and the zero isobar omitted.

pðy,zÞ ¼

1 X

An expðbn yÞZn ðzÞ,

ð12Þ

n¼0

where sloping shelf of width of order L (Fig. 1). Take Cartesian axes with 0x along the shelf, 0y out to sea and 0z vertical, with the depth profile depending on y alone. Let the flow be Boussinesq with the total density r0 ðzÞ þ rðx,y,z,tÞ, pressure p0 ðzÞ þ pðx,y,z,tÞ and equilibrium values r0 and p0 in hydrostatic balance, i.e. @p0 =@z ¼ r0 g. For simplicity of exposition the perturbation flow is also taken to be hydrostatic although the extension to non-hydrostatic flow is immediate (WH). Introduce the buoyancy frequency N 2 ¼ ðg=r0 Þ dr0 =dz, and scalings 0

0

0

ðx ,y ,z Þ ¼ ðx=L,y=L,z=HÞ, p0 ¼ p=r0 ð0ÞfUL,

0

0

ðu ,v ,w Þ ¼ ðu=U,v=U,wL=UHÞ,

r0 ¼ rgH=r0 ð0ÞfUL, N ðzÞ ¼ N 0 bðzÞ,

ð1Þ

iouv ¼ ikp,

ð2Þ

iovþ u ¼ py ,

ð3Þ

pz r ¼ 0,

ð4Þ

iku þ vy þwz ¼ 0,

ð5Þ

iorb2 B2 w ¼ 0:

ð6Þ

The non-dimensional parameter B ¼ N 0 H=fL measures the dynamical importance of stratification relative to rotation, and is often described as a Burger number, with small B corresponding to weak stratification effects and large B to strong effects. Eqs. (2), (3), (5) and (6) combine to give the velocity components in terms of the pressure alone as

ð7Þ

The vanishing of the normal component of velocity on the bottom and coast gives ðy ¼ 0Þ,

ð8Þ

ðz ¼ hðyÞ, y r 1Þ:

ð9Þ

2 2

a pz þB b ð0Þp ¼ 0,

ð14aÞ

Zn0 ¼ 0,

ð14bÞ

z ¼ 1,

a2 Zn0 þ B2 b2 ð0ÞZn ¼ 0,

z ¼ 0,

ðz ¼ 0Þ,

ð10Þ

1=2

where a ¼ ðgHÞ =fL is the Rossby radius of deformation. Substituting for u,v,w in (5), (8) and (9) gives the standard problem (Wang and Mooers, 1976) pyy þ ð1o2 ÞB2 ðb2 pz Þz ¼ k2 p,

ð14cÞ

and the constants An give the amplitude of mode n. All modes on the shelf project onto the modes in (12), and all solutions of (11) decay exponentially at infinity.1 The rate of decay in the flat outer ocean is given by b0 , a monotone increasing function of the wavenumber k and the smallest eigenvalue l0 of (14). In general this decay is sufficiently fast to guarantee high accuracy in the method presented below. The outer ocean decay is slow for extremely long waves (kI0) in rigid lid (aI1) flows. A novel formulation for this limit is presented in Section 4.3.

3. Reduction to an algebraic eigenvalue problem For fixed o system (11) is a nonlinear differential eigenvalue problem for k. The system can be recast into a linear eigenvalue problem by writing ð15Þ

and the flow domain mapped to a rectangular domain by introducing the vertical stretching

Z ¼ y, x ¼ zmðyÞ where mðyÞ ¼ 1=hðyÞ,

ð11aÞ

ð16Þ

to give pZZ þ ð2m0 x=mÞpZx þ ðm00 x=mÞpx þðm0 x=mÞ2 pxx ð17aÞ

oð1o2 ÞB2 b2 mpx þ oh0 ½pZ þðm0 x=mÞpx  ¼ h0 kp, ðx ¼ 1Þ, ð17bÞ

o½pZ þ ðm0 x=mÞpx  ¼ kp, ðZ ¼ 0Þ,

ð17cÞ

mpx þ ½Bbð0Þ=a2 p ¼ 0, ðx ¼ 0Þ,

ð17dÞ

p-0,

The boundary condition at the free surface is 2

2

B2 ðb2 Zn0 Þ0 þ ln Zn ¼ 0,

þ ð1o2 ÞB2 mðb2 mpx Þx ¼ kq,

v ¼ iðkp þ opy Þ=ð1o2 Þ,

w ¼ ðio=b2 B2 Þpz :

w ¼ vhy ,

the Zn are the eigenfunctions (normalised to have maximum value unity) of the standard Sturm–Liouville problem

q ¼ kp,

u ¼ ðpy þ kopÞ=ð1o2 Þ,

ð13Þ

0

where L is the shelf width, H is the open-ocean depth, U is a typical horizontal velocity, f is the Coriolis parameter and N 0 is the maximum buoyancy frequency while b(z) contains the z dependence. The non-dimensional equations of motion (dropping the primes) in the body of the fluid are then, omitting a common factor expfiðoftkxÞg,

v ¼ 0,

bn ¼ ½ð1o2 Þl2n þk2 1=2 ,

Z-1,

ð17eÞ

where b ¼ bðxhÞ. Transformation (16) requires the depth at the coast to be nonzero however in practice hð0Þ can be taken to be arbitrarily small. 1 Huthnance (1978) uses a finite-difference implementation of a matching to (12) to provide the open-ocean condition for the flow on the shelf.

E.R. Johnson, J.T. Rodney / Continental Shelf Research 31 (2011) 1481–1489

Eq. (11a) shows that motions occupying the whole fluid depth have a natural horizontal scale of order B. Thus for strong stratification, where B is large, the sidewall boundary becomes dynamically equivalent to an almost vertical wall and the vertical stretching (16) is no longer an efficient mapping to a rectangle. It then becomes effective to follow Johnson (1991), write the depth profile as y ¼ dðzÞ, and introduce the transformation

Z ¼ ydðzÞ, x ¼ z,

ð18Þ

which is increasingly efficient as B-1 and also allows hð0Þ ¼ 0. Examples using this transformation are given in Section 4.4. The non-constant coefficient system (17) is approximated pseudo-spectrally following Trefethen (2000), Weideman and Reddy (2000), and Shen (2001), by forming a two-dimensional grid (tensor product grid) of N Chebyshev points, xi , in the vertical and M Laguerre points, Zj (with Z1 ¼ 0), in the horizontal. It is important to note that choosing a Laguerre grid in the offshore direction is equivalent to choosing the offshore basis functions to be Laguerre functions thus automatically satisfying the offshore boundary condition. Boyd (2008) describes this as a behavioral boundary condition at infinity and notes that exponential accuracy is usually achieved when, as here, the desired solution is the only solution bounded at infinity. The numerical results for the present problem, discussed in detail later, bear this out. The unknown values of p on the grid, taken row by row (i.e. in lexicographical, or reading, order), then form a vector p of length NM. It is shown in Appendix B that the spectral approximation to (17) can then be written as the 2NM  2NM generalised linear algebraic eigenvalue problem FðoÞr ¼ kGr,

ð19Þ

where FðoÞ ¼



AðoÞ 0

 0 , I

G¼



B I

 C , 0

r¼

p q

! ,

I is the NM  NM identity matrix and 0 is the NM  NM null matrix. Problem (19) can be solved directly by many packages however in practice it has proved more efficient to convert (19) to the more familiar regular eigenvalue problem by inverting F. Note as F comes from a discrete form of an elliptic operator it is invertible, whereas G is not. Then (19) becomes 1

HðoÞr ¼ k

r,

ð20Þ

where HðoÞ ¼ F1 ðoÞG ¼

AðoÞ1

0

0

I

!

B I

C 0

 ¼

AðoÞ1 B

AðoÞ1 C

I

0

! :

Boyd (2008) notes that the spectral discretisation of regular eigenvalue problems like (19) can introduce eigenvalues that do not correspond to physical solutions but are simply products of the discrete formulation and suggests two methods for separating ‘‘good’’ eigenvalues from ‘‘bad’’. The first relies on the position of the eigenvalues, when ordered by increasing absolute value, remaining unchanged at different resolutions: the ‘‘ordinal’’ difference

dn,ordinal ¼ jknðN1 ,M1 Þ knðN2 ,M2 Þ j=minðjkn j, sn Þ,

ð21Þ

knðNi ,Mi Þ

is the nth eigenvalue computed at resolution N ¼ Ni , where M ¼ Mi and the intermodal separation, computed using the finer resolution, is defined as

s1 ¼ jk1 k2 j, sn ¼ 12ðjkn kn1 jþ jkn þ 1 kn jÞ, n 4 1:

ð22Þ

For some problems, including (19) and (20) here, good and bad eigenvalues become interlaced at different resolutions. Therefore comparison needs to be made between the nth low resolution eigenvalue and whichever high resolution eigenvalue agrees most

1483

closely with it. This is achieved by introducing

dn,nearest ¼

min

l A ½1,N2 ,M2 

1 ,M1 Þ jkðN klðN2 ,M2 Þ j=minðjkn j, sn Þ: n

ð23Þ

Even this selection criterion can fail for complex eigenvalues if the real/imaginary part has a large/small relative change, with the real part close to zero (or vice versa). Therefore in the determination of good eigenvalues here, the denominator in (21) and (23) is changed to minðjkn j,jReðkn Þj,jImðkn Þj, sn Þ,

ð24Þ

with good eigenvalues determined by large values of 1/dn,nearest . The accuracy of the solution can be significantly improved by choosing the furthest offshore collocation point to lie where the solution is smaller than some tolerance (Weideman and Reddy, 2000; Shen, 2001). For Laguerre polynomials this corresponds to the change of variable Z ¼ aM Z~ , mapping the interval ½0,1Þ to itself. Introducing the threshold e, where jpðZ, xÞj r e for Z 4 LZ , gives the scaling factor

aM ¼ ZM =LZ ,

ð25Þ

where ZM is the largest root of the Laguerre polynomial of degree M1. Scaling also requires the horizontal differentiation matrices in the Appendix to be multiplied by the appropriate power of aM . It is important to note that neither the solution nor any of its derivatives is constrained to be zero at ZM ¼ LZ . The solution continues to decay exponentially for Z 4LZ but is simply required to be smaller than e there. The eigenvalue problems (19) and (20) can be solved by any standard algorithm. The commands eig(FðoÞ,G) and eig(HðoÞÞ in MATLAB use the QZ algorithm, making it possible to calculate all modes (both propagating and evanescent) simultaneously without initial guesses. The QZ algorithm requires OðN3 Þ operations to find all eigenvalues and eigenvectors of an N  N matrix and so if only the lowest modes are of interest it can be more efficient computationally to use inverse iteration to improve the accuracy of individual eigenmodes, i.e. solving the linear system ½HðoÞlIr^ j þ 1 ¼ r j ,

r j þ 1 ¼ r^ j þ 1 =sj þ 1 ,

ð26Þ

for iterate r j þ 1 , where l is a fixed estimate for k1 and sj þ 1 is a normalising factor for r^ j þ 1 , chosen here to be the element of r^ j þ 1 of largest modulus. The matrix on the left side of (26) is fixed and so it can be LU decomposed once-and-for-all reducing each solution of (26) to a rapid back substitution. The scale factors converge to an eigenvalue, s1 (say), of ½HðoÞlI1 , giving the eigenvalue k1 ¼ ðls1 þ1Þ=s1 of H. Given a reasonable estimate, inverse iteration is both accurate and fast. The regular eigenvalue problem (19) can thus be solved initially using the QZ algorithm, with eigenvalues sensitive to changes in (N,M) discarded. The remaining approximations to good eigenvalues and corresponding eigenvectors can then be used as initial estimates and their accuracy significantly improved by inverse iteration at increased resolution, with optimised scale parameter aM . The discussion to date has considered finding the form of modes, both propagating and evanescent, that exist at a given frequency (as needed, for example, in discussions of wavescattering). However it is equally feasible to restrict attention to propagating modes of a given wavenumber and determine the corresponding frequencies, as needed, for example, in discussions of the initial value problem of estimating the wavefield forced by a disturbance of given spatial structure. Setting q ¼ op and r ¼ oq gives a discrete representation of system (11), analogous to (19), as F1 ðkÞr 1 ¼ oG1 r 1 ,

ð27Þ

1484

E.R. Johnson, J.T. Rodney / Continental Shelf Research 31 (2011) 1481–1489

4.2. Moderate to weak stratification (B t 2)

where 0

A1 ðkÞ B F1 ðkÞ ¼ @ 0 0

0 I

1 0 C 0 A,

0

I

0

B1 B G1 ¼ @ I 0

C1 0 I

1 D1 C 0 A

and

0 1 p B C r 1 ¼ @ q A,

0

r

and the matrices A1 , B1 , C1 , and D1 are simple recastings of the matrices A, B, and C in (19). As the matrices have dimensions 3NM  3NM, the time taken by eig is increased by ð3=2Þ3  3:4. Although it is more expensive computationally to fix k and find o there can be computational advantages, as the number of good eigenvalues may increase. In particular for B small, the real part of the wavenumber k controls the decay rate in the open ocean, expðRe kÞ, and so the scaling parameter aM can be determined immediately. Since Huthnance (1978) shows that for given B and k frequencies decrease monotonically with increasing mode number fixing k and then computing o allows modes to be identified immediately. Section 4.2 gives examples of dispersion relations calculated from (27).

4. Applications 4.1. Internal Kelvin waves If the surface is taken to be rigid and both hðzÞ  1 and bðzÞ  1, the modes are baroclinic Kelvin wave with wavenumbers km ¼ omp=B

ðm ¼ 1,2,3, . . .Þ:

ð28Þ

hðyÞ ¼ 1ð1H0 Þ expfg1 ½1ð1y2 Þ2 g:

Table 1 Absolute error for selected eigenvalues for weak (B ¼0.1) and strong (B ¼2) stratification for o ¼ 0:2, hðzÞ  1, bðzÞ  1, LZ ¼ 0:3 ðB ¼ 0:1Þ and LZ ¼ 4 ðB ¼ 2Þ. Here N ¼20 and M ¼35. km error

m

1.441 109

1.734 107

2.046 1011

5.893 1012

4

7.182 1011

7.420 1011

6

5.422 1011

5.223 1011

8

2.944 109

3.824 1010

10

1.447 106

7.247 108

12

1.602 104

8.012 106

13

1.117 103

5.845 105

25

25

10

B¼ 2

1

30

15

B ¼0.1

2

30

20

ð29Þ

Here H0 is the depth at the coast, taken to be H0 ¼ 0:01 in the examples presented here, and the parameter g controls the steepness of the shelf, taken to be g ¼ 2 here. The pressure fields for the CTWs (only that for mode 2 shown here, included in Fig. 1) behave as those in Huthnance (1978): the flow is almost barotropic over the shallow shelf and upper shelf, the number of zero crossings across the continental shelf equals the mode number, and the largest displacement in the pressure field is found near to the coast. Fig. 3(a) shows the convergence of the mode 1 frequency (for k¼0.1, B ¼0.01, bðzÞ  1 and N ¼8) with increasing number M of Laguerre points at various fixed values of offshore scale factor LZ , calculated using formulation (27). When the solution is chosen to be negligible for Z 410, i.e. LZ ¼ 10, the frequency becomes resolution-independent at about MI40 Laguerre points. Increasing LZ shows however that LZ has been chosen too small. Taking the domain scale LZ ¼ 30 requires MI60 Laguerre points for resolution-independence. However this eigenvalue is now independent of LZ , as subsequently increasing LZ to 40 does not alter the eigenvalue. Note that the innate compression of the Laguerre grid near the coast and stretching over the outer ocean means that the optimum value of M increases only slowly with LZ . In all subsequent results presented here, unless otherwise stated, the values of LZ and M have been chosen sufficiently large so that the

log (1/δn,nearest)

log (1/δn,nearest)

Fig. 2 shows logðdn,nearest Þ for resolutions ðN1 ,M1 Þ ¼ ð30,25Þ and ðN2 ,M2 Þ ¼ ð35,30Þ for the first 200 numerical eigenvalues with hðzÞ  1 and bðzÞ  1. The ‘‘good’’ eigenvalues are clearly distinguishable, with values of logðdn,nearest Þ exceeding 10. Table 1 compares the good eigenvalues and analytical values given by (28), with agreement, in general, to within 1010 for the lowest modes. To demonstrate the effect of the scaling (25) the values of LZ have been held constant for each B. This affects the accuracy of the lowest eigenvalue, corresponding to the mode with slowest offshore decay, and in practice a larger value of LZ would be used for this mode if higher accuracy was required, as described in more detail below. Similarly N has been fixed to demonstrate the drop-off in accuracy for mode numbers greater than N=2, as happens with the horizontal resolution of barotropic modes in Kaoullas and Johnson (2010) (KJ here). In practice, for moderate stratifications and slopes, N should be taken to be more than twice the mode number of the highest mode required. For smaller stratifications and slopes modes have less vertical structure and N can be smaller.

Consider the shelf profile, constructed to be flat at y¼0, 1,

20 15 10 5

5 0

0 0

50

100 kn

150

200

0

50

100 kn

150

200

Fig. 2. The ratio 1=dn,nearest as a function of eigenvalue kn for o ¼ 0:2, bðzÞ  1 and hðzÞ  1. The eigenvalues are ordered in increasing absolute value and the two resolutions are ðN1 ,M1 Þ ¼ ð30,25Þ and ðN2 ,M2 Þ ¼ ð35,30Þ. (a) B¼0.1, LZ ¼ 0:3. (b) B¼2, LZ ¼ 4.

E.R. Johnson, J.T. Rodney / Continental Shelf Research 31 (2011) 1481–1489

0.279

1485

0.28 0.277

0.277 c = ω/k

c = ω/k

0.274 0.275

0.271 0.268

Lη = 10

0.273

Lη = 30

Offshore pts = 55 Offshore pts = 110 Offshore pts = 220 Offshore pts = 440

0.265

Lη = 40 0.271

0.262 20 30 40 50 60 70 80 90 100 110 120

2

4

6

8

10

12

14

16

Offshore boundary

Laguerre points (M)

Fig. 3. The mode 1 frequency, plotted as the phase speed c ¼ o=k as k is fixed, for depth profile (29) with B¼0.01, bðzÞ  1 and k¼ 0.1. (a) The phase speed as a function of offshore resolution for formulation (27) with N ¼8 and offshore scales LZ ¼ 10,30,40: (b) The phase speed as a function of distance to the offshore boundary for BL4 with 55 vertical grid points and 55, 110, 220 and 440 offshore grid points. The dashed line in (b) is the converged eigenvalue from (a).

0.8 B = 0.8

0.8

B = 0.4 0.6 B = 0.2

0.6

ω

ω

B = 0.01 0.4

0.4

0.2

0.2

0

0 0

2

4

6

8

10

k

0

2

4

6

8

10

k

Fig. 4. Dispersion curves (solid lines) for varying stratification strength B over the depth profile (29). The crosses give frequencies calculated using BL4. (a) Mode 1 for uniform stratification bðzÞ  1. The dots give wavenumbers calculated for barotropic flow using KJ. (b) Modes 1, 2 and 3 at B¼ 0.4 for bðzÞ ¼ ez .

eigenvalues are independent of LZ and M to the stated accuracy. Fig. 3(b) shows the behaviour of the same eigenvalue calculated using the CTW program ‘‘BIGLOAD4’’ of BC (BL4 here) with increasing domain size, at various fixed resolutions. The dashed line gives the accurate eigenvalue from Fig. 3(a). When the openocean boundary condition is applied at Z ¼ 2, all resolutions give approximately the same value, in error by approximately 1% but also showing that for this position of the offshore boundary condition the standard 55 offshore points is a sensible resolution. However if the offshore boundary is moved outwards to obtain a more accurate estimate of the semi-infinite domain eigenvalue then more grid points are required before a resolution-independent eigenvalue is obtained, as shown by the divergence with increasing offshore boundary distance of the lines in Fig. 3. In particular it is not clear that even 440 offshore points are sufficient to match the accuracy of the Laguerre method. Fig. 4 shows the dispersion curves for shelf profile (29) calculated using formulation (27) with values of LZ and M sufficiently large for the frequencies to be resolution- and scaling-independent. Fig. 4(a) shows curves for mode 1 with uniform stratification bðzÞ  1 and Fig. 4(b) modes 1, 2 and 3 with the surface-intensified stratification profile given by bðzÞ ¼ ez . The crosses give the frequencies from BL4 (with the standard

resolution of 55 offshore grid points, 50 vertical grid points and offshore boundary at Z ¼ 2) which are typically accurate to within 1.5%. Fig. 4(b) illustrates a further difficulty for searching methods like BL4. At any fixed k, the increment in o in moving from one mode to the next higher mode becomes smaller and smaller as the mode number increases. This requires a finer and finer search in o and increases both the computational difficulty and time required to find these modes. The direct spectral method requires no such searching and so calculates even high modes rapidly (having chosen the number of vertical grid points to be at least twice the largest vertical mode number required). In fact, to construct Fig. 6 efficiently the values from the direct Laguerre method were needed as the initial guesses for the searches in BL4. For sufficiently weak stratification the flow is essentially barotropic and so for small B the present method approaches ‘‘shelfwaves_exactbc:m’’ of KJ. Fig. 4(a) includes a comparison between dispersion curves calculated using the present method with B ¼0.01, and the barotropic model wavenumbers (the dots) solved using KJ for depth profile (29). The fractional difference between the calculated wavenumbers is less than 1%, i.e. of order B, as expected. Table 2 shows the convergence with increasing resolution of the mode 1 wavenumber for profile (29) with B ¼0.2, o ¼ 0:3 and

1486

E.R. Johnson, J.T. Rodney / Continental Shelf Research 31 (2011) 1481–1489

Table 2 Convergence ofthe wavenumber with increasing resolution for modes 1 and 3 for depth profile (29) for B¼0.2, 0.01 (mode 3), o ¼ 0:3 and LZ ¼ 5 (mode 1) and LZ ¼ 2 (mode 3). M

(a) Mode 1

40 50 60 70 80 85

(b) Mode 3

N ¼8

N ¼12

N ¼8

1.227933885 1.227885224 1.227868534 1.227864163 1.227863052 1.227863051

1.227933882 1.227885178 1.227868480 1.227864117 1.227863015 1.227863015

8.5202979  7.0098335 8.5201601  7.0102695 8.5201997  7.0102617 8.5201972  7.0102584 8.5201959  7.0102591 8.5201958  7.0102593

LZ ¼ 5, using the wavenumber-eigenvalue formulation of (20). On a 2 GHz PC a resolution of N ¼8 and M ¼80 gives eight significant figures in 1.8 s using inverse iteration. Ten eigenvalues, with accuracy ranging from three to six digits, are found in 9.8 s. BL4 takes approximately 3 s to compute a single eigenvalue to three digit accuracy. For moderate to weak stratifications at moderate frequencies, higher modes are evanescent and, as in KJ, above an order one cut-off frequency at weak stratification, all modes are evanescent. As in WH, formulation (20) computes these modes automatically. In particular, for profile (29) with B ¼0.01 and o ¼ 0:3, all modes above mode 2 are evanescent. Table 2 also shows the convergence for the evanescent mode 3. The fractional difference compared with the evanescent mode 3 wavenumber calculated using KJ is well below 1%, i.e. of order B, as expected. Note that the wavenumber at the lowest resolution differs from the high resolution result by less than 0.1%. This contrasts with the difficulty in the finite-difference method of WH (and so also BC) where the lowest resolution mode 3 was computed to be propagating and it was only by doubling resolution that the mode was correctly determined as evanescent. Even at this resolution the decay rate differed from the analytical result by over 30%.

N¼ 12 i i i i i i

8.5202980 7.0098335 8.5201601 7.0102696 8.5201997  7.0102617 8.5201972  7.0102584 8.5201963  7.0102594 8.5201958  7.0102593

^ o ^ to the jth estimates pðjÞ and oðjÞ , by requiring that corrections p, the ðj þ1Þth estimates ^ pðj þ 1Þ ¼ pðjÞ þ pðy,zÞ,

oðj þ 1Þ ¼ oðjÞ þ o^ ,

2

2 ^ ½@yy þ ð1oðjÞ ÞB2 ðb2 @zz k2 Þp2 oðjÞ o^ B2 ðb2 pzðjÞ Þz ¼ ½  pðjÞ , z @z þ b

ð33aÞ 2

^ ð12oðjÞ Þb2 B2 pzðjÞ ½oðjÞ ð1oðjÞ Þb2 B2 @z þ oðjÞ hy þkhy p^ þ o ^ hy pðjÞ ¼ ½  pðjÞ , þo

ðz ¼ hðyÞ, y r 1Þ,

^ pyðjÞ ¼ ½  pðjÞ , ½oðjÞ @y þkp^ þ o p^ z ¼ pzðjÞ

ðy ¼ 0Þ,

ðz ¼ 1, y 4 1Þ,

½a2 @z þB2 b2 ð0Þp^ ¼ ½  pðiÞ ,

pðy,zÞ-A0 expðb0 yÞZ0 ðzÞ,

as y-1:

ð30Þ

Thus ðlog pÞyy -0 as y-1, giving the nonlinear-boundary condition 2

ppyy ðpy Þ -0

as y-1:

ð31Þ

To avoid excessively large computational domains when k is extremely small (31) can be applied at some finite offshore position. It might be expected that this would need to be large but, as shown below, in some cases applying (31) at the shelf edge is sufficient. The error in applying (31) at some finite distance is given by the first neglected term in (12) which has offshore decay rate b1 which is of order unity or greater for moderate and weak stratification (strong stratification is treated in Section 4.4). Thus the error in applying (31) at some finite value yL is of order expðb1 yL Þ. In particular, since b1 is of order B1 for B of order unity or smaller it is sufficient to choose yL to be of order B from the shelf-ocean boundary. The nonlinear eigenvalue problems (11) and (31) can be reduced to a linear iteration by the Newton–Kantorovich method, following Boyd (2008). System (11) and (31) are solved for the

ð33bÞ ð33cÞ ð33dÞ

ðz ¼ 0Þ, 2

In the rigid lid limit (a-1), the smallest eigenvalue of (14) is close to zero, giving a barotropic (depth-independent) mode in the far-field. Then for extremely long waves (kI0) in the rigid-lid limit the offshore decay rate b0 from (12) of this barotropic mode is very small and the mode decays only slowly over the open ocean. For these waves at sufficiently large

ð32Þ

satisfy (11) and (31) exactly. Substituting (32) into (11) and (31) and neglecting products of small quantities gives the linear system,

ðjÞ ^ ðjÞ ðjÞ ðjÞ ^ , pðjÞ p^ yy þ pðjÞ yy p2py p y ¼ p pyy p

4.3. Long waves

i i i i i i

ð33eÞ ðy ¼ yL Þ,

ð33fÞ

where (following the notation of Adamou et al., 2007) in each equation the symbol ½   denotes the bracketed differential operator appearing on the left hand side of the same equation (and so represents a different operator in each equation). Even with the nonlinear boundary condition (31), any constant multiple of a solution p of (11) remains a solution. This arbitrary multiplicative constant can be fixed by requiring the iterate remain unchanged at one fixed point. Choosing this point to be ðy0 ,z0 Þ gives the additional equation ^ 0 ,z0 Þ ¼ 0: pðy

ð34Þ

Discretising system (11), (31) and (34) on a N  M Chebyshev grid, as in Section 3, gives an ðNM þ1Þ  ðNM þ 1Þ linear system of equations of the form ^ A^ p^ ¼ b,

ð35Þ

^ M ,zN Þ, o ^ 0 ,z0 Þ, . . . , pðy ^ ÞT . This linear system can be where p^ ¼ ðpðy solved by any standard method. As the coefficient matrix A^ changes slightly with each iteration, the iterations are slightly more computationally intensive than the inverse iterations of (26) where the coefficient matrix can be factorised once-and-for-all. Since, like the Newton method, the Newton–Kantorovich iteration ^j converges quadratically, the iteration can be terminated when jo falls below some acceptable level of error, with the final correction providing an upper bound on the error. With the accurate initial guess from the Laguerre scheme convergence is rapid, typically requiring five iterations for an accuracy of 108 to 1010 . For extremely long waves the computational cost of obtaining an initial guess from formulation (27) can be reduced

E.R. Johnson, J.T. Rodney / Continental Shelf Research 31 (2011) 1481–1489

1487

0.28 0.2815

0.279

0.281

0.278 c = ω/k

c = ω/k

0.2805 0.28 0.2795

0.277 0.276 0.275

0.279 0.2785

0.274

0.278

0.273

0.2775

Lη = 10 Lη = 30 Lη = 40 Lη = 60

0.272 0

10 20 30 40 50 60 70 80 90 100

20 30 40 50 60 70 80 90 100 110 120

Chebyshev points

M

Fig. 5. The mode 1 frequency, plotted as the phase speed c ¼ o=k as k is fixed, as a function of offshore resolution for depth profile (29) with B¼ 0.01, bðzÞ  1, N¼ 8 and k¼ 0.01. (a) The nonlinear-boundary-condition method. (b) Formulation (27) for offshore scales LZ ¼ 10,30,40,60. The dashed line in (b) is the converged eigenvalue from (a).

linear portion of the dispersion relations of Fig. 4. For higher accuracy the nonlinear-boundary-condition method should be preferred.

0.3

0.25

4.4. Moderate to strong stratification (B \ 0:5)

ω/B

0.2

As noted above, with increasing stratification it becomes more efficient to use the horizontal translation (18) rather than the vertical stretching (16). Fig. 6 shows the dispersion relations for modes 1–4 using (18) for the depth profile

0.15

hðyÞ ¼ ð1e2y Þ=ð1e2 Þ,

0.1

i:e:

dðzÞ ¼ ð1=2Þ logfz½1expð2Þþ 1g, ð36Þ

0.05

0 0

1

2 k

3

4

Fig. 6. Dispersion curves for the modes1–4 for uniform stratification bðzÞ  1 with B¼ 3 over the depth profile (36). The crosses give frequencies obtained using BL4.

by a factor of 27 by considering the long-wave limit (k-0) directly. This formulation is given in Appendix A. Fig. 5(a) shows the convergence of the mode 1 frequency with increasing offshore resolution for depth profile (29) with B ¼0.01 and k¼0.01 calculated using the nonlinear-boundary-condition method with boundary condition (31) applied at the shelf-ocean boundary (which gave the same accuracy for this example as when the boundary was moved further outwards). Convergence is rapid with seven digit accuracy achieved when M ¼45 and N ¼8. Fig. 5(b) shows the convergence of the same eigenvalue calculated using the Laguerre scheme for various values of the offshore scaling LZ . The dashed line shows the converged eigenvalue from the nonlinear-boundary-condition method. As for Fig. 3, the innate stretching of the Laguerre grid means that M¼80 is sufficient to give a resolution-independent eigenvalue for all the LZ here, accurate to within 0.1% when LZ ¼ 60. Since the phase speed c ¼ o=k for smaller k differs from the value here by order k2 (as o is odd in k), this value at k¼ 0.01 from the general Laguerre method gives the phase speed to within 0.1% even for k¼0. This is equivalent to noting that k¼0.01 lies well within the

and uniform stratification bðzÞ  1. The curves presented here were calculated using formulation (20) as this makes identifying modes easier however formulation (27) is equally as effective and accurate. BL4 failed to converge for k \0:1 for mode 3 and at all wavenumbers for mode 4 even with accurate initial guesses from the spectral method. For stratifications in the range 0:5r B r 2 both mappings are efficient and accurate. Only outside this range do their differing approaches begin to give significant efficiency gains however with sufficient resolution each method converges.

5. Discussion CTWs have been shown to have a highly economical Chebyshev–Laguerre representation that automatically incorporates the infinite offshore extent of the flow domain. This allows the transformed linear eigenvalue problem (20) to be solved rapidly and directly by standard packages for both real and evanescent modes. In particular, unlike WH and BC, no difficulties appear at low resolution for evanescent modes, in identifying higher modes at moderate stratification or in computing moderate modenumber waves at stronger stratification. WH note that the linearisation described here allows more general CTW problems to be addressed including the scattering of baroclinic waves and determining the structure of baroclinic CTWs in the presence of a baroclinic mean flow. In particular the formulation here allows the adjoint modes required for forced wave calculations to be obtained immediately from the eigenvalues and eigenmodes of the conjugate transpose of H in (20). A MATLAB version of these routines will be available on the ‘‘MATLAB Central: File Exchange’’ website.

1488

E.R. Johnson, J.T. Rodney / Continental Shelf Research 31 (2011) 1481–1489

Acknowledgement J.T.R. was supported by the NERC research studentship NE/ F008260/1 from the National Oceanography Centre, Liverpool. We are also grateful to an anonymous referee whose keen interest in the long-wave problem encouraged us to derive the novel method for slowly decaying modes presented here.

The present method can be modified to address the long-wave limit of (11) directly by taking the limit o-0, k-0 with wave speed c ¼ o=k fixed to give pyy þ B

2

ðb

pz Þz ¼ 0,

ðy ¼ 0Þ,

a2 pz þB2 b2 ð0Þp ¼ 0, pz ¼ 0,

ðz ¼ hðyÞ, y r1Þ,

ðB:2Þ

DðnÞ lag

where is the one-dimensional Laguerre differentiation matrix and the factor anM arises from the scaling (25). The matrix representation of the differential operator in (17a)

þ ð1o2 ÞB2 m2 ððb2 Þx @x þb2 @xx Þ,

ðB:3Þ

is then L1 ðoÞ ¼ Dð2Þ Z 2

ðA:1aÞ

hy p ¼ cðb2 B2 pz þ hy py Þ, p ¼ cpy ,

ðnÞ n @n =@Zn -DðnÞ Z ¼ aM IN  Dlag ,

L1 ¼ @ZZ þ ð2m0 x=mÞ@Zx þðm00 x=mÞ@x þ ðm0 x=mÞ2 @xx

Appendix A. The long-wave limit ðk-0Þ

2

arises from mapping the Chebyshev interpolation points on the standard interval ½1,1 to the interval ½1,0. Similarly

2 ð2Þ  diagðn  h1 ÞDZx þ diagðn  h01 ÞDð1Þ x þdiagðn  h1 Þ Dx

þ ð1o2 ÞB2 IN  diagðh

ðA:1bÞ

2

2

Þ½diagðDð1Þ x b

ÞDð1Þ x þ diagðb

2

ÞDð2Þ x , ðB:4Þ

ðA:1cÞ where ðz ¼ 0Þ,

ðz ¼ 1, y4 1Þ:

ðA:1dÞ ðA:1eÞ

System (A.1) is a linear eigenvalue problem for c, which, following Section 3, can be written as a NM  NM generalised eigenvalue problem (one quarter the size of (19) and one ninth the size of (27)) H1 p ¼ cE1 p:

ðA:2Þ

To resolve the structure of the mode over the shelf and cater for the slower decay scale, it is advisable to choose larger values for offshore scale LZ in (25) and the number of offshore grid points, M. The first four modes for the depth profile (29) with B ¼0.01 can be found to within 1% with LZ ¼ 10, M ¼50 and N ¼8. If higher accuracy is required then these values can be used as initial guesses in the nonlinear-boundary-condition method of Section 4.3 or in the reduced form it takes when applied to equations (A.1). It may be computationally more efficient to use the method of Schmidt and Johnson (1993), specifically designed for arbitrary stratification in the long-wave limit, where the problem is reduced to a one-dimensional integral equation across the shelf, allowing arbitrarily high resolution there. Alternatively, for barotropic flow, the programs of KJ allow all the resolutions to be placed over the shelf.

 Dð1Þ ÞaM , DZx ¼ ð2Dð1Þ cheb lag 2

b

ðB:5Þ

¼ b2 ðn  hÞ,

ðB:6Þ 2

and n, h, h1 , and h are vectors with elements xi , hðZj Þ, h0 ðZj Þ=hðZj Þ, 2½h ðZj Þ=hðZj Þ2 h00 ðZj Þ=hðZj Þ, and h2 ðZj Þ, respectively. For any l-vector x, diagðxÞ is the l  l matrix with the elements of x on the main diagonal and zeroes elsewhere. Similarly the differential operator in (17b) h10 0

L2 ¼ oð1o2 ÞB2 b2 mpx þ oh0 ½pZ þðm0 x=mÞpx 

ðB:7Þ

becomes the matrix L2 ðoÞ ¼ oð1o2 ÞB2 diagðb

2

1

Þ½IN  diagðh

ÞDð1Þ x

0

ð1Þ þ o½IN  diagðh Þ½Dð1Þ Z diagðn  h1 ÞDx , 1

ðB:8Þ 0

where h is the vector with elements h1 ðZj Þ and h is the vector with elements h0 ðZj Þ. The differential operator in (17c) L3 ¼ o½pZ þ ðm0 x=mÞpx 

ðB:9Þ

is represented by the matrix ð1Þ L3 ðoÞ ¼ o½Dð1Þ Z diagðn  h1 ÞDx :

ðB:10Þ

Thus equations (15) and (17) can be written as the matrix eigenvalue problem Appendix B. Construction of the differential operator matrices The spectral approximation to any derivative of p on the tensor product grid is given by multiplying p by an NM  NM differentiation matrix. Weideman and Reddy (2000) give the form of these matrices for the one-dimensional (1D) case for both Chebyshev and Laguerre points. On the two-dimensional grid differentiation with respect to the offshore direction operates row by row, with partial differentiation with respect to the vertical direction operating column by column. The two-dimensional versions on the reordered grid follow by forming the Kronecker tensor product (Trefethen, 2000). The Kronecker tensor product of two matrices A and B of dimension p  q and r  s, respectively, is denoted A  B and forms a new p  q block matrix of dimension pr  qs, where the i,j block is A i,j B. So, for example n

n ðnÞ @n =@x -DðnÞ x ¼ 2 Dcheb  IM ,

DðnÞ cheb

AðoÞp ¼ kðBp þCqÞ,

ðB:11aÞ

q ¼ kp,

ðB:11bÞ

where A,B and C are NM  NM matrices. Matrix A is constructed by replacing those rows of the field equation operator matrix L1 ðoÞ that represent points lying on the boundary with the corresponding rows of the boundary operator matrices L2 ðoÞ, L3 ðoÞ and Dð1Þ x . To satisfy the surface boundary condition the first M rows of L1 ðoÞ are replaced by the first M rows of 1 2 ½IN  diagðh ÞDð1Þ x þ½Bbð0Þ=a INM . To satisfy the bottom boundary condition the last M rows of L1 ðoÞ are replaced by the last M rows of L2 ðoÞ. To satisfy the sidewall boundary condition rows ðn1ÞM þ 1 (n ¼ 2, . . . ,N) of L1 ðoÞ are replaced2 by the corresponding rows of L3 ðoÞ. It remains to incorporate the right side of system (17) into the matrices B and C. Since q appears only on the right side of the governing equation (17a), matrix C is

ðB:1Þ

where is the N  N one-dimensional Chebyshev differentiation matrix, IM is the M  M identity matrix and the factor 2n

2 This is why, following Weideman and Reddy (2000), the point Z1 ¼ 0 is added to the Laguerre points.

E.R. Johnson, J.T. Rodney / Continental Shelf Research 31 (2011) 1481–1489

constructed by replacing rows ðn1ÞM þ 1, (n ¼ 2, . . . ,N) and the first and last M rows, of INM with zeros. Matrix B is assembled by replacing rows ðn1ÞM þ 1 (n ¼ 2, . . . ,N) of the NM  NM null matrix with the corresponding rows of INM . Finally, to incorporate the right side of (17b) into B, the last M rows are replaced 0 with the corresponding rows of the matrix IN  diagðh Þ. The matrix representation for the transformation (18) follow analogously to those for vertical stretching and so the details are omitted here. References Adamou, A.T.I., Craster, R.V., Llewellyn Smith, S.G., 2007. Trapped edge waves in stratified rotating fluids: numerical and asymptotic results. Journal of Fluid Mechanics 592, 195–220. Boyd, J.P., 2008. Chebyshev and Fourier Spectral Methods, second ed. Dover Publication INC., New York. Brink, K.H., Chapman, D.C., 1987. Programs for computing properties of coastal trapped waves and wind-driven motions over the continental shelf and slope, second ed. Technical Report, Woods Hole Oceanographic Institution, Woods Hole, MA. WHOI Technical Report 87-24, 119 pp.

1489

Huthnance, J.M., 1978. On coastal trapped waves-analysis and numerical calculation by inverse iteration. Journal of Physical Oceanography 8 (1), 74–92. Johnson, E.R., 1991. The scattering at low-frequencies of coastally trapped waves. Journal of Physical Oceanography 21 (7), 913–932. Kaoullas, G., Johnson, E.R., 2010. Fast accurate computation of shelf waves for arbitrary depth profiles. Continental Shelf Research 30 (7), 833–836. Martinez, J.A., Allen, J.S., 2004. A modelling study of coastal trapped wave propagation in the Gulf of California. Part I: response to remote forcing. Journal of Physical Oceanography 34, 1313–1331. Schmidt, G.A., Johnson, E.R., 1993. Direct calculation of low-frequency coastally trapped waves and their scattering. Journal of Atmospheric and Oceanic Technology 10 (3), 368–380. Shen, J., 2001. Stable and efficient spectral methods in unbounded domains using Laguerre functions. SIAM Journal on Numerical Analysis 38 (4), 1113–1133. Trefethen, L.N., 2000. Spectral Methods in MATLAB. SIAM, Philadelphia. Wang, D.P., Mooers, C.N.K., 1976. Coastal trapped waves in a continuously stratified ocean. Journal of Physical Oceanography 6, 853–863. Webster, I., Holland, D., 1987. A numerical method for solving the forced baroclinic coastal-trapped wave problem of general form. Journal of Atmosphere and Oceanic Technology 4 (1), 220–226. Weideman, J.A.C., Reddy, S.C., 2000. A MATLAB differentiation matrix suite. ACM Transactions of Mathematical Software 6, 465–519.