A numerical study of wind-induced motions in shallow coastal seas: model and basic experiments

A numerical study of wind-induced motions in shallow coastal seas: model and basic experiments Milivoj KuzmiC Center for Marine

Research

Zagreb,

Rudjer

BoSkoviC

Institute,

Zagreb,

Yugoslavia

A three-dimensional linear hydrodynamical numerical model of wind-induced motions in a homogeneous, shallow sea has been developed. The model is based upon the spectral, eigenfunction approach originated by N. S. Heaps. It accommodates arbitrary variation of vertical eddy viscosity by phaseamplitude transformation of eigenfunctions and numerical solution of the vertical, eigenvalue problem. The horizontal problem is solved using a forward-time, staggered-space finite difference scheme. Comparisons to analytical solutions show that the applied method yields, after a finite inverse transformation, meaningful results of satisfactory accuracy. After verification, the model was used to simulate several cases of vertically variable eddy viscosity grouped around four scenarios: constant eddy viscosity, eddy viscosity decreasing throughout the water column, eddy viscosity influenced by the surface and bottom boundary layers, and eddy viscosity derived from observational data. The most pronounced dtfferences due to vertically variable eddy viscosity have been observed within the third scenario. Tentative comparison to appropriately processed current data from the Northern Adriatic has shown that vertically variable eddy viscosity can contribute to a more successful match of model predictions and empirical data. Keywords:

numerical model, eigenfunction method, Priifer transformation, shallow seas, Northern

Adriatic

Introduction A number of papers have appeared in the literature reporting numerical modeling studies of wind-induced motions in the sea. The models and approaches differ considerably in both mathematical formulation (finite difference, finite element, semianalytical, etc.) and physical sophistication (simple linear, nonlinear threedimensional, turbulence resolving, etc.). An important point of physical necessity and computational consequence has been the vertical resolution of such models. The problem has been attacked by solving directly the three-dimensional equations of motion,‘v2 by vertically layering the computational domain,3,4 by combining numerical and analytical techniques5*‘j and by using vertical expansion. 7,8 Current profiles have been also extracted from a vertically integrated model by combining the Galerkin and convolution methods.’ A particularly appealing approach was pioneered by HeapslO who identified an eigenvalue problem which formed a basis upon which the solution through the vertical could be expanded. He later elaborated on the method.“~i2 A model of wind-driven motions in the Northern

Address reprint requests to Dr. KuzmiC at the Center for Marine Research Zagreb, Rudjer BoSkoviC Institute, P.O. Box 1016, 41001 Zagreb, Yugoslavia. Received February 1988; revised August 1988.

178

Appl. Math. Modelling,

1989, Vol. 13, March

Adriatic, based on the Heaps approach, has been recently developed by KuzmiC et a1.13 Our modeling philosophy has been based on two major premises: advance through stepwise refinement and verification against observational data. Such a philosophy led to a simple, linear, three-dimensional dynamical model and the study of influence of the bottom topography and magnitude of vertically constant eddy viscosity coefficient.i3 The spectral eigenfunction method seemed suitable for such a less aspiring start, but also very promising for future advancements. We later used the same model to study the influence of heterogeneity in the wind field over the Northern Adriatic.i4,r5 The magnitude of the eddy viscosity coefficient proved to be of considerable importance while the wind curl emerged as the most energetic source of variability in the fields considered (the sea surface and velocity at various depths). Comparison of model predictions to available and appropriately processed field data has shown significant similarity, i3,i5 but also left room for improvements. Observed discrepancies provided further incentive for the present work. In keeping with the declared philosophy, the goal of this study is threefold: (1) to explore the feasibility of a particular numerical treatment of the equations of motion within the framework of the Heaps spectral method; (2) to simulate and comparatively analyze several characteristic depth-dependent eddy viscosity functions; (3) to assess the preliminary importance of the vertical eddy viscosity variability for explaining the

0 1989 Butterworth

Publishers

A numerical

study of wind-induced

observed winter, barotropic dynamics of the Northern Adriatic. To that end the paper is organized as follows. The model-governing equations are described in the second section. Numerical solution of the model is the subject of the third section. Various formulations of the vertical eddy viscosity variability are briefly presented in the fourth section. Characteristic cases are selected for computer experiments. Verification of the model and presentation of computational experiments constitute the fifth and sixth section, respectively. In the sixth section modeling results for the Northern Adriatic are also tentatively examined in the light of observational evidence from the modeled area. Finally, the results are summarized and further research steps suggested in the seventh section. Governing

equations

Aiming to model mathematically the interaction of heterogeneous wind stress at the surface, variable vertical eddy viscosity, and real bottom topography in a homogeneous sea, one can start from linear equations of motion and continuity. Assuming incompressible water, hydrostatic motion, constant Coriolis parameter, and neglecting the nonlinear interactions and lateral shear, we have

(2) (3)

P

= 0

N(u) --=

au _-

N(a) au -h aa

rxb

h&

p’

= -Tats

= 1 (5)

The stresses in equations (4) and (5) are defined as TX,= GPcddz-x,

(6)

rys = cDp,v,~ rxb

=

kpu(l),

qb

=

h(l)

(7)

The symbols u, and v, denote the wind velocity components, pa is the density of air, p is the density of seawater, CD is the nondimensional drag coefficient, and k is the coefficient of bottom friction. The use of linear slip at the bottom, although practiced,13J8 may seem somewhat unjustified. However, results of the correlation and regression analysis of the data for a site in the Northern Adriatici suggest that the linear and quadratic laws of bottom friction could be of comparable validity. Furthermore, by increasing (decreasing) the value of k one can, within the same conceptual framework, reasonably well approximate no-slip (freeslip) conditions at the bottom. Forristall et a1.,18 for example, found that k = 0.005 m-s-’ represented very limited slip, while k = lop5 mess’ produced almost pure slip. The quadratic law is most often used in modeling the wind stress, although other relations have been proposed. Some problems with the wind stress parameterization have been studied by KuzmiC and OrliC. l9 Along the solid boundary no normal flow is allowed: (U,V), = w, = 0

(8)

En+@=0

Equations (l)-(3) are written in the coordinate system (x,y,a) in which the horizontal coordinates are laid at the undisturbed sea surface and the vertical coordinate is a-transformed (o = zlh(x,y)) and directed downward. The transformation is similar to the one originally proposed by PhillipsI and is often used when both the bottom topography and vertical resolution are important (e.g., Paul and Lick”). The initial conditions are satisfied, assuming a state of rest. Conditions prescribed at the surface and bottom boundaries are N(cr) au -- 7.x.T -,o’ h am = _-7YS at fl

and

The problem is closed by postulating a radiation condition of the form

