- Email: [email protected]
Modelling currents in highly sheared surface and bed boundary layers
ContinentalShelfResearch,Vol. 12, No. 1, pp. 189-211,1992.
0278-4343/92 $5.00+ 0.00 ~) I991PergamonPressplc
Printedin Great Britain.
Modelling currents in highly sheared surface and bed boundary layers ALAN M. DAVIES* (Received 31 June 1990; in revised form 20 February 1991; accepted 13 March 1991) Abstract--A method is presented for accurately determining currents in the high shear surface and bed boundary layers of a homogeneous sea region under conditions of meteorological and tidal forcing. The Galerkin technique is used in the vertical, with a mixed basis set composed of eigenfunctions (modes) of the eddy viscosity profile and two additional functions to account for shear in the two boundary layers. Calculations of wind induced currents in a closed rectangular basin, and oscillatory pressure driven flow at a point clearly demonstrate that this new approach is significantly more accurate than the "classical" eigenfunction method. For wind induced flow, in which the surface stress is specified, there is no restriction on the form of the additional function used at the surface. However in the case of the additional function at the seabed, analysis and calculations show that there is a stability condition which must be satisfied.
1. I N T R O D U C T I O N
WITH the development of new instruments to measure wind driven currents in the near surface layer and their associated shear (GRWFITHS, 1989), there is increasing demand to model these near surface currents. Similarly with improved methods of measuring current in the bottom boundary layer (LOHRMANet al., 1990) and the associated bed stress, there is a need to develop models which can accurately reproduce these bed stresses, and associated near bed currents. The conventional method of modelling current profiles is to use a finite difference grid (e.g. JAMES, 1990) or finite elements (WERNER, 1987) in the vertical. Generally this grid is uniform, although it has been shown (DAVIESand STEPHENS,1983), that by using finer grids in the near surface and near bed layers, accuracy can be improved in these regions. However in general the order of over 30 grid points in the vertical is required to accurately reproduce currents in these boundary layers, particularly when eddy viscosity is low. An alternative approach is to use a functional expansion in the vertical (CRAI6 1989; FURNES, 1983; FURNES and MORK, 1987; GORDON and SPAULDING, 1987). DAVIES a n d OWEN (1979) showed that highly accurate solutions could be obtained using as few as four Chebyshev polynomials, due to the highly sheared form of these functions in the surface and bed boundary layers. However, the major difficulty with these functions is that they have a high computational overhead and their rapid rate of convergence is offset by their
*Proudman Oceanographic Laboratory, Bidston Observatory, Birkenhead, Merseyside L43 7RA, U.K. 189
190
A.M. DAVIES
computational cost (DAVIES and STEPHENS, 1983). An alternative approach is to use an expansion of eigenfunctions of the eddy viscosity profile (a modal method). The modal technique is computationally very economic, also it is ideal for the new generation of multi-processing computers, in that each mode can be integrated in parallel on a separate processor (DAVIESet al., 1992). The modal method has been shown to be ideal for tidal calculations in continental shelf seas (DAVIES and FURNES, 1980; DAVIES 1986). However for problems involving wind driven flows, where there is a high shear surface layer, the modal expansion converges very slowly (DAVIES, 1983). This problem has recently been examined (DAVIES, 1991b; LARDNER,1990; ZITMAN, 1992), using a mixed basis set, in a similar, though less general approach to that developed here. At the sea bed, the bottom shear boundary layer can be accurately reproduced using a modal method with a no-slip bottom boundary condition (DAVIES, 1991a). Similarly the shear layer associated with a linear slip condition can also be reproduced by ensuring that each mode satisfies this condition exactly. However this method is only appropriate for linear bottom friction with constant viscosity. The application of a quadratic friction law cannot be included by this means, and this can lead to a slower rate of convergence in the near bed region (DAVIES, 1983, 1987). In this paper we present a method for improving the rate of convergence and hence the accuracy of the modal approach in the near surface and near bed highly sheared boundary layers. The technique is developed within the framework of a modal basis set which can be used with a quadratic friction law. The mathematical formulation is presented in Section 2, with comparisons of surface and bed currents computed with the new approach in subsequent sections. Calculations show a significant improvement in the accuracy of boundary layer currents computed with the new approach. However in the case of the bottom boundary layer the method can be restricted by a stability condition which is discussed in the Appendix. 2. MATHEMATICAL DEVELOPMENT (i) Hydrodynamic equations The linear hydrodynamic equations using sigma coordinates, cr = z/h, are given by 0 h
Ou Ot
yv = - g
+
~xx
h
+
vdcr = 0
~ ~oo(/~ ~ )
1 o{ O--v-rot+ y u = - g ~y -v -~ -~ (,u -~ I
(1)
(2) (3)
where t denotes time, x,y,z Cartesian coordinates with u,v, the x,y components of velocity. Acceleration due to gravity g, and geostrophic coefficient 9' are taken as constant, with p vertical eddy viscosity, and h water depth. The linear equations are used here to illustrate the major steps in the mathematical analysis using the Galerkin approach. The Galerkin method can also be readily applied to the non-linear terms (DAVIES, 1980, 1986).
Currents in highlysheared surface and bed boundarylayers
191
Boundary conditions at sea surface and sea bed are given by
(e °ut
- p \h Oa]° =F,,
-p~-~-OO)o = Gs
(4)
--P
(5)
and -Pfh~a),
FB'
1 = GB.
With Fs, Gs the x,y components of surface wind stress and FB, GB components of bed stress. (ii) Expansion method Here we consider a mixed basis set function approach, by expanding the u and v components of current as U = ~ , + ~uB + Ue,
V = ~t'~ + ~
+ lie
(6)
with Ue and Ve given by expressions of time and spatial dependent coefficients Ar(x,y,t), and Br(X,y,t) and functions fr(a) of the form v__m
Ue = / , mrfr(O),
Ve =
r=l
2 Brfr(O). r=l
(7)
The functions ~ s , ~Sbeing specifically chosen to improve convergence in the near surface layer, with Wg, ~Bv for the bed layer. In theory the basis functions fr in equation (7) are arbitrary, but here we will choose eigenfunctions of the eddy viscosity profile. Thus, expressing the eddy viscosity as a time dependent part a(x,y,t) and a function ~(a), gives
lz = a(x,y,t,)cb(o).
(8)
The associated eigenvalue problem is given by d (~ ~--~fo)= - s f
(9)
with e eigenvalues and f eigenfunctions of equation (9), solved subject to boundary conditions all = 0
~-~fal=0.
(10)
The boundary conditions (10) leading in a linear model to an uncoupled set of equations (DAWES, 1983, 1987). Also the first mode is constant in the vertical and is the only mode which contributes to the mean flow. By this means a highly computationally efficient time splitting method can be derived (DAVIES, 1987). Calculation of the eigenfunctions of equation (9) is accomplished using the Galerkin method. Thus (9) is multipled by fk, and the viscosity term is integrated by parts, yielding
192
A . M . DAVIES
rb ~afk l -dP ~ fklo -Ilo dPdfr dfk da da = -er )'i
(11)
Equation (11) can be readily solved for an arbitrary profile of eddy viscosity. Methods are given in DAVIES(1983), DAWESand FURNES(1986) and will not be repeated here. Solving (9) subject to boundary conditions (10), then equation (11) reduces to
o da do (iii)
(12)
ofrfkda"
Galerkin form of the hydrodynamic equations
Substituting (6) into equation (1), and using equation (7), and expressing
Ar = Ur~r, with
B r Vr~)r
(13)
=
fl
~r= 1/ fZda
(14)
o
gives OA = -- L 0t
lh
0X (
2 Urerar __ 0 h r=l
1
Vrdprar + Ce r=l
(15)
where
Ce = _ OOx[[h)o(1.Sdaj _ ~y [h Ii *Sda] - ~xO[h ~o*~da] - O-~y(h ~0*vSda] with ar =
;1
(16)
f~da.
0
Applying the Galerkin method to the solution of equation (2), multiplying by fr, integrating from sea surface to seabed, and integrating the term involving viscosity by parts, followed by substituting surface and bed boundary conditions into (2), gives
fo ~fr 0u da =Y ~oVfrda-g~x ~ofrda-~hhfr(1) + ~fr(O)-~ ~o#c)udfr ~aadaa do.
(17)
[Only the briefest details of deriving (17) have been presented here, and the interested reader is referred to DAvIEs, 1983, 1987, for full details of the steps involved.] Substituting (6) into (17) and rearranging, gives
IloOUe~da --~-'r _- q}"lo Vefrda -g~xO~~ofrda- -~-~ f~(1) + ~h fr(O) - -~1 ~olu OUedfr da da + (18)
Currents in highly sheared surface and bed boundary layers
193
where
De = - o
frda + 7 o attSvfrda- -~ J0/t ~
dfrdOdo
(1 0xi),B [1 1 ~ ° O~2,Bdfr - Jo z-~-frdo + Y o~Bvfrda---~ tz--~-a -~ad°"
(19)
Substituting expansion (7) into (18), using expressions (8), (12), (13) and (14), we obtain from (18), the Galerkin form of the U equation of motion, namely
OUr -~ -
O~ ar -4- ~h fr(0) -7V~ - g ~x
~hhf~(1)
Urer + D e. -- a --h-7--
(20)
In a similar manner from equation (3), we obtain
OVr_ Ot --
-~/
Ur - - g ' ~ O~ y ar +-~fr(O)--~fr(1 ) - - ~
VrEr
--'~- -4- E e
(21)
where
(I oxltSv do
ee=-Jo--g-f -
([ *Sfrda
1(1 o*S dfr
-y__
IlOTJrOxxFBfd(loXI/BfrdO gr 1(1 -
d°
O*~Jdf~do.
(22)
Comparison of equations (15), (20) and (21) with the Galerkin form of these equations derived in earlier papers (DAvIES, 1983, 1987) shows them to be similar with the addition of terms Ce, De and Ee arising from the additional specified functions in expansion (6). In the next section we will show the importance of these specified functions in the accurate determination of boundary layer currents.
(iv) Form of function Previous calculations (DAWES, 1983), have shown that eigenfunctions computed by solving (9), subject to zero derivative boundary conditions (10), converge slowly in high shear regions. Also wind driven currents computed with a small number of modes (of the order of less than 15), exhibit spurious oscillations in the vertical. However when the external shearing force is zero, e.g. the wind stress at the sea surface is zero, then the modal expansion method is very accurate in the surface region. Consequently it is essential to make ~s, ~s proportional to the externally applied wind stress, thus
~ s _ hFs fl~s(o )
(23)
hGs fl~S(o )
(24)
-
p
(o)
194
A.M. DAVIES
with fl a free parameter and ~ ( a ) determining the profile of ~s, ~s. Similarly at the seabed, we can determine a similar type of function, thus • ~-
hFB /~vB(a)
(25)
PkL(1) x[rB v -- hGB /~xItB(a).
(26)
p/~(1) In order for the two functions ~ff,v and ~S,v to be independent of each other, it is necessary to chose functions for U s and ~B that are only non zero over the near surface and near bed regions with ~Fs and its derivatives vanishing at the bed, and xFB and its derivatives vanishing at the surface. Here we choose ~s(o), to be a piecewise function given by: (a) Function (1)
¢)
.s(o ) = 1(o 2_ qa+ o< q\2 -2q (27) =0.0
~->q.
For/3 = 1.0, the function ~s(o) satisfies the surface boundary condition. However for cr - q, ~s(a) and its vertical derivative are zero. Similarly at the seabed 't'B(o) = 0.0
~ -< q
(28)
=(q 1
-qo+
~>q.
Again when/~ = 1.0, equation (28) satisfies the bottom boundary condition. In essence by using piecewise functions such as (27) and (28) in expansion (6), to represent the sheared surface and bed layers, a much smoother function remains that can be readily approximated with the modal expansion (7). Also by an appropriate choice of q and ~ the extent of these piecewise functions can be limited to the near surface and near bed shear regions where modal convergence is poor without influencing the interior regions where the modal expansion is accurate (DAVIES, 1983). In the case of oscillatory flows (tidal flows) in deep water, with low eddy viscosity, the sheared region is confined to a very thin bottom boundary layer, with an inviscid constant current above this (see later). For this particular problem a very effective function is given by a trigonometric expansion of the form: (b) Function (2) • B(ff) = ~-- [(3ff2 -- 6 f f . 2)~2 -- Z I COSrTrff] r=l
(29)
r2 j
where~=l-o. This function is ideal in oscillatory flow problems, as its vertical derivative at the surface is zero, also the trigonometric expansion is very effective in reducing spurious "ripples" in
Currents in highly sheared surface and bed boundary layers
195
the region above the bottom sheared boundary layer (see later). This function was first suggested by HEAPS (1972) as a means of completing the expansion to infinity in the case of a surface wind stress. 3. IDEALIZED CALCULATIONS (a) North Sea rectangle problem To gain some insight into the improvement in representing shear in the sea surface and seabed boundary layers and the associated changes in surface and bed currents, it is instructive to consider some idealized problems. To compare results with previous published calculations (HEAPS, 1972; DAVIES and OWEN, 1979; DAVIES and STEPHENS, 1983; DAVIES, 1983) we consider the problem of wind induced flow in a closed rectangular basin, ("the Heaps rectangle") used by HEAPS (1972) to develop three-dimensional numerical models. The rectangular basin has a water depth h = 65 m, and dimensions of 400 km in the x direction, 800 km in the y direction, with grid spacing Ax = 400/9 Am, Ay -- 800/17 km. Details of the grid and the basin dimensions can be found in HEAPS (1972), DAVIES (1983). Other parameters used in the calculation were g = 9.81 m 2 s -1, p = 1025 kg m -3 and y -1.2 x 10 -4 S-1. Motion was started from rest by a uniform northerly wind stress with values Fs = 0.0, Gs = - 1 . 5 N m -2. To be consistent with previous calculations a linear slip condition was applied at the bed with k = 0.002 m s -1, and eddy viscosity constant in the vertical. A series of calculations were performed with eddy viscosity values of 0.0650 m 2 s -1 and 0.0130 m 2 s -1. These values were used previously by HEAPS (1972). (i) Wind induced profiles,/~ = 0.0650 m 2 s -1. In an initial series of calculations eddy viscosity was fixed at 0.0650 m 2 s - i. Values of surface and bottom current in the centre of the rectangle, 30 h after the onset of the wind field are given in Table 1, together with values computed using the original eigenfunction approach of HEAPS (1972). In the method of HEAPS (1972) the bottom boundary condition was satisfied by each eigenfunction. Consequently each eigenfunction was sheared in the near bed region (see Fig. 1, profile b), although not at the sea surface. The fact that these eigenfunctions satisified exactly the bottom boundary condition means that they converged rapidly in the near bed region (see Table 1). However this method was restricted to a linear slip condition (HEAPS, 1972). Also at the sea surface these modes were shear free and hence converged slowly for the V component of current (the component in the wind direction) (see Table 1) and exhibited some spurious ripples in the computed profile (see Fig. 2a), although for the u component of current (no external wind stress) convergence was rapid (see Table 1). The modes computed using zero derivative surface and bed boundary conditions [equation (10)] can however be used with a quadratic bottom friction law and also can be integrated with a time splitting method with significant computational saving (DAVIES, 1983, 1987). The principal difficulty with these modes is that because they have zero derivative bottom boundary conditions they converge slowly in the near bed region (see Table 1, calculation 2). A consequence of this slow convergence, is that bottom shear is underestimated, together with its associated retarding force on the bottom current, which is larger than that
196
A.M.
DAVIES
Table 1. U and V components of surface and bed current (cm s - l ) , 30 h after the onset of the wind field at the centre of the basin. A value o f p = 0 . 0 6 5 0 m 2 s - 1 o r 0 . 0 1 3 0 m Es - 1 was used in the calculation, with k = 0 . 0 0 2 m s - 1 Calculation 1. O r i g i n a l eigenfunction approach of HEAPS (1972) Number of eigenfunctions m 10
20 /~ = 0 . 0 6 5 0 m z s -1
30
10
20 /~ = 0 . 0 1 3 0 m 2 s -1
30
-15.06 -32.14 7.12 11.09
-15.07 -33.73 7.12 10.97
-15.07 -34.25 7.12 10.95
-66.71 -65.44 9.06 4.95
-66.90 -73.26 9.04 4.46
-66.92 -75.82 9.04 4.36
Current Surface Bed
U V U V
Calculation 2. Additional function ~ s with q = 0 . 2 5 , fl = 1.0
Surface Bed
U V U V
6
10 /~ = 0 . 0 6 5 0 m 2 s -1
20
6
10 /~ = 0 . 0 1 3 0 m 2 s -1
20
-14.71 -35.56 7.52 11.35
-14.85 -35.30 7.35 11.24
-14.93 -35.31 7.26 11.08
-65.21 -83.23 10.86 5.28
-66.41 -81.40 9.89 5.00
-66.75 -81.24 9.44 4.59
Calculation 3. Additional functions ~ s (q = 0 . 2 5 , fl = 1.0) a n d ~ a ( q = 0 . 9 5 , / J = 1.0)
Surface Bed
U V U V
6
10 kt = 0 . 0 6 5 0 m 2 s -1
20
6
10 # = 0 . 0 1 3 0 m 2 s -1
20
-14.81 -35.44 7.34 11.09
-14.93 -35.19 7.21 10.99
-15.00 -35.22 7.14 10.95
-65.24 -82.20 9.95 4.55
-66.38 -80.61 9.20 4.46
-66.68 -80.69 9.03 4.29
computed with the method of HEAPS (1972) (see Table 1). Although as the number of terms m is increased the bottom current converges towards that obtained with the HEAPS (1972) method. Considering surface currents, by adding in a function Ws with q = 0.25, fl = 1.0 [equation (27)] the V component of current converges much more rapidly and there are no significant oscillations (Fig. 2b). An improved rate of convergence in the bottom currents are obtained by including the additional function ~ a , see calculation 3 in Table 1. It is apparent from Table 1, that adding xIta reduces both the u and v components of current with the values computed with m = 20 being in excellent agreement with the method of HEAPS (1972) using m = 30. (ii) Wind induced profiles,/~ = 0.0130 m 2 s -1. As a further test of the accuracy of the method the calculations were repeated with/.t = 0.0130 m 2 s-1 (Table 1). With this lower value of eddy viscosity, shear in the surface layer is increased significantly. Comparing surface currents computed with the addition of ~ s with the original method of Heaps, it is clear that without ~ s the series converges slowly (Table 1, calculation 1), and even with m = 30 the series has not converged and the surface current is the order of 0.06 m s -1 (6 cm s -1) below its converged values (compare calculations I and 2). The addition of ~ s significantly improves the rate of convergence of the surface current and the current no longer exhibits spurious ripples (see Fig. 2c,d). However, it is apparent from Table 1, that
Currents in highly sheared surface and bed boundary layers 0
2
012
012
012
197 012
BED r=l
r=2
r=5
r=4
r=5
Fig. 1. Vertical profiles of the first five modes computed (a) satisfyinga zero derivative boundary condition, (b) a linear slip bottom boundary condition with/~ = 0.0130 m2s - 1and k = 0.002 m s - 1.
the u and v bed currents in calculation 2 converge m o r e slowly, and exceed the values given in calculation 1. This rate of convergence is improved by including the function ~B (compare calculations 1, 2 and 3), again with m = 20 giving comparable values to those obtained using H e a p s ' method with m = 30. (b)
Oscillatory motion
To examine in m o r e detail the influence of ~]3 upon the computed current profile and bed stress, it is instructive to consider flow induced by an oscillatory pressure gradient. In this case current structure is generated by the retarding forces at the sea bed. To examine this p r o b l e m in detail, it is instructive to consider the simple model
Ouot- hxw cos cot +
10u{ 0~)
~-~~ ~/~
(30)
with hx external forcing, chosen to given an inviscid (free stream) current of 1.0 m s -1 and period co, taken as 30 ° h -1 (a semi-diurnal period). In the calculations considered later h = 100 m, with initially/z = 0.0130 m 2 s -1 and subsequently 0.0010 m 2 s -1. B o t t o m friction coefficient was fixed at 0.002 m s -1. Motion was started from rest, and once a periodic state was reached, a harmonic analysis was p e r f o r m e d to yield the amplitude and phase of the current through the water column. (i) Oscillatory profiles, kt = 0.0130 m 2 s -1. Amplitude and phase of the current, in the near bed region, where current shear is a m a x i m u m (see Fig. 3a) computed with an increasing n u m b e r of eigenfunctions are given in Table 2, calculation 1. These eigenfunctions were c o m p u t e d subject to a zero derivative at the bed [equation (10)] and consequently are not sheared in the near bed region (compare profiles a and b in Fig. 1). The fact that they do not satisfy a linear slip condition means that they can be readily used in a m o r e computationally efficient time split approach (DaviEs et al., 1992) with the m o r e physically appropriate quadratic form of b o t t o m friction (DAviES, 1986). H o w e v e r it is clear from Fig. 3a, that very close to the bed, current shear is reduced, with the current
198
A.M.
(a) -0.5
-0.4
o. . . . . . . . . .
I
V (ms-') -0.2
-0.3
.......
~,[4.~
DAVIES
.......
-0.1
I .........
o. I .........
I .........
0.1 I .....
0.2--
0.4-
O0.6--
0.8--
r
.
- -
(b)
V(ms-') -0.4
-0.5
-0.3
-0.2
-0-
1
O.
O. I
I
O-
,
,
,
,
,
,
,
,
,
I
.
.
.
.
.
.
.
.
.
,
.
.
.
.
.
.
.
.
.
,
.
.
.
.
.
.
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
0.2--
0;4
-
O0.6-
0.8-
P
.
- -
Fig. 2. (a) V c o m p o n e n t of velocity, at the centre of the basin, 30 h after the onset of the w i n d field computed using an expansion of eigenfunctions with ~ = 0.0650 rn2 s- 1 w i t h m = 6 ( ) and m = 30 ( - - - ) . (b) A s (a) but with a surface function U s for m = 4 ( ) and m = 10 ( - . - ) .
showing little vertical change, as a consequence of the bottom boundary condition [equation (10)]. Also above the bottom boundary, the profile contains physically unrealisitic "ripples" (Fig. 3a, m = 6, profile A) although these are reduced, as m is increased (Fig. 3a, m = 20, profile B). It is evident from Table 2, calculation 1, that as m increases,
Currents in highly sheared surface and bed boundary layers
199
V(cmls) (C)
-100.
o.
-80.
-60.
I
I
-40.
-20.
[
I
".<__
o.
f
20.
40.
I
[
0.2--
0.4--
(7" 0.6--
0.8--
] , --
Vlcmls) (d)
-100.
-80.
1
o.
-so.
[
-40.
]
-20.
I
o.
20.
I
I
0.2--
0.4--
o" 0.6--
0.8--
[.
--
Fig. 2. (c) As (a) using# = 0.0130 m 2 s- 1 with m = 6 eigenfunctions. (d) As (a) using# = 0.0130 m2 s-1 with m = 6 eigenfunctions and a surface function 't~s.
b o t t o m s h e a r i n c r e a s e s , a n d b o t t o m c u r r e n t is r e d u c e d , a l t h o u g h t h e r a t e o f c o n v e r g e n c e is slow. T h e a d d i t i o n o f ~ B ( q = 0.95, fl = 1.00), i n c r e a s e s t h e r a t e o f c o n v e r g e n c e ( T a b l e 2, c a l c u l a t i o n 2), with t h e c u r r e n t profile c o m p u t e d with m = 6, Fig. 3b, (profile A ) b e i n g significantly s h e a r e d in t h e b e d r e g i o n . H o w e v e r a b o v e this r e g i o n , s t r o n g r i p p l e s a r e
200
A . M . DAVIES
hu (cm s-J)
(a ) 60.
80.
1
o.
100.
,
,
120-
140.
I
I
0,2--
0.4--
(A)
o-
0.6-0.8-I.--
S
m=6
hu(cms -I) 20 o.
40. f
60. r
80. I
~oo.
T2o. I
~ 40. I
0,2--
0,4
- -
o-
(B)
0.6--
~
0.8--
~
m=20 Fig. 3a. Profile of current amplitude induced by an oscillatory flow of 30 ° per hour, with /x = 0.0130 m 2 s - 1computed with m = 6 (Profile A) and m = 20 (Profile 13) eigenfunctions satisfying a zero derivative bottom boundary condition.
201
Currents in highly sheared surface and bed boundary layers
hu(cms-') (b)
20
40.
o.
I
60.
.
80.
.
.
100.
.
120.
140.
I
I
0.2
0.4-
o-
(A)
0.6--
0.8--
m--6
hu(cms -I) 20.
40. I
o.
60. I
80. I
~ 00.
~ 20. I
,40. I
0.2---
0.4
--
(B)
0.6--
1
0.8--
~
1o - -
m=20 Fig. 3b.
As Fig. 3a, but with the addition of bed function ~B (piecewise).
202
A.M. DAVIES
hu(cms-')
(c)
20. o.
0.2
40. I
60. I
so. 1
~oo.
~2o. I
140. I
-
0.4
O0.6-
/
0.8J
m=6 Fig. 3c. As Fig. 3a, but only rn = 6 with the addition of bed function WB(trigonometric).
evident, which are only removed with m the order of 20 (Fig. 3b, m = 20, profile B). It is interesting to note that the profile computed with m = 20 with ~B, is significantly smoother (Fig. 3b) in the upper part of the water column than that computed with the modes alone (Fig. 3a). In theory the choice of q and fl should be arbitrary. However unlike the surface boundary condition where the stress that multiplies XI~B is specified externally, with the bottom boundary condition, this stress is computed from the bed current at the earlier time step. This introduces a stability condition (see Appendix and Table A1) which restricts the choice of q, fl depending upon the coefficient of bottom friction. Using the trigonometric function [equation (29)], even scaled down by setting fl = 0.5, (Table 2, calculation 3) significantly improves the rate of convergence of the series. When fl is increased to fl = 1.0, the series converges very rapidly, with profiles computed with (m = 6) indistinguishable from those computed with (m = 30) (Fig. 3c and Table 2, calculation 4). It is evident from Fig. 3c, that even with m -- 6, a smooth flow free from spurious "ripples" is obtained above the bottom shear layer. (ii) Oscillatory profiles, /~ = 0.0010 m 2 s -1. In a final series of calculations, eddy viscosity was reduced to 0.0010 m e s -1, giving a thin, highly sheared bottom boundary layer (see Fig. 4a). In this case, the modal expansion converged very slowly (Table 3, calculation l), with m -- 6 giving very strong ripples in the vertical, which are still present with m = 30 (Fig. 4a, profiles A,B). Including the additional function ~B [equation (28)] with a range of q and fl values, improves the rate of convergence of the series, with values ofq = 0.99 andfl = 0.75 giving a
203
Currents in highly sheared surface and bed boundary layers
Table 2. Amplitude and phase of the U component of near bed current, induced by a sinusoidalforcing of period 30° per hour computed with/~ = 0.0130 m 2 s- l, k = 0.002 m s- l
Calculation 1. Eigenfunctions only 10 20
6
30
o
h(cm s l)
gO
h(cm s -1)
gO
h(cm s -1)
gO
h(cm s -1)
gO
0.728 0.92 1.00
112.7 67.7 57.0
264 244 236
104.7 74.4 50.6
266 247 238
106.7 83.1 47.2
265 247 239
106.3 81.8 46.0
266 248 241
6 0.728 0.92 1.00
110.9 76.4 49.3
Calculation 2. Additional function ~B (q = 0.95, fi = 1.00) 10 20 265 247 240
105.1 80.0 45.2
266 249 241
106.1 83.4 43.9
266 250 242
30 106.1 83.6 44.0
Calculation 3. Additional function ~B (trigonometric) (fl = 0.5) M 6 10 20
266 249 242
30
ty
h(cm s 1)
gO
h(cm s -1)
gO
h(cm S-1)
gO
h(cm s -l)
gO
0.728 0.92 1.00
109.3 77.2 50.0
265 246 239
105.6 79.4 47.3
265 247 239
106.4 83.4 45.6
265 248 240
106.2 82.8 45.1
265 248 240
Calculation 4. Additional function ~B (trigonometric) (/3 = 1.0) M 6 10 20
0.728 0.92 1.00
30
h(cm s -1)
gO
h(cm s -1)
gO
h(cm S-1)
gO
h(cm s -1 )
gO
106.4 84.2 44.2
266 248 242
106.1 83.6 44.1
265 249 241
106.1 83.6 44.1
265 249 241
106.1 83.6 44.1
265 249 241
s o l u t i o n c o n t a i n i n g s o m e r i p p l e s with m -- 6, a n d a fairly s m o o t h profile with m = 20 (Fig. 4b profiles A , B ) . C o m p a r i n g Fig. 4b with Fig. 4a shows t h a t n e a r b e d s h e a r is i n c r e a s e d a n d t h e s p u r i o u s r i p p l e s a r e r e d u c e d b y i n c l u d i n g ~ B in t h e e x p a n s i o n . T h e r e s t r i c t i o n o f t h e s t a b i l i t y c o n d i t i o n (see A p p e n d i x a n d T a b l e A 1 ) m e a n t h a t a v a l u e o f q = 0.85, c o u l d n o t b e u s e d with a v a l u e o f f l a b o v e a b o u t 0.07, w i t h o u t p r o d u c i n g instabilities (see T a b l e A 1 ) . I d e a l l y a fl v a l u e o f 1.0 w o u l d h a v e b e e n o p t i m a l . U s i n g t h e t r i g o n o m e t r i c f u n c t i o n [ e q u a t i o n (29)] with fl = 0.25, gave t h e m o s t significant i m p r o v e m e n t in t h e r a t e of c o n v e r g e n c e o f t h e series ( T a b l e 3, c a l c u l a t i o n 5). C u r r e n t profiles (Fig. 4c) c o m p u t e d with m = 6, c l e a r l y s h o w a thin highly s h e a r e d b o t t o m b o u n d a r y l a y e r , a l t h o u g h s p u r i o u s " r i p p l e s " a r e a g a i n e v i d e n t a b o v e this layer. I n c r e a s i n g m to rn = 20 r e d u c e s t h e s e ripples.
204
A . M . DAvIEs
(a)
hu(cms-i) 60. f
o.
80. t
loo. I
12o.
~4o. I
0.2--
0.4--
o
(A)
0.6--
0.8--
1 . --
S " ' "
m=6
hu(cms -I) 20.
o.
40.
I
60.
80.
F
J
TO0.
~/
120.
J
140.
I
) 0.2--
)
) 0.4--
)
o-
I
0.6--
(B)
0.8-
r. --
m=20 Fig, 4a. Profile of current amplitude induced by an oscillatory flow of 30 ° per hour, with /~ = 0.0010 m 2 s -1 computed with m = 6 (Profile A) and m = 20 (Profile B) eigenfunctions satisfying a zero derivative bottom boundary condition.
205
Currents in highly sheared surface and bed boundary layers
hu(cms-')
(b) 20. o.
,o. 1
60. I
80. I
,00.
,20. I
,40. L
0.2
0.4
o-
(A)
0.6-
0.8-
m--6
hu(cms-') 2o. o.
40. I
60. I
80. I
0.2--
100. I
<
12o. I
~4 0 . I
<
(
0.4--
(
(
(B)
0.6 -
0.8 --~
1.
m--20 Fig. 41). As Fig. 4a, but with the addition of bed function ~ n (piecewise).
206
A . M . DAVIES
(C)
hu(cms-') 20. o.
40. 1
80. f
80. I
,oo.
12o. I
14o. I
0.2--
0.4--
(A)
o-
0.6--
0.8--
m=6
hu(cms-j) 20. o.
40. I
60. [
80. 1
~00. I/
~20. I
,40. I
0.2--
0.4
--
(B)
00.6--
0.8--
1.
- -
m=20 Fig. 4c.
As Fig. 4a, but with the addition of bed function ~B (trigonometric) with ,3 = 0.25.
Currents in highly sheared surface and bed boundary layers
207
hu(cms-') (d )
o.
~o.
4o.
6o.
8o.
I
I
I
I
1 oo.
120.
I
140.
1
0.2--
0.4
--
O0.6--
0.8--
I.
m=30 Fig. 4d. As Fig. 4a, but only for rn = 30 ( ) with the addition of bed function xltB (trigonometric) with ~ = 1.0. The nature of the trigonometric function is such that its value ~ B at the seabed is larger for m = 6 than m = 30. A consequence of this is that when fl was increased to 0.5, and finally to fl = 1.0, the solutions c o m p u t e d with m less than 30 were unstable. Profiles c o m p u t e d with fl = 0.5, m = 30 exhibit some oscillations, although when fl is increased to 1.0 these oscillations are r e m o v e d (Fig. 4d). These calculations clearly show that for oscillatory flow, when eddy viscosity is reduced, the b o t t o m boundary layer becomes m o r e highly sheared, and occupies a small fraction of the water column. To maintain accuracy under these conditions it is essential to use an enhanced basis set with the trigonometric function appearing preferable to the other functions. 4. CONCLUDING REMARKS This p a p e r has dealt with the mathematical development and numerical solution of the linear hydrodynamic equations using the Galerkin method in the vertical. The solution has been in terms of an expansion of modes, with the addition of extra functions to take account of highly sheared surface and b o t t o m boundary layers. By this means the computational efficiency of the modes in terms of a time split solution algorithm which can be efficiently integrated on multi-processor vector computers (DAVIES et al., 1992) is retained. Also the m e t h o d can readily incorporate quadratic b o t t o m friction. Calculations of wind induced flow in a closed rectangular basin clearly show that the rate of convergence of the modes can be significantly increased, particularly in the near surface region, by including functions to take account of near surface and near bed shear regions.
208
A . M . DAVIES
Tab& 3. Amp~tude and phase of the U component of near bed current, induced by a sinusoidalforcing of period 30°per hour compu~d with ~ = 0.0010 m 2 s -1, k = 0.002 m s -1 Calculation 1. Eigenfunctions only 10 20
6 a 0.884 0.92 0.96 1.00
h(cm s -1) 73.2 61.1 55.3
54.5
6
30
g°
h(cm s - l )
gO
h(cm s -1)
gO
h(cm s -1)
gO
251 236 220 214
109.0 78.8 44.7 35.3
265 256 232 208
102.6 114.4 70.4 22.9
269 260 247 215
106.5 104.6 88.5 20.3
268 263 247 220
Calculation 2. Additionalfunction~B (~ = 0.85, fl = 0.05) 10 20
30
a
h(cm s-1)
gO
h(cm s -1)
gO
h(cm s -1)
gO
h(cm s -1)
gO
0.884 0.92 0.96 1.00
86.1 76.0 63.8 47.7
250 238 228 222
109.3 82.7 51.1 33.8
265 255 232 222
102.6 114.1 72.0 22.6
269 260 247 216
106.3 104.7 88.8 20.1
268 263 247 220
6 a 0.884 0.92 0.96 1.00
h(cm s - j ) 85.5 79.7 76.0 44.0
6
Calculation 3. Additionalfunction~B (~ = 0.95, fl = 0.15) 10 20
30
gO
h(cm s -1)
gO
h(cm s -~)
gO
h(cm s -1)
gO
252 241 231 226
108.5 85.4 60.8 31.8
266 255 233 216
102.8 113.3 75.6 21.8
269 261 248 218
106.0 105.0 89.6 19.7
268 264 248 221
Calculation 4. Additionalfunction~B(~ = 0.99, fl = 0.75) 10 20
30
o
h(cm s -1)
gO
h(cm s -x)
gO
h(cm s -a)
gO
h(cm s -1)
gO
0.884 0.92 0.96 1.00
84.8 78.1 74.9 42.4
253 242 233 228
108.6 85.4 60.6 30.6
266 256 235 219
102.4 114.0 76.6 20.8
269 261 248 221
105.9 105.0 91.3 18.7
268 264 249 224
Calculation 5. Additional function WB (trigonometric) (/5 = 0.25) 6 10 20
30
~r
h(cm s -1 )
gO
h(cm s -1 )
g°
h(cm s -1 )
gO
h(cm s -1)
gO
0.92 0.96 1.00
102.6 77.6 30.0
249 244 242
95.4 70.6 28.2
254 236 224
112.8 77.9 21.2
262 246 221
105.2 90.1 19.3
263 248 223
Currents in highly sheared surface and bed boundary layers
209
I n the case of c u r r e n t s d r i v e n b y a n oscillatory p r e s s u r e g r a d i e n t , w h e r e c u r r e n t s t r u c t u r e is g e n e r a t e d b y b o t t o m drag, the i n c l u s i o n of a n a d d i t i o n a l f u n c t i o n is particularly v a l u a b l e . T h e a p p l i c a t i o n of a t r i g o n o m e t r i c f u n c t i o n in the vertical a p p e a r s to be o p t i m a l , giving a n i n c r e a s i n g rate of c o n v e r g e n c e in the n e a r b e d region. This f u n c t i o n also m i n i m i z e s s p u r i o u s oscillations in the u p p e r part of the w a t e r c o l u m n . U n f o r t u n a t e l y the stability c o n d i t i o n c o n s i d e r e d in the A p p e n d i x m e a n s that the e x t e r n a l f u n c t i o n s c a n n o t always b e i n c l u d e d with fl = 1.0. A t p r e s e n t t h e r e does n o t a p p e a r to be a n y m e a n s of o v e r c o m i n g this. T h e a p p l i c a t i o n of the m e t h o d d e v e l o p e d in this p a p e r to m o r e physically realistic r e g i o n s is p r e s e n t l y in progress a n d results will b e r e p o r t e d s u b s e q u e n t l y . Acknowledgements--This paper is dedicated to Bruno M. Jamart who died on 17th November 1990. Bruno was an excellent scientist, who made a major contribution to numerical modelling using both functional and finite difference methods. He was also a good friend, who will be sorely missed. The care taken by Mr R. A. Smith in preparing the diagrams and Mrs L. Parry in typing the paper is very much appreciated. REFERENCES CRAIG P. D. (1989) Constant eddy viscosity models of vertical current structure forced by periodic winds. Continental Shelf Research, 9, 343-358. DAVIESA. M. (1980) Application of the Galerkin method to the formulation of a three-dimensional non-linear hydrodynamic numerical sea model. Applied Mathematical Modelling, 4,245-256. DAVIESA. M. (1983) Formulation of a linear three-dimensional hydrodynamic sea model using a GalerkinEigenfunction method. International Journal for Numerical Methods in Fluids, 3, 33-60. DAvmsA. M. (1986) A three-dimensional model of the northwest European ContinentalShelf, with application to the M4 tide. Journal of Physical Oceanography, 16,797-813. DAVIESA. M. (1987) Spectral models in continental shelf sea oceanography. In: Three-dimensionalcoastal ocean models, N. S. HEAPS,editor, A.G.U., Washington, pp. 71-106. DAVIESA. M. (1991a) On using turbulence energy models to develop spectral viscositymodels. ContinentalShelf Research, 11, 1313-1353. DAXaESA. M. (1991b) Solution of the three dimensional linear hydrodynamic equations using an enhanced eigenfunction approach. International Journal for Numerical Methods in Fluids, 13, 235-250. DAVIESA. M. and A. OWEN(1979) Three-dimensional numerical sea model using the Galerkin method with a polynomial basis set. Applied Mathematical Modelling, 3,421-428. DAVIESA. M. and G. K. FURNES(1980) Observed and computed M2 tidal currents in the North Sea. Journal of Physical Oceanography, 10, 237-257. DAVIESA. M. and C. V. STEPHENS(1983) Comparison of the finite difference and Galerkin methods as applied to the solution of the hydrodynamic equations. Applied Mathematical Modelling, 7,226-240. DAWESA. M. and G. K. FURNES(1986) On the determination of vertical structure functions for time-dependent flow problems. Tellus, 38A, 462-477. DAVIESA. M., R. B. GRZONKAand C. V. STEPnENS(1992) Implementation of a three-dimensional hydrodynamic numerical model using parallel processing on a CRAY X-MP series computer. Advances in Parallel Computing, to appear. FURNESG. K. (1983) A three-dimensional numerical sea model with eddy viscosityvarying piecewise linearly in the vertical. Continental Shelf Research, 2,231-242. FURNESG. K. and M. MORK(1987) Formulation of a continuously stratified sea model with three-dimensional representation of the upper layer. Coastal Engineering, 11,415-445. GORDONR. B. and M. L. SPAULDING(1987) Numerical simulations of the tidal- and wind-driven circulation in Narragansett Bay. Estuarine, Coastaland Shelf Science, 24,611-636. GRIFFITHSG. (1989) Estimates of eddy diffusion in the upper ocean based on the backscatter strength measurements of an acoustic Doppler current profiler. In: Advances in water modelling and measurement, M. H. PALMER,editor, BHRA, Cranfield, pp. 259-270. HEAPSN. S. (1972) On the numerical solution of the three-dimensional hydrodynamic equations for tides and storm surges. M~moires de la Soci~tid des Sciences de Lidge, Series, 6, 2, 143-180.
210
A . M . DAVIES
JAMES I. D. (1990) Numerical modelling of density-driven circulation in shelf seas. In: Modeling marine systems, Vol. 2, A. M. DAVIES, editor, CRC Press, Florida, U.S.A., pp. 345-372. LARDNER R. W. (1990) Numerical solution of tile linearized three-dimensional tidal equations using eddy viscosity eigenfunctions. Journal of Geophysical Research (Oceans), 95, 22,269-22,274. LOHRMANN, A., B. HACKE~ and L. P. ROED (1990) High resolution measurements of turbulence, velocity and stress using a pulse-to-pulse coherent sonar. Journal of Atmospheric and Oceanic Technology, 7, 19-37. WERNER F. E. (1987) A numerical study of secondary flows over continental shelf edges. Continental Shelf Research, 7, 379--409. Z1TMANT. J. (1992) Quasi three-dimensional current modelling based on a modified version of Davies' shape function approach. Continental Shelf Research, 12, 143-158.
APPENDIX In order to understand why the choice offl and ~ affects the stability of the solution it is necessary to consider a simpler model than that given by equations (1)-(3) and examine the role of bottom friction in the solution. By this means a simple stability condition can be readily derived without resorting to complex mathematical analysis. An analysis of bed currents from the tidally forced model prior to the model becoming unstable showed time step oscillation of increasing amplitude before the model became unstable. To understand this and why enhancing the bottom stress produces the problem, it is useful to consider the simple model
dUb - G - FB dt ph
(A1)
Fn = kpU h
(A2)
where
with G the external specified forcing term (e.g. wind or tidal forcing which need not be constant). Physically the forcing term G causes the flow to accelerate. The current Uh increases with time and would continue to increase, except that the term F n which is proportional to the flow, also increases with Uh, and retards the flow, with a steady state being reached when FB/ph = G. Expressing Uh as
Uh = Ue(1 ) + W~(1),
(A3)
*~(1) = 2*n(1)Fn
(A4)
then from equation (25) we have
with 2-
hfl q/~(1)'
(A5)
Substituting (A4) into (A3) and using (A2), gives for (A1)
dUhdt - G - ~h (kpUe(1) + kp2~B(1)FB)
(A6)
dUhdt = G - p-~(FB - kP2~B(1)FB)
(A7)
which can be expressed as
where FB = k0U~(1). In order to satisfy the bottom boundary condition and to make physical sense ~ (A1) must be negative and must be scaled with FB (Fig. 5). It is evident on physical grounds from the simple flow picture shown in Fig. A1 that when Uh is positive, FB is also positive, consequently the addition of FB ~ n only retards the flow in the near bed region if qtB is negative. By scaling ~B by F B the additonal function ~ n also enhances bed shear when flow is in the opposite direction.
Currents in highly sheared surface and bed boundary layers
Us
Us
Uh
0
uh Fig. A1.
211
+
Uh
Ks ql s
=
u
Schematic figure showing the influence of adding a piecewise bottom layer function, scaled by F B, on to the current profile.
Table A1.
Values of K = [kp~h[ and influence upon stability for a range of k, q and fl values
k
kp
~
fl
~h
K
0.002 0.002
2.05 2.05
0.99 0.99
0.95 1.00
-0.464 -0.488
0.9504 1.000
Stable Weak instability
0.002 0.002
2.05 2.05
0.85 0.85
0.05 0.10
-0.366 -0.732
0.7503 1.5106
Stable Unstable
0.02
20.5
0.99
0.5
-0.049
1.005
Weak instability
For physically realistic k, p, 2. which are all positive, a consequence of ~B(1) being negative is that ~ and FB appear with opposite signs. In order to proceed to a simple stability condition it is necessary to assume that FB and F a are not substantially different, which is borne out by the fact that F a is formed by correcting FB, consequently (AT) can be written as d U h _ G - FB + K FB
dt
ph
(g8)
ph
where
K = [kp~hl
(A9)
and ~ h = 2 US(l) is the value of the additional function at the bed. The physical reason for the numerical instability is clearly evident from (A8). If K is unity then there is no damping in the model and the flow continues to accelerate. For K greater than one, the bed stress rather than acting as a damping term continues to force the solution which rapidly grows. Consequently in order to obtain a stable solution
K = Ikp~hl < 1.
(A10)
Namely the value of the additional function ~B at the bed is constrained by (A10) and ~ and fl must be adjusted to ensure that K < 1. Although we have only considered a very simple model here, calculations using two values of bottom friction coefficient with a range ~ and fl values appear to confirm that K < 1 (see Table A1) is required to maintain stability.
0278-4343/92 $5.00+ 0.00 ~) I991PergamonPressplc
Printedin Great Britain.
Modelling currents in highly sheared surface and bed boundary layers ALAN M. DAVIES* (Received 31 June 1990; in revised form 20 February 1991; accepted 13 March 1991) Abstract--A method is presented for accurately determining currents in the high shear surface and bed boundary layers of a homogeneous sea region under conditions of meteorological and tidal forcing. The Galerkin technique is used in the vertical, with a mixed basis set composed of eigenfunctions (modes) of the eddy viscosity profile and two additional functions to account for shear in the two boundary layers. Calculations of wind induced currents in a closed rectangular basin, and oscillatory pressure driven flow at a point clearly demonstrate that this new approach is significantly more accurate than the "classical" eigenfunction method. For wind induced flow, in which the surface stress is specified, there is no restriction on the form of the additional function used at the surface. However in the case of the additional function at the seabed, analysis and calculations show that there is a stability condition which must be satisfied.
1. I N T R O D U C T I O N
WITH the development of new instruments to measure wind driven currents in the near surface layer and their associated shear (GRWFITHS, 1989), there is increasing demand to model these near surface currents. Similarly with improved methods of measuring current in the bottom boundary layer (LOHRMANet al., 1990) and the associated bed stress, there is a need to develop models which can accurately reproduce these bed stresses, and associated near bed currents. The conventional method of modelling current profiles is to use a finite difference grid (e.g. JAMES, 1990) or finite elements (WERNER, 1987) in the vertical. Generally this grid is uniform, although it has been shown (DAVIESand STEPHENS,1983), that by using finer grids in the near surface and near bed layers, accuracy can be improved in these regions. However in general the order of over 30 grid points in the vertical is required to accurately reproduce currents in these boundary layers, particularly when eddy viscosity is low. An alternative approach is to use a functional expansion in the vertical (CRAI6 1989; FURNES, 1983; FURNES and MORK, 1987; GORDON and SPAULDING, 1987). DAVIES a n d OWEN (1979) showed that highly accurate solutions could be obtained using as few as four Chebyshev polynomials, due to the highly sheared form of these functions in the surface and bed boundary layers. However, the major difficulty with these functions is that they have a high computational overhead and their rapid rate of convergence is offset by their
*Proudman Oceanographic Laboratory, Bidston Observatory, Birkenhead, Merseyside L43 7RA, U.K. 189
190
A.M. DAVIES
computational cost (DAVIES and STEPHENS, 1983). An alternative approach is to use an expansion of eigenfunctions of the eddy viscosity profile (a modal method). The modal technique is computationally very economic, also it is ideal for the new generation of multi-processing computers, in that each mode can be integrated in parallel on a separate processor (DAVIESet al., 1992). The modal method has been shown to be ideal for tidal calculations in continental shelf seas (DAVIES and FURNES, 1980; DAVIES 1986). However for problems involving wind driven flows, where there is a high shear surface layer, the modal expansion converges very slowly (DAVIES, 1983). This problem has recently been examined (DAVIES, 1991b; LARDNER,1990; ZITMAN, 1992), using a mixed basis set, in a similar, though less general approach to that developed here. At the sea bed, the bottom shear boundary layer can be accurately reproduced using a modal method with a no-slip bottom boundary condition (DAVIES, 1991a). Similarly the shear layer associated with a linear slip condition can also be reproduced by ensuring that each mode satisfies this condition exactly. However this method is only appropriate for linear bottom friction with constant viscosity. The application of a quadratic friction law cannot be included by this means, and this can lead to a slower rate of convergence in the near bed region (DAVIES, 1983, 1987). In this paper we present a method for improving the rate of convergence and hence the accuracy of the modal approach in the near surface and near bed highly sheared boundary layers. The technique is developed within the framework of a modal basis set which can be used with a quadratic friction law. The mathematical formulation is presented in Section 2, with comparisons of surface and bed currents computed with the new approach in subsequent sections. Calculations show a significant improvement in the accuracy of boundary layer currents computed with the new approach. However in the case of the bottom boundary layer the method can be restricted by a stability condition which is discussed in the Appendix. 2. MATHEMATICAL DEVELOPMENT (i) Hydrodynamic equations The linear hydrodynamic equations using sigma coordinates, cr = z/h, are given by 0 h
Ou Ot
yv = - g
+
~xx
h
+
vdcr = 0
~ ~oo(/~ ~ )
1 o{ O--v-rot+ y u = - g ~y -v -~ -~ (,u -~ I
(1)
(2) (3)
where t denotes time, x,y,z Cartesian coordinates with u,v, the x,y components of velocity. Acceleration due to gravity g, and geostrophic coefficient 9' are taken as constant, with p vertical eddy viscosity, and h water depth. The linear equations are used here to illustrate the major steps in the mathematical analysis using the Galerkin approach. The Galerkin method can also be readily applied to the non-linear terms (DAVIES, 1980, 1986).
Currents in highlysheared surface and bed boundarylayers
191
Boundary conditions at sea surface and sea bed are given by
(e °ut
- p \h Oa]° =F,,
-p~-~-OO)o = Gs
(4)
--P
(5)
and -Pfh~a),
FB'
1 = GB.
With Fs, Gs the x,y components of surface wind stress and FB, GB components of bed stress. (ii) Expansion method Here we consider a mixed basis set function approach, by expanding the u and v components of current as U = ~ , + ~uB + Ue,
V = ~t'~ + ~
+ lie
(6)
with Ue and Ve given by expressions of time and spatial dependent coefficients Ar(x,y,t), and Br(X,y,t) and functions fr(a) of the form v__m
Ue = / , mrfr(O),
Ve =
r=l
2 Brfr(O). r=l
(7)
The functions ~ s , ~Sbeing specifically chosen to improve convergence in the near surface layer, with Wg, ~Bv for the bed layer. In theory the basis functions fr in equation (7) are arbitrary, but here we will choose eigenfunctions of the eddy viscosity profile. Thus, expressing the eddy viscosity as a time dependent part a(x,y,t) and a function ~(a), gives
lz = a(x,y,t,)cb(o).
(8)
The associated eigenvalue problem is given by d (~ ~--~fo)= - s f
(9)
with e eigenvalues and f eigenfunctions of equation (9), solved subject to boundary conditions all = 0
~-~fal=0.
(10)
The boundary conditions (10) leading in a linear model to an uncoupled set of equations (DAWES, 1983, 1987). Also the first mode is constant in the vertical and is the only mode which contributes to the mean flow. By this means a highly computationally efficient time splitting method can be derived (DAVIES, 1987). Calculation of the eigenfunctions of equation (9) is accomplished using the Galerkin method. Thus (9) is multipled by fk, and the viscosity term is integrated by parts, yielding
192
A . M . DAVIES
rb ~afk l -dP ~ fklo -Ilo dPdfr dfk da da = -er )'i
(11)
Equation (11) can be readily solved for an arbitrary profile of eddy viscosity. Methods are given in DAVIES(1983), DAWESand FURNES(1986) and will not be repeated here. Solving (9) subject to boundary conditions (10), then equation (11) reduces to
o da do (iii)
(12)
ofrfkda"
Galerkin form of the hydrodynamic equations
Substituting (6) into equation (1), and using equation (7), and expressing
Ar = Ur~r, with
B r Vr~)r
(13)
=
fl
~r= 1/ fZda
(14)
o
gives OA = -- L 0t
lh
0X (
2 Urerar __ 0 h r=l
1
Vrdprar + Ce r=l
(15)
where
Ce = _ OOx[[h)o(1.Sdaj _ ~y [h Ii *Sda] - ~xO[h ~o*~da] - O-~y(h ~0*vSda] with ar =
;1
(16)
f~da.
0
Applying the Galerkin method to the solution of equation (2), multiplying by fr, integrating from sea surface to seabed, and integrating the term involving viscosity by parts, followed by substituting surface and bed boundary conditions into (2), gives
fo ~fr 0u da =Y ~oVfrda-g~x ~ofrda-~hhfr(1) + ~fr(O)-~ ~o#c)udfr ~aadaa do.
(17)
[Only the briefest details of deriving (17) have been presented here, and the interested reader is referred to DAvIEs, 1983, 1987, for full details of the steps involved.] Substituting (6) into (17) and rearranging, gives
IloOUe~da --~-'r _- q}"lo Vefrda -g~xO~~ofrda- -~-~ f~(1) + ~h fr(O) - -~1 ~olu OUedfr da da + (18)
Currents in highly sheared surface and bed boundary layers
193
where
De = - o
frda + 7 o attSvfrda- -~ J0/t ~
dfrdOdo
(1 0xi),B [1 1 ~ ° O~2,Bdfr - Jo z-~-frdo + Y o~Bvfrda---~ tz--~-a -~ad°"
(19)
Substituting expansion (7) into (18), using expressions (8), (12), (13) and (14), we obtain from (18), the Galerkin form of the U equation of motion, namely
OUr -~ -
O~ ar -4- ~h fr(0) -7V~ - g ~x
~hhf~(1)
Urer + D e. -- a --h-7--
(20)
In a similar manner from equation (3), we obtain
OVr_ Ot --
-~/
Ur - - g ' ~ O~ y ar +-~fr(O)--~fr(1 ) - - ~
VrEr
--'~- -4- E e
(21)
where
(I oxltSv do
ee=-Jo--g-f -
([ *Sfrda
1(1 o*S dfr
-y__
IlOTJrOxxFBfd(loXI/BfrdO gr 1(1 -
d°
O*~Jdf~do.
(22)
Comparison of equations (15), (20) and (21) with the Galerkin form of these equations derived in earlier papers (DAvIES, 1983, 1987) shows them to be similar with the addition of terms Ce, De and Ee arising from the additional specified functions in expansion (6). In the next section we will show the importance of these specified functions in the accurate determination of boundary layer currents.
(iv) Form of function Previous calculations (DAWES, 1983), have shown that eigenfunctions computed by solving (9), subject to zero derivative boundary conditions (10), converge slowly in high shear regions. Also wind driven currents computed with a small number of modes (of the order of less than 15), exhibit spurious oscillations in the vertical. However when the external shearing force is zero, e.g. the wind stress at the sea surface is zero, then the modal expansion method is very accurate in the surface region. Consequently it is essential to make ~s, ~s proportional to the externally applied wind stress, thus
~ s _ hFs fl~s(o )
(23)
hGs fl~S(o )
(24)
-
p
(o)
194
A.M. DAVIES
with fl a free parameter and ~ ( a ) determining the profile of ~s, ~s. Similarly at the seabed, we can determine a similar type of function, thus • ~-
hFB /~vB(a)
(25)
PkL(1) x[rB v -- hGB /~xItB(a).
(26)
p/~(1) In order for the two functions ~ff,v and ~S,v to be independent of each other, it is necessary to chose functions for U s and ~B that are only non zero over the near surface and near bed regions with ~Fs and its derivatives vanishing at the bed, and xFB and its derivatives vanishing at the surface. Here we choose ~s(o), to be a piecewise function given by: (a) Function (1)
¢)
.s(o ) = 1(o 2_ qa+ o< q\2 -2q (27) =0.0
~->q.
For/3 = 1.0, the function ~s(o) satisfies the surface boundary condition. However for cr - q, ~s(a) and its vertical derivative are zero. Similarly at the seabed 't'B(o) = 0.0
~ -< q
(28)
=(q 1
-qo+
~>q.
Again when/~ = 1.0, equation (28) satisfies the bottom boundary condition. In essence by using piecewise functions such as (27) and (28) in expansion (6), to represent the sheared surface and bed layers, a much smoother function remains that can be readily approximated with the modal expansion (7). Also by an appropriate choice of q and ~ the extent of these piecewise functions can be limited to the near surface and near bed shear regions where modal convergence is poor without influencing the interior regions where the modal expansion is accurate (DAVIES, 1983). In the case of oscillatory flows (tidal flows) in deep water, with low eddy viscosity, the sheared region is confined to a very thin bottom boundary layer, with an inviscid constant current above this (see later). For this particular problem a very effective function is given by a trigonometric expansion of the form: (b) Function (2) • B(ff) = ~-- [(3ff2 -- 6 f f . 2)~2 -- Z I COSrTrff] r=l
(29)
r2 j
where~=l-o. This function is ideal in oscillatory flow problems, as its vertical derivative at the surface is zero, also the trigonometric expansion is very effective in reducing spurious "ripples" in
Currents in highly sheared surface and bed boundary layers
195
the region above the bottom sheared boundary layer (see later). This function was first suggested by HEAPS (1972) as a means of completing the expansion to infinity in the case of a surface wind stress. 3. IDEALIZED CALCULATIONS (a) North Sea rectangle problem To gain some insight into the improvement in representing shear in the sea surface and seabed boundary layers and the associated changes in surface and bed currents, it is instructive to consider some idealized problems. To compare results with previous published calculations (HEAPS, 1972; DAVIES and OWEN, 1979; DAVIES and STEPHENS, 1983; DAVIES, 1983) we consider the problem of wind induced flow in a closed rectangular basin, ("the Heaps rectangle") used by HEAPS (1972) to develop three-dimensional numerical models. The rectangular basin has a water depth h = 65 m, and dimensions of 400 km in the x direction, 800 km in the y direction, with grid spacing Ax = 400/9 Am, Ay -- 800/17 km. Details of the grid and the basin dimensions can be found in HEAPS (1972), DAVIES (1983). Other parameters used in the calculation were g = 9.81 m 2 s -1, p = 1025 kg m -3 and y -1.2 x 10 -4 S-1. Motion was started from rest by a uniform northerly wind stress with values Fs = 0.0, Gs = - 1 . 5 N m -2. To be consistent with previous calculations a linear slip condition was applied at the bed with k = 0.002 m s -1, and eddy viscosity constant in the vertical. A series of calculations were performed with eddy viscosity values of 0.0650 m 2 s -1 and 0.0130 m 2 s -1. These values were used previously by HEAPS (1972). (i) Wind induced profiles,/~ = 0.0650 m 2 s -1. In an initial series of calculations eddy viscosity was fixed at 0.0650 m 2 s - i. Values of surface and bottom current in the centre of the rectangle, 30 h after the onset of the wind field are given in Table 1, together with values computed using the original eigenfunction approach of HEAPS (1972). In the method of HEAPS (1972) the bottom boundary condition was satisfied by each eigenfunction. Consequently each eigenfunction was sheared in the near bed region (see Fig. 1, profile b), although not at the sea surface. The fact that these eigenfunctions satisified exactly the bottom boundary condition means that they converged rapidly in the near bed region (see Table 1). However this method was restricted to a linear slip condition (HEAPS, 1972). Also at the sea surface these modes were shear free and hence converged slowly for the V component of current (the component in the wind direction) (see Table 1) and exhibited some spurious ripples in the computed profile (see Fig. 2a), although for the u component of current (no external wind stress) convergence was rapid (see Table 1). The modes computed using zero derivative surface and bed boundary conditions [equation (10)] can however be used with a quadratic bottom friction law and also can be integrated with a time splitting method with significant computational saving (DAVIES, 1983, 1987). The principal difficulty with these modes is that because they have zero derivative bottom boundary conditions they converge slowly in the near bed region (see Table 1, calculation 2). A consequence of this slow convergence, is that bottom shear is underestimated, together with its associated retarding force on the bottom current, which is larger than that
196
A.M.
DAVIES
Table 1. U and V components of surface and bed current (cm s - l ) , 30 h after the onset of the wind field at the centre of the basin. A value o f p = 0 . 0 6 5 0 m 2 s - 1 o r 0 . 0 1 3 0 m Es - 1 was used in the calculation, with k = 0 . 0 0 2 m s - 1 Calculation 1. O r i g i n a l eigenfunction approach of HEAPS (1972) Number of eigenfunctions m 10
20 /~ = 0 . 0 6 5 0 m z s -1
30
10
20 /~ = 0 . 0 1 3 0 m 2 s -1
30
-15.06 -32.14 7.12 11.09
-15.07 -33.73 7.12 10.97
-15.07 -34.25 7.12 10.95
-66.71 -65.44 9.06 4.95
-66.90 -73.26 9.04 4.46
-66.92 -75.82 9.04 4.36
Current Surface Bed
U V U V
Calculation 2. Additional function ~ s with q = 0 . 2 5 , fl = 1.0
Surface Bed
U V U V
6
10 /~ = 0 . 0 6 5 0 m 2 s -1
20
6
10 /~ = 0 . 0 1 3 0 m 2 s -1
20
-14.71 -35.56 7.52 11.35
-14.85 -35.30 7.35 11.24
-14.93 -35.31 7.26 11.08
-65.21 -83.23 10.86 5.28
-66.41 -81.40 9.89 5.00
-66.75 -81.24 9.44 4.59
Calculation 3. Additional functions ~ s (q = 0 . 2 5 , fl = 1.0) a n d ~ a ( q = 0 . 9 5 , / J = 1.0)
Surface Bed
U V U V
6
10 kt = 0 . 0 6 5 0 m 2 s -1
20
6
10 # = 0 . 0 1 3 0 m 2 s -1
20
-14.81 -35.44 7.34 11.09
-14.93 -35.19 7.21 10.99
-15.00 -35.22 7.14 10.95
-65.24 -82.20 9.95 4.55
-66.38 -80.61 9.20 4.46
-66.68 -80.69 9.03 4.29
computed with the method of HEAPS (1972) (see Table 1). Although as the number of terms m is increased the bottom current converges towards that obtained with the HEAPS (1972) method. Considering surface currents, by adding in a function Ws with q = 0.25, fl = 1.0 [equation (27)] the V component of current converges much more rapidly and there are no significant oscillations (Fig. 2b). An improved rate of convergence in the bottom currents are obtained by including the additional function ~ a , see calculation 3 in Table 1. It is apparent from Table 1, that adding xIta reduces both the u and v components of current with the values computed with m = 20 being in excellent agreement with the method of HEAPS (1972) using m = 30. (ii) Wind induced profiles,/~ = 0.0130 m 2 s -1. As a further test of the accuracy of the method the calculations were repeated with/.t = 0.0130 m 2 s-1 (Table 1). With this lower value of eddy viscosity, shear in the surface layer is increased significantly. Comparing surface currents computed with the addition of ~ s with the original method of Heaps, it is clear that without ~ s the series converges slowly (Table 1, calculation 1), and even with m = 30 the series has not converged and the surface current is the order of 0.06 m s -1 (6 cm s -1) below its converged values (compare calculations I and 2). The addition of ~ s significantly improves the rate of convergence of the surface current and the current no longer exhibits spurious ripples (see Fig. 2c,d). However, it is apparent from Table 1, that
Currents in highly sheared surface and bed boundary layers 0
2
012
012
012
197 012
BED r=l
r=2
r=5
r=4
r=5
Fig. 1. Vertical profiles of the first five modes computed (a) satisfyinga zero derivative boundary condition, (b) a linear slip bottom boundary condition with/~ = 0.0130 m2s - 1and k = 0.002 m s - 1.
the u and v bed currents in calculation 2 converge m o r e slowly, and exceed the values given in calculation 1. This rate of convergence is improved by including the function ~B (compare calculations 1, 2 and 3), again with m = 20 giving comparable values to those obtained using H e a p s ' method with m = 30. (b)
Oscillatory motion
To examine in m o r e detail the influence of ~]3 upon the computed current profile and bed stress, it is instructive to consider flow induced by an oscillatory pressure gradient. In this case current structure is generated by the retarding forces at the sea bed. To examine this p r o b l e m in detail, it is instructive to consider the simple model
Ouot- hxw cos cot +
10u{ 0~)
~-~~ ~/~
(30)
with hx external forcing, chosen to given an inviscid (free stream) current of 1.0 m s -1 and period co, taken as 30 ° h -1 (a semi-diurnal period). In the calculations considered later h = 100 m, with initially/z = 0.0130 m 2 s -1 and subsequently 0.0010 m 2 s -1. B o t t o m friction coefficient was fixed at 0.002 m s -1. Motion was started from rest, and once a periodic state was reached, a harmonic analysis was p e r f o r m e d to yield the amplitude and phase of the current through the water column. (i) Oscillatory profiles, kt = 0.0130 m 2 s -1. Amplitude and phase of the current, in the near bed region, where current shear is a m a x i m u m (see Fig. 3a) computed with an increasing n u m b e r of eigenfunctions are given in Table 2, calculation 1. These eigenfunctions were c o m p u t e d subject to a zero derivative at the bed [equation (10)] and consequently are not sheared in the near bed region (compare profiles a and b in Fig. 1). The fact that they do not satisfy a linear slip condition means that they can be readily used in a m o r e computationally efficient time split approach (DaviEs et al., 1992) with the m o r e physically appropriate quadratic form of b o t t o m friction (DAviES, 1986). H o w e v e r it is clear from Fig. 3a, that very close to the bed, current shear is reduced, with the current
198
A.M.
(a) -0.5
-0.4
o. . . . . . . . . .
I
V (ms-') -0.2
-0.3
.......
~,[4.~
DAVIES
.......
-0.1
I .........
o. I .........
I .........
0.1 I .....
0.2--
0.4-
O0.6--
0.8--
r
.
- -
(b)
V(ms-') -0.4
-0.5
-0.3
-0.2
-0-
1
O.
O. I
I
O-
,
,
,
,
,
,
,
,
,
I
.
.
.
.
.
.
.
.
.
,
.
.
.
.
.
.
.
.
.
,
.
.
.
.
.
.
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
,
0.2--
0;4
-
O0.6-
0.8-
P
.
- -
Fig. 2. (a) V c o m p o n e n t of velocity, at the centre of the basin, 30 h after the onset of the w i n d field computed using an expansion of eigenfunctions with ~ = 0.0650 rn2 s- 1 w i t h m = 6 ( ) and m = 30 ( - - - ) . (b) A s (a) but with a surface function U s for m = 4 ( ) and m = 10 ( - . - ) .
showing little vertical change, as a consequence of the bottom boundary condition [equation (10)]. Also above the bottom boundary, the profile contains physically unrealisitic "ripples" (Fig. 3a, m = 6, profile A) although these are reduced, as m is increased (Fig. 3a, m = 20, profile B). It is evident from Table 2, calculation 1, that as m increases,
Currents in highly sheared surface and bed boundary layers
199
V(cmls) (C)
-100.
o.
-80.
-60.
I
I
-40.
-20.
[
I
".<__
o.
f
20.
40.
I
[
0.2--
0.4--
(7" 0.6--
0.8--
] , --
Vlcmls) (d)
-100.
-80.
1
o.
-so.
[
-40.
]
-20.
I
o.
20.
I
I
0.2--
0.4--
o" 0.6--
0.8--
[.
--
Fig. 2. (c) As (a) using# = 0.0130 m 2 s- 1 with m = 6 eigenfunctions. (d) As (a) using# = 0.0130 m2 s-1 with m = 6 eigenfunctions and a surface function 't~s.
b o t t o m s h e a r i n c r e a s e s , a n d b o t t o m c u r r e n t is r e d u c e d , a l t h o u g h t h e r a t e o f c o n v e r g e n c e is slow. T h e a d d i t i o n o f ~ B ( q = 0.95, fl = 1.00), i n c r e a s e s t h e r a t e o f c o n v e r g e n c e ( T a b l e 2, c a l c u l a t i o n 2), with t h e c u r r e n t profile c o m p u t e d with m = 6, Fig. 3b, (profile A ) b e i n g significantly s h e a r e d in t h e b e d r e g i o n . H o w e v e r a b o v e this r e g i o n , s t r o n g r i p p l e s a r e
200
A . M . DAVIES
hu (cm s-J)
(a ) 60.
80.
1
o.
100.
,
,
120-
140.
I
I
0,2--
0.4--
(A)
o-
0.6-0.8-I.--
S
m=6
hu(cms -I) 20 o.
40. f
60. r
80. I
~oo.
T2o. I
~ 40. I
0,2--
0,4
- -
o-
(B)
0.6--
~
0.8--
~
m=20 Fig. 3a. Profile of current amplitude induced by an oscillatory flow of 30 ° per hour, with /x = 0.0130 m 2 s - 1computed with m = 6 (Profile A) and m = 20 (Profile 13) eigenfunctions satisfying a zero derivative bottom boundary condition.
201
Currents in highly sheared surface and bed boundary layers
hu(cms-') (b)
20
40.
o.
I
60.
.
80.
.
.
100.
.
120.
140.
I
I
0.2
0.4-
o-
(A)
0.6--
0.8--
m--6
hu(cms -I) 20.
40. I
o.
60. I
80. I
~ 00.
~ 20. I
,40. I
0.2---
0.4
--
(B)
0.6--
1
0.8--
~
1o - -
m=20 Fig. 3b.
As Fig. 3a, but with the addition of bed function ~B (piecewise).
202
A.M. DAVIES
hu(cms-')
(c)
20. o.
0.2
40. I
60. I
so. 1
~oo.
~2o. I
140. I
-
0.4
O0.6-
/
0.8J
m=6 Fig. 3c. As Fig. 3a, but only rn = 6 with the addition of bed function WB(trigonometric).
evident, which are only removed with m the order of 20 (Fig. 3b, m = 20, profile B). It is interesting to note that the profile computed with m = 20 with ~B, is significantly smoother (Fig. 3b) in the upper part of the water column than that computed with the modes alone (Fig. 3a). In theory the choice of q and fl should be arbitrary. However unlike the surface boundary condition where the stress that multiplies XI~B is specified externally, with the bottom boundary condition, this stress is computed from the bed current at the earlier time step. This introduces a stability condition (see Appendix and Table A1) which restricts the choice of q, fl depending upon the coefficient of bottom friction. Using the trigonometric function [equation (29)], even scaled down by setting fl = 0.5, (Table 2, calculation 3) significantly improves the rate of convergence of the series. When fl is increased to fl = 1.0, the series converges very rapidly, with profiles computed with (m = 6) indistinguishable from those computed with (m = 30) (Fig. 3c and Table 2, calculation 4). It is evident from Fig. 3c, that even with m -- 6, a smooth flow free from spurious "ripples" is obtained above the bottom shear layer. (ii) Oscillatory profiles, /~ = 0.0010 m 2 s -1. In a final series of calculations, eddy viscosity was reduced to 0.0010 m e s -1, giving a thin, highly sheared bottom boundary layer (see Fig. 4a). In this case, the modal expansion converged very slowly (Table 3, calculation l), with m -- 6 giving very strong ripples in the vertical, which are still present with m = 30 (Fig. 4a, profiles A,B). Including the additional function ~B [equation (28)] with a range of q and fl values, improves the rate of convergence of the series, with values ofq = 0.99 andfl = 0.75 giving a
203
Currents in highly sheared surface and bed boundary layers
Table 2. Amplitude and phase of the U component of near bed current, induced by a sinusoidalforcing of period 30° per hour computed with/~ = 0.0130 m 2 s- l, k = 0.002 m s- l
Calculation 1. Eigenfunctions only 10 20
6
30
o
h(cm s l)
gO
h(cm s -1)
gO
h(cm s -1)
gO
h(cm s -1)
gO
0.728 0.92 1.00
112.7 67.7 57.0
264 244 236
104.7 74.4 50.6
266 247 238
106.7 83.1 47.2
265 247 239
106.3 81.8 46.0
266 248 241
6 0.728 0.92 1.00
110.9 76.4 49.3
Calculation 2. Additional function ~B (q = 0.95, fi = 1.00) 10 20 265 247 240
105.1 80.0 45.2
266 249 241
106.1 83.4 43.9
266 250 242
30 106.1 83.6 44.0
Calculation 3. Additional function ~B (trigonometric) (fl = 0.5) M 6 10 20
266 249 242
30
ty
h(cm s 1)
gO
h(cm s -1)
gO
h(cm S-1)
gO
h(cm s -l)
gO
0.728 0.92 1.00
109.3 77.2 50.0
265 246 239
105.6 79.4 47.3
265 247 239
106.4 83.4 45.6
265 248 240
106.2 82.8 45.1
265 248 240
Calculation 4. Additional function ~B (trigonometric) (/3 = 1.0) M 6 10 20
0.728 0.92 1.00
30
h(cm s -1)
gO
h(cm s -1)
gO
h(cm S-1)
gO
h(cm s -1 )
gO
106.4 84.2 44.2
266 248 242
106.1 83.6 44.1
265 249 241
106.1 83.6 44.1
265 249 241
106.1 83.6 44.1
265 249 241
s o l u t i o n c o n t a i n i n g s o m e r i p p l e s with m -- 6, a n d a fairly s m o o t h profile with m = 20 (Fig. 4b profiles A , B ) . C o m p a r i n g Fig. 4b with Fig. 4a shows t h a t n e a r b e d s h e a r is i n c r e a s e d a n d t h e s p u r i o u s r i p p l e s a r e r e d u c e d b y i n c l u d i n g ~ B in t h e e x p a n s i o n . T h e r e s t r i c t i o n o f t h e s t a b i l i t y c o n d i t i o n (see A p p e n d i x a n d T a b l e A 1 ) m e a n t h a t a v a l u e o f q = 0.85, c o u l d n o t b e u s e d with a v a l u e o f f l a b o v e a b o u t 0.07, w i t h o u t p r o d u c i n g instabilities (see T a b l e A 1 ) . I d e a l l y a fl v a l u e o f 1.0 w o u l d h a v e b e e n o p t i m a l . U s i n g t h e t r i g o n o m e t r i c f u n c t i o n [ e q u a t i o n (29)] with fl = 0.25, gave t h e m o s t significant i m p r o v e m e n t in t h e r a t e of c o n v e r g e n c e o f t h e series ( T a b l e 3, c a l c u l a t i o n 5). C u r r e n t profiles (Fig. 4c) c o m p u t e d with m = 6, c l e a r l y s h o w a thin highly s h e a r e d b o t t o m b o u n d a r y l a y e r , a l t h o u g h s p u r i o u s " r i p p l e s " a r e a g a i n e v i d e n t a b o v e this layer. I n c r e a s i n g m to rn = 20 r e d u c e s t h e s e ripples.
204
A . M . DAvIEs
(a)
hu(cms-i) 60. f
o.
80. t
loo. I
12o.
~4o. I
0.2--
0.4--
o
(A)
0.6--
0.8--
1 . --
S " ' "
m=6
hu(cms -I) 20.
o.
40.
I
60.
80.
F
J
TO0.
~/
120.
J
140.
I
) 0.2--
)
) 0.4--
)
o-
I
0.6--
(B)
0.8-
r. --
m=20 Fig, 4a. Profile of current amplitude induced by an oscillatory flow of 30 ° per hour, with /~ = 0.0010 m 2 s -1 computed with m = 6 (Profile A) and m = 20 (Profile B) eigenfunctions satisfying a zero derivative bottom boundary condition.
205
Currents in highly sheared surface and bed boundary layers
hu(cms-')
(b) 20. o.
,o. 1
60. I
80. I
,00.
,20. I
,40. L
0.2
0.4
o-
(A)
0.6-
0.8-
m--6
hu(cms-') 2o. o.
40. I
60. I
80. I
0.2--
100. I
<
12o. I
~4 0 . I
<
(
0.4--
(
(
(B)
0.6 -
0.8 --~
1.
m--20 Fig. 41). As Fig. 4a, but with the addition of bed function ~ n (piecewise).
206
A . M . DAVIES
(C)
hu(cms-') 20. o.
40. 1
80. f
80. I
,oo.
12o. I
14o. I
0.2--
0.4--
(A)
o-
0.6--
0.8--
m=6
hu(cms-j) 20. o.
40. I
60. [
80. 1
~00. I/
~20. I
,40. I
0.2--
0.4
--
(B)
00.6--
0.8--
1.
- -
m=20 Fig. 4c.
As Fig. 4a, but with the addition of bed function ~B (trigonometric) with ,3 = 0.25.
Currents in highly sheared surface and bed boundary layers
207
hu(cms-') (d )
o.
~o.
4o.
6o.
8o.
I
I
I
I
1 oo.
120.
I
140.
1
0.2--
0.4
--
O0.6--
0.8--
I.
m=30 Fig. 4d. As Fig. 4a, but only for rn = 30 ( ) with the addition of bed function xltB (trigonometric) with ~ = 1.0. The nature of the trigonometric function is such that its value ~ B at the seabed is larger for m = 6 than m = 30. A consequence of this is that when fl was increased to 0.5, and finally to fl = 1.0, the solutions c o m p u t e d with m less than 30 were unstable. Profiles c o m p u t e d with fl = 0.5, m = 30 exhibit some oscillations, although when fl is increased to 1.0 these oscillations are r e m o v e d (Fig. 4d). These calculations clearly show that for oscillatory flow, when eddy viscosity is reduced, the b o t t o m boundary layer becomes m o r e highly sheared, and occupies a small fraction of the water column. To maintain accuracy under these conditions it is essential to use an enhanced basis set with the trigonometric function appearing preferable to the other functions. 4. CONCLUDING REMARKS This p a p e r has dealt with the mathematical development and numerical solution of the linear hydrodynamic equations using the Galerkin method in the vertical. The solution has been in terms of an expansion of modes, with the addition of extra functions to take account of highly sheared surface and b o t t o m boundary layers. By this means the computational efficiency of the modes in terms of a time split solution algorithm which can be efficiently integrated on multi-processor vector computers (DAVIES et al., 1992) is retained. Also the m e t h o d can readily incorporate quadratic b o t t o m friction. Calculations of wind induced flow in a closed rectangular basin clearly show that the rate of convergence of the modes can be significantly increased, particularly in the near surface region, by including functions to take account of near surface and near bed shear regions.
208
A . M . DAVIES
Tab& 3. Amp~tude and phase of the U component of near bed current, induced by a sinusoidalforcing of period 30°per hour compu~d with ~ = 0.0010 m 2 s -1, k = 0.002 m s -1 Calculation 1. Eigenfunctions only 10 20
6 a 0.884 0.92 0.96 1.00
h(cm s -1) 73.2 61.1 55.3
54.5
6
30
g°
h(cm s - l )
gO
h(cm s -1)
gO
h(cm s -1)
gO
251 236 220 214
109.0 78.8 44.7 35.3
265 256 232 208
102.6 114.4 70.4 22.9
269 260 247 215
106.5 104.6 88.5 20.3
268 263 247 220
Calculation 2. Additionalfunction~B (~ = 0.85, fl = 0.05) 10 20
30
a
h(cm s-1)
gO
h(cm s -1)
gO
h(cm s -1)
gO
h(cm s -1)
gO
0.884 0.92 0.96 1.00
86.1 76.0 63.8 47.7
250 238 228 222
109.3 82.7 51.1 33.8
265 255 232 222
102.6 114.1 72.0 22.6
269 260 247 216
106.3 104.7 88.8 20.1
268 263 247 220
6 a 0.884 0.92 0.96 1.00
h(cm s - j ) 85.5 79.7 76.0 44.0
6
Calculation 3. Additionalfunction~B (~ = 0.95, fl = 0.15) 10 20
30
gO
h(cm s -1)
gO
h(cm s -~)
gO
h(cm s -1)
gO
252 241 231 226
108.5 85.4 60.8 31.8
266 255 233 216
102.8 113.3 75.6 21.8
269 261 248 218
106.0 105.0 89.6 19.7
268 264 248 221
Calculation 4. Additionalfunction~B(~ = 0.99, fl = 0.75) 10 20
30
o
h(cm s -1)
gO
h(cm s -x)
gO
h(cm s -a)
gO
h(cm s -1)
gO
0.884 0.92 0.96 1.00
84.8 78.1 74.9 42.4
253 242 233 228
108.6 85.4 60.6 30.6
266 256 235 219
102.4 114.0 76.6 20.8
269 261 248 221
105.9 105.0 91.3 18.7
268 264 249 224
Calculation 5. Additional function WB (trigonometric) (/5 = 0.25) 6 10 20
30
~r
h(cm s -1 )
gO
h(cm s -1 )
g°
h(cm s -1 )
gO
h(cm s -1)
gO
0.92 0.96 1.00
102.6 77.6 30.0
249 244 242
95.4 70.6 28.2
254 236 224
112.8 77.9 21.2
262 246 221
105.2 90.1 19.3
263 248 223
Currents in highly sheared surface and bed boundary layers
209
I n the case of c u r r e n t s d r i v e n b y a n oscillatory p r e s s u r e g r a d i e n t , w h e r e c u r r e n t s t r u c t u r e is g e n e r a t e d b y b o t t o m drag, the i n c l u s i o n of a n a d d i t i o n a l f u n c t i o n is particularly v a l u a b l e . T h e a p p l i c a t i o n of a t r i g o n o m e t r i c f u n c t i o n in the vertical a p p e a r s to be o p t i m a l , giving a n i n c r e a s i n g rate of c o n v e r g e n c e in the n e a r b e d region. This f u n c t i o n also m i n i m i z e s s p u r i o u s oscillations in the u p p e r part of the w a t e r c o l u m n . U n f o r t u n a t e l y the stability c o n d i t i o n c o n s i d e r e d in the A p p e n d i x m e a n s that the e x t e r n a l f u n c t i o n s c a n n o t always b e i n c l u d e d with fl = 1.0. A t p r e s e n t t h e r e does n o t a p p e a r to be a n y m e a n s of o v e r c o m i n g this. T h e a p p l i c a t i o n of the m e t h o d d e v e l o p e d in this p a p e r to m o r e physically realistic r e g i o n s is p r e s e n t l y in progress a n d results will b e r e p o r t e d s u b s e q u e n t l y . Acknowledgements--This paper is dedicated to Bruno M. Jamart who died on 17th November 1990. Bruno was an excellent scientist, who made a major contribution to numerical modelling using both functional and finite difference methods. He was also a good friend, who will be sorely missed. The care taken by Mr R. A. Smith in preparing the diagrams and Mrs L. Parry in typing the paper is very much appreciated. REFERENCES CRAIG P. D. (1989) Constant eddy viscosity models of vertical current structure forced by periodic winds. Continental Shelf Research, 9, 343-358. DAVIESA. M. (1980) Application of the Galerkin method to the formulation of a three-dimensional non-linear hydrodynamic numerical sea model. Applied Mathematical Modelling, 4,245-256. DAVIESA. M. (1983) Formulation of a linear three-dimensional hydrodynamic sea model using a GalerkinEigenfunction method. International Journal for Numerical Methods in Fluids, 3, 33-60. DAvmsA. M. (1986) A three-dimensional model of the northwest European ContinentalShelf, with application to the M4 tide. Journal of Physical Oceanography, 16,797-813. DAVIESA. M. (1987) Spectral models in continental shelf sea oceanography. In: Three-dimensionalcoastal ocean models, N. S. HEAPS,editor, A.G.U., Washington, pp. 71-106. DAVIESA. M. (1991a) On using turbulence energy models to develop spectral viscositymodels. ContinentalShelf Research, 11, 1313-1353. DAXaESA. M. (1991b) Solution of the three dimensional linear hydrodynamic equations using an enhanced eigenfunction approach. International Journal for Numerical Methods in Fluids, 13, 235-250. DAVIESA. M. and A. OWEN(1979) Three-dimensional numerical sea model using the Galerkin method with a polynomial basis set. Applied Mathematical Modelling, 3,421-428. DAVIESA. M. and G. K. FURNES(1980) Observed and computed M2 tidal currents in the North Sea. Journal of Physical Oceanography, 10, 237-257. DAVIESA. M. and C. V. STEPHENS(1983) Comparison of the finite difference and Galerkin methods as applied to the solution of the hydrodynamic equations. Applied Mathematical Modelling, 7,226-240. DAWESA. M. and G. K. FURNES(1986) On the determination of vertical structure functions for time-dependent flow problems. Tellus, 38A, 462-477. DAVIESA. M., R. B. GRZONKAand C. V. STEPnENS(1992) Implementation of a three-dimensional hydrodynamic numerical model using parallel processing on a CRAY X-MP series computer. Advances in Parallel Computing, to appear. FURNESG. K. (1983) A three-dimensional numerical sea model with eddy viscosityvarying piecewise linearly in the vertical. Continental Shelf Research, 2,231-242. FURNESG. K. and M. MORK(1987) Formulation of a continuously stratified sea model with three-dimensional representation of the upper layer. Coastal Engineering, 11,415-445. GORDONR. B. and M. L. SPAULDING(1987) Numerical simulations of the tidal- and wind-driven circulation in Narragansett Bay. Estuarine, Coastaland Shelf Science, 24,611-636. GRIFFITHSG. (1989) Estimates of eddy diffusion in the upper ocean based on the backscatter strength measurements of an acoustic Doppler current profiler. In: Advances in water modelling and measurement, M. H. PALMER,editor, BHRA, Cranfield, pp. 259-270. HEAPSN. S. (1972) On the numerical solution of the three-dimensional hydrodynamic equations for tides and storm surges. M~moires de la Soci~tid des Sciences de Lidge, Series, 6, 2, 143-180.
210
A . M . DAVIES
JAMES I. D. (1990) Numerical modelling of density-driven circulation in shelf seas. In: Modeling marine systems, Vol. 2, A. M. DAVIES, editor, CRC Press, Florida, U.S.A., pp. 345-372. LARDNER R. W. (1990) Numerical solution of tile linearized three-dimensional tidal equations using eddy viscosity eigenfunctions. Journal of Geophysical Research (Oceans), 95, 22,269-22,274. LOHRMANN, A., B. HACKE~ and L. P. ROED (1990) High resolution measurements of turbulence, velocity and stress using a pulse-to-pulse coherent sonar. Journal of Atmospheric and Oceanic Technology, 7, 19-37. WERNER F. E. (1987) A numerical study of secondary flows over continental shelf edges. Continental Shelf Research, 7, 379--409. Z1TMANT. J. (1992) Quasi three-dimensional current modelling based on a modified version of Davies' shape function approach. Continental Shelf Research, 12, 143-158.
APPENDIX In order to understand why the choice offl and ~ affects the stability of the solution it is necessary to consider a simpler model than that given by equations (1)-(3) and examine the role of bottom friction in the solution. By this means a simple stability condition can be readily derived without resorting to complex mathematical analysis. An analysis of bed currents from the tidally forced model prior to the model becoming unstable showed time step oscillation of increasing amplitude before the model became unstable. To understand this and why enhancing the bottom stress produces the problem, it is useful to consider the simple model
dUb - G - FB dt ph
(A1)
Fn = kpU h
(A2)
where
with G the external specified forcing term (e.g. wind or tidal forcing which need not be constant). Physically the forcing term G causes the flow to accelerate. The current Uh increases with time and would continue to increase, except that the term F n which is proportional to the flow, also increases with Uh, and retards the flow, with a steady state being reached when FB/ph = G. Expressing Uh as
Uh = Ue(1 ) + W~(1),
(A3)
*~(1) = 2*n(1)Fn
(A4)
then from equation (25) we have
with 2-
hfl q/~(1)'
(A5)
Substituting (A4) into (A3) and using (A2), gives for (A1)
dUhdt - G - ~h (kpUe(1) + kp2~B(1)FB)
(A6)
dUhdt = G - p-~(FB - kP2~B(1)FB)
(A7)
which can be expressed as
where FB = k0U~(1). In order to satisfy the bottom boundary condition and to make physical sense ~ (A1) must be negative and must be scaled with FB (Fig. 5). It is evident on physical grounds from the simple flow picture shown in Fig. A1 that when Uh is positive, FB is also positive, consequently the addition of FB ~ n only retards the flow in the near bed region if qtB is negative. By scaling ~B by F B the additonal function ~ n also enhances bed shear when flow is in the opposite direction.
Currents in highly sheared surface and bed boundary layers
Us
Us
Uh
0
uh Fig. A1.
211
+
Uh
Ks ql s
=
u
Schematic figure showing the influence of adding a piecewise bottom layer function, scaled by F B, on to the current profile.
Table A1.
Values of K = [kp~h[ and influence upon stability for a range of k, q and fl values
k
kp
~
fl
~h
K
0.002 0.002
2.05 2.05
0.99 0.99
0.95 1.00
-0.464 -0.488
0.9504 1.000
Stable Weak instability
0.002 0.002
2.05 2.05
0.85 0.85
0.05 0.10
-0.366 -0.732
0.7503 1.5106
Stable Unstable
0.02
20.5
0.99
0.5
-0.049
1.005
Weak instability
For physically realistic k, p, 2. which are all positive, a consequence of ~B(1) being negative is that ~ and FB appear with opposite signs. In order to proceed to a simple stability condition it is necessary to assume that FB and F a are not substantially different, which is borne out by the fact that F a is formed by correcting FB, consequently (AT) can be written as d U h _ G - FB + K FB
dt
ph
(g8)
ph
where
K = [kp~hl
(A9)
and ~ h = 2 US(l) is the value of the additional function at the bed. The physical reason for the numerical instability is clearly evident from (A8). If K is unity then there is no damping in the model and the flow continues to accelerate. For K greater than one, the bed stress rather than acting as a damping term continues to force the solution which rapidly grows. Consequently in order to obtain a stable solution
K = Ikp~hl < 1.
(A10)
Namely the value of the additional function ~B at the bed is constrained by (A10) and ~ and fl must be adjusted to ensure that K < 1. Although we have only considered a very simple model here, calculations using two values of bottom friction coefficient with a range ~ and fl values appear to confirm that K < 1 (see Table A1) is required to maintain stability.