1. Introduction

Coastal-trapped waves (CTWs) are a class of subinertial waves within the weather frequency band, which have elegant mathematical formulation; in the Northern Hemisphere, they propagate unidirectionally such that the coast lies on the right. In coastal hazard study, the role of CTWs is rarely touched, though its importance in exchanging mass and energy has been well recognized (e.g., [1]). For example, the wind-driven variabilities over the continental margin can be largely explained in terms of CTW properties. These waves play a key role in coastal oceans’ response to atmospheric weather changes; they make a major mechanism in the coastal ocean that spreads the atmospheric energy in the alongshore direction, leading to the formation of the remotely driven upwelling, and hence regulating coastal processes from far away. Historically, the low-frequency CTW wave motions were first identified along the East Australian coast [2], and soon studied in terms of continental shelf waves (CSW) by Robinson [3], reflecting the early effort in this field at Harvard University ([4], cf. [5]). Systematic CTW studies began with Allen [6], Wang and Mooers [7], Clarke [8], Huthnance [9], to name several. A comprehensive review is referred to Brink [10].

The East China Sea (ECS), including the Yellow Sea (YS), is one of the major marginal seas in the Northwestern Pacific. Despite its long coastline, CTWs are rarely investigated for this region. It seems that most of the CTW studies are in 80s to early 90s, though recently there is a renewed interest with the aid of Hilbert-Huang transform [11]. For example, Feng [12] examined the dissipative trapped waves on a broad shelf; Liu and Qin [13] investigated the instability of CTWs; Kong et al. [14] studied the impact of the ECS CTWs on Kuroshio; and Chen and Su [15] gave a systematic investigation of continental shelf waves (CSW) along the East China coast. This is in contrast to the South China Sea (SCS), which has attracted much more interest for CTW research; a recent study is referred to Zheng et al. [16].

As the northern part of the ECS, the YS is even less studied in the light of CTWs. This is because the YS is a nearly closed basin with a broad, shallow double shelf, a deep trough lying in between (Figure 1). Strictly speaking, CTWs in the traditional sense may not be significant, and, indeed, the YS wave motion is dominated by tides. But previous studies did confirm the existence of CTWs in the YS. Particularly, Chen and Su [15] found they have frequencies around 0.21–0.32 cycle per day (cpd), and a speed of roughly 15 m/s; Wang et al. [17] identified two free CTWs, with phase speeds of 14–19 m/s. A careful and detailed study of the YS CTW was conducted by Hsueh and his group (cf. Hsueh and Pang [18]); it is particularly for the wave motion along the Korea coast.

Figure 1.
(a) Bathymetry of the Yellow Sea. Abutted to the left is China, while on its right lies Korea. The solid line across Subei Bank marks the cross-section which we will study. Also marked are the points A and B where time series will be extracted from the HYbrid Coordinate Ocean Model (HYCOM) dataset for analysis. (b) The topography along the section marked in (a). (c) A profile of the buoyancy frequency squared N 2 (for 1995).

So far there is almost no study on the CTWs along the Jiangsu coast (particularly off Subei Bank) in the YS. This may be due to the lack of a deep ocean to its east. But, if we focus on some motions local enough, it is conceivable that the trapping condition, that is, vanishing far field, that admits the CTWs may still hold. Besides, such subinertial signals do have been identified (e.g., [15, 17]). Considering the potential importance of CTWs on the processes off Subei Bank, we henceforth conduct a preliminary study for this region. To fit for the scope of this book, we focus on the slow CTW modes and examine their climatic variabilities (see below). On the other hand, we take this opportunity to present a brief but self-contained introduction of the CTW theory, in order to provide environmentalists a reference for this important topic, which is yet to be explored in coastal hazard research.

2. A brief introduction of the coastal-trapped wave theory

2.1. Governing equations

Take a right-hand coordinate system x y z with x directing offshore and z upward. With the assumption of hydrostaticity and incompressibility, the linearized Boussinesq equations for the free-wave problem on an f -plane in the absence of a basic flow are