where t = time 5 = elevation of the surface U,V = horizontal components of current h = undisturbed depth of water (T = transformed vertical coordinate iV(cr) = coefficient of vertical eddy viscosity f = Coriolis parameter g = acceleration of gravity

N(c) au --= haa

motions in shallow coastal seas: M. Kuzmid

(4)

(9)

at the open boundary. The overbar denotes vertical averaging. Appropriateness of this condition for the Northern Adriatic was investigated by KuzmiC and 0rliC.20 Numerical solution As pointed out in the introduction, one way to obtain better vertical resolution in numerical circulation models is to use some kind of expansion of the horizontal current components along the vertical. Heaps’O was apparently the first to suggest such an expansion based on the solution of an eigenvalue problem. He assumed a constant vertical eddy viscosity coefficient and solved the eigenvalue problem analytically. In a later paper” Heaps extended his method to a case in which vertical eddy viscosity was separately prescribed in two layers, explicitly resolving the bottom boundary layer. The vertical eigenvalue problem was again treated analytically. Heaps’* further extended the method to three

Appl. Math. Modelling,

1989, Vol. 13, March

179

A numerical

study of wind-induced

motions in shallow coastal seas: M. KuzmiC

layers prescribing a constant, but different, eddy viscosity coefficient in each. As in previous cases the eigenvalue problem was treated analytically, providing required coefficients and functions for the horizontal problem solved numerically, by integrating in time the transformed equations of motion and continuity. An alternative formulation of the problem was proposed by Davies,’ who obtained detailed vertical representation of the horizontal current components using the B-spline method. Davie? generalized his method to a case of vertically variable eddy viscosity using the Galerkin method and a basis set of cosine functions. In the same paper he showed the correspondence between his formulation and the one proposed by Heaps. Using the Galerkin method, Davies and Owen** explored the feasibility of the Chebyshev and Legendre polynomials as basis functions. Davies23,24 later applied the Galerkin method with eigenfunctions as the basis set. In doing so he essentially developed a twostep procedure in which the horizontal current components are expanded in terms of eigenfunctions, which in turn are expanded employing the B-splines. Davies and Furnes2s have recently developed a procedure of solving the three-dimensional hydrodynamic equations using the Galerkin technique and the Runge-KuttaMerson (RKM) method. The RKM method was used as an alternative to combining the Galerkin technique with an expansion of B-splines in order to solve the vertical eigenvalue problem. In the present work the integral transform framework proposed by HeapslO is retained, but the eigenvalue problem is solved numerically, after applying the Prtifer transformation.26 Detining the transformed current components as Tlu(a),

‘d

u(u)

+ 0

$+fu,=

N(u)a+r

--

h&

[

h

da

(14)

da 1

-ga,%+BV

ay

where BU

_

h(O)Tsx ; ; a;r) N;’ h

u(o)

P

*r(l) - -+u(l)

1 aw) wu -41) CT h

- j--

&(I)

1 am)

- hkU(l)

- j--

(16)

N(1) h 41)

(17)

Equations (14) and (15) could be further simplified if I,!J~were eigenfunctions of the appropriate SturmLiouville problem. The problem arises in many practical applications and can be generally stated as

