shallow water equations for channels with power law geometry∗
2
Geir Pedersen †
3 4
Abstract. The present study was originally motivated by some intriguing exact solutions for waves propagating 5
in nonuniform media. In particular, for special depth profiles reflected waves did not appear and ray 6
theory became exact. Herein, geometrical optics is employed to obtain asymptotic series for waves of 7
general shapes in non-uniform narrow channels, within the framework of linear shallow water theory.
8
While being kept simple, the series incorporate higher order contributions that may describe the 9
evolution of waves with high accuracy. The higher orders also contain a secondary wave system.
10
For selected classes of geometries and wave shapes explicit solutions are calculated and compared 11
to numerical solutions. Apart from the vicinity of shorelines, say, higher order expansions generally 12
may provide very accurate approximations to the full solutions. The asymptotic series are analyzed 13
for different wave shapes and are found to be convergent for cases where the basic wave profiles have 14
compact support (finite length). A number of new, closed form, exact solutions are also found. The 15
asymptotic expansion is put into a context by employing it for the transmission of waves from a 16
uniform channel section into a nonuniform one. Additional results and side topics are presented in 17
a supplement.
18
Key words. Ray theory, asymptotic expansions, shallow water waves, channels 19
AMS subject classifications. 68Q25, 68R10, 68U05 20
1. Introduction. For wave propagation in uniform media solutions are available as uni-
21
form, harmonic wave trains. When the the medium is slowly varying such solutions may be
22
generalized to a slowly varying wave train by means of ray theory. Variants of ray theory are
23
employed in a vast number of articles and are included in most textbooks on waves. Such the-
24
ories may be developed by direct invocation of dispersion relation and energy expressions, by
25
perturbation schemes, most notably the WKB method, [19,13,2,4] or a variational approach
26
[35]. Ray theory is a useful vehicle for the conception of wave phenomena such as focusing,
27
amplification, wave trapping and shadowing. It is important for understanding the behaviour
28
of water waves in a bathymetry and propagation of acoustic and seismic waves in the layered
29
atmosphere and earth, respectively [27,1]. On the other hand, standard ray theory has severe
30
limitations concerning quantitative description of waves in even moderately complex condi-
31
tions due to singularities at caustics and lack of diffraction. Ray theory is mostly applied to
32
linear harmonic modes which then may be combined in a Fourier integral to describe waves
33
of more general shapes. There are a few examples of similar approaches applied to specific
34
nonlinear wave forms, such as shocks [34,35,3] and solitary waves [14,20,16,17,22].
35
In shallow water theory waves of normal incidence on laterally uniform beaches lead to
36
mathematical formulations with one spatial dimension in addition to time. The depth may
37
then be written h(x) where x is the coordinate along an axis parallel to the direction of
38
∗February 3, 2021
†Department of Mathematics, University of Oslo, PO box 1053, 0316 Oslo, Norway ([email protected],https://
www.mn.uio.no/math/english/).
1
wave advance. Some depth profiles then yield equations of standard types with analytic
39
solutions. For an inclined plane (h=αx) periodic waves may be expressed in terms of Bessel
40
functions. In this case the standard WKB solution, equivalent to asymptotic expansion of
41
the Bessel functions, predicts that the amplitude is proportional to x−14. While the inclined
42
plane generally is conceived as the most basic one, other configurations are simpler and allow
43
for a more transparent analysis, in particular for the linear case. Simple exact solutions are
44
known for theh=αx43 profile for which the surface elevation evolves according to ray theory.
45
[33] exploited this to obtain explicit solutions for waves generated by slides on a slope, in a
46
corresponding form as in constant depth. Later, [10] retrieved theh=αx43 profile as the only
47
one where ray theory for periodic waves gave an exact result for the surface elevation. Then
48
the exact solution for the surface elevation for a wave of any form followed.
49
Linearized shallow water theory for long waves in narrow channels also leads to equations
50
in one spatial dimension for the surface elevation and the longitudinal velocity, while the
51
geometry may be described in terms of the channel width and the effective depth, which will
52
be defined later. By inclusion of finite amplitude and non-hydrostatic effects these may be
53
generalized to Boussinesq or KdV type equations [26,25,28,29].
54
Within linear shallow water theory there is a particularly simple exact solution for which
55
changing widths and depths of a channel counterbalances to allow wave propagation without
56
amplification. Inclined channels with parabolic cross sections allow for other exact solutions
57
[5, 6, 7]. For this case the linear solution has a similar structure as for the plane, normal
58
incidence, case with the h= αx43 profile. Furthermore, [8] applied the WKB method to the
59
class of channels for which both the transverse shape and lengthwise slope were defined by
60
power functions in x. They found a set of combinations of shapes and slopes, including those
61
mentioned above, that made the expression for the surface elevation exact. Then, surface
62
excursions with finite wavelength can propagate without the evolution of a secondary system
63
of surface elevations. This may seem to indicate absence of any form of diffraction or reflection.
64
However, these surface elevations come with velocity fields that have modified shapes relative
65
to the ray solution and a trailing region of constant volume flux. The cases where ray theory is
66
exact can be conceived as cases for which the variable coefficient equations may be transformed
67
to a standard wave equation with constant coefficients. A slightly more general approach
68
transforms the shallow water equations to a constant coefficient Klein-Gordon equation [15,
69
9, 24]. However, this and a few other exact solutions of the shallow water equations with
70
non-constant depth, such as the oscillations in parabolic basins [32], are less closely related to
71
the present work.
72
The first objective of the present article is to obtain a wider perspective on the exact ray
73
solutions from the literature. Are they unique with properties that are distinct in relation
74
to solutions which cannot be exactly expressed in this manner? Higher order expansions are
75
designed for investigation of the cases where low orders do not provide exact solutions. Since
76
the medium is non-dispersive we employ a direct formulation of ray theory for general wave
77
shapes; notably single-crested pulses. This is preferred over application through the Fourier
78
transform, which, in the present case, is a cumbersome detour at best. To obtain solutions
79
that are quite explicit, and hence open to detailed analysis and interpretation, a selection of
80
particular channel geometries and wave shapes is made. This include also cases where uniform
81
and non-uniform channel geometries are joined and where the transmission and reflection at
82
the junction may be analyzed by employing the asymptotic expansions. The validity and
83
accuracy of the asymptotic solutions are checked through comparison with numerical solutions
84
and by the derivation of a theoretical error constraint. In the supplement additional details,
85
solutions and plots are presented, as well as secondary topics such as conservation properties
86
and implications of the transmission properties for calculation of runup on a beach.
87
2. Formulation. Exempting curved channels from this study we introduce the coordinates
88
x, y and z in the longitudinal, transverse and vertical directions respectively. A channel
89
geometry is defined in figure1. The depth is expressed ash(x, y), while the equilibrium width
90
at the surface is B. The surface elevation and longitudinal velocity are denoted by η and
91
u, respectively. Dimensionless expressions are used throughout this article. The depth scale
92
is hc and a characteristic length xc, associated with the change of the geometry, is used for
93
scalingx in the hydrodynamic equations. Withg as the acceleration of gravity the time scale
94
is xc/√
ghc and the ratio between the surface and velocity becomes ηc/uc = p
hc/g. When
95
convenient, the scaling may give a reference positionx0 = 1 at t= 0 withB(x0) = ¯h(x0) = 1
96
and, say, a unitary amplitude. A shorter spatial scale,xc/κ, whereκ is large, will be used for
97
wave lengths.
98
2.1. Linear shallow water waves in a channel. The equations are described in, for in-
99
stance [18]§185. Within the shallow water approximationηanduare uniform in the span-wise
100
directions. In the linear approximation the surface width, B, is regarded as fixed in time. It
101
is convenient to introduce the effective depth, ¯h, which is defined as the equilibrium cross
102
section area, divided toB. For the dynamic cross-section area we then find ¯hB+Bη and the
103
momentum and continuity equations may be expressed according to
104
(2.1) ut=−ηx, ηt=−1
B(B¯hu)x,
105
where indices indicate differentiation. Elimination of u yields
106
(2.2) ηtt= 1
B B¯hηx
x= Bx
B hη¯ x+ (¯hηx)x.
107
WhenBis constant we retrieve the standard shallow water equations for the plane case. For a
108
rectangular channelB and ¯hare the width and equilibrium depth of the channel, respectively.
109
It is noteworthy that channels of any shape correspond to a rectangular one, with depth
110
equal to ¯h and width B, by means of (2.1) or (2.2). Numerical solution of (2.2) is quite
111
straightforward and is briefly discussed in sectionSM5.1. However, a proper comparison with
112
the asymptotic solutions may require a high resolution.
113
3. The ray expansion. Usually ray and transport equations are formulated for harmonic
114
modes. The evolution of a more general wave form may then be found by inversion of Fourier
115
type integrals. However, such inversions are cumbersome. Herein we exploit the lack of
116
dispersion for shallow water waves to employ a method that is directly applicable to waves of
117
any shape, while it also extends to higher orders. In general, such an approach may not be
118
applied to dispersive and nonlinear waves.
119
y z
h(x, y) B
h
Figure 1. Definition sketch. Cross section of a generic channel. The side walls and bottom of the equivalent rectangular channel is indicated by the thin vertical and horizontal lines.
3.1. Leading order asymptotic approximation in a channel. Waves in slowly varying
120
channels may approximately retain their shape, while the amplitude and length change. From
121
[18]§ 185 a leading order approximation is adapted as
122
(3.1) η=A(x)F(Θ), u=∓h¯−12A(x)F(Θ),
123
whereF is an arbitrary shape function and
124
(3.2) Θ =κ(t±τ), τ =
x
Z
x0
¯h−12 dˆx, A=σB−12¯h−14.
125
Here,τ is the propagation time from a reference positionx0with the local wave speedc0=√¯h,
126
while the identical appearance of t in Θ is due to the invariance of the wave period. The
127
relation forAin terms ofB and ¯hfollows from energy conservation and is often referred to as
128
Green’s law. WhenB¯h12 is decreasing in the direction of wave advance the wave is amplified.
129
In the opposite case the amplitude attenuates. The constant σ can be made equal to unity
130
by the appropriate scaling (see section 2), but is kept in the relations for now. The factor κ
131
governs the wavelength and a large value makes the wave short in comparison to the length
132
scales of the geometry.
133
The relations (3.1) will be reproduced in the derivation of higher order theory. Then,
134
AF, renamed to A0F0, is the leading order approximation forη and will be referred to as the
135
principal wave.
136
3.2. Higher order asymptotic approximations. With A and F renamed to A0 and F0,
137
respectively, a higher order approximation may be expressed as an asymptotic series on the
138
form
139
(3.3) η∼A0(x)F0(Θ) +κ−1A1(x)F1(Θ) +κ−2A2(x)F2(Θ)...
140
Here κis used as an ordering parameter, but the magnitudes of terms could have been made
141
equally transparent by observing the number of derivatives throughFj in relation to derivatives
142
of ¯h, B,Aj etc. The expansion (3.3) is a special case of geometrical optics as defined in [34]
143
§ 7.7. If harmonic shapes,Fj = (−i)jexp(iΘ), are used (3.3) will become the standard WKB
144
technique. However, herein we will study other shapes. The asymptotic series may then have
145
quite different properties than that of the WKB expansion.
146
Substituting (3.3) into the single PDE forη, (2.2), and sorting by magnitude we obtain
147
0 =κh
∓(A0B)−1
B¯h12A20
xF0′i +κ0h
−B−1 BhA¯ 0,x
xF0 ∓(A1B)−1
Bh¯12A21
xF1′i +...
148
where indices after commas indicate differentiation. The O(κ2) part vanishes due to τx =
149
±¯h−12, wheras the O(κ) balance reproduces the leading order approximation
150
(3.4) A0=σB−12h¯−14.
151
From O(κ0) we have a balance between terms with factors from two subsequent orders in
152
(3.3). Since Θ and x must be treated as independent variables these orders are annihilated
153
by the recursion (j = 0,1...)
154
(3.5) Fj+1(Θ) = Z Θ
−∞
Fj( ˆΘ)dΘ,ˆ A−1j+1
B¯h12A2j+1
x=∓ B¯hAj,x
x.
155
When wave propagation to the left (Θ =κ(τ+t)) is assumed for simplicity (3.4) may be used
156
with (3.5) to give the explicit amplitude recursion
157
(3.6) Aj+1=−1
2¯h12Aj,x+1 2σ−2A0
x
Z
xr
¯hBA0,xAj,xdˆx.
158
The reference position, xr, may be chosen freely for each j. Changing the value of an xr 159
corresponds to adding a constant timesA0toAj. This ambiguity corresponds to a redefinition
160
of the principal wave at each level in the expansion and is exploited in section 5.
161
Employing both equations in (2.1) we find the velocity
162
(3.7) u∼ −h
h¯−12A0F0(Θ) +κ−1(¯h−12A1+A0,x)F1(Θ) +κ−2(¯h−12A2+A1,x)F2(Θ)...i
+Qf(B¯h)−1,
163
where Qf is a constant. On its own, the last term, combined with η = 0, is a solution of
164
the linear shallow water equations. Hence, it is unrelated to the principal wave and yields
165
a volume flux, of constant strength, along the channel. Accordingly, we set Qf = 0. The
166
omitted term should not be confused with the similar, but local, velocity field in section 4.2.
167
It is now assumed that the principal wave, defined through F0 has compact support in
168
the sense that F0 is zero for Θ <Θ0 and Θ >ΘL. If the total integral ofF0 is nonzero, as
169
for a principal wave of elevation, then F1 is constant for Θ >ΘL, while Fj is a polynomial
170
of degree j−1. When F0 instead decays exponentially for large Θ the Fj functions still
171
have a dominant polynomial behaviour for large Θ. Hence, there is a a trailing system with
172
surface excursion and fluid velocity behind the principal wave. The trailing system is needed
173
to balance the volume and momentum in the principal wave and does also drain energy from
174
it (details in sections SM7, SM2.4). Similar trailing systems are described for solitary waves
175
that shed volume due to geometry or lateral non-uniformity of the crest [16,17,22].
176
3.3. Waves propagating in the positivex-direction. When waves are propagating toward
177
increasingxan asymptotic expression may be readily deduced from (3.2), (3.3), (3.5) and the
178
momentum equation in (2.1)
179
(3.8) η∼
∞
X
j=0
(−1)jκ−jAjFj( ˆΘ), u∼¯h−12A0F0( ˆΘ) +
∞
X
j=1
(−1)jκ−j
¯h−12Aj −Aj,x Fj( ˆΘ),
180
where all quantities, except the modified phase ˆΘ =κ(t−τ), are defined as previously.
181
4. Solutions for selected profiles. Explicit expressions for the terms in (3.3) may only
182
be obtained for selected classes of geometries and principal wave shapes. Geometries defined
183
by single power functions or exponentials yield simple recursions for the amplitudes A0, A1 184
etc. When ¯h is monotonic it can be used as spatial variable and (3.6) rewritten accordingly.
185
Then, certain relations between ddx¯h and ¯h yield new explicit recursions for the amplitudes,
186
in addition to those mentioned above. For brevity, only the power function geometries are
187
analyzed herein. The other options are briefly outlined in the sectionsSM3.4and SM6of the
188
supplement.
189
Two types of principal wave shapes are employed; one whereηtapers off exponentially at
190
the outskirts and another one with compact support (strictly finite length). There are some
191
striking differences concerning the asymptotic series for these two groups.
192
4.1. Amplitude recursions for geometries defined through power functions. Non-planar
193
geometries defined through power functions have attracted attention in the literature, partic-
194
ularly in relation to exact solutions of shallow water equations [8].
195
When both ¯h and B are simple power functions it may be assumed that asymptotic
196
expansions do exist with amplitudes as simple power functions as well.
197
(4.1) ¯h=h0xα, B =B0xβ, Aj =Cjxqj.
198
Insertion of these expressions in (3.6), withxras 0 or∞according to the power in the integral,
199
provides us with recursion formulas forCj and qj. The different choices of xr are required to
200
obtain Aj as a single power, but will lead to a jump in the solution when α crosses 2. For
201
α6= 2, we may write
202
(4.2) qj =−p−µj, Cj+1 =νjbjCj,
203
wherep= 14α+ 12β,µ= 1− 12α and
204
νj = 1 2h
1
02µ(j+ 1), bj =
1− µ−p
(j+ 1)µ 1− p (j+ 1)µ
.
205
Forj → ∞the factorbj approaches unity, whereas the factorνj causesCj to grow fast withj.
206
However, the discussion of the implications for the series expansion forη awaits specification
207
of the shape functionsFj.
208
Expressions can be made more compact by introducing τ(x) = h−
1
0 2µ−1xµ, which is the
209
travel time from 0 to x if α < 2, and the travel time from ∞ to x if α > 2. For waves
210
propagating toward decreasing x the amplitudes Aj and the phase function, defined in (3.2),
211
may be written
212
(4.3) Aj =j! (2τ(x))−j
j
Y
i=1
bj
!
A0, Θ =κ(τ(x)−τ(x0) +t),
213
wherex0 is the reference position for which Θ = 0 at t= 0.
214
A change in the properties of the asymptotic expansion occurs at α = 2. For smaller
215
values ofαthe shoreline,x= 0, is reached in a finite time and theAj are expressed as inverse
216
powers of x. For α > 2 the shoreline is not reached (in linear theory) and the Aj inherit
217
positive powers of x for sufficiently large j. For α = 2 an alternative recursion, involving
218
logarithmic terms, can be devised. This, and a discussion of solutions forα values close to 2,
219
are found in section SM2.
220
4.2. Exact solutions. The most apparent example of a ray solution that becomes exact is
221
B¯h12 = const.Then A0 in (3.4) is constant, such that the series (3.3) and (3.7) truncate after
222
the first term, and η=A0F0,u=−¯h−12A0F0 form an exact solution. Correspondingly, when
223
the right hand side of the amplitude recursion in (3.5) vanishes for j= 0 the amplitudesA1,
224
A2 etc. become zero andη =A0F0 is still an exact solution, even thoughA0 is non-constant.
225
A more general group of exact solutions is found whenA1 is required to be constant, but not
226
necessarily zero. Then, A2 etc, becomes zero and η =A0F0+κ−1A1F1 is an exact solution.
227
For constantA1 and j= 0 the rightmost relation in (3.5) is integrated to give
228
(4.4) A0,x=−
A1¯h−12 +Dσ2(B¯h)−1
=−h¯−12 A1+DA20 ,
229
whereDis a arbitrary constant and (3.4) has been invoked. ForA1 = 0 equation (4.4) yields
230
(4.5) A0(C+Dτ(x)) = 1 or, equivalently B= ¯h−12
C+D
Z ¯h−12dx 2
,
231
where C is another constant. For D = 0 the case B¯h12 = const. is reproduced, while ¯h =
232
h0 = const.yields quadratic width B = (C+Dh−
1
02x)2. For channels of constant width (4.5)
233
is readily integrated again to ¯h=h0(x+ const.)43, which is discussed in [33,10].
234
When A1 = 0 the surface elevation is completely described by the principal wave and
235
η approaches zero behind it, provided F0 approaches, or becomes, zero for large arguments.
236
Still, combining (4.4) with (3.7) we find u = −¯h−12A0F0−Dσ2(κ¯hB)−1F1 which implies a
237
trailing system with a volume flux of constant strength q = Dσ2F1(∞)/κ. The lack of a
238
surface elevation in the trailing system has been interpreted as lack of reflection. Through
239
comparison with solutions for a uniform channel we may well construe another interpretation.
240
In a uniform channel the solution may be written η = A0F(t+ ¯h−12x) +B0G(t−¯h−12x),
241
u = −¯h12A0F(t+ ¯h−12x) + ¯h12B0G(t−h¯−12x), where ¯h, A0 and B0 now are constant. A
242
volume flux of constant strength q may then be obtained by choosing G = F = const. and
243
−A0 =B0 = 12q/(¯h12F). In view of this we may describe the trailing system as combination
244
two progressive waves; one following the principal wave and one being a reflection. In a
245
nonuniform channel this argument is debateable, not least since there is no obvious definition
246
of pure unidirectional wave. Be that as it may, it is more clearly important that the presence
247
of the trailing system leads to an irreversible energy loss in the principal wave (more details in
248
sectionSM7). Moreover, as will be demonstrated in section5.2.2, the trailing volume flux will
249
couple to reflections from an apex combining an uniform and a non-uniform channel section.
250
For nonzeroA1 the equation (4.4) integrates to expressions
251
(4.6) A0 =C−A1τ, A0 =K−1tan(C−A1Kτ), A0=K−1tanh(C−A1Kτ), K = s
D A1 ,
252
for D = 0, D/A1 > 0 and D/A1 < 0, respectively. The above relations may readily give
253
B expressed in terms of ¯h. If F0 has compact support there is a constant surface elevation
254
in the trailing system. When D = 0 (3.7) and (4.4) imply u = −h¯−12A0F0. The water is
255
then quiescent in the trailing system. On the other hand, there is volume flux of constant
256
strength there ifDis nonzero. When the width is constant the leftmost relation in (4.6) gives
257
h=h0(x+ const.)4.
258
4.2.1. Exact solutions for geometries defined by power functions. [8] studied channels
259
where the geometries were defined by power functions and discovered cases where the surface
260
elevations from the basic ray solution (3.1) were exact. Below we will show that these are
261
only the first in a sequence of exact solutions of increasing complexity. (4.2) and (3.5) imply
262
that the asymptotic series will truncate after termn(Cn+1 = 0), providedbn= 0. This gives
263
two families of exact solutions,
264
(4.7) α(i)n = 4n+ 4−2β
2n+ 3 , α(ii)n = 4n+ 2β 2n−1 .
265
In both casesαnapproaches 2 asn→ ∞. Moreover,α(i)n <2 andα(ii)n >2 for alln. Naturally,
266
a truncated series means that the remaining finite sum is an exact solution.
267
For the first family we have cancellation after a single term in (3.3) when α = 43 −23β.
268
Forβ = 0 the case α= 43 is reproduced. Another interesting case is α= 1, β = 12 which may
269
correspond to an inclined channel with parabolic cross-sections [5]. The series (3.7) foru will
270
contain two terms in these cases. For β= 0 we find α(i)1 = 85. Then the exact expressions for
271
η anduhave two and three terms, respectively. For largernthe solutions may contain strong
272
modifications of the principal wave (examples will be given in sec.4.4).
273
For the second group there is cancellation aftern= 0 only for the trivial case ¯h−14B−12 =
274
const. Forβ = 0 we findα(ii)1 = 4. While (3.3) now yields a two term expression forη, only a
275
single term remains for u. For all αn(ii) the factor ¯h−12An+An−1,x in (3.7) vanishes, making
276
the series foruone term shorter than that ofη. Forn= 1 and anF0 that either has compact
277
support, or is tapering off exponentially at the outskirts, the surface elevation approaches a
278
constant for large Θ.
279
It is stressed that the truncation of the series (3.3) does not imply that these solutions
280
are of a different nature from those with a full series. In particular, energy is not better
281
preserved in the principal wave, as shown in section SM7.1. Moreover, the variation of the
282
wave characteristics withα and β is continuous.
283
Y−1
Y0
Y1
Y2
Y8
Θ Y
P0
P1 P2
P3
P8
Θ Pj(2)
Figure 2. The shape functions. Left panel: hyperbolic function shapes. Right panel: piece-wise polynomials forM = 2. Do observe the different extent (inΘ) of the functions in the panels.
Exact solutions, similar to those described above, are described in section SM6.2 for
284
geometries defined through exponentials.
285
4.3. Shape functions for the surface elevation. As mentioned in section3.2 our expan-
286
sion (3.3) becomes equivalent to the WKB method when, instead of a pulse, we invoke a
287
periodic shape Fj = (−i)jexp(iΘ). Then |AjFj| behaves like const.×(√
h0|µ|κ−1)jj!. This
288
gives the typical asymptotic and divergent series generally associated with WKB expansions.
289
Below two different shapes are investigated. Firstly, F0 is given by hyperbolic functions
290
that vanish exponentially at either side. Secondly, F0 are defined by piece-wise polynomials
291
and has compact support. The two options give series with quite different properties. Also
292
other shapes have been used, but these gave little additional insight and are omitted.
293
4.3.1. Thesech2 shape. A natural choice for the shape of the principal wave isF0(Θ) =
294
aY0(Θ) (single crest) or F0(Θ) =aY−1(Θ) (so called N-wave, [31]), wherea is an amplitude
295
and Yj are hyperbolic functions defined as
296
(4.8) Y−1(Θ) =−2 cosh−2Θ tanh Θ, Y0(Θ) = cosh−2Θ,
Y1(Θ) = tanh(Θ) + 1, Y2(Θ) = log(2 cosh(Θ)) + Θ,
297
where each is the integral of the former. The subsequent functions,Y3etc., cannot be expressed
298
that simply, but may be obtained by numerical integration (see appendixSM5.2 for details).
299
The option F0 =aY0 implies Fj = aYj, while F0 =aY−1 leads to Fj = aYj−1. In the latter
300
case a trailing system will appear first to order κ−2 in (3.3). In constant depth aY0 is the
301
leading order approximation to a solitary wave whenβ = 0 andκ=q
3a 4¯h3
xc
hc. However, we are
302
working with linear shallow water equations and no relation between κ and the wave height
303
needs to be imposed. Shapes are depicted in figure 2. For Θ < 0 the Yj may be expanded
304
jmin
js
x j
0
5
10
x η
Figure 3. Properties of the asymptotic series for κ= 3,h0 = 1, α= 1, x0 = 1 andt= 0.75. Left panel:
Number of minimum term, as defined in the text. Right panel: partial sums for η with different number of terms. When 10 terms are retained the divergence is becoming visible nearx= 0.
according to
305
(4.9) Yj = 22−j
∞
X
n=1
(−1)j+1n1−je2nΘ, j ≥0.
306
For Θ ≪ 1 we have |Yj+1/Yj| ≈ 12 and the series for η again diverges since Aj grows as j!
307
(see (4.3)). The smallest term in the series (3.3) for Θ≪1 is then found for js ∼4κ|τ(x)|=
308
4κxµ/(√
h0|µ|). For illustration the caseκ= 3,h0 = 1,α= 1,β = 0, andx0 = 1 is presented
309
in figure3 for the time when the wave is about to run up on the beach. The number of the
310
minimum term, jmin, increases with x in the front of the wave, as indicated by the formula
311
forjs. Moreover,jmin continues to increase at the rear of the principal wave and in the tail.
312
Depth profiles withα >2 yields ajminthat decreases withx(results not shown). The different
313
behaviours ofjmin forα <2 andα >2 comply with the discussion in section SM2.1.
314
4.3.2. Piece-wise polynomial shape functions. An interval J, corresponding to Θ0 ≤
315
Θ≤ΘL, is divided inL sub-intervals. In each sub-interval F0 may be represented as a poly-
316
nomial in Θ. This will define a principal wave which is truly confined to J, and which may
317
inherit continuous derivatives of any order provided the local polynomials are chosen accord-
318
ingly. Moreover, the integration to obtainFj is straightforward, although the expressions may
319
become lengthy. Naturally, Fj, from some positivej and on, will be nonzero for Θ>ΘL.
320
In appendix Ait is shown that if F0 is any piece-wise polynomial with compact support,
321
and Aj is given through (4.2), the series (3.3) converges if (still µ= 1−12α)
322
(4.10) K≡ Θ−Θ0
2κ|τ(x)| = |µ|
2κh012x−µ(Θ−Θ0)<1,
323
and diverges if K >1.
324
When the wave approaches the shore (x = 0) the front of the wave, xf(t) is given by
325
Θ(xf, t) = Θ0. By means of (4.3) Θ−Θ0 is then expressed byxandxf and (4.10) is rewritten
326
as
327
(4.11) K = 1
2
1−τ(xf) τ(x)
= 1 2
1−xf
x µ
<1.
328
Forα <2 the series (3.3) thus converges everywhere until the wave front reaches the shoreline.
329
Stopping the summation of the series after AnFn we find a residue −κ−nB−1(B¯hAn,x)xFn 330
which goes to zero as n → ∞. Provided that the wave equation on the form (2.2) is well
331
posed, in the sense that vanishing residuals imply vanishing errors (demonstrated for the L2
332
norm in section SM8), the converged series is then an exact and complete solution. It is also
333
noteworthy that the convergence is independent of κ, meaning that the solutions are valid
334
regardless of wave length and hence go beyond what is normally conceived as ray theory. Since
335
all continuous functions with compact support may be approximated to arbitrary accuracy by
336
piece-wise polynomials, it follows that these properties apply to allF0 with compact support.
337
The above results are obtained with the amplitude recursion formula (3.6) for which wave
338
propagation to the left, interpreted as Θ =κ(τ+t), was assumed. However, this interpretation
339
is strictly valid only in the limitκ→ ∞. For finite κthe variables xand Θ in the convergent
340
series (3.3) are intrinsically mixed and the sum is η(x,Θ), which is equivalent to η being a
341
function ofxandt. Hence, (3.3) may contain a reflected wave. What may be seen as surprising
342
is the fact that the possible reflected wave is inherent in an expansion of geometrical optics
343
like (3.3). On the other hand, the notion of no reflection can be questioned even for the exact,
344
closed form solutions in4.2.
345
The front of the incident wave reaches the shore at time tc = τ(x0) + Θ0/κ. For larger
346
times the front of the wave that is reflected from the shore is at xR, where τ(xR) = t−tc.
347
When Θ−Θ0 is expressed in terms ofx and xR the criterion (4.10) becomes
348
(4.12) K = 1
2 +1 2
τ(xR) τ(x) = 1
2+1 2
xR x
µ
<1.
349
Hence, the series (3.3) will diverge in the region, adjacent to the shoreline, which is affected
350
by a reflection. Outside this region, the series still converges.
351
For α > 2 equation (4.11) implies that (3.3) converges when x < 3α−22 xf. This region
352
decreases withα and diminishes as the wave approaches the shoreline.
353
In some subsequent examples F0(Θ) = P0(M) will be used, where P0(M) is a simple, nor-
354
malized, piece-wise polynomial, with continuous derivatives through orderM −1:
355
(4.13) P0(M)= 4M(Θ−1
2)M(Θ + 1
2)M for −1
2 <Θ< 1 2,
356
and P0(M) = 0 elsewhere. The integrals Fj =Pj(M) are given by (A.3) combined with (A.8).
357
A selection of Pj(M) is depicted in figure 2.
358
4.4. Evolution of wave on a slope; examples. The partial sum of (3.3) where the terms
359
up to, and including,O(κ)−nare retained is now denoted byηn. In a boundary value problem
360
on the interval xa ≤ x ≤ xb the approximation ηn may then be used for initial conditions
361
at t =t0, as well as boundary conditions at x = xa and x = xb. The computed solution is
362
denoted as ηn:numerical. Then, the temporal growth of the error ∆ηn = ηn:numerical−ηn will
363
give an indication of the accuracy ofηn. Naturally, this requires that the discretization error
364
of the numerical solution is much less than η−ηn. The geometry is normalized as to give
365
h0 = 1 and x0 = 1. The initial time,t0, may be chosen.
366
For the inclined plane (α = 1,β = 0) the first few terms of (3.3) read
367
(4.14) η
C0 ∼x−14
F0(Θ) + 1
16κx12F1(Θ) + 9
512κ2xF2(Θ) + 147
8192κ3x32F3(Θ) +...
,
368
where Θ =κ(2√
x−2 +t) and errors of order κ−4 are implicit.
369
For n = 0 and κ = 4 figure 4 shows errors up to 2.3% at a time when the wave has
370
reached the left boundary xa = 0.05. In fact, the trailing system, which appears in the
371
higher orders in (3.3), gives an elevation that is larger than the errors observed. However,
372
the trailing system is modified when η0 is used for the boundary condition (more details in
373
section SM3.3). For n= 6 and x >0.05, we have errors as small as |∆η6|<2.2·10−5. The
374
norm L2 ={(xb−xa)−1Rxb
xa(∆η)2dx}12 of ∆η6is 1.5·10−5, while the corresponding theoretical
375
L2 bound from section SM8is 4·10−4. Optimization of the number of terms, in the sense of
376
stopping the summation at smallest one, reduces these errors only slightly. More examples
377
are given in section SM2.5.
378
As shown in figure 5 the trailing systems of the exact solutions, with α given by (4.7),
379
become increasingly strong with increasing n (number of terms in exact solution). Also the
380
heights of the leading pulses are strongly modified.
381
The expansion provides accurate, or sometimes even exact, mathematical solutions also
382
for cases with α close to and above 2. However, the results look strange and have strayed
383
far from the principal wave shape. When an α near 2 is employed in the geometry with an
384
apex (section5) none of the strange features appear due to a modification of the wave shape
385
(sectionSM3.2).
386
5. Waves conveyed to a non-uniform channel. The simple power function solutions
387
described in section 4.1 display some puzzling properties. First, no reflected wave is clearly
388
discernible in the solutions. Then, as shown in figure 5, the substantial changes of the wave
389
shapes due to the higher order terms in the expansions become visible. However, the wave
390
forms, such as (4.14), are merely mathematical solutions on a special class of slopes, as the
391
waves are presented without an origin or a generation mechanism. To put the expansions
392
into perspective propagation of waves from uniform channel section into a non-uniform one,
393
through a apex, is described. The expansion (3.3) is employed in the variable part, while
394
the solution in the uniform region is described by waves of permanent form. Then the local
395
solutions are patched at the apex by requirements of continuity. Relevant examples from the
396
literature, where such patching is used, are [30,12,11,23]. Two problems related to that of
397
transmission at the apex are outlined in section SM3.
398
5.1. Transmission at an apex. An example sketch of a channel transect is shown in figure
399
6. Forx > x0 bothB and ¯h are constant, while the channel is non-uniform forx < x0. Atx0
400
nu. 6 0
x
η C0
-0.29 0.35
0.99
0 1 2 6
x
∆η C0
Figure 4.Normalized surfaces on an inclined plane forF0=Y0,κ= 4,xa= 0.05,xb= 1.6andt0=−1.57.
Left panel: numerical surface elevation and selected ηn at times as given above the crests. Right panel: ∆ηn, marked byn, as explained in the text, for the last time displayed in the left panel. The thin solid line marks zero.
i,0,1.33 i,1,1.60 i,5,1.85 ii,1,4.00 ii,3,2.40 ii,5,2.22
x
η C0
i,0,1.33 i,1,1.60 i,4,1.82 i,7,1.88 i,10,1.91
x
η C0
Figure 5. Surfaces from selected exact solutions forF0=Y0,κ= 3,β= 0 andt0= 0. Curves are marked by family, value ofnand value ofαaccording to (4.7).
we have continuous B and ¯h, but there may be jumps in their derivatives. Still,x0 is also the
401
location where Θ = 0 att= 0, as appearing in (3.2). Moreover, it is also crucial to choosexr 402
in (3.6) equal tox0. Scales are chosen as to giveB = ¯h= 1 for x > x0.
403
At x0 the fluxes of mass, momentum and energy must be continuous, which gives con-
404
tinuous η and u (see, for instance, sectionSM7). The surface elevation of the incident wave,
405
I, is specified while the reflected wave, represented by the elevation R, and the fields of the
406
transmitted wave, denoted by η(t) and u(t), must be found. For the transmitted wave the
407
shape (F0) and the amplitude (A0(x0)) of the principal wave must be determined. The rest
408
x=x0
←η
←I R→
x z
Figure 6. Sketch of geometry inspired by wave tanks.
of (3.3) and (3.7) then follows from (3.4), and the recursion formulas (3.5) and (3.6).
409
Forx > x0 it is convenient to write the solution according to
410
(5.1) η=I+R, I =I(Θi), R=R(Θr), Θi =κ(t+ (x−x0)), Θr =κ(t−(x−x0)).
411
At the apex all the phases then become equal Θi = Θr = Θ = κt. Continuity of η and u,
412
respectively, yield the following relation atx=x0
413
(5.2)
η(t)(x0, t) ∼
∞
X
i=0
κ−iAiFi(Θ) ∼ I(Θ) +R(Θ), u(t)(x0, t) ∼ −
∞
X
i=0
κ−iAiFi(Θ) +
∞
X
i=1
κ−iAi−1,xFi(Θ)
!
∼ −I(Θ) +R(Θ),
414
Combining the two we find
415
(5.3) I(Θ)∼A0F0(Θ) +
∞
X
i=1
κ−i
Ai+1 2Ai−1,x
Fi(Θ).
416
Since (3.6), with xr=x0, yieldsAj+1(x0) =−12Aj,x(x0) the sum vanishes and we obtain
417
(5.4) A0(x0)F0=I, R(Θ)∼ −1 2
∞
X
i=1
κ−iAi−1,x(x0)Fi(Θ) =O(κ−1).
418
When the integral of the incident wave, V = R∞
−∞I dΘ = A0(x0)F1(∞), is non-zero the
419
leading order reflection will be a surface elevation of heightV Hr, whereHr =κ−1(18ddx¯h+14dBdx).
420
This is a long wave that persists after the passing of the incident wave. The incident shape
421
I =aY−1 yieldsV = 0 and the leading surface reflection isaHrY0. If the transition at x0 is
422
smooth, in the sense that the first derivatives ofB and ¯h are zero, the leading order reflected
423
wave becomes−A0(x0)κ−2(161 ddx2¯h2 +18ddx2B2)F2.
424
ProvidedV >0 and that the wave is amplifying (A0,x<0) forx < x0 the equation (3.6)
425
yieldsA1(x0)>0. This implies that the wave height is increased in the transmission, which
426
num.
6 0
x
η a
I=Y−1, κ= 6 I=Y0, κ= 6 I=Y0, κ= 4
x
η a
Figure 7. Reflection at an apex. Left panel: Normalized surfaces for F0 = I = aY0, α = 1, β = 0 and κ= 6. Right panel: The reflected waves for a few selected cases. The bold dashed curves in the legend correspond toη6, numerical solutions are depicted with thin solid lines andη1,η2 with thin dashes. The latter two may be distinguished byη1 approaching a constant asx→1+
may seem counter intuitive as the energy of the incident wave must support both transmission
427
and reflection. However, this is resolved through the modified energy densities on the slope
428
as elaborated in SM7.4.
429
5.2. Examples of transmitted waves.
430
5.2.1. Transmission to an inclined plane. For constant channel width and ¯h=xwe have
431
x0 = 1 and the amplitude recursion yields
432
(5.5) η∼x−14F0+ 1 16κ
x−34 +x−14
F1+ 1 512κ2
9x−54 + 2x−34 + 5x−14
F2+O(κ−3),
433
whereF0 =I and the phase is Θ =κ(2√
x−2 +t). Comparing (5.5) with (4.14) we observe
434
an extra term inA1 with the samex-dependence asA0 (const.×x−14). This is a modification
435
of the wave shape which, in turn, leads to the mid term in A2. The third term in A2 is a
436
shape modification again, and in this manner the series expansion continues. The reflected
437
wave becomes
438
(5.6) R= α
8κF1(Θr) + α2
32κ2F2(Θr) +O(κ−3).
439
Relating the coefficient of (5.6) and (5.5) in view of the continuity forη, we observe that half
440
of the leading order reflection (theF1 term) is linked to the shape change in the transmission,
441
while the other half is linked to the higher order terms in the expansion on the slope.
442
Now the shape of the incident wave is chosen as eitherI =aY0, orI =aY−1 whereais an
443
amplitude factor andYj is as defined in (4.8). Figure7, right panel, shows that the asymptotic
444
expansion reproduces the qualitative features of the transmission with a few terms retained.
445
On the other hand, for the seven term approximation (η6) the error is of order 10−5 in finite
446
depth (x >0.05, say) for the time shown. In the reflection from the N-wave (I =aY−1) there
447
is a long O(κ−2) tail in addition to the single (κ−1) pulse, with shape Y0.
448