u t − fv = − 1 ρ 0 P x + 1 ρ 0 τ z x E1

v t + fu = − 1 ρ 0 P y + 1 ρ 0 τ z y E2

ρ = − 1 g P z E3

u x + v y + w z = 0 E4

ρ t − ρ 0 g N 2 w = 0 E5

where u v w , ρ , and P are the perturbation velocity, density, and pressure, respectively. The turbulent stress τ x τ y will be needed when bottom drag and wind forcing are applied; otherwise, it is usually neglected in the ocean interior. Here, the density field is decomposed into three parts: ρ 0 + ρ ¯ z + ρ x y z t with the perturbation field ρ ≪ ρ ¯ ≪ ρ 0 . ρ 0 (usually taken as 1000 or 1025 km/m³) and g = 9.8 m/s² are constants, and N 2 = − g ρ 0 ∂ ρ ¯ ∂ z is the Brunt-Väisälä or buoyancy frequency; N 2 = N 2 z is given. Suppose that the bottom depth has a dependence on x only, h = h x , as schematized in Figure 2. The boundary conditions are that the perturbation vanishes as x → + ∞ , and that no flux crosses the solid boundary. The former is the very costally trapped condition. In the figure, the location x = b is pretended to be where the coast is. In practice, it has been chosen such that h b = 22.5 m (Clarke and Brink, 1985), and such that h b is equal to three times the Ekman layer e -folding decay scale δ [19].

Figure 2.
A schematic of the topographic configuration and its discretization ( k and j are the indices). The location x = b is pretended to be where the coast lies.

In the absence of turbulent stress, (1)–(5) can be reduced to one single equation for P . Let ζ = ∂ v ∂ x − ∂ u ∂ y and D = ∂ u ∂ x + ∂ v ∂ y , respectively, be vorticity and divergence, then it is easy to obtain

ζ t + f D = 0 , D t − f ζ = − 1 ρ 0 ∇ 2 P .

But D = − W z = ∂ ∂ z 1 ρ 0 N 2 P zt by the continuity and hydrostatic equations (4) and (5). These yield

f 2 + ∂ 2 ∂ t 2 ∂ ∂ z P zt N 2 + ∇ 2 P t = 0 . E6

In the “long-wave limit”

∂ 2 / ∂ t 2 ≪ f 2 , ∂ 2 / ∂ y 2 ≪ ∂ 2 / ∂ x 2 ,

the equation is simplified to

P xxt f 2 + P zt N 2 z = 0 . E7

With the same long-wave assumption, the boundary conditions are (Clarke and Van Gorder [20]):

at x = b , P xt + rP x h + f P y = f τ y h E8

at z = 0 , P z = 0 rigid lid , or P z + N 2 g P = 0 free surface E9

at z = − h , f 2 P zt N 2 + h x P xt + fP y + rP x x − h x rP xz = 0 , E10

as x → ∞ , P x → 0 , E11

where τ y is the alongshore component of the wind stress, and r is the coefficient of the y -bottom drag in a linear form: τ B y = ρ 0 rv B . Note Eq. (11) is the far field condition; in practice, it may be replaced by another convenient condition ( u x = 0 at x = L ) (see [10]). Eq. (8) is the no-flux condition: ∫ − h b 0 udz = 0 .

Eq. (9) is the boundary condition at z = 0 for a rigid-lid or a free surface. In this study, we will adopt the rigid-lid assumption, which implies a nondivergent horizontal flow, and effectively filters out the fast external waves. In the YS/ECS case, there have been strong arguments against using it because the number Γ = fL 2 gH , where L is the width of continental shelf and H the typical depth, is not small (e.g., [11, 15, 18]). Here for the purpose of this monograph, we will focus on the climatic variability of the YS CTWs, which are only associated with slow modes. As such, a rigid lid on top of the ocean makes a more convenient configuration. For discussions on the free surface condition, see, e.g., Brink [19].

Motivated by the form of (7), Clarke and Van Gorder [20] introduced a vector

M = P x f 2 i + P z N 2 k , E12

and the equation is then simply ∇ ⋅ M t = 0 .

2.2. Eigenvalue problem

In order to avoid getting into too much technical details, we only consider the frictionless situation: r = 0 . For free CTWs, τ y = 0 . We seek a wave solution in the form

P x y z t = F x z ϕ y − ct . E13

Substitution of it back into (7) and the boundary conditions (8)–(11) gives

F xx f 2 + F z N 2 z = 0 , E14

at z = 0 , F z = 0 , E15

at x = b , c f F x − F = 0 , E16

at z = − h x , c h x F x + f 2 N 2 F z − fh x F = 0 , E17

as x → ∞ , F x → 0 . E18

Note the minus signs in the conditions at x = b and z = − h x ; they are positive in Brink [19] and Clarke [8] because a left-hand coordinate system is adopted therein. In terms of the vector M introduced by Clarke and Van Gorder [20]), the solution form is M = B ϕ y − ct , where

B = F x f 2 i + F z N 2 k . E19

Let the whole domain for the configuration (cf. Figure 2) be Ω and denote the boundary that encloses it as ∂Ω . Introduce an inner product space H over Ω . Its inner product F m F n H is defined such that, for any F , G ∈ H ,

F G H = ∫ − h b 0 FG x = b dz + ∫ b ∞ h x FG z = − h x dx . E20

Considering the vanishing far field condition and the rigid-lid surface condition, this is actually an inner product evaluated at the boundary ∂Ω . This is quite different from the inner product space in the usual sense. It is easy to show that F m F n H satisfies the needed axioms.

Let F m be a set of the eigenfunctions of (14)–(18). Clearly F m ∈ H for all integers m ; they possess the following important property [7, 8]:

Theorem 2.1 F m forms an orthogonal set in H ; that is to say, F m F n H = 0 if m ≠ n .

PROOF. The orthogonality can be proved in different ways. In Clarke and Van Gorder [20], it is easily obtained through manipulating the identity

∇ ⋅ F n B m − F m B n = 0 ,

which integrated over Ω must also vanish. By Gauss’ theorem,

∫ ∂Ω n ⋅ F n B m − F m B n ds = 0 ,

where ds is the arc length along ∂Ω . With the boundary conditions taken into account, this becomes

1 c n − 1 c m ∫ x = b F m F n ds + ∫ z = − h x h x F m F n 1 + h x 2 ds = 0 .

For m ≠ n , this is

∫ x = b F m F n ds + ∫ z = − h x h x 1 + h x 2 F m F n ds = 0 ,

or

∫ − h b 0 F m 0 z F n 0 z dz + ∫ x = b ∞ F m x − h x F n x − h x h x dx = 0 . E21

Q.E.D.

This theorem states that the eigenvectors are orthogonal in the inner product space H . They can always be normalized by the induced norm ∥ F m ∥ = F m F m H , and then form an orthonormal set. [ ∥ F m ∥ is shown to be f ; see Clarke [8]]. Caution should be used, however, in the normalization, which must be made such that energy is conserved as waves propagate through a slowly varying medium [21].

From the orthogonality condition Wang and Mooers [7] connected CTWs to internal Kelvin waves (flat bottom), shelf waves, and topographic Rossy waves ( b = 0 ). They also showed that there exists another orthogonality condition over Ω rather than along ∂Ω ; that is:

∫ b ∞ ∫ − h x 0 ∂ F n ∂ x ∂ F m ∂ x + 1 N 2 ∂ F n ∂ z ∂ F m ∂ z + ℓ 2 F n F m dxdz = 0 ,

if c n ≠ c m , where ℓ is the wavenumber in y . But whether this satisfies the inner product axioms is yet to be checked.

In addition to orthogonality, it is also possible to prove that F m forms a complete set (see Courant and Hilbert, pp. 424–426). That is to say, the normalized F m makes a basis of H .

2.3. Excitation of CTWs

By the forgoing argument, F n makes a basis of H . Thus, any function of H , P in particular, can be written as a linear combination of F n , n = 1 , 2 , … ,

P = ∑ n = 1 ∞ F n x z ϕ n y t . E22

When a forcing τ y is applied at the boundary, we can then solve for the amplitudes ϕ n .

Usually, the amplitude equation can be obtained through taking the inner product of F n with the governing Eq. (7), or ∇ ⋅ M t = 0 . In this particular case, it is not this straightforward because the inner product is evaluated along the boundary ∂Ω , rather than in Ω . Besides multiplying F n with ∇ ⋅ M t , Clarke and Van Gorder [20] considered a duality P t B n . They found that

∇ ⋅ F n M t − P t B n = 0 ,

which of course vanishes if integrated over Ω . By Gauss’ theorem,

∮ ∂Ω n ⋅ F n M t − P t B n ds = 0 .

Expanding and taking into account of the boundary conditions, this is eventually reduced to (see Clarke and Van Gorder [20])

∫ − h b 0 F n P y − P t c n x = b dz + ∫ − h b 0 F n rP x fh x = b dz + ∫ b ∞ F n P y − P t c n z = − h x dx + ∫ b ∞ F n x − h d dx rP x x − h x f dx = ∫ − h b 0 F n b z dz τ y y t h b .

Substitution of (22) for P , and use of the orthogonality property: F m F n H = 0 if m ≠ n , yields

− 1 c n ∂ ϕ n ∂ t + ∂ ϕ n ∂ y + ∑ m = 1 ∞ a mn ϕ m = b n τ y y t ,

where

b n = 1 D n h b ∫ − h b 0 F n b z dz ,

a mn = 1 D n ∫ b ∞ F n x − h 1 f d dx rF mx x − h dx + ∫ − h b 0 F n rF mx fh x = b dz

D n = F n F n H = ∥ F n ∥ 2 .

In the case, only one mode n is excited, there is no coupling: ∑ 1 ∞ a mn ϕ m = a nn ϕ n (e.g., Battisti and Hickey [22]), the amplitude equation can be solved analytically in terms of characteristic coordinates [19].

s = t − ∫ y 0 c n − 1 d y ′ , ξ = y .

In this case,

ϕ n y t = ϕ n 0 t − ∫ y 0 c n − 1 dξ ⋅ e − ∫ y 0 a nn dξ + ∫ y 0 b n ξ τ y ξ t − ∫ y ξ c n − 1 d ξ ′ ⋅ e − ∫ y ξ a nn d ξ ′ dξ ,

where y = 0 signifies the starting point of the calculation. See Brink [19] for details.

Once ϕ is solved, P may be obtained through reconstruction with modes. From (1)–(2), it is easy to obtain v tt + f 2 v = − 1 ρ 0 P yt + f ρ 0 P x , which in the long-wave limit simply shows the geostrophic balance: f 2 v = fP x / ρ 0 . The fields then can be obtained accordingly:

u = − 1 ρ 0 f 2 P xt + f P y , v = 1 ρ 0 f P x , w = − P zt ρ 0 N 2 , ρ = − 1 g P z .

3. The CTWs off the Subei Bank in the Yellow Sea

3.1. Background fields

We choose a cross-section outside the Subei Bank as shown in Figure 1. It is chosen in the middle since in the north the shelf is too shallow, while in the south the Changjiang diluted water complicates the circulation system very much and hence may disguise the wave propagation. The HYCOM data are used to set up the model; they are advantageous in that tides are excluded in the reanalysis fields. The horizontal resolution is 0.1° in longitude, by 0.0702° in latitude, which results in a horizontal spacing of 12 km along the section. In the vertical direction, the levels are at 0, 2, 4, 6, 8, 10, 12, 15, 20, 25, 30, 35, 40, 45, 50, 60, 70, 80 m, respectively. The data are available for 1993 through 2013. After 2013, the vertical resolution is too coarse for our purpose, and hence they are not considered in this study.