$ [Pdq +[q(a) + Wu)I+(u)

= 0 (18)

with boundary conditions: &I

4(+)1= [4,~,1

NC+)

+

~*P(d

W(u)

7

o

atu = (+I

=

(19)

and

where &(c+) are kernels of the transformation and r is the modal index, and applying the transformation to the equations of motion, (2) and (3), one obtains

In the above equations Cc,is an eigenfunction of the Sturm-Liouville problem, and A is the corresponding eigenvalue; (Y], a2, pi, and p2 are arbitrary constants (in the nonsingular case considered here). The coefficient functions p(a), q(a), and s(u) are expected to be real and continuous, with p(u) and s(u) also being positive. In our case equations (18)-(20) take the form

$wy] cy, =0 krCl,(u)

where

+

J

ICl,(u)du

(13)

After twice integrating by parts and applying (4) and (5), we have

Appl. Math. Modelling,

1989, Vol. 13, March

atu=O

(22)

a*#0

N(u) W(u) du

o

atu=l

(23)

assuming N(u) = Non(u)

0

(21)

-_=

h

a, =

180

and

+h&(u)=0

(24)

In order to simplify expressions (16) and (17), we apply the condition &(O) = 1. Condition (23), based on the

A numerical

study of wind-induced

assumption of linear slip at the bottom, detined by (7), further simplifies expressions (16) and (17). It is computationally advantageous to solve this problem after applying the Prtifer transformation:

h(fl) = Adsin

%(a),

-drCl,(u)= p,(u)cos

n(u)

do

(25)

O,(u)

where e,.(u) and pL,(o) denote the phase and amplitude of &(c+), respectively. When (25) are applied to (21), one obtains two first-order equations:

W(u)

-

do

d/.44

-

da

1

&gsin’&(u)

= -cos2 n(c)

&(a) +

= p,(a)sin

1 B,(c~)cos O,(C) n(g)

0

(26)

1

(27) The phase function satisfies the boundary

conditions

(Y~sin 0,(O) + ff2 cos 0,(O) = 0

(28)

pi sin 0,( 1) + p2 cos I%(1) = 0

(29)

which are equivalent

MN = ff,

(Y =

to the simpler ones arctan

- -

,

ffl

O%a
h(l) = P +jK

- -P2

p = arctan (

o
PI

ing point, a numerical integrator with error estimation capability is required. The SLEIGN error assessment is based on the numerical integrator GERK,32 a Fortran code that uses Runge-Kutta-Fehlberg formulae and provides global error estimates by producing a full-step solution and two half-step solutions at each integration step. After applying the Priifer transformation, one can obtain the eigenvalues by solving only equation (26) with boundary conditions (30) and (31); to calculate the eigenfunction, one needs the amplitude and must solve equation (27). The SLEIGN code actually solves the amplitude equation after an exponential change of variable to avoid instability, stiffness, and scaling difficulties. Once the phase and amplitude functions are known, one can calculate (13) using (25). With the eigenvalue problem solved, expression (13) calculated, and expressions (16) and (17) simplified, we can return to the equations of motion, (14) and (15), and write them as follows: ah-

z

-

fv, =

-g&Z- hJ.4, + 5

(32)

+s

(33)

2 +fu, = -ga,-

( > a2

motions in shallow coastal seas: M. Kuzmid

(30)

,

-

ay

h,v,

=

r

where [[I# = $A$$(a) du. Replacing llll~# by cprand truncating the sum to M terms (an approximation to be justified in the verification section), we can write (34) as

[4uMa)l

= 5 [whlcp,~r(~)

(35)

i-=1

Finally, using (35), we obtain the transformed of continuity:

equation

Equations (32), (33), and (36) constitute the horizontal problem, which is solved numerically by using a forward-time, staggered-space, finite difference scheme originally proposed by HeapslO and used in our previous calculations. 13-15Notation for the grid is given in Figure 1. This particular schematization then yields the following 1 + 2M difference equations: S,(5), = - 5 [~,(h”&%&)i r=l

and integrating the corresponding variational equation. Since success of the root-finding process depends on accurate calculation of the phase function difference and phase function derivative difference at the match-

5 ‘ur7;$u) (34)

r= 1

(31)

The amplitude function satisfies the condition of never being zero in the interval 0 I (+ 5 1. Several authors have solved the Sturm-Liouville problem employing the phase/amplitude transformation. Godartz7 and Banks and Kurowskiz8 used the standard Prtifer transformation for solution of eigenvalue problems. Baileyz9 used the phase function to calculate Sturm-Liouville eigenvalues. The same author later proposed a slightly modified Prtifer transformation30 that offered computational advantages. The transformation is built into the SLEIGN code,31 which was used in solving the present problem. In the code the phase function equation is solved by shooting from the end values toward the interior. The discrepancy observed at the matching point measures the success in finding the eigenvalue and vanishes when the right one is found. A combination of bisection and Newton’s method is used to find that value. Newton’s method requires differentiating the phase function

ph

Using the complete set of orthogonal functions I/J,., we can apply the inverse Sturm-Liouville transform33 to obtain the current components: T- ‘[u,,v,] = [ u(u),v(u)]

)

0,1,2

al

&(uJ, = - XiA&,)r

+ fA&), I

Gs,i

phr ’

Appl. Math. Modelling,

+ $W’&cp:‘4,1

(37)

- ga,u,i&@)i+1 r=

l,...,M

1989, Vol. 13, March

(38)

181

A numerical

study of wind-induced

motions

in shallow

coastal

seas:

M. KuzmiC

bility of the scheme is assured by obeying the CFL criterion. Vertical eddy viscosity

AX _/

/

’

Oi+n /

lkYYTLX

Figure 1.

Finite difference grid used in the model. Symbols 0, 0 and x respectively denote places where the sea level, u component and u component of velocity are calculated. Solid, dash and dash-dot lines connect points used in calculating respective variables. Number of computational elements along x axis is n=O

+7ys.i

r= l,...,M

phy ’

(39)

The difference and averaging operators used above are S,(X), = 6

(xl,

x

Xi(t + At) - Xi(t) At ’ xi(t)

=

(x)

Y

xi-

ltt)

Ax

sx(x)i+I = s

-

I

=

Xi+I(t + At) - Xi(t + At) Ax

xi-m -

n

xi(t)

AY

Xi(t + At) - Xi+,(t + At) ayW>i+n

9

=

AY

9

AZ(X), = t[Xi(t + At) + X;(t)] Ad(

= $[Xi_l(t + At) + X<(t + At) + Xi+,_I(t + At) + Xi+,@ + At)1

A4(X)n = a[Xi-n(t) + Xi-,+ i(t) + Xi(t) + Xi+ I(t)1 Superscripts u and u denote, respectively, values taken at the u and u points on the computational grid. Sta182

Appt. Math.

Modelling,

1989, Vol. 13, March

In mathematical modeling of the mean fluid flow a necessity arises to represent turbulent shearing stresses due to the cross-correlation of fluctuating velocities. A simple concept, attributed to Boussinesq, relates the turbulent stress to the mean velocity gradient and a coefficient of proportionality called virtual, eddy, or turbulent viscosity. An obvious first choice is to keep this coefficient isotropic and time-independent-an approach followed by Ekman in his famous paper.34 Inadequacy of the constant-eddy-viscosity coefficient, particularly in conjunction with the no-slip bottom boundary condition, has long been recognized. The concept of slip bottom velocity was introduced as a remedy allowing an additional degree of freedom. Rapid advances in the field of computing technology have made feasible development of numerical models that incorporate various aspects of eddy viscosity variability (e.g., Heaps” or Davies35). Advanced numerical models of recent development incorporate vertical mixing employing various turbulence closure schemes (e.g., Blumberg and Mellot-’ or Svensson36). In the present paper, however, eddy viscosity is treated in a less aspiring way by considering just one additional degree of freedom within the Boussinesqian framework-an explicit depth dependence. As an immediate improvement over the fully isotropic and time-independent case, the eddy viscosity coefficient is modeled as a product of a constant part and vertically variable function (expression (24)). This is a choice in keeping with the stepwise refinement approach, rather than a suggestion that the spatial and temporal variability do not exist or are not important. Various functions have been advanced for modeling the vertical variability of the eddy viscosity coefficient. In this section some of those functions are briefly presented grouped around four general scenarios with a view to test the model capabilities and comparatively assess the consequences of a particular choice. The comparison was planned to obtain observational evidence from the Northern Adriatic. To stress this relation, we kept the values of No and k equal to those previously chosen.13 The following scenarios have been considered: Scenario A. Eddy viscosity coefficient constant throughout the water column. Only one case, from our previous work,i3 is included for comparison: Al

n(a) = 1

Scenario B. Eddy viscosity coefficient decreasing throughout the water column. With such a distribution the surface is recognized as a source of turbulence and (in shallow water) the bottom as an element of constraint. Empirical studies that substantiate such a distribution in shallow water seem to be rare. Apparently

A numerical

study of wind-induced

some support has been drawn from studies of similar problems. Klug,37 for example, derived from radon measurements a decreasing eddy viscosity profile in an upper oceanic layer. This scenario translated into simple functional form has often been treated analytically. Three cases are presented here: Bl

n(a) = ClU + c2

Linear change was used by Thomas38 and related to the idea of mixing length proportional to the distance from the bottom. B2

n(cr) = exp(c3cr)

Witten and Thomas39 used exponentially decreasing viscosity, assuming that it “reflects the expected physical features of turbulent momentum transport, without introducing severe mathematical complications.” Similar reasoning led Dobroklonsky40 to use it in a study of the drift currents in the surface layer of a deep sea. 12((T)= (1 + c4 (1 +

B3

ctj(+

+

c7

1

c4y5

forOIusu, foru,%u%

1

The n(u) grows linearly from a small value at the surface to a constant value below the surface layer (of depth u,). A distribution of this kind has been used by Pearce and Cooper* and by Davies.35 c2

n(u) =

csu + {1

c9

forO%uIu, foru,susub forub(crI

CIOU + Cl1

c3

n(u) =

1 +

Cl3

1

uYs

Scenario C. Constraining effect of the surface (airsea interface) and bottom (sea bed) boundaries is recognized in an analogy with the atmospheric boundary layer. Accordingly, eddy viscosity is assumed to grow with increasing distance from boundaries. Results of more sophisticated turbulence closure models of Svensson36 (for surface layer) and Mofjeld and Lavelle4* (for bottom layer) support such a conceptualization. Five cases are presented for this scenario. n(c+) =

to allow a decrease from a larger value at the surface to a constant value below the surface layer. The case is evidently an exception to this scenario; it has been included since it represents a counter concept to the previous one. It has been used by Davies.2 As in previous cases several values for ci were tried. We will take a closer look at the cases of two orders of magnitude smaller (larger) surface viscosity increasing (decreasing) to a constant value. The increase (decrease) was chosen symmetrically and offered a reasonable basis for comparison. Of course, more careful choice is needed in any attempt to approximate reality. Required coefficient values were cg = 4.95, c7 = 0.01, cs = -495, and c9 = 100. In both cases the depth of surface layer was set equal to 0.2, after Pearce and Cooper.8

ClZ(+ -

Fjeldstad4* arrived at the power law distribution by analyzing some North Siberian shelf data. For each of the above functions several values of the coefficients ci were tried. In this paper the coefficients were chosen to yield at the bottom viscosity coefficient an order of magnitude smaller than the surface value. The choice of such a change was deemed appropriate for the purpose of comparison. Consequently, everything but the form of n(o) was kept unchanged, allowing better assessment of the flow sensitivity to the form of n(g). The ci values were c, = -0.9, c2 = 1.0, c3 = -2.3, c4 = 0.005, and c5 = 0.434 (the c5 coefficient was set to Fjeldstad’s original value of 0.75 in the fourth experiment when the B3 case was compared to the D cases).

Cl

motions in shallow coastal seas: M. KuzmiC

forO%uIu, foru,sus

1

The formal description of n(u) is the same as in the previous case, but values of the coefficients are chosen

The n(u) increases linearly from the boundaries toward the interior, staying constant in the interior region (a, and (Tb mark the layers). This shape of n(u) has been also used previously (e.g., by Davies43). C4

n(a) =

c14(T2

+

c15(+

+

c1,j

forO117%

1

The layers at both ends are modeled without explicating their extent and assuming a parabolic shape of n(u). Etling et ~1.~ obtained a similar profile of N(a), decreasing toward the ends from a maximum near the middle, using the mixing length approach. C5

n(u) = -i

Cl’(+

+

Cl8

C19@

+

c20

forOIus(+, forumIu%

1

In this case the eddy viscosity increases linearly from both ends to some point in the interior (u,) without forming a layer of constant viscosity. The case is based on the eddy viscosity model put forward by Madsen,45 but without iterating to obtain the bottom stress. After experimenting with the ci values, we chose three cases for presentation in this subgroup. To facilitate comparison the coefficients for the selected cases read: cl0 = 4.95, cl1 = 0.01, ~12 = -4.95, ~13 = 4.96, -3.96, cl5 = 3.96, c]6 = 0.01, cl7 = 1.98, Cl4 = Cl8 = 0.01, C]9 = - 1.98, and c20 = 1.99. As before us = 0.2 and SymmetriCally wb = 0.8, urn was Set equal to 0.5 for simplicity. The ci values were selected to assure a symmetrical reduction of n(u) to 0.01 at both ends. The forms of n(u) described in the three scenarios are summarized in Figure 2. Another scenario has been also formulated following somewhat different logic. It can incorporate cases similar to the ones already described although no explicit spatial schematization is applied. Scenario D. Vertical variability is derived directly from observational data. Variability can be rather general and is commonly presented in discrete, tabular form. Modem instruments like the Doppler current profiler can supply data for this scenario. Appl. Math. Modelling,

1989, Vol. 13, March

183

A numerical study of wind-induced

motions in shallow coastal seas: M. KuzmiC

periments tion.

0.07

we will briefly address the model verifica-

02 1

Model verification

04.

G

A mathematical model is an approximate representation of a natural or artificial system. Its symbolic nature enables replication of the system dynamics on the computer and an assessment of its behavior under various conditions. To trust the results of such a simulation, however, one must first ensure that the computerized model indeed exhibits the anticipated behavior with acceptable accuracy, i.e., verify the model. At later stages the model responses must be compared with observational data to substantiate its validity. Model validation is not a subject of this paper. It is only tentatively considered in a later section. Several tests were performed, however, to verify the model. First, to verify solutions of the vertical eigenvalue problem, we compared numerical values of h,, er, a, and qo, with their respective values obtained analytically for the first 20 eigenfunctions of case Al. Results of the comparison are summarized in Table 2. One can observe satisfactory agreement between numerical and analytical values, although some differences do appear. When used in the model, to obtain the sea-level and current predictions, numerical and analytical values of those parameters gave identical results within practical limits of submillimeter and submillimeter per second, respectively. It should be noted that all calculations were performed on a 16-bit minicomputer with single precision (six significant digits). To verify predictions for variable N a reduced computational setup (straight boundaries, flat bottom, and uniform wind stress) was further simplified by turning off the rotational effect (f = 0). The analytical solution of the problem was obtained by further assuming a

06,

W 10 I

Al

61

82

83

Cl

Figure 2. Vertical distributions of eddy viscosity coefficient considered in numerical experiments. Dash line indicates the B3 distribution used in the fourth experiment (comparison with the D cases)

With no appropriate Northern Adriatic data at hand some other historical data are used to illustrate this kind of eddy viscosity formulation. Analyzing some North Siberian shelf current data, Sverdrup46 found an example of wind-induced currents in shallow (22 m), practically homogeneous water and concluded that “the change of the current with depth must be an effect of the eddy viscosity.” Rearranging equations of motion and using a graphical method, he computed the vertically variable eddy viscosity; his assistant independently repeated the calculation. We therefore consider two cases here: Dl Sverdrup’s derivation of vertically variable eddy viscosity from a set of current measurements46 D2 Mosby’s derivation of vertically variable eddy viscosity from the same data set46 The data for both cases are presented in Table I. All described cases were-used to perform several computational experiments. Before presenting the ex-

Table 1. Vertical distribution of eddy viscosity derived by Sverdrup Mosby (D2) from observed wind current Depth (m) Dl viscosity (m*.s-‘) D2 viscosity (m*.s-‘)

0

5

10

15

20

0.0335 0.0196

0.0223 0.0213

0.0198 0.0213

0.0181 0.0178

0.0074 0.0080

for convenience,

water

Original cgs values are expressed density of 1000 kg.m-3.

Table 2.

Comparison

of analytically

r 2 4 6 8 10 12 14 16 18 20

x x x x x x x x x x

1O-4 10m4 1O-3 1om3 1O-3 10m3 10-Z lo-* 10m2 IO-2

For each parameter

184

and numerically

in mks units assuming,

obtained values of some model parameters

A, 1.1587 6.5028 1.6523 3.1390 5.1160 7.5851 1.0547 1.4002 1.7950 2.2392

Appl.

Math.

a, 1.1587 6.5027 1.6523 3.1389 5.1159 7.5848 1.0547 1.4002 1.7950 2.2392

x x x x x x x x x x

lo-“ lO-4 1O-3 1O-3 10-z 1O-3 10-z 1O-2 1o-2 1o-2

-0.2133 - 0.0686 - 0.0322 -0.0182 -0.0115 - 0.0079 - 0.0058 - 0.0044 - 0.0034 - 0.0028

analytical solution is in the left and numerical

Modelling,

1989,

(Dl) and

Vol. 13, March

for the Al case *All

cpr -0.2134 - 0.0686 - 0.0322 -0.0182 -0.0115 - 0.0079 - 0.0058 - 0.0043 - 0.0034 - 0.0027

1.8444 1.9065 1.9466 1.9673 1.9785 1.9849 1.9889 1.9915 1.9933 1.9946

in the right column.

1.8442 1.9067 1.9468 1.9676 1.9785 1.9846 1.9889 1.9910 1.9939 1.9791

- 0.3955 -0.7141 - 0.8518 -0.9132 - 0.9440 - 0.9612 -0.9716 - 0.9784 - 0.9830 - 0.9863

- 0.3955 -0.7140 - 0.8517 -0.9131 - 0.9437 -0.9610 -0.9713 - 0.9782 - 0.9924 - 0.9858

A numerical

0

study of wind-induced

motions in shallow coastal seas: M. KuzmiC

tive equal zero, imposed by (22), contradicts the assumption of nonzero surface stress without ruining the solution but slowing it down near the surface boundary. Heaps circumvented the problem for the case of constant N by introducing a corrective term in the finite inverse transformation, thereby approximating the completion to infinity. Davies43 encountered the same problem near the surface when expanding the solution by the Galerkin method over otherwise rapidly converging eigenfunctions. He proposed23 a way to overcome this difficulty by assuming the nonzero vertical derivative proportional to the eigenfunction itself-a condition that corresponds to assuming nonzero (Y~in the present formulation. A nonzero value of this coefficient allows for a tuning which should improve the convergence in the way reported by Davies. In all simulations presented in a later section, (Y~was kept equal to zero. For the purpose of cross-comparison of the model predictions for various N(a) and tentative comparison with some Northern Adriatic current data, it was deemed sufficient to run the model with M = 20 bearing in mind the slower convergence near the surface. The choice of larger M enabled model decomposition of the predicted current at selected points (into first five modes, next five, etc.) and thereby offered an insight into the role the higher modes play when eddy viscosity varies with depth. It also helped develop a feeling for the convergence problem at the surface. For the purpose of verification, the present model was also used to obtain predictions for the Northern Adriatic with the parameters set as in a previously published case.i3 Successful comparison with predictions obtained earlier with another, constant-N model and the results presented in the Tables 2 and 3 were taken as an adequate evidence to warrant the next step, production runs.

b

Figure 3. Extended (a) and reduced (b) computational domains used in numerical experiments. Letters J and K indicate, respectively, computational elements in reduced and extended domains used in comparisons

steady state. Numerical values of various parameters needed in simulations are given in the next section, before the experiments are presented. Let us only mention here that the verification runs were done assuming 7”s = 0.15 Nom-*, T,, = 0.0 Nom-*, with default values for other parameters. Verification calculations were performed for several n(cr) functions-a particular case of C 1 was chosen for presentation. Analytical and modelgenerated values of the sea-level displacement along the longitudinal (northwestward, downwind) axis and velocity distribution along the vertical at point J (Figure 3) are comnared in Table 3. One can see from the table’that the sea-level displacements agree to within a few tenths of a millimeter. Predictions of the velocity distribution are also in rather good agreement, differences being about 1 mm*s’. There is one exception, however: the behavior of the solution near the surface. The surface, as well as bottom, boundary condition is satislied exactlv onlv with the inverse transformation completed to infinity (equation (34)). Satisfactory accuracy can be achieved with a finite transformation (equation (35)), provided convergence of the solution is fast. It does prove to be rather fast except near the surface, where even 20 modes are not enough to match the analytically obtained values. As pointed out by Heaps4’ the condition that the eigenfunction deriva-

Computational

experiments

After verifying the model, we carried out computational experiments to analyze several characteristic cases of vertically variable eddy viscosity and to further test model capabilities. Of particular interest was possible relevance of those simulations for explaining the winter, barotropic dynamics of the Northern Adriatic. Five numerical experiments are presented in the section. One of them was performed using a computational domain that resembles the geometry of the Northern

Table 3. Comparison of analytically and numerically obtained values of the sea-level displacement along the longitudinal velocitv distribution alona vertical axis in the comoutational element J of the reduced model for case Cl Sea-level displacement (cm) Computational element, 1. downwind Analytical solution 0.00 Numerical solution 0.00 Vertical distribution of current (cm s’) Normalized depth 0.0 Analytical solution 57.79 Numerical solution 50.58

3.

5.

7.

9.

0.83 0.83

1.66 1.66

2.49 2.48

3.31 3.31

4.14 4.14

0.6 - 5.98 - 5.90

0.8 - 6.02 - 5.96

1.0 - 2.58 - 2.68

0.2 4.50 4.56

0.4 -2.48 -2.51

Appl.

11.

13.

15.

4.97 4.97

Math. Modelling,

5.80 5.80

axis and

17. 6.63 6.62

1989, Vol. 13, March

19. 7.46 7.45

185

A numerical

study of wind-induced

motions in shallow coastal seas: M. KuzmiC

Adriatic (this domain will be referred to as extended); another three were done on a reduced one; in one experiment only the eigenvalue problem was considered. Boundaries of the Northern Adriatic were schematized, as in previous studies,r3-I5 to fit a field of 29 x 22 rectangular, computational boxes of 7.5 km in the x (northeastward) and y (northwestward) directions. To speed up the execution, most of the comparison runs were done on the reduced domain of 9 x 19 boxes of the same size (see middle of the Figure 3). The Northern Adriatic bottom topography is indicated in the figure for completeness, but in the experiments we assumed a flat bottom with average depth 40 m in both domain models. We have modeled the bottom slope contribution to the Northern Adriatic current tield previously;13 comparison shows that its inclusion is not critical for the tentative assessment to be presented later. With selected depth and grid size the CFL stability criterion was satisfied with a time step of 2 min. Values of 1.04 x 10e4 s-i, 1025 kg*mp3, and 9.81 rn.s’ were respectively selected for f, p, and g in both domain models. Also, both domain models were forced with suddenly applied wind stress kept constant in time. The stress was calculated for CD = 2.5 x 1O-3 and pU = 1.247 kg.mp3 with a NE wind (locally known as bura), as in previous Northern Adriatic studies. In the reduced domain experiments, spatial homogeneity of the stress field was maintained. With a wind speed of 5 mess’, T,~,= -0.078 N*m-* and 7sY= 0. Since the importance of the Northern Adriatic wind field heterogeneity has been demonstrated previously,‘4~‘5 the homogeneity assumption was abandoned in the extended domain runs. Two heterogeneities were considered: one based on lo-year monthly mean statistics for six stations along the Northern Adriatic coast (see Refs. 14, 15), and the other synthesized to reflect the empirical evidence of bura being weaker near Rovinj than near Pula and Trieste (see Figure 3). Details of the applied stress heterogeneities are summarized in Table 4. Note that in both cases the wind direction is that of bura. The reference value of eddy viscosity, No, was kept equal to 10e2 rn*.s-l in accord with our previous experiments. The only exception was the fourth experiment, in which a non-Adriatic situation was considered. In cases B3, Dl, and D2, we assumed NO = 0.0385, 0.0335, and 0.0196 m2.ss1, respectively. The

heterogeneities

value k = 2.5 x lo-3 m.s-r was also retained from our previous modeling studies. Reduced domain calculations were carried out for a simulated 24 h-long enough to obtain major response features, and 72 h was used to obtain representative response in the Northern Adriatic simulations. In all runs M = 20. For each run the fields of sea-level displacement, vertically averaged current, surface current, and bottom current were plotted. For selected computational boxes vertical distribution of current was also calculated and displayed in several forms. In the rest of this section part of that information is used to present and analyze the experiments. In the first, reduced domain, experiment the B cases were used to assess the sensitivity of the flow to the shape of n(cr) within the framework of continuously decreasing eddy viscosity. The first five eigenfunctions that correspond to those cases are displayed in Figure 4(a). An almost identical first mode is readily observed. The higher modes exhibit rather consistent pattern: zeros of the Bl functions lie between the B2 and B3 zeros. Vertical distributions of current (in the computational element J, Figure 3) for these three cases are given in Figure 4(b) in terms of magnitude and deflection angle. The magnitude is expressed as a percentage of the wind speed, and deflection angle is given relative to the wind direction. One can observe significant qualitative (and also quantitative) similarity among these cases of linear (Bl), “concave” (B2), and “convex” (B3) eddy viscosity decrease. Still, the most pronounced decrease of the eddy viscosity coeflicient (B2) can be related to the most pronounced magnitude (particularly in the upper part of the water column). The Coriolis effect is also most notable in this case. The purpose of the next experiment was to assess the flow alternation due to the changes of N in the surface layer. In particular, the case of constant N (Al) has been compared with the cases of increasing (Cl) and decreasing (C2) N. The first five eigenfunctions for the three cases are contrasted in Figure 5(a). Some difference can be traced in agood part of the first mode, and higher modes exhibit pronounced discrepancy. The C2 eigenfunctions, in particular, are almost constant in the upper 20% of the water column (assumed depth of the surface layer)-a feature responsible for the uniformity of the predicted current response in the same segment. The current response for these cases (as before, in computational element J) are given in Figure

Table 4.

Simulated

in the wind stress field over the Northern Adriatic

i

Slfi)

Wj)

i

Slfi)

S2(il

i

Sl(i)

S2(i)

i

Sl(i)

S2(i)

1 2 3 4 5 6

1.19 1.00 0.85 0.73 0.65 0.58

1.43 1.32 1.19 1.06 0.90 0.71

7 8 9 10 11 12

0.56 0.51 0.69 0.88 1.15 1.35

0.45 0.32 0.45 0.55 0.71 0.78

13 14 15 16 17 18

1.43 1.47 1.47 1.42 1.34 1.23

0.86 0.96 1 .Ol 1.08 1.15 1.19

19 20 21 22

1.08 0.90 0.69 0.42

1.25 1.30 1.38 1.43

Sl : as abstracted from climatological data; S2: as synthesized to reflect empirical evidence of bura being weaker near Rovinj, then near Pula and Trieste. The stress applied along the x axis and variable along the y axis is calculated using the relation rXs = -0.15 S2( j) while 7Ys = 0. Points j = 2, 6, and 12 correspond respectively to the cities Pula, Rovinj, and Trieste marked in Figure 3.

186

Appl.

Math.

Modelling,

1989,

Vol. 13, March

A numerical

study of wind-induced

motions in shallow coastal seas: M. Kuzmid Elgenfunctions

Elgenfunctlons -2 -

r=l

2

r=2

r=3

r-4

I!

/,

3

02

4\ ‘\ ’ \\

I\\ (\ ‘\

I’ II:
0.6

Q

Ii 1.0

’ \

1

‘\ 1:

PIY7?!kI ----

\ ‘\ ./ /:

‘\

‘I

/’ /,’

0.L

0.6

,, I’

.\ ,/:

,I ‘\ _\

01 02 03

Figure 4(a). Eigenfunctions corresponding to the Bl, B2, and 83 distributions of the vertical eddy viscosity

Magnitude

r=5

r=3

r-2

r=5

Deflection

(%I

Figure 5(a). Eigenfunctions corresponding to the Al, C2 distributions of the vertical eddy viscosity

Magnitude

angle (degrees)

(%)

Cl, and

Deflect ion angle (degrees)

O.L-

G 0.6-

-

81

---

82

-.-

83

1 1.0

Figure 4(b). Vertical distributions of current for the Bl, B2, and 83 cases, in computational element J, given in terms of magnitude, expressed as a percentage of the wind speed and deflection angle relative to the wind direction

Figure 5(b). Vertical distributions of current for the Al, Cl, and C2 cases, in computational element J, given in terms of magnitude, expressed as a percentage of the wind speed and deflection angle relative to the wind direction

5(b). Limited vertical exchange of momentum (Cl) allows for larger velocities and, consequently, more pronounced veering than with constant N. On the other hand, an eddy viscosity coefficient decreasing from a large value at the surface permits pronounced exchange of momentum and consequent development of a slab of uniform velocity and “synchronous” veering. We note that in the Cl case the near-surface magnitude is considerably larger than in the other two cases despite the convergence problem which tends to diminish the near-surface currents. In the lower part of the water column the three predictions converge under the influence of the same values of N and k. These responses are influenced by selected numerical values of various parameters. An example of such a parameter is the surface layer depth. Its influence on the surface current has been studied, for example, by Davies.48 In the third experiment the three cases in which both the surface and bottom layers are recognized (C3, C4, and C5) are compared. Remember that the n(a) for the three cases differ in the rate of change within the layers and in the depth over which the change takes place. All three functions were set symmetrically to the middle. The three cases have rather similar first eigen-

function. At higher modes, the similarity between the C3 and C4 is more pronounced. Vertical distributions of current (computational element J) for these cases are compared in Figure 6. With overall appearance much the same, one can observe progressively higher velocities and more pronounced veering when going from C3 to CS. The similarity observed at the eigenfunction level (not shown) is readily observed in this figure. The fourth experiment differs somewhat from the previous three. It compares, at the eigenfunction level only, three vertical eddy viscosity formulations derived by different authors (Fjeldstad (B3), Sverdrup (Dl), and Mosby (D2)) from the same empirical data. Detailed modeling reconstruction of the region where the data were taken (the North Siberian shelf), although representing an interesting problem, was deemed out of the scope of the present study. Since the D cases were defined in tabular form the experiment was also used to test such a formulation of eddy viscosity variability. When solving the eigenvalue problem for the D cases, we used a cubic spline interpolation procedure49 to calculate appropriate n(o) values. It was also necessary to obtain the bottom n(c) values (between 20

Appl.

Math. Modelling,

1989, Vol. 13, March

187

A numerical

study of wind-induced

motions in shallow coastal seas: M. Kuzmid

Y

4

NE ---NW

Bm

n

b

Figure 6. Vertical distributions of current for the C3, C4, and C5 cases in the computational element J: (a) perspective view; (b) horizontal projection 19-27 3 1982

and 22 m), which was accomplished using the method of local procedures. 5oThe results of these calculations are presented in Figure 7. As in previous comparisons the first five eigenfunctions are given. The observed discrepancies reflect the differences exercised in constructing the eddy variability functions. The Fjeldstad formula has two “free” coefficients allowing for tuning, but that was not attempted in the experiment. The figure also testifies to successful treatment of tabular formulation of eddy viscosity variability. Spline-interpolated discrete values presented no computational difficulties. The final experiment was the only one performed on the extended domain in an attempt to assess the relevance of previous numerical experiments for the Northern Adriatic. Namely, in several field studies conducted in the area a consistent pattern in the sea response to bura forcing has been observed at the station offshore of the city of Rovinj (the station location corresponds to the computational element K marked in Figure 3). This observational evidence is summarized in Figure 8, in which four characteristic bura episodes and sea current responses are presented. The wind data were collected at the Pula Airport at the Eigenfunctions r-2

r=3

r=4

Figure 7. Eigenfunctions corresponding to the B3, Dl, and D2 distributions of the vertical eddy viscosity

188

Appl.

Math. Modelling,

1989, Vol. 13, March

11-I8 4 1982 5 - 1011985. 2712 1984-L 11985

Figure 8. Characteristic intervals of wind and currents (registered at 8 m and 35 m) taken from low-pass filtered time series measured at Pula Airport (wind) and the station off Rovinj (currents) during the MEDALPEX (1982) and ASCOP (1984185) cruises

height of 73 m (10 m above ground), and the currents were registered at depths of 8 and 35 m in water 40 m deep. Both wind and current data were resolved into NE and NW components and low-pass filtered to remove tidal and high-frequency components. The current intervals were taken from longer series of measurements carried out within the framework of the Mediterranean Alpine Experiment (MEDALPEX) and Italo-Yugoslavian Adriatic Scientitic Cooperation Programme (ASCOP). In the described figure one can readily observe rather prompt current response directed upwind at both measurement levels. This dynamics has been explainedi in terms of frictionally controlled flow. A curl in the wind field over the area, also explored before,15 can be taken to explain the buildup of water along the Italian coast and the consequent prominence of the slope current in the observed response. The influence that vertically variable viscosity might have on that response is tentatively considered here. To that end the current data from Figure 8 were further decimated and then vector averaged to obtain a pair of vectors (one for each depth) suitable for comparison. A similar procedure was employed previously. i3,i5 An additional requirement was imposed on decimated current magnitudes: only those larger than 10 cm-s-’ were used. This pair of empirical vectors is presented at the top of Figure 9. Keeping in mind the frictional adjustment time for the case considered, these vectors can be taken to represent the steady-state response. One aggregated pair was produced instead of four separate ones because the purpose of this comparison was not a detailed explanation of particular episodes in the Northern Adriatic dynamics but consideration of the role of variable eddy viscosity. The rest of the figure is used to present mod-

A numerical

study of wind-induced

motions in shallow coastal seas: M. KuzmiC

heterogeneity) in hydrodynamic models of the Northern Adriatic and justifies further efforts along this line.

b

Figure 9. Comparison of model predictions of current vectors at two depths (8 m and 36 m) in computational element K with the ones obtained empirically (depths of 8 m and 35 m): (a) empirical data; (b) model predictions after 24 h; (c) model predictions after 48 h; (d) model predictions after 72 h. Three eddy viscosity distributions (Al, Bl, and C5) are combined with two wind stress fields (Sl, S2)

eling results (for computational element K in Figure 3). The adequacy of four combinations of eddy viscosity variability and wind heterogeneity was tested. The first one, Al-Sl, was a combination of constant viscosity and “climatological” curl (S 1, Table 4) used in previous studies. The same stress field was also combined with the Bl and C5 cases in the next two runs. In the fourth run the CS viscosity case was combined with another stress field (S2 in Table 4) to indicate significance of the curl. For each case a progression of states is shown spanning a period of three simulated days. It is readily apparent that all cases are in some way inadequate. Cases with Sl stress mispredict general direction of the current response. Furthermore, as observed in our previous study,is constant eddy viscosity leads to an underprediction of current magnitudes at both levels. Simultaneously, the upper-to-lower current ratio is overpredicted as is the angle between the current vectors. Assumption of linearly decreasing viscosity improves the current ratio while overpredicting the magnitudes and somewhat affecting the angular difference. Recognition of linear growth of viscosity away from both the surface and bottom boundaries causes further overprediction of magnitudes, but reduces the angular difference. It is quite obvious from the figure that it takes a change in the stress field to affect the general direction of the current vectors. As in reduced domain studies, we should remember that our results were obtained with a particular set of numerical values of various parameters; the experiment did not include selection of the most appropriate parameter values for any of the combinations. We should also remember that a flat bottom was assumed in all simulations. The experiment nevertheless does demonstrate the significance of incorporating vertically variable eddy viscosity (in addition to the wind stress

Conclusions A three-dimensional numerical model of wind-induced motions in a homogeneous sea has been constructed by using the spectral eigenfunction decomposition pioneered by N. S. Heaps. Besides accommodating realbottom topography and wind stress heterogeneity, the model allows rather general formulation of vertically variable eddy viscosity. The feasibility of a particular numerical treatment of the equations of motion, within the spectral decomposition framework, has been explored. It has been demonstrated, and verified by comparison with analytical solutions, that the approach based on a phase-amplitude transformation of eigenfunctions yields meaningful results of satisfactory accurac y . The approach is similar to the Galerkin-RKM formulation of Davis and Furnes. The authors used the Galerkin method to formulate a transformed problem and eigenfunctions as the basis set. The RKM method was used to obtain the eigenvalues, eigenfunctions, and related integrals (a, and cpJ. In the present approach an integral transformation is applied to the equations of motion with eigenfunctions as kernels of the transformation-a step that in the Galerkin formulation corresponds to setting the inner product of the weighted residual to zero. Thus the generated eigenvalue problem is solved after applying the Prtifer transformation. Combined with the shooting method, the transformation offers a reliable procedure to locate the eigenvalues and subsequently to calculate the eigenfunctions. The same numerical integrator used to solve the Prtifer equations is also used to obtain a, and cpr. Both approaches require solution of an eigenvalue problem and use of a numerical integrator, but differ in the choice of equations to integrate, in selection of Runge-Kutta integrator, in the eigenvalue finding procedure, and in error control. In both approaches, the solution of the eigenvalue problem contributes to accurate calculation of eddy-viscosity-influenced vertical modal structure of coastal sea flows. Solution procedures most likely differ regarding computational complexity and/or accuracy. To obtain quantitative evidence of such differences, one should compare the models at the execution level for various physically representative cases. Generally, with ever increasing computing power one can afford more elaborate procedures if accuracy or some other desirable feature is gained. The present model allows for conceptual and computational improvements. Implementation of the nonzero vertical derivative, as discussed in the numerical solution section, is an option offering improved surface current prediction. Unfortunately, general lack of reliable measurements in the near-surface zone impairs verification of such an effort. As in all previous North-

Appl.

Math.

Modelling,

1989,

Vol. 13, March

189

A numerical

study of wind-induced

motions in shallow coastal seas: M. Kuzmid

ern Adriatic studies, a linear bottom slip was also assumed in the present work and offered good results. It would be interesting, however, to explore the differences that quadratic slip introduces. To give the comparison proper weight, one should have adequate empirical evidence; in the Northern Adriatic a bottom boundary layer current experiment is still something to be desired. An important further consideration within the barotropic framework is treatment of the coastal boundary layer. After verification, the model was used to analyze several formulations of vertically variable eddy viscosity grouped around four scenarios: constant eddy viscosity, eddy viscosity, decreasing throughout the water column, eddy viscosity influenced by the surface and bottom layers, and eddy viscosity derived from observational data. Simulations within the second scenario have shown significant qualitative and quantitative similarity. The most pronounced difference in the current response has been observed in experiments within the third scenario when the constraining effect of the air-sea interface was first recognized and then ignored. An additional object of numerical experimenting was comparative assessment of vertical viscosity variability having in mind some observational evidence from the Northern Adriatic. In the light of that evidence the present model offers a new dimension to interpretation of available empirical data. In particular, tentative analysis suggests that, with an appropriate choice of numerical values, vertically variable eddy viscosity function can contribute to a more successful match of model predictions and empirical data. Numerical experiments also demonstrated the importance of the heterogeneous wind stress field for the general direction of the current response in the studied area. It will take further numerical experimenting as well as field studies to obtain more conclusive results concerning the Northern Adriatic. The work in progress along both lines will hopefully contribute to that end.

3 4 5 6 7

8 9 10 11 12

13

14 15

16 17 18 19

Acknowledgments The author wishes to thank Mr. Mirko OrliC for his careful reading of the typescript and many discussions in the course of the study. The Hydrometeorological Office of the Croatia SR kindly provided the wind data. The MEDALPEX and ASCOP current data were collected in cooperation with the Rovinj Center for Marine Research. One of the reviewers kindly drew the author’s attention to the Davies and Furnes paper. This work was supported by the Self-Management Community of Interest for Scientific Research of the Croatia SR.

20

Rapp. Comm. Int. Mer. Medit. 1985, 29(3), 75-78

21 22

1 2

190

25

Blumberg, A. F. and Mellor, G. L. Diagnostic and prognostic numerical circulation studies of the South Atlantic Bight. .I. Geophys. Res. 1983, 88, 4579-4592 Davies, A. M. Application of sigma coordinate sea model to

Appl. Math. Modelling,

1989, Vol. 13, March

Davies, A. M. On formulating a three-dimensional hydrodynamic sea model with an arbitrary variation of vertical eddy viscosity. Compur. Mefh. Appl. Mech. Engng. 1980,22, 187-211 Davies, A. M. and Owen, A. Three dimensional numerical sea model using the Galerkin method with a polynomial basis set. Appl. Math. Modelling

23

24

References

the calculation of wind-induced currents. Continental ShelfRes. 1985, 4, 389-423 Simons, T. J. Development of three-dimensional numerical models of the Great Lakes. Sci. Ser. No 12, 1973, Canada Centre for Inland Waters Leendertse, J. J. and Liu, S.-K. A three-dimensional model for estuaries and coastal seas. Vol. II. Aspects of computation. RAND Corp. Rep. R-1764-OWRT, 1975, Santa Monica Jelesnianski, C. P. Bottom stress time-history in linearized equations of motion for storm surges. Mon. WearherRev. 1970, 98, 462-478 Forristall, G. Z. Three-dimensional structure of storm-generated currents. J. Geophys. Res. 1974,19, 2721-2729 Davies, A. M. The numerical solution of the three-dimensional hydrodynamic equations using a B-spline representation of the vertical current profile. In Bottom Turbulence ed. J. C. J. Nihoul. Elsevier,Amsterdam, The Netherlands, 1977, pp. l-25 Pearce, B. R. and Cooper, C. K. Numerical circulation model for wind induced flow. J. Hydraulics Div. ASCE 1981, 107, 285-302 Davies, A. M. On extracting current profiles from vertically integrated numerical models. Coastal Engng. 1987,11,445-477 Heaps, N. S. On the numerical solution of the three-dimensional hydrodynamic equations for tides and storm surges. Mem. Sot. R. Sci. Liege 1972, I(6), 143-180 Heaps, N. S. Three-dimensional model for tides and surges with vertical eddy viscosity prescribed in two layers. I. Mathematical formulation. Geophys. J.R. Astr. Sot. 1981,64,291-302 Heaps, N. S. Development of a three-layered spectral model for the motion of a stratified shelf sea. I. Basic equations. In Physical Oceanography of Coastal and ShelfSeas ed. B. Johns. Elsevier, Amsterdam, The Netherlands, 1983, pp. 387-400 KuzmiC, M., OrliC, M., Karabeg, M., and Jeftic, Lj. An investigation of wind-driven topographically controlled motions in the Northern Adriatic. Estuarine Coasral Shelf Sci. 1985, 21, 481-499 OrliC, M., KuzmiC, M., and VuCak, Z. Wind-curl currents in the Northern Adriatic and formulation of bottom friction. Oceanol. Acta 1986, 9, 425-431 KuzmiC, M. and OrliC, M. Wind-induced vertical shearing: ALPEX/MED-ALPEX data and modeline. exercise. Annal. Geophys. 1987, 5B, 103-112 Phillips, N. A. A coordinate system having some special advantages for numerical forecasting. J. Meteorology 14,184-185 Paul, J. F. and Lick, W. J. A numerical model for thermal plumes and river discharges. In Proceedings of the 17th Conference on Great Lakes Research, 1974, pp. 445-455 Forristall, G. Z., Hamilton, R. C., and Cardone, V. J. Continental shelf currents in tropical storm Delia: Observations and theory. J. Phys. Oceanogr. 1977, 7, 532-546 Kuzmic, M. and OrliC, M. Predicting the currents in the Northern Adriatic and the problem of ill-defined wind forcing. Geofizika 1987, 4, 137-145 KuzmiC, M. and OrliC, M. A study of the influence of open boundary conditions on the predictions of a wind-driven model.

26

1979, 3, 421-428

Davies, A. M. Application of a Galerkin-eigenfunction method to cornouting currents in homogeneous and stratified seas. In Proceedings-of the IMA Conference on Computational Fluid Dvnamics. Readina Univ., 1982. pp. 287-301 Davies, A: M. Formulation of a linear three-dimensional hydrodynamic sea model using a Galerkin-eigenfunction method. Inr. J. Num. Methods Fluids 1983, 3, 33-60 Davies, A. M. and Fumes, G. K. On the determination of vertical structure functions for time-dependent flow problems. Tellus 1986, 38A, 462-477 Prtlfer, H. Neue Harleitung der Sturm-Liouvilleschen Reihenentwicklung Stetiger Funktionen. Marh. Ann. 1926,95,499-518

A numerical 27 28 29 30 31 32 33 34 35 36 37 38 39

study of wind-induced

Godart, M. An iterative method for the solution of eigenvalue problem. Math. Comput. 1966, 20, 399-406 Banks, D. 0. and Kurowski, G. J. Computation of eigenvalues of singular Sturm-Liouville systems. Math. Comput. 1968, 22, 304-310 Bailey, P. B. Sturm-Liouville eigenvalues via a phase function. J. SIAM Appl. Math. 1966, 14, 242-249 Bailey, P. B. An slightly modified Prtifer transformation useful for calculating Sturm-Liouville eigenvalues. J. Computut. Phys. 1978, 29, 306-310 Bailey, P. B. SLEIGN An eigenvalue-eigenfunction code for Sturm-Liouville problems. SAND77-2044, Sandia Laboratories, Albuquerque, N.M., 1978 Shampine, L. F. and Watts, H. A. GERK: Global error estimation for ordinary differential equations. ACM Trans. Math. Soffwnre 1976, 2, 200-203 Johnson, D. E. and Johnson, J. R. Mathematical Methods in Engineering and Physics. Prentice-Hall, Englewood Cliffs, N.J., 1982, pp. 229-241 Ekman, V. W. On the influence of the earth’s rotation on ocean-currents. Ark. Mat. Astr. och. Fys. 1905, 2, 1-53 Davies, A. M. A three dimensional modal model of wind induced flow in a sea region. Prog. Oceunogr. 1985, 15, 71-128 Svensson, U. The structure of the turbulent Ekman layer. Tellus 1979, 31, 340-350 Klug, W. Eddy viscosity profile in upper oceanic layers derived from radon measurements. Tel/us 1974, 26, 36-38 Thomas, J. H. A theory of steady wind-driven currents in shallow water with variable eddy viscosity. J. Phys. Oceunogr. 1975, 5, 136-142 Witten, A. J. and Thomas, J. H. Steady wind-driven currents in a large lake with depth-dependent eddy viscosity. J. Phys. Oceunogr.

40

41 42

43

44

motions in shallow coastal seas: M. KuzmiC Dobroklonsky, S. V. Drift currents in the sea with an exponentially decaying viscosity coefficient (in Russian). Okeunologiyu 1969, 9, 26-32 Fjeldstad, J.E. Ein Beitrag zur Theorie der Winderzeugten Meeresstromungen. Gerlunds Beitr. Geophysik 1929,23,237-247 Motjeld, H. 0. and Lavelle, J. M. Setting the length scale in a second-order closure model of the unstratified bottom boundary layer. J. Phys. Oceunogr. 1984, 14, 833-839 Davies, A. M. Comparison of the finite difference and Galerkin methods as applied to the solution of the hydrodynamic equations. Appl. Math. Modelling 1983, I, 226-240 Etling, D., Detering, H. W., and Theunert, F. On the simulation of wind-driven currents in shallow water. Arch. Met. Geouh. Biocl. 1985, 33A. 355-363

45 46

47

48 49

50

Madsen, 0. S. A realistic model of the wind-induced Ekman boundary layer. J. Phys. Oceunogr. 1977, 7, 248-255 Sverdrun. H. U. The waters on the North-Siberian Shelf. the Norwegian North Polar expedition with the MAUD 191881925. Scientific Results 1929, 4(2), 107-113 Heaps, N. S. Spectral method for the numerical solution of the three-dimensional hydrodynamic equations for tides and suges. In Mathematical Modelling of Estuurine Physics eds. J. Sundermann and K.-P. Holz. Springer-Verlag, pp. 75-90 Davies, A. M. On determining the profile of steady wind-induced currents. Appl. Mufh. Modelling 1985, 9, 409-418 Forsythe, G. E., Malcolm, M. A., and Moler, C. B. Computer Methods for Mathematical Compurutions. Prentice-Hall, Englewood Cliffs, N.J., 1977, pp. 86-94 Akima, H. Algorithm 433. Interpolation and smooth curve litting based on local procedures (E2). Comm. ACM 1972, 15, 914-918

1976, 6, 85-92

Appl.

Math.

Modelling,

1989,

Vol. 13, March

191