Coastal Engineering, 8 (1984) 219232
219
Elsevier Science Publishers B.V., Amsterdam  Printed in The Netherlands
VERIFICATION OF A PARABOLIC EQUATION FOR PROPAGATION OF WEAKLYNONLINEAR WAVES
JAMES T. KIRBY 1and ROBERT A. DALRYMPLE 2
1Marine Sciences Research Center, State University of New York at Stony Brook, Stony Brook, N Y 11794 (U.S.A.) 2Department of Civil Engineering, University of Delaware, Newark, DE 19711 (U.S.A.) (Received May 11, 1983; revised and accepted November 2, 1983)
ABSTRACT Kirby, J.T. and Dalrymple, R.A., 1984. Verification of a parabolic equation for propagation of weaklynonlinear waves. Coastal Eng., 8: 219232. Verification tests of linear wave propagation models for cases of waves propagating through a refractive focus typically overpredict peak amplitudes in the vicinity of the focus. Using a parabolic equation model for the combined refractiondiffraction of weaklynonlinear waves [valid for Ursell number Ur < O(1)], we show that, in a particular experiment, much of the discrepancy between laboratory measurements and linear wave model predictions can be accounted for by including the effect of nonlinearity, which is important in the region of the focus.
INTRODUCTION
Since the introduction of the elliptic mildslope model equation of Berkhoff (1972), great strides have been made in the modelling of propagating surface waves under the simultaneous influence of refraction and diffraction. The introduction of the splitting matrix method by Radder (1979), who derived parabolic equations governing forwardscattered wave motion for waves travelling principally in one direction (here denoted by x), has allowed for the development of models which can encompass large enough physical domains to be of practical importance in engineering applications. The development of new computational techniques for predicting the characteristics of waves in inhomogeneous physical domains has led to the concurrent need for experimental data sets designed to test the predictions of these models. Lozano and Liu (1980) have tested the predictions of a linear parabolic equation against data from the physical experiment of Whalin (1971), where the topography consisted of a series of semicircular steps approximating a tilted cylindrical surface; waves incident on this surface were focused into a cusped caustic. Recently, Berkhoff et al. (1982} have presented data obtained in the vicinity of a focus and cusped caustic
03783839/84/$03.00
© 1984 Elsevier Science Publishers B.V.
220
created b y an elliptic shoal resting on a plane slope, and have compared the data to the predictions of three linear wave models: a refraction scheme involving averaging over bundles o f adjacent rays (Bouws and Battjes, 1982), a parabolic equation model for the scattered incident wave, and an elliptic b o u n d a r y value model for the entire wave field. While the results of each computational model differ in particulars, an important conclusion may be drawn with respect to comparison between linear model predictions and the data set. Linear models uniformly tend to overpredict m a x i m u m wave amplitudes in the vicinity o f focused waves, where wave steepness may become large and nonlinear effects become important. This tendency is clear in the comparison between Berkhoff et al.'s parabolic model results and data. The conclusions to be drawn from a comparison of data and their elliptic m o d e l results are less clear. The elliptic model tends to predict lower amplitudes than the parabolic model in the region of the focused waves; however, the model results are significantly contaminated with waves reflected from the downwave boundary, as evidenced b y the rapid oscillation in amplitude along x
In a recent paper, Kirby and Dalrymple (1983) derived a parabolic equation, in the form of a cubic SchrSdinger equation, governing the complex amplitude A of the fundamental frequency c o m p o n e n t of a Stokes wave. The equation is given by:
2ikCCgAx + 2k(k  ko)CCgA + i(kCCg)xA + (CCgAy)y k(CCg)K' fAI2A = 0
(1)
A simple derivation o f eq. (1) is presented in the Appendix. The corresponding velocity potential for a wave in depth h(x,y) is given by:
221 ¢(x,y,z,t) 
ig
2co o
e A ( x , y ) cosh k ( h + z ) ei(kox_~o0t ) cosh k h
3ko0 cosh 2k(h + z) + complex conjugate  ~ • e2A 2 e 2i(kox  c°°t) 16 sinh 4 k h + complex conjugate + O(e 3)
(2)
where e is the Stokes exPansion parameter. The local wave number k(x,y) and wave frequency co o are related b y the linear dispersion relation: ¢Oo2 = g k tanh k h ( x , y )
(3)
where g is the gravitational acceleration. Note that this formulation neglects interaction with the steady return flow at O(e2). In eqs. (1) and (2), k0 is a reference wavenumber related to the initial conditions of the incident wave; consequently, A must absorb any difference between the reference phase koX and the actual phase. The coordinate x represents the principal propagation direction as in the formulation of Radder (1979). The remaining coefficients are given by: C = ¢o o/k = phase speed
(4)
Cg = ~¢o o/~k = group velocity
(5)
K ' = k 3 (C/Cg)
cosh 4 k h + 8  2 tanh2kh 8 sinh 4 kh
(6)
The linearized version of the parabolic equation (1) is equivalent to the linear mildslope formulation of Radder (1979), which can be verified by substituting the O(e) portion of ~ into Radder's equation. The added nonlinearity provides for the effect of amplitude dispersion and for weak interaction between intersecting waves. Kirby and Dalrymple {1983) obtained eq. (1) b y following the multiple scale expansion of Yue and Mei (1980} b u t allowing for slow [O(e2)] variations in the water depth h(x,y); the reader is referred to those papers for the mathematical details. The nonlinear equation is valid under the same conditions as the Stokes theory, with the principal condition given by: Ur = ( ~  ~ ) / ( k h ) 2 <: O(1)
(7)
where Ur is the UrseU number. For conditions where Ur approaches O i l ) , wave models based on the Boussinesq equations, such as that of A b b o t t et al. (1978) b e c o m e more appropriate. Also, the multiple scale expansion used in obtaining eq. (1) technically requires t h a t v h / ( k J A I ) b e O(e), although this has n o t been found to be a necessary restriction in applications
222 of the linear model. L a b o r a t o r y experiments designed to test the predictions of linear wave models typically satisfy the condition of small Ur in the vicinity of the wave maker. However, as waves propagate into shallower water and are focussed b y b o t t o m irregularities, Ur may increase rapidly, and care must be taken in verifying that Stokes wave theory is valid for the problem being considered. The parabolic model may be solved by a forwardmarching (in x) procedure in analogy to the heattransfer equation. After discretizing a domain ( 0 ~ x ~ Xmax, 0 ~ y ~< Y m a x ) b Y a finitedifference grid, initial conditions A ( x = 0, y) are specified, and the problem then requires lateral boundary conditions at y = 0, Yrnax for each value of x; these are discussed in the following section. Solutions are obtained efficiently by the CrankNicolson implicit scheme. Iteration is required to estimate the factor IAI 2 in the nonlinear term at the updated positions of x; here, we use an explicit intermediate step as in Yue and Mei (1980) in order to obtain an estimate of A at the updated x. In the event that linear results are to be obtained, the cubic term is d r o p p e d and the CrankNicolson scheme proceeds without iteration. EXPERIMENTAL TOPOGRAPHY
AND COMPUTATIONAL
DOMAIN
For the purpose of testing the predictions of the nonlinear model, eq. (1), and of determining the qualitative differences to be found between the nonlinear model and its linear counterpart, we have chosen to reinvestigate the experiment described recently by Berkhoff et al. (1982). This data set is well suited to the present purpose due to the great detail in which data on wave amplitude was obtained over the entire vicinity of a refractive focus, where diffractive and nonlinear effects both become significant. In addition, the Ursell parameter remains of a reasonably small size over the entire domain of interest, as will be shown below, thus indicating that the Stokes theory should be a valid representation of the experiment. Details of the experimental procedure and setup may be obtained from Berkhoff et al. (1982) or the companion report Berkhoff (1982), which includes a photograph o f the experimental wave field; we summarize the important points here. The experimental topography consists of an elliptic shoal resting on a plane sloping b o t t o m with a slope o f 1:50. The plane slope rises from a region of constant depth h = 0.45 m, and the entire slope is turned at an angle of 20 ° to a straight wave paddle. B o t t o m contours are shown in Fig. 1 along with the chosen computational domain, which is indicated by the dashed line surrounding the contours. The offshore boundary of the computational domain is chosen so that water depth is constant along x  0. The initial condition for the wave then corresponds to the uniform wave train generated at the wave paddle; we set: A ( X = 0, y) = A o
(8)
where A0 = 0.0232 m is the amplitude of the incident wave. The wave period
223
T = 1 sec. R e f e r r i n g t o Fig. 1, we establish s l o p e  o r i e n t e d c o o r d i n a t e s {x',y' } which are r e l a t e d t o t h e c o m p u t a t i o n a l c o o r d i a t e s {x,y} a c c o r d i n g t o : x ' = (x  1 0 . 5 ) c o s 20 °  (y  10)sin 20 °
(9a)
y ' = (x  10.5)sin 20 ° + (y  1 0 ) c o s 20 °
(9b)
0 
I
m _ _ ~0.45
I n c i d e n t wove T = sec
x(m
20
0
I 5
I I0
I 15
I 20
y(m)
Fig. 1. Topography and computation domain for experiment of Berkhoff et al. (1982). Dashed lines indicate transects labelled as sections 18 in Fig. 4.
T h e origin {x',y' } = {0,0} c o r r e s p o n d s t o the c e n t e r o f t h e shoal. T h e slope is d e s c r i b e d by: t 0 . 4 5 m,
x' < 5.82 m
(10a)
x ' 1>   5 . 8 2 m
(10b)
h 0.45  0.02 (5.82 + x ' ) m,
T h e b o u n d a r y o f t h e elliptic shoal is given by:
(T
+
=1
and t h e d e p t h in t h e shoal region is m o d i f i e d a c c o r d i n g t o :
(11)
224
h =hslope0.5
[(1
3.75
+0.3
(12)
resulting in a depth h(x' = 0, y' = 0) = 0.1332 m. The lateral boundaries at y = 0, 20 m are open, and transmitting conditions are specified according to: /
b_A/ = ik sin0 A OY~ ~0, 20
(13)
where we have assumed that IAI varies slowly in the ydirection at the boundaries. The wavenumber c o m p o n e n t k sinO at each grid row is evaluated numerically based o n values of A at the previous grid row. Preliminary studies of the relevant parameters in the wave field were performed b y determining Ur at several locations in the wave field corresponding to the initial condition, the value at the shoal crest, and the value at the measured amplitude maximum in the focus, yielding values of Ur = 0.014, 0.290 and 0.213 at each position. Stokes theory should thus be valid t h r o u g h o u t the region of principal interest. COMPARISON OF COMPUTATIONAL MENTS
RESULTS AND EXPERIMENTAL MEASURE
Data from the laboratory experiment of Berkhoff et al. (1982) are available for the labelled transects 1 through 8 indicated in Fig. 1. The computational domain was discretized into square grids (Ax = Ay = grid spacing), and the grid scheme was established so that grid rows coincided with the measurement transects. Grid size was decreased until the point was reached where further reduction did not affect model predictions significantly. The final numerical calculations were performed using a spacing of ~ x = 0.25 m, corresponding to a grid of 87 X 81 rows. Calculations were performed using b o t h the nonlinear model eq. (1) and a linearized version obtained by dropping the nonlinear term, yielding a model equivalent to Radder's (1979) approximation. C o n t o u r lines of wave amplitude relative to incident wave amplitude for the linear and nonlinear results are shown in Figs. 2 and 3, respectively. A comparison o f the figures shows that the nonlinear model predicts a broadening of the region o f focused waves and a resulting lateral displacement of the diffraction fringes away from the centerline of the focused region. Also, the area containing waves with amplitude in excess o f twice the incident wave amplitude is greatly reduced in the nonlinear results in comparison to the linear results. A comparison of Fig. 3 to the measured contours reported in Fig. 2 of Berkhoff et al. (1982) points o u t the improvement in qualitative similarity between c o m p u t e d and measured waves provided by the nonlinear model.
225
/ /
/
/I
.
\ \\
/
0.3
~02~
i
Fig. 2. Amplitude contours relative to incident wave amplitude in region of focused waves, using the linear model.
Results of the nonlinear and linear models are shown in comparison to the experimental data in Figs. 4ah for the labelled transects 18 respectively. In each figure, nonlinear and linear model results are indicated by solid and dashed lines, respectively, while open circles indicate experiment data points. We first consider the results on each individual transect. On section 1, focusing effects have only begun to be apparent, and only a slight difference exists between the linear and nonlinear models, which both agree well With data. Sections 2 and 3 describe the region of the development of the cusped caustic. For both sections (Figs. 4b, c), data typically fall between the predictions o f the two models in the region of maximum amplitude, with the nonlinear model underpredicting the maximum amplitude by about 8% on section 3. However, on sections 4 and 5, where the wave has passed through
226
"0.3 .~
// //
/ z/
~
~ \ \\\
0.2~
Fig. 3. Amplitude contours relative to incident wave amplitude in region of focused waves, using the nonlinear model.
the caustic cusp, the nonlinear model predictions and the data points are in striking agreement, with both the height and width of the central focused region and the shape of side lobes in the diffraction pattern being very well predicted. Agreement between the nonlinear model and the data is also very good along the longitudinal sections 68. On section 7, which is just off center of the axis of the focused region, the nonlinear model predicts both the drop in amplitude in the region of the cusp (with respect to linear model predictions) and the slower decay o f amplitude towards the shoreward end of the transect (large x). On section 6, the nonlinear model quite accurately predicts a sharp drop in amplitude which is apparent in the data but totally absent in the results of the linear model. Berkhoff et al. (1982) interpreted the
227
discrepancy between the linear model and the data in this region as being the result of the proximity of an amphidromic point in the wave phase surface and the linear model's inability to resolve this point. Here, we see that the discrepancy is entirely explained in terms of nonlinear effects and the resulting change in the geometry of the focused region. 2.5
I
,
i
i
i
i
,
(o} Section I
2.0
1.5
1.0
0
0
0
0
0
0
0
O0 

O
~
0
0.5
0.0 5
i
I
i
I
i
6
?
8
9
I0
ill
i
I
!
12
13
14
i
i
i
15
y(m) 2.5
i
i
i
l
i
!
(b) Section 2
2.0
1.5
1.0
o o o
o
o
0.5 0.0
5
I
i
i
i
6
7
8
9
[ I0
I
I
I
I
II
I;~
13
14
i
i
i
15
y(m)
2.5
i
l
i
!
i
i
/o\ 2.0
/
/
~
(c) Section 3
1.5
1,0 o
0.0
5
[
I
I
I
6
7
8
9
I I0
y(m)
Fig. 4.
For
explanation
see
p. 229.
i II
oo
o
I
I
I
12
13
14
15
LD
o
o
~
b
b~
b i
b~
p 0
0 b,
/
t/

J
~1 ~. I
....I
0
D
m
P 0
o
s ~'
P
o
I
o~l i
~ _ i
L~ O0
229 2.5
i
i
i
i
i
i
i
i
i
i
2.0
i.5 1.0 0.5
(g) S e c t i o n 7
0.0 0
I
[
I
I
I
l
I
I
I
I
I
2
:3
4
5
6
7
8
9
I0
(X
2.5
1
i
i
i
10.5)(
i
m}
i
i
i
i
i
(h) Section 8
2.0
1.5 I.O o
o
o
.*"
o
o
0.5 0.0
i
I
I
I
2
3
4
i
I
5 6 (X10.5)(m)
I
I
i
I
7
8
9
I0
Fig. 4. C o m p a r i s o n o f linear and nonlinear m o d e l results to e x p e r i m e n t a l data, sections 1   8 ( a   h , respectively). Solid line = nonlinear results; dashed line = linear results; o p e n circles = data.
Differences between the linear and nonlinear results are less striking along section 8; however, the nonlinear model does reproduce the relatively high amplitudes in the region 69 m downwave of the shoal, in comparison to the linear model. Taken as a whole, it is apparent that the results of the nonlinear model exhibit closer agreement with the experimental data than do the results of the linear model, with the only obvious discrepancy occurring in the region of initial development of the wave focus. It is possible that the discrepancy in this initial region of growing amplitude, seen in sections 2 and 3, is due in part to a failure of the approximate equations to allow waves to focus rapidly enough to obtain a realistic amplitude. This discrepancy is seen to be shortlived, as evidenced by the remainder of the data set.
230 DISCUSSION W e have shown that a parabolic equation governing the fundamental mode of a forwardscattered Stokes wave is capable of producing qualitatively and quantitatively accurate results in comparison to laboratory data in a situation where the conditions for the validity of Stokes theory are approximately satisfied.Further, we have shown that the discrepancies remaining between data and linear wave theory are largely due to the neglect of nonlinear effects rather than to inaccuracies in the modelling technique. It is not known at present h o w the Stokes wave model will perform under conditions beyond the range of its strict validity. This point should be examined in the future with the objective of delineating the limits of applicability of such models and determining the possible range of overlap with models based on the Boussinesq equations. ACKNOWLEDGEMENT This study was sponsored by the Office of Naval Research, Sciences Program through a grant to the University of Delaware.
Coastal
APPENDIX The nonlinear equation (1) governing the complex amplitude A ( x , y ) may be obtained through a simple twostep procedure. A more formal derivation may be found in Kirby and Dalrymple (1983). A hyperbolic equation governing the time
~ • ( C C g ~ 7 ) + COo2
Cg (1  C ) ~ = 0
(14)
where ~0 is given by eq. (3) and ~ is the horizontal gradient { ~ / a x , a / ~ y }. For the case of constant depth and, without loss of generality, propagation in the x~lirection, eq. (14) is satisfied by linear plane waves of the form: ~? = A e i(kx  ~o ot) ;
A
constant
(15)
We now assume a modified form of eq. (14) containing an unknown nonlinear term ~(~), given by: a 2~ ~t 2
Cg
v(VCgv'V) +co02 (1 C)V +~(~) =0
(16)
In the constant depth case, this equation is to be satisfied exactly by the leading component of a nonlinear Stokes wave given by:
231
= A e i(kx  ~ t) ,"
A constant
(17)
where co and coo are related t o O(klAI 2) by (Whitham, 1967): co2 = COo2 (1 + CgK']AI2/COo)
(18)
Substituting eq. (17) into eq. (16) leads t o the result: ~(~) = (co:  ¢o0:)~ = ¢o0CgK'[Al:~
(19)
where K' is given by eq. (6). We n o w take the alternate viewpoint t ha t t he wave is allowed to travel at small angles to the x
IA'I = IAI
(20)
where A' must absorb t he effects o f angle o f propagation and nonlinearity. Substituting eq. (20) into eq. (16) with [(~) given by eq. (19) leads to a parabolicelliptic model including diffraction effects in the xdirection, given by: 2 i k C C g A ' x + i ( k C C g ) x A ' + v " (CCg v A ' )  k C C g K ' IA'12A ' = 0
(21)
Assuming t h a t A' varies m or e slowly in the xdirection than the y~lirection, since the principal variation appears in the phase ( f X k d x  coot), and introducing th e phase shift: A' = Ae i(kox  fXk(x)dx)
(22)
gives eq. (1). A somewhat different derivation is given in Liu and Tsay (1984).
REFERENCES Abbott, M.B., Petersen, H.M. and Skovgaard, O., 1978. On the numerical modelling of short waves in shallow water. J. Hydraulic Res., 16(3): 173204. Berkhoff, J.C.W., 1972. Computation of combined refractiondiffraction. Proc. 13th Int. Conf. Coastal Eng., Vancouver. Berkhoff, J.C.W., 1982. Refraction and diffraction of water waves: wave deformation by a shoal. Rept. W. 154VIII, Delft Hydraulics Laboratory. Berkhoff, J.C.W., Booij, N. and Radder, A.C., 1982. Verification of numerical wave propagation models for simple harmonic linear waves. Coastal Eng., 6: 255279. Booij, N., 1981. Gravity waves on water with nonuniform depth and current. Rept. 811, Dept. of Civ. Eng., Delft Univ. Technology.
Bouws, E. and Battjes, J.A., 1982. A Monte Carlo approach to the computation of refraction of water waves. J. Geophys. Res., 87(C8): 57185722. Kirby, J.T. and Dalrymple, R.A., 1983. A parabolic equation for the combined refraction
232 Liu, P. LF. and Tsay, TK., 1984. Refractiondiffraction model for weaklynonlinear waterwaves. J. Fluid Mech., in press. Lozano, C. and Liu, P. LF., 1980. Refractiondiffraction model for linear surface water waves. J. Fluid Mech., 101(4): 705720. Radder, A.C., 1979. On the parabolic equation method for waterwave propagation. J. Fluid Mech., 95(1): 159176. Whalin, R.W., 1971. The limit of applicability of linear wave refraction theory in a convergence zone. R.R. H713, U.S. A r m y Corps of Engineers, Waterways Experiment Station, Vicksburg. Whitham, G.B., 1967. Nonlinear dispersion of water waves. J. Fluid Mech., 27(2): 399412. Yue, D. KP. and Mei, C.C., 1980. Forward diffraction of Stokes waves by a thin wedge. J. Fluid Mech., 99: 3352.