The topography along the section is displayed in Figure 1b; it is slightly modified to avoid the bumps (the dashed line). The basic buoyancy frequencies differ year by year. They are calculated from the annual mean density anomaly profiles averaged over along the section. Shown in Figure 1c is such an example.

3.2. Eigenvalues and eigenmodes

With the background field, we solve the eigenvalue problem. The solution method is referred to the appendix; it is essentially the same as that by Brink [10], Clarke and Gorder [20], and Wilkin and Chapman [23]. Considering the broad shallow shelf, we artificially set the coastline at the isobath H = 10 m. By computation, four typical trapped eigenmodes have been identified. For the year of 2000, the celerities are, in order of magnitude, − 0.12 , − 0.19 , − 0.21 , − 0.29 m/s; the corresponding modes will be henceforth refer to as M1, M2, M3, and M4, respectively. Besides, many slower modes such c = − 0.0046 m/s have also been found, because the model is too coarsely resolved that we do not consider them here. Similarly, a few higher modes, particularly the one with c = − 0.57 m/s, are evident. But these higher modes are not trapped to the western topography or coast and hence should be discarded. An observation is that all these celerities are negative; that is to say, all the corresponding waves propagate to the south. This is the very important property of CTW: in the Northern Hemisphere, the wave propagates with the coast on the right hand side.

In order to see whether these eigenvalues indeed exist, we examine the time series at two points A 122.6 o E 34.825 o N and B 122.9 o E 34.035 o N (see Figure 1). They are roughly along the isobath H ≈ 64 m, located ∼ 92 km apart. Again, we use the HYCOM reanalysis data to form the time series. Since too many data are missing for the flow field, consider only the density anomaly. The series are from January 01, 1993 through December 31, 2012, with an interval of 1 day.

The dominant modes of the series have periods of 365 and 183 days, respectively; these are the annual and biannual modes (not shown). However, in the weather frequency band, a local peak is also significant. We use the wavelet method to filter out the signal within this band. Choose the series between October 12, 2001 to December 31, 2012 to ensure that they have a length of 4096 = 2 12 steps (a length of power of 2 is required). We use the orthonormalized cubic spline wavelets constructed in Liang and Anderson (2007) to fulfill the filtering. The scale window is chosen to be bounded by periods 2 and 16 days, respectively. This roughly corresponds to the weather frequency band. We then perform a time-delayed correlation analysis with the two band-pass filtered series, written ρ A t and ρ B t . Given a delay τ (in days), we then compute the correlation coefficient γ τ between ρ A t and ρ B t + τ . ( τ should not be too large, in accordance with the weather frequency band requirement.) Suppose that γ τ takes its maximum at T , then the phase speed can be estimated to be d / T . Notice that it is not a good idea to compute γ for all the 4096 days, since the CTWs are coastal wind-driven and within the weather frequency band (1 week or 2 weeks or so). Different modes may be excited during different periods. We hence choose to do the correlation analysis every 30 days. The results may differ, but the CTW modes consistent with our calculation frequently show up. For example, in Figure 3a, the maximum γ takes place at a delay of 8 days for ρ B , so the modal phase velocity is 92000 / 8 × 24 × 3600 = 0.13 m/s. Similarly, in Figure 3b, γ takes its maximum at 5 days, which corresponds to a speed of 0.21 m/s. Considering there may be a slight change in magnitude with time (see below), these celerities agree very well with the celerities of M1 ( c = 0.12 ) and M3 ( c = 0.21 ) in 2000. We actually find many cases with a speed of 0.53, a high mode that have been discarded, and many other slower modes, which we do not consider here. Note the time resolution (1 day) of the series cannot distinguish M2 and M3.

Figure 3.
The correlation coefficients between ρ A and the delayed ρ B for the periods (a) 12 Aug 2005–2011 Sep 2005 and (b) 04 Feb 2011–2106 Mar 2011. The locations A and B are marked in Figure 1.

The eigenmodes corresponding to the four celerities are shown in Figure 4. From it we see that, the slower the wave, the more topographically trapped the structure. Particularly, in the first subplot, we see a bottom-trapped mode ( c = − 0.12 m/s).

Figure 4.
The cross-section structures of the four coastal trapped wave modes M1, M2, M3, and M4, which have celerities of c = − 0.12 , − 0.19 , − 0.21 , − 0.29 m/s, respectively.

The modal structures for other years are similar, but the propagating speeds differ. In the next section, we will see more details about this.

4. The change in CTW celerities during the past 20 years

The buoyancy frequency profiles change with time, which will inevitably affect the CTW propagating properties. We have computed the eigenvalues for the six identified modes for all the past 20 years. The results are shown in Figure 5.

Figure 5.
Time variation of the celerities of the four modes M1–M4 as shown in Figure 4.

Interestingly, from the figure, all the identified modes seem to have the same evolutionary trend. On the whole, the computed speeds are all on decline (note the negative sign) during the past 20 years. This might be related to global warming, but what makes it puzzling is that, simply warming the ocean cannot change the CTW wave propagating; rather, it is the alteration of the vertical stratification that may cause the change. For such a shallow sea, it is still unclear how this can be fulfilled. Even though it is related to global warming, it must be through some third mechanism(s). In the Yellow Sea, a conspicuous hydrographic feature is its cold water mass in summer; maybe that functions to mediate the events–we will explore this in future studies.

Another interesting feature of Figure 5 is the maxima (or minima if magnitude is considered) around 1997 and 2009. This reminds one of El Niño, as around those are the El Niño years; particularly, the 97–98 El Niño is very strong, as shown in Figure 6. However, from the figure one may argue that 02–03 is also a moderately strong El Niño year, while in Figure 5, the c ’s reach their respective minima. We hence need to pick more quantity to compare. As shown in Figure 6, the La Niña intensity may be such a candidate. Notice that immediately after the 97–98 El Niño, there is a strong La Niña; same happens after the 09–10 one, while after 2003 no La Niña occurs at all. Therefore, it seems that the CTWs are greatly slowed down when a strong El Niño followed by a strong La Niña occurs.

Figure 6.
ENSO events and their intensities since 1990. The figure is based on the Oceanic Niño Index; see http://www.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ensoyears.shtml.

We emphasize that we do not claim that ENSO has caused the slow-down of the CTWs. What we have shown is just a remarkable correlation between the wave property change and the ENSO occurrence—In fact, even correlation may not be an appropriate term here. We should have computed the correlation coefficients between the c ‘s and some ENSO index, say, the Oceanic Niño Index as shown in Figure 6, but the time resolution of the c series are very low; there is only one value per year. We have also tried CTW configurations with monthly mean fields, but the dominating seasonal cycle completely disguises the signal, which just worsens the situation.

5. Concluding remarks

Coastal trapped waves (CTW), which flourished during 60–80s, thanks to the efforts from investigators such as Allen [6], Wang and Mooers [7], Clarke [8], Huthnance [9], Mysak [24], Brink [19], Chapman [25], among many others, play an important role in the processes of coastal and continental oceans (cf. [10]). However, this important topic seems to be missing in the heated debate in recent years on coastal environmental change and protection. In this chapter, we have given a brief but self-contained introduction of the CTW theory, in the hope of gaining renewed interest in this field. The nice properties such as the mutual orthogonality between the eigenmodes make this class of waves not just important in practice, but also theoretically appealing. For reference, a numerical scheme is supplied in the appendix to solve the resulting eigenvalue problem.

As a demonstration, we have conducted a preliminary study of the CTWs for a section across Subei Bank in the Yellow Sea, a region where traditionally the role of CTWs has been overlooked due to its shallow depth. Considering the scope of this monograph, we have particularly focused on the slow wave modes and their variabilities during the past decades. Interestingly, all the identified prominent modes, including a bottom-intensified one, seem to be slowing down during the past 20 years. Particularly, it seems that they are greatly slowed down during the events when a strong El Niño is followed by a strong La Niña, such as the 97–98 and 09–10 El Niño events.

However, caution should be used in interpreting these results. Apart from the poor time resolution (one value per year), the obtained c -series are not long enough for climate variability studies. It would be better if the data can be extended so that decadal variability can be examined; particularly one would like to see how they behave during the two super El Niño years: 1982–1983 and 2015–2016 (currently data are not available). We will investigate this in future studies.

Acknowledgments

Section 2 and Appendix are extracted from an unpublished technique report (Liang [26]) by the author in the former Harvard Ocean Group, which was prepared for the project AOSN-II then (cf. [27]). Thanks are due to Kenneth Brink for the help in drafting that report. The HYCOM group is acknowledged for proving the reanalysis data. Yineng Rong prepared the background field of the cross-section and the point series for analysis and drew Figure 1. This study was partially supported by the 2015 Jiangsu Program for Innovation Research and Entrepreneurship Groups, by the National Program on Global Change and Air-Sea Interaction (GASI-IPOVAI-06), and the National Science Foundation of China under Grant 41,276,032 [28, 29, 30].

A.

The x – z domain is discretized into a mesh grid, as sketched in Figure 2. The spacings Δ x and Δ z may vary with the respective indices j and k . We use j , instead of i , to indicate the x index, in order not to confuse with the complex number − 1 . In the grid, the eigenfunction F is defined at integer points k , while the buoyancy frequency N 2 is at half integer points (Figure 7).

Figure 7.
Discretization in the z direction.

For j , let j b and j m be the boundary limits, which correspond to x = b and x = L ( L some large number). Note for each horizontal level k , the limits for j should be j cst k and j m , with j cst dependent on k . Similarly denote the bottom index as k m , then k runs from 1 to k m , and the rigid-lid and bottom correspond to k = 1 2 and k = k m + 1 2 , respectively. The limit k m = k m j is a function j , which is determined by the function h x .

Let

δ x = 1 , x = 0 ; 0 , x ≠ 0 ,

and

δ sea ≡ δ j − j ∞ ; c sea ≡ 1 − δ sea δ cst ≡ δ j − j cst k ; c cst ≡ 1 − δ cst δ top ≡ δ k − 1 ; c top ≡ 1 − δ top δ bot ≡ δ k − k m j ; c bot ≡ 1 − δ bot .

Further define

Δ x + ≡ Δ x j + 1 + Δ x j 2 , Δ z + ≡ Δ z k + 1 + Δ z k 2 , Δ x − ≡ Δ x j − 1 + Δ x j 2 , Δ z − ≡ Δ z k − 1 + Δ z k 2 , γ x + ≡ Δ x j − Δ x j + 1 2 Δ x + + Δ x − , γ z + ≡ Δ z k − Δ z k + 1 2 Δ z + + Δ z − , γ x − ≡ Δ x j − 1 − Δ x j 2 Δ x + + Δ x − , γ z − ≡ Δ z k − 1 − Δ z k 2 Δ z + + Δ z − .

Then

∂ 2 F ∂ x 2 = 1 Δ x j ⋅ c sea ⋅ F j + 1 , k − F j , k Δ x + + γ x + c cst F j + 1 , k − F j , k Δ x + − F j , k − F j − 1 , k Δ x − + δ sea ⋅ ∂ F ∂ x ∞ k − c cst ⋅ F j , k − F j − 1 , k Δ x − + γ x − c sea F j + 1 , k − F j , k Δ x + − F j , k − F j − 1 , k Δ x − − δ cst ⋅ ∂ F ∂ x j cst k k ,

and

f 2 F z N 2 z = f 2 Δ z k ⋅ c top N k − 1 / 2 2 ⋅ F j , k − 1 − F j , k Δ z − + γ z − c bot F j , k − 1 − F j , k Δ z − − F j , k − F j , k + 1 Δ z + + δ top N k − 1 / 2 2 ⋅ ∂ F ∂ z j 1 − c bot N k + 1 / 2 2 ⋅ F j , k − F j , k + 1 Δ z + + γ z + c top F j , k − 1 − F j , k Δ z − − F j , k − F j , k + 1 Δ z + × − δ bot N k + 1 / 2 2 ⋅ ∂ F ∂ z ( j , k m j ,

where the boundary conditions are

∂ F ∂ x ∞ k = 0 , ∂ F ∂ x j cst k k = f c F j cst k , k , ∂ F ∂ z j 1 = 0 , ∂ F ∂ z j k m j = 0 .

Remarks

Note here the bottom boundary condition is no longer

c h x F x + f 2 N 2 F z − fh x F = 0 , at z = − h x E23

as shown before. Because we adopt a staircase-like topography in our discretized domain, (23) is replaced by a horizontal wall condition and a flat-bottom condition.

In the equation, γ x + , γ x − , γ z + , and γ z − have been multiplied respectively by c cst , c sea , c bot , and c top to handle those points near the boundary. In this case, the discretization at these points have only precision of order 1 if the mesh is not uniform.

With the above F xx and F zz , the governing equation

F xx f 2 + F z N 2 z = 0

is discretized. Collect the coefficients of F j − 1 , k , F j , k − 1 , F j , k , F j , k + 1 , and F j + 1 , k , as correspond to the points 1, 2, 3, 4, and 5 in Figure 2. They are:

F j − 1 , k : 1 Δ x j 1 Δ x − γ x + c cst c sea − c cst − 1 + γ x − c sea ≡ a 1 j k , E24

F j , k − 1 : f 0 2 Δ z k 1 Δ z − c top N k − 1 / 2 2 1 + γ z − c bot − c bot N k + 1 / 2 2 γ z + c top ≡ a 2 j k , E25

F j , k + 1 : f 0 2 Δ z k 1 Δ z + c top N k − 1 / 2 2 γ z − c bot − c bot N k + 1 / 2 2 − 1 + γ z + c top ≡ a 4 j k , E26

F j + 1 , k : 1 Δ x j 1 Δ x + c sea 1 + γ x + c cst − γ x − c cst c sea ≡ a 5 j k , E27

F j , k : a 3 − η j , k ⋅ 1 c , E28

where

a 3 j k = 1 Δ x j c sea ⋅ − 1 Δ x + + γ x + c cst − 1 Δ x + − 1 Δ x − − c cst ⋅ 1 Δ x − + γ x − c sea − 1 Δ x + − 1 Δ x − + f 2 Δ z k c top N k − 1 / 2 2 ⋅ 1 Δ z − + γ z − c bot − 1 Δ z − − 1 Δ z + − c bot N k + 1 / 2 2 ⋅ 1 Δ z + + γ z + c top − 1 Δ z − − 1 Δ z + ,

and

η j , k = 1 Δ x j δ cst f . E29

We would like to make (24)–(28) into a matrix with eigen value 1 c . But here the problem is unusual, since 1 c take values only on boundaries [see the η in (28)], reflecting the fact that the inner product is defined along boundaries. To make it valid throughout the domain, let ε be some small number reflecting the machine precision, say, 10 − 7 , and define a new η such that

η j , k = ε ⋅ max l = 1 5 a l if η j , k = 0 .

Now, the equation with the new η should be the same as the old one up to the machine precision. Write

a ˜ l = a l η , l = 1 , 2 , … , 5 . E30

The discretized equation is then

a ˜ 5 F j + 1 , k + a ˜ 4 F j , k + 1 + a ˜ 3 − 1 c F j , k + a ˜ 2 F j , k − 1 + a ˜ 1 F j − 1 , k = 0 . E31

And the whole problem is now translated into finding the eigenvalue 1 c for a pentadiagonal matrix A , with entries along the five diagonal being a ˜ i , i = 1,2,3,4,5 .

Riad Benzaid

Chapters authored

Related collaborators

Long Yu

Martin Masuelli

Mônica Garcez

Theodoros Rousakis

George Manos

Manal Zaki

Leila Menegthetti

Luiz Carlos Pinto Silva Filho

Habib-Abdelhak Mesbah

Kostas V. Katakalos