A long-time-scale, density-stratified shelf circulation model

Continental Shelf Research, Vol. 12, No. 1, pp. 115-141, 1992.

0278-4343/92 $5.00 + 0.00 © 1991 Pergamon Press plc

Printed in Great Britain.

A long-time-scale, density-stratified shelf circulation model RICHARD J. GREATBATCH* a n d ALLAN GOULDING*

(Received 2 October 1990; in revised form 12 March 1991; accepted 2 April 1991) Abstract--We describe a computationally efficient, finite-difference numerical model for use in

understanding circulation and variability in a density-stratified shelf and slope region on seasonal time scales and longer. The model achieves its efficiency by reducing the horizontal momentum equations to a balance between the Coriolis, pressure gradient and friction terms but leaving the temperature and salinity equations in their full generality. The model uses a sigma coordinate in the vertical, with a = 0 at the ocean surface and a = - 1 at the ocean bottom. There is no restriction on the form of the parameterizations that can be used for mixing in the temperature, salinity and momentum equations. The most important limitation of the model is the requirement to use a linear parameterization of bottom stress. This can be in terms of either bottom velocity, vertically averaged velocity or bottom geostrophic velocity. The model also lacks information on coastal, baroclinic Kelvin waves. This restricts its use to shelf regions that are wide compared to the internal radius of deformation. Fortunately, many shelf regions (e.g. the Newfoundland/Labrador shelf, the Middle Atlantic Bight, the West Florida shelf and the South Australia shelf) satisfy this requirement. We illustrate the behaviour of the model by describing an experiment showing the development of a shelf-break jet such as the Labrador Current. The model ocean has a uniform depth of 50 m along the western and northern boundaries and has an adjoining shelf and slope region. Initially, the density is everywhere uniform. Lighter water is then flushed in through the northern boundary and allowed to exit through the southern boundary. The inflow water enters the model domain at its shallowest depth, has uniform vertically averaged velocity and has a density which matches the initial value at the eastern boundary and then decreases linearly westwards. As time progresses the flow is increasingly concentrated at the shelf-break. One year into the integration, a well-defined shelf-break jet has developed. At this time, the shelf-break appears as a convergence zone in the bottom level of the model, with offshore Ekman transport on the shelf being opposed by weak upslope flow offshore. Subsequently, the J E B A R term plays a dramatic role in the splitting of the jet, creating a new jet further offshore. After 2000 days, the model exhibits two jets with roughly half the inflow transport feeding directly into the offshore jet. The development of the split jet is sensitive to the parameterization used for the bottom stress.

1. I N T R O D U C T I O N

IN this p a p e r we describe a computationally efficient, finite-difference numerical model for

use in understanding circulation and variability in a wide, density-stratified shelf region on seasonal time scales or longer. The model is an extension, to allow for a three-dimensional d e n s i t y f i e l d , o f t h a t d e s c r i b e d b y HENDERSHOTr a n d PaZZOLI ( 1 9 7 6 ) a n d is a m o d i f i c a t i o n a p p r o p r i a t e t o a s h e l f r e g i o n o f t h o s e p r o p o s e d b y HASSELMANN ( 1 9 8 2 ) a n d KILLWORTH

*Department of Physics, Memorial University of Newfoundland, St. John's, Newfoundland, Canada A1B 3X7. 115

116

R.J. GREATBATCHand A. GOULDING

(1985). It achieves its efficiency by reducing the horizontal momentum equations to a balance between the Coriolis, pressure gradient and friction terms but leaving the temperature and salinity equations in their full generality. Time-stepping is limited only to the temperature and salinity equations with the velocity field being diagnosed at each time step from the calculated density field. Since coastal trapped waves have been filtered from the governing equations, the only important limitation on the time-step is imposed by the advective time scale. Time-steps of the order of 1 day are then possible, allowing integrations covering many years model time using only modest computer resources. The principal difficulty in applying this technique to a shelf region is in the treatment of the bottom stress. Bottom stress not only depends on the velocity field but is also instrumental in shaping that field. However, the velocity field must be diagnosed at each time-step which means that the parameterization of bottom stress must be amenable to a semiimplicit treatment. This is possible if a linear parameterization is employed. Fortunately, on the time scales for which the model equations are appropriate (about a month or longer) a linear parameterization is probably quite reasonable (CSANADY, 1982). This, however, is the most important limitation of our model. The restriction to a "wide" shelf arises from the fact that the model does not contain information about coastal Kelvin waves. This is discussed in detail at the end of Section 2 (by "wide" we mean that the crossshelf scale must be large compared to the internal radius of deformation). Our original motivation for this work was to find a computationally efficient method of modelling the long-time-scale circulation and variability on the Newfoundland and Labrador shelf and slope (for an example of interannual variability in this region see PETRIEet al., 1988). We took as our starting point the work of HtJKtJDAet al. (1989). These authors used a computationally efficient method to solve for the three-dimensional, quasisteady circulation off Newfoundland. The efficiency of the model, as with the one to be described, came from neglecting the local time derivative and non-linear terms in the momentum equations. The circulation pattern computed using the model is very similar to that obtained by GREENBER~and PETRIE(1988). These authors used a vertically integrated, uniform density model which included the terms neglected by HtJI~UDAet al. (1989) and also allowed for quadratic bottom friction. Hukuda et al.'s model is restricted, however, by the assumption of uniform density. The desirability of including density stratification can be understood from Fig. 1. The density front at the shelf-break is associated with the southward flowing Labrador Current. Clearly, to correctly model the three-dimensional flow field the geostrophic shear associated with this front must be included. This is important if we are to understand the movement of icebergs onto the shelf (e.g. PETRIEand ISENOR, 1985) or the flux of nutrients across the shelf-break. Indeed, as already pointed out by Hukuda et al., a feature of their calculations is that the transport across the shelf-break can be in different directions in different depth ranges. Including density stratification will place constraints on this transport not present in their uniform density model. In developing stratified, shelf circulation models, the usual approach is to apply timestepping to the momentum equations with the local time-derivative and non-linear terms retained (e.g. B LUMBERGand MELLOR,1987). When the interest is in instability processes, as in JAMES (1987), retaining these terms is essential. However, on time-scales long compared to the adjustment time by coastal trapped waves, it should be possible to neglect these terms in many shelf regions. IKEDA (1986, 1987, 1988) has made extensive use of a model (and variants thereon) in which the zonal momentum equation is reduced to a geostrophic balance between the

A long-time-scale, density-stratified shelf circulation model THETA (DEG. C)

SALINITY (PPT)

STATIONS

~¢

117

|t

20

le tn

IT

Is

STATIONS i | 14 t~

tl

*

It

!

l( 41

i.,

-

f

200

~

300

I

o ~, ~,0

5oc

SIGMA- TH (kg/m~)

.

';-i:3

BRUHT - VAISALA FREQ (CYCLES/HR)

STATIONS

.'t~'°~'/,".,-,~r,~.,',-z~'~,~-~" "

"

"/"-'__f~.,.ZJ

i '°°

" > " ~ . . - , 4 , , ' ,...". .

STATIONS

?"

~ " ' - "

"

"

'° " "

"

"""

"'"

:

""" ,

"

...

-,..

Fig. 1. CTD data from a section taken across Hamilton Bank, Labrador, during October 1985 (after WRIGHTet al., 1988). Note the change in vertical scale below 50 m depth.

n o r t h - s o u t h flow and the east-west pressure gradient. This simplification also leads to a computationally efficient model, although not as efficient or flexible as the model we shall describe. The plan of this paper is as follows. In Section 2 we set up the model equations and discuss their validity. Section 3 discusses the formulation of the model in terms of finite differences and in Section 4 we describe a model experiment which illustrates the development of a shelf-break current (such as the Labrador Current). In Section 4 we also discuss in detail the handling of the boundary conditions for both "open" and "closed" boundaries. Section 5 provides a summary and conclusions.

118

R . J . GREATBATCHand A. GOULDING

2. THE MODEL EQUATIONS AND THEIR VALIDITY We use spherical coordinates (2, ~b, a) where 2 is longitude, ~blatitude and a = z/His the vertical coordinate and equals 0 at the ocean surface, z = 0, and - 1 at ocean bottom, z = - H . Here H --- H(2, ¢) is the ocean depth. In this coordinate system, the governing equations are

-fv =

1

Op

acos~bp0O2

fu=

10p apo O~

1

OH + 1 02 -~(vuo)o

acos~b°

l b a O H + 1 (vvo)o a Odp -H2

0 = -po - pobH

(2) (3)

( H u ) + (Hvcose) +

acose

(1)

(Ha)=0

(4)

where =~

w-

acoscp 02 + a ~ / j

(5)

and

atO(bH) + a cosl {0~b (bHu) + O----(bHvc°s~)} + -~a O (bHY~) =

(6)

(1) and (2) are the horizontal momentum equations, (3) is the hydrostatic equation, (4) is the equation of continuity and (6) is the buoyancy equation. For convenience, all the symbols are defined in Appendix 1. In (1) and (2), the balance is between the Coriolis force, the horizontal pressure gradient force and the vertical mixing of momentum (note that the horizontal pressure gradient terms include a part involving gradients of H that arises from the fact that derivatives with respect to 2 and q~ are along constant a and not constant z surfaces). We also have the following boundary conditions on the velocities at the ocean surface (a = 0) and bottom (a = - 1) 1

rz

~vUo-po,S" vuo

_

~

Z'b.

1

r~

HVV°=--po 1

at

a=0

(7)

at

a=-I

(8)

r~

Po' ~vvo=--po

and fl=0

at

a=0

and

a=-l.

(9)

Here (r~, r~) is the surface wind stress and (r~, r~) is the bottom stress. Putting ~ = 0 at a = 0 corresponds to the rigid-lid approximation and is valid provided horizontal length scales are small compared to the barotropic Rossby radius of deformation. On the other hand, putting fl = 0 at a = - 1 corresponds to the kinematic condition w=

aH

vaH

a cos ~p a2

u

a d¢~

at

z = -H.

119

A long-time-scale, density-stratified shelf circulation model

It should be noted that the buoyancy equation (6) can be replaced by separate equations for temperature, T, and salinity, S. The density, p (and hence the buoyancy b), is then calculated from T and S using the equation of state for sea water (e.g. BRYANAND COX, 1972; FlUEDRICHand LEvrrus, 1972). As outlined in the introduction, the method of solution involves time-stepping the buoyancy equation and then diagnosing the velocities u, v, f/from the new buoyancy field, b. This process involves two steps: (i) calculating the vertically integrated flow field; and (ii) calculating the vertical structure of u and v and the corresponding f~ field. In order to carry out (i) we first vertically integrate (4) and use (9) to give +

0-~

8

0

(10)

vdo.

(11)

(~H cos q~)

where

~=~_ uda

and

~=~_

1

1

(10) enables us to define a volume transport streamfunction W with

affH =

-~0~ and

av-H cos ¢ = ~px.

(12)

The equation satisfied by ~ is obtained by first vertically integrating (1) and (2) to obtain _f~

=

f~=

1

{of °

a COS ~ P 0 -

~-~ --1

1 {0__~I0

- -

apo

p da

OHm_ pda+po-~~ 1 bada} + p ~1( % - z +

-1

Po0__~__ ba da} 1

+

1

(r~

-

r~)

~)

(13)

(14)

~ 0H

where use has been made of (7) and (8).

now gives 0 f 0 [ ~ (~,t ~) - ~-~(~0~ f ) ] = ~oo a 1_~ [ 0 ~--------H~ (rs~- r~l1 H

JEBAR is the Joint Effect of Baroclinicity And Relief (SARKISVANand IVANOV,1971) and is given by JEBAR

=

-

O {O_~I~l ba da} + ~-~ O {OH -~- ~1 _ ba da} •

~-~

(16)

In order to solve (15) for ~0, all the terms on the right hand side must involve either ~0itself or known quantities which are either given at each time step [e.g. the surface wind stress field _r s = ( r ~s, rs~)] or can be calculated from the buoyancy, b. JEBAR clearly falls into the

120

R.J. GREATBATCHand A. GOULDING

latter category and therefore does not present a problem. In general, however, it is not possible to express the bottom stress in such a form. In this paper, we shall consider three parameterizations of _.rbfor which this is possible; namely a linear parameterization in terms of either bottom velocity (Ub, Vb) (this case is discussed in Appendix 2), the vertically averaged velocity (~, ~) or bottom geostrophic velocity (Ugh, "Ugb).The simplest case is that of vertically averaged velocity. Putting

r__= por(u, v)

(17)

(15) becomes

coZs;], + rrw',>c°s+

_a[O(7) P0 g

° - ~

+ JEBAR.

(18)

This equation is now easily solved given appropriate boundary conditions for q~ (these are discussed in more detail later). We next consider the case of bottom geostrophic velocity for which

[b = pOr(Ugb, "Ogb)'

(19)

[Note that for the method to be described to work, _rb can be any linear function of (/Agb, Dgb). In particular, the direction of r__bcan be rotated through some angle from that of (Ugb, Vgb).] In order to express (Ugb, l)gb) in terms of ~pand the density field, we first rewrite (1) and (2) in the form

1 t

- f v = -fVg + - ~ (vuo) o

(20)

fu lug + ~ (vvo)o where (Ug, Vg) is the geostrophic velocity given by 1{ 1 0 p + l b c r O H l Ug = - 7

apo-~

a

-~J

and

(21)

1{

1

Op

boOHl

1

Vg = ~ [a cos qOpo02 t- a cos ~p

02 J

Vertically integrating (20) and using (19) gives

- f g = -fgg -t- Tzs poll f~ = f~

+ ~__L _ rv~b

poll where

rugb ] H H "

(22)

A long-time-scale,density-stratifiedshelfcirculationmodel

ifo

Ug= ~ _HUgdZ and

Vg

=~

121

H vg dz.

We now put and /)g = ~)g -- ~-g

ag = Ug -- Ug

(23)

so that and "Ogb = ~g + 0gb (24) where the subscript b denotes that the quantity is being evaluated at the ocean bottom. Substituting for ~g and ~g from (22), (24) can be rearranged to give Ugb = ~g -t- /lgb

--

Ugb Ugb

_

1 (1+~

tigb+~

1

(1 nt- " ~

{[

_

r s~ ]+/?/)gb+ pofHJ

r~ ]_ [

/)gb "1- ~ -t- ~ ]

--.q_ uS U

E //gb + ~

(25)

_ r~ ]l pofH]j

where e = rlfH. In (25), agb and Ogb can be calculated from the buoyancy using (3), (21) and (23); u, v are expressed in terms of ~0using (12) and (z~, v~) is the known surface wind stress field. Substituting for (Ugh, Vgb)into (15) using (19) now gives the following equation for ~f which can be solved at each time step { ( r~2 ~)+ 0 {r~q~cosdpl 0 (_ re lp¢) 0 H2( 1 + e2) cos 0-~ ~ / H 2 ( 1 + e2)) + ~-~ + e2) H 2

re

0

f

0 ( ( 1 e2) H~---@}+ {0-~(~Paf ) - ~-~(IP~ ~)} 04~ + a

He_

a

_0 ircos" u/l

+ e2)H] + e2)H] a{O ( - a

0-2 (1 + {0-~( (1 +

rrs~

O~ \(1 +

Otp \ (1 + ]

O(

eZ)pofH2] + -~

0

e2)H]J

ea)H]J

r cosq~s ]~ (1 + e2)oofH2]J

re~ ] 0 ( recosq~s ]~. eZ)pofH2] - 0-~ (1 + e2)pofH2]J

(26)

Both (18) and (26) are elliptic equations for % They are solved using boundary conditions of two forms: (i) on inflow boundaries or where f/H contours enter the model domain (in the sense of long topographic wave propagation) v2 is specified; on other boundaries Ovd/On= 0 where O/Ondenotes derivative normal to the boundary. The latter is a form of open boundary condition and forces the flow to be normal to the boundary. Other forms of open boundary condition are also possible as long as (26) can still be solved. The condition that 7) be specified where f/H contours enter the model domain is important. Experience

122

R.J. GREATBATCHand A. GOULDING

has shown that when f/H contours connect two open boundaries, spurious flows can be generated. Once ~0 has been calculated using either (18) or (26), u, v can be calculated using (12). To obtain the full three-dimensional (u, v) field, we first subtract (22) from (20) to give 1 x + V~b l --fY = --fl)g + ~

(Yaa)cr -- T~s

t"

(27)

poll - - ~ where in (22), (rUgb, rVgb) has been replaced by Vb/Po in order to accommodate both parameterizations (17) and (19). (27) is now easily solved for (a, 7)) to which (if, Y) is then added to given (u, v). Finally, 1~ is obtained from equation (4). We next discuss the validity of the model equations. We noted in the introduction that neglecting the local time derivative terms from the momentum equations filters out coastal trapped waves. It follows, therefore, that our model equations should be valid on time scales long compared to the adjustment time by these waves, assuming that the Rossby number is sufficiently small that the non-linear terms can also be neglected. This time scale is typically of the order of a few weeks. We have also neglected the lateral mixing terms from the momentum equations. This is not essential for our model to work, although retaining them would greatly complicate the solution procedure. In particular, the equation for the streamfunction [replacing either (18), (26) or (50)--see Appendix 2] would be fourth order and in addition, an elliptic operator would have to be inverted on each o-level and at each time step to obtain the three-dimensional flow field. The model is also strictly valid only in "wide" shelf/slope regions, where by "wide" we mean S=

<< 1

(28)

where L is a length scale characterizing the width of the shelf and R is the internal Rossby radius of deformation. This restriction arises from the fact that coastal Kelvin waves and, in particular, any knowledge of coastal Kelvin wave propagation, is completely missing from our model equations. Information propagates through the model domain along f/H contours (in the sense of long, topographic wave propagation). This means that the model does have a "memory" of topographic wave propagation. A similar "memory" of coastal Kelvin waves is absent. This means that our model equations are not strictly valid within an internal deformation radius of the coast. Fortunately, as noted by MIDDLETONand WRIGHT (1990) for many shelf regions the inequality in (28) is satisfied. They give a value of S = 0.02 for the Labrador shelf, 0.01 for the West Florida shelf, South Australian shelf and Middle Atlantic Bight (CLARKE and BRINK, 1985) and 0.04 for the Weddell Sea shelf (MIDDLETON et al., 1987). Our model should, therefore, prove useful for modelling longtime (seasonal and longer) circulation and variability in these shelf regions. 3. T H E M O D E L F O R M U L A T I O N

The grid arrangement of the dependent variables u, v, 12 and b in both the horizontal and vertical is shown in Fig. 2. In the horizontal, the Arakawa B-grid is employed. Since both u and v are stored at the same grid points, this greatly facilitates the calculation of ti and 73using (27). Streamfunction ~ is stored at the same points in the horizontal as b, as is

A long-time-scale,density-stratifiedshelfcirculation model

Z~

X

X

x

x

123

• ~,b, 9 x

u,v

AX

(o) Horizontal grid arrangement • x

o'=0 IlZ~)l

x

l (Z~)Z

* 'f/' x u,v, b

X

[b) Vertical grid arrangement Fig. 2. The grid arrangement of the dependent variablesu, v, fl and b in both the horizontal and vertical. the surface wind stress _rs and the ocean depth H. In order to obtain values of ~ and H at (u, v) points [necessary for the calculation of t~, ~3 from (27)], four point averaging is employed. We now describe the finite differencing, beginning with the buoyancy equation (6). The advective flux terms on the left hand side are each separately treated, following JAMES (1986), as a weighted average of"forward in time, upstream differencing (FU)" and "leapfrog in time, centred in space differencing (LC)". It should be noted that the form of (6) ensures that the resulting finite-difference operator for these terms is conservative. We illustrate the method for the u-advection term (O/al)(bHu). This is approximated at grid point i, j, k (where i numbers the grid points in the direction of t increasing, j is the direction of ~ increasing and k is the direction of tr increasing) by

(~2 (bHu))

= (CU)uk(FU} + (1- (CU)ijk){LC} • i,j,k

(CU) is chosen at each time step and at each grid point so that F U differencing is important only in regions where gradients of b are changing rapidly. This has the effect of removing the computational dispersion associated with LC differencing while at the same time avoiding the strong numerical dissipation associated with FU differencing. Specifically, (CU)i,j,k is defined, following JAMES (1986), by bl (CU),j,k = I IA,+x/2,j, [Ai+l/2,/,kb I + iA,_t/2j,kbl I

P

(29)

124

for

R . J . GREATBATCHand A. GOULDING

[Ai+UEj,kb[ + [Ai-1/Ej,kb[ > e and (CU)i,j,k

=

0 for

[Ai+l/Ej,kb[ + Ai_l/E,Lkb[ <--e where

Ai+ l/2,Lkb = bi+ l,Lk -- bi,Lg

p = 1 and e = 10 -4. A similar procedure is adopted for each of the (O/Ocp)(bHv cos ~) and (O/Oa)(bHl~) terms but with the testing by equation (29) this time being for changes of gradient in each of the q~and a directions, respectively. To control time splitting associated with the leap frog scheme, an ASSELIN (1972) time filter is applied every five time-steps. The filter replaces the value b n of b at time level n by the value x , ' ~In - 1 b nI = b" + g'g~o 2b n + b n + l ) after the new value b "+a has been calculated. A value of v = 0.5 is used. Next we consider the "MIXING" terms in equation (6). This is a combination of vertical mixing and lateral mixing with the latter being along constant tx surfaces, following MELLOR and BLUMBERG(1985). These authors argue that in the presence of topography, this is a physically more realistic method of parameterizing lateral mixing than the usual horizontal mixing often employed in numerical models (e.g. BRYAN, 1969). The finite difference form used for "MIXING" in equation (6) is, therefore, an approximation to --

l

la_._(AI4H(gblO

1 0 /

Ob\

+ ~o(KH~o)"

(30)

The vertical mixing term is treated implicitly in time in order to avoid numerical instability in shallow water, otherwise forward-in-time differencing is employed. The spatial derivatives are approximated using centred-differencing. AH is the lateral mixing coefficient and KH is the vertical mixing coefficient. Next, we consider equations (27) and (4) used to diagnose the three-dimensional velocity field. The usual centred space differencing associated with the B-grid is used throughout (see MESINeER and ARAKAWA, 1976). Care, however, must be taken in the treatment of the fig, 0g terms defined using (21) and (23). This is because the horizontal pressure gradient terms in (21) have two parts: one involving derivatives of p along constant o surfaces and the other, derivatives of H. Error can result when these parts have opposite signs but comparable magnitudes (see, for example, BAa'rEEN, 1988). For example, if the buoyancy, b, has no variation in constant depth (as distinct from or) surfaces, (ug, Vg) defined by (21) should be depth independent. This means that the ((1/a cos dp)bH~, (1/a)bH¢) term in (21) must be exactly cancelled by the depth varying part of ((1/a cos ~po)Op/O)~, 1/apo, Op/Odp).Unfortunately, when these terms are written in finite-difference form, this is not, in general, guaranteed. In order to reduce this error, we reformulate these terms following BLUMBER~and MELLOR(1987). We therefore define two new functions A(2, q~, o) and B(2, q~, (7) by A= and

do

1

(31)

A long-time-scale,density-stratifiedshelfcirculation model

125

Using integration by parts, it is easy to see that B = -ab - A.

(32)

P = Ps + PoH A

(33)

Also, integrating (3), we see that

where Ps is the surface pressure (at tr = 0). Using (32) and (33), we now see that Op Op OH Ops OJ,z = ~-o + O o b a - - ~ = p o ( H A ~ - H ~ B } + ~-~

and

ap 0p + ae~z a~o

I poba #tl = Ro{HA o - HoB} +

a~

(34)

09

where z is vertical distance measured upwards (i.e. z = trH). It turns out that when the horizontal pressure gradients are approximated by finite differencing the expressions involving A and B on the right hand side of (34), then they are exactly zero if b depends only on z and varies linearly with z (the finite differencing being the usual centred differencing associated with the B-grid--see MESlNGERand ARAKAWA,1976). In this way, the error associated with these terms when there is steep topography is reduced. Finally, we discuss the streamfunction equation [either (18), (26) or (50)--see Appendix 2]. Once again the usual centred space differencing associated with the B-grid is employed except that a five-point Laplacian is used for the first two terms on the left-hand side of each equation. Deriving these equations in finite differences from the finite difference form of (13) and (14) will lead to a nine-point Laplacian. However, experience has shown that nine-point Laplacians are somewhat less stable in their behaviour than five-point ones, with a tendency for spurious computational modes to be excited. For this reason, we prefer the five-point Laplacian. 4. A N E X P E R I M E N T I L L U S T R A T I N G T H E D E V E L O P M E N T OF A SHELF-BREAK JET

In this section we describe a model experiment which illustrates the behaviour of the model and also the form of the boundary conditions which can be applied when solving the buoyancy equation (6). The experiment illustrates the development of a shelf-break current such as the Labrador current (see Fig. 1). In order to verify the model, we have also reproduced the solutions obtained by SHAWand CSANADY(1983) and shown in their Figs 4, 5, 8 and 9. The model domain used for the experiment is shown in Fig. 3. It consists of a rectangular domain in ;t - q~space bounded by latitudes 40°N and 60°N and of width 10° longitude. The bottom topography slopes down from a uniform depth of 50 m around the northern and western boundaries to a maximum depth of over 1000 m in the southeast corner. There are three distinct regions within the model domain: the gently sloping shelf in the north and west, a slope region in which the depth changes much more rapidly and a deep water region in the southeast. There is no wind forcing or surface buoyancy flux in this experiment. Rather, water is flushed in through the northern boundary and allowed to exit through the southern boundary. No transport is allowed through the western and eastern boundaries. Initially the density of the water is everywhere uniform. Lighter water is then fed in

126

R.J. GREATBATCHand A. GOULDING

0

T o p o g r a p h y5 [

50N~

I

I I1|111

45N~

I

I IIIIII

40 N

0

~

5E

10 [

4,50N

I

445N

~

,

4

0

IOE

N

Interval 100 m Fig. 3. The model domain. The bottom topography is contoured with interval 100 m with the shallowest water along the western and northern boundaries.

through the northern boundary with question to be addressed is whether concentrated in a narrow current at divided into two parts, the first dealing the second with the model results.

the lightest water being on the western side. The or not the southward flow becomes progressively the shelf break. The remainder of this section is with the numerical treatment of the boundaries and

(a) The numerical treatment of the boundaries As previously noted, the horizontal grid arrangement is that of the A r a k a w a B-grid shown in Fig. 2. For the model domain shown in Fig. 3, we should like the eastern and western boundaries to be closed with the normal velocity, u, set to zero there. Also, the condition of zero transport through these boundaries requires that ~p be constant along them. It turns out that a convenient way to implement these conditions is for the boundary itself to run through ~0 grid points with u, v grid points on either side. T h e requirement that u = 0 on the boundary is then achieved by putting u at the grid point immediately outside the boundary equal to minus its value at the first grid point inside and at the same latitude. This is illustrated in Fig. 4a for the eastern boundary. This choice for u implies, from

A long-time-scale,density-stratifiedshelfcirculationmodel

~+~UhV XI

x ~_~

•

•

Uz,V2 x

• t~tVt

127

grid points

XU

"U2,V2 x

(o) boundary

Ul~Vl

Uh° VI

X

X

u~,v~

i u2,-v2

X

X

(hi bo~dory Fig. 4. (a) The implementation of the condition u = 0 on the eastern boundary. (b) The modification used to eliminate the Killworth boundary instability.

equation (12), that the streamfunction ~ must have the values shown at the first set of ~-grid points beyond the boundary. It then follows from (12) that v immediately outside the boundary must be set equal to v at the first grid point inside, as shown. It should be noted that in running the model, it is only the values of u and v (not ~) that need to be set outside the boundary. These values are used in the calculation of 12 on the boundary from the finite difference form of (4). KILLWORXn(1985) has shown that the condition of zero normal velocity at a boundary is associated in his model with a computational wave which propagates along the boundary. The phase speed of this wave depends on the density stratification and inversely on the grid spacing normal to the boundary. If the CFL criterion associated with this wave is violated, an instability develops near the boundary and can destroy the entire calculation. As noted in the Introduction, our model is very similar to Killworth's. Indeed, both models combine diagnostic momentum equations with time-stepping of the density equation. It is not surprising, therefore, that we have also encountered this instability. It is especially a problem in high resolution calculations because of the inverse dependence of the phase speed on grid spacing. To overcome this we have used the alternative boundary formulation shown in Fig. 4b. In this formulation, the values of u on either side of the boundary are set equal, with v outside the boundary being set equal to minus v inside. Although it is no longer guaranteed that the velocity normal to the boundary will be zero at all depths, the condition that ~pis constant along the

128

R.J. GREATBATCHand A. GOULDING

boundary still ensures that there is no vertically integrated transport through the boundary. This form of boundary condition has been found to be very useful. The removal of the requirement that u = 0 at the boundary filters out the computational mode responsible for the instability, allowing high grid resolutions without the need to use a prohibitively small time step. One point that should be noted is that with u and v chosen as in Fig. 4b, calculated along the boundary from the finite difference form of (4) is identically zero at all depths. In a flat-bottomed model, this means that the vertical velocity w is zero at all depths. This should be borne in mind when using the boundary formulation in Fig. 4b since -- 0 along the boundary may not be appropriate. KILLWORTH(1983) has discussed the question as to what are appropriate boundary conditions to use with the thermocline equations. This difficult problem will not be pursued here. In the experiment to be described, the formulation in Fig. 4b is used on both the eastern and western boundaries. As it happens, this particular experiment also runs with the boundary formulation shown in Fig. 4a, the results using either boundary formulation being essentially the same. A condition of no normal flux of buoyancy is also applied. Next we turn to the northern and southern boundaries. Along the former, ~ is specified according to = -Wo2/10 °

(35)

and therefore varies linearly from a value of 0 on the western boundary to -~p0 on the eastern boundary. In the experiment to be described ~0 has value 2 Sv. Since the depth of the model ocean is a uniform 50 m along the northern boundary, this corresponds to specifying the vertically integrated velocity normal to the boundary to be a uniform 0.07 m s -1 into the model domain. At the southern boundary 0__~_~= 0.

(36)

This is an open boundary condition which forces the vertically integrated flow to leave or enter the domain normal to the boundary. Next we turn to the boundary conditions applied to the buoyancy equation (6) along these boundaries. The grid arrangement at the two boundaries is shown in Fig. 5, with the northern boundary located at j = ny and the southern boundary at j = 0. At each level in the vertical, we first test for inflow or outflow at the boundary using the most recently updated velocities. In the case of the northern boundary, this means that when considering how to update b at grid points (i, ny) on a given level we form the average (1)i_l/2.ny_l/2 -t- Vi+l/2,ny_l/2)/2. If this quantity is positive or zero we apply an "outflow" boundary condition, if it is negative an "inflow" boundary condition. For each of these, we solve, respectively, the following equations along the boundary: outflow

~t (bH) + Hv Ob _ 1 0 (AHH~_x?_. Obl 2~._gS _~g_2~ 1 O / _~] Ob\ a O, a2cosdpO2\,,u~wo~,]+,,o~,\Knoo/

(37)

inflow

O (bH) 0t

1

0 (AHH ObI

1 0 /

Ob\

= a 2 co-~~ 02 \cos ~ 021 + H-~a ~KH ~-~) - aH(b - bsPEC).

(38)

A long-time-scale, density-stratified shelf circulation model j=ny

--.

....

j =ny-~

• .... x i-V2

j=ny-I

• i-I

(o)

northern

j

=~

x i+~

(b)

Fig. 5.

i-I

• i+l

boundary

•

--o

•---

• i

•

i-v2

j--O

....

southern

129

•

b, lk

X

UtV

•

i+V~ • .... i

o--i+l

boundary

The grid arrangement at the northern and southern boundaries.

In both cases, lateral mixing along the boundary and vertical mixing are retained [compare (37) and (38) with (30)] but we do not attempt mixing across the boundary. In the outflow case, buoyancy is advected out of the model domain by the (v/a)(ab/a~) term. This term is finite-differenced so that the combination (ab/at) + (v/a)(ab/O~) uses forward in time, upstream differencing (this was found to be desirable for reasons of computational stability). In the inflow case, a term - a H ( b - bsPEC) is included. This term relaxes b to a specified buoyancy field bsPEc on an e-folding time scale of 1/a. For the experiment to be described, 1/a = 10 days and bsPEC is given on the northern boundary by bSPEC(~,, O) = bw + b~_~ 10°

(39)

where bo = 2g/po, bw = -g/Po and ,~ is measured in degrees from the western boundary where b has value bw. bsPEC is essentially the (specified) buoyancy associated with the water flowing in through the northern boundary. For the experiment to be described a = 0 on the southern boundary, implying that the density of the inflowing water is not specified but is allowed to evolve and may depend on the past history of the model experiment (for example if previously there had been outflow there). These conditions are similar to those used by BLUMBERGand MELLOR(1987). They differ primarily in the inclusion of the lateral and vertical mixing terms on the boundary. We found, using text experiments with a more complex density structure specified for the inflowing water, that if these mixing terms were not included, spurious density gradients could develop along the topographic slope near the inflow boundary. These in turn lead to spurious values for the J E B A R term in (15) with consequences for the flow field downstream and in the interior of the model domain. We concluded, therefore, that it is important to include mixing terms along the boundary in order to ensure that the inflowing water is compatible with that already in the interior of the model domain.

130 (b)

R.J. GREATBATCHand A. GOULDING

Model results

The parameter values used for the model experiment are as follows. For simplicity, the mixing coefficients A n and KH in equation (30), and also the vertical eddy viscosity, v, in equations (1) and (2), are taken to be uniform constants with values 100m2s -1, 5 × 10 -4 m 2 s -1, and 2 x 10 -3 m 2 s -1, respectively. The values chosen are probably on the small side for a continental shelf region but are adequate for demonstrating the behaviour of the model. Bottom stress is parameterized in terms of bottom geostrophic velocity using equation (19) with r = 0.001 m s -a. The horizontal grid resolution is ~° × ~° and there are 20 equally spaced levels in the vertical extending from a = 0 at the surface to o -- - 1 at the bottom. The total grid size is, therefore, 30 × 60 × 20. The computational efficiency of the model is demonstrated by the fact that all the model results to be presented were obtained on a V A X 8800 computer. Initially, the perturbation density, p', has the value + 1 kg m -3 everywhere in the model domain. As can be seen from equation (39), the perturbation density of the inflowing water is independent of depth and varies linearly from this value on the eastern boundary to - 1 kg m -3 on the western boundary. Furthermore, the depth of the model ocean is uniform along the inflow boundary and the vertically integrated velocity associated with this inflow also has a uniform value (0.07 m s-l). Given the "unbiased" nature of this inflow, the question arises as to whether or not the presence of the shelf-break will lead to the development of a shelf-break jet. Figure 6 shows the volume transport streamfunction 72, the perturbation density at the top level of the model and the horizontal velocity vectors at the top and bottom levels 53 days into the model integration. At this time, the light, inflow water has spread only a short distance from the northern boundary. As a consequence, the surface velocity vectors (Fig. 6) show no concentration of the current at the shelf-break but, instead, rather a uniform flow over the shelf. The flow in the bottom level of the model looks quite different and is dominated by the Ekman transport associated with the bottom friction. Figure 7 shows the same output fields 388 days into the integration. By this time, the light, inflow water has penetrated to the southern boundary. However, it has also been carried towards the shelf-break by the offshore Ekman transport. As a consequence, the largest surface velocities (of the order of 0.2 m s -1) are now found at the shelf-break. WRIGHT (1989) has shown how the frictional stress associated with bottom geostrophic velocities can result in the cross-shelf migration of an idealized density front. Our results suggest that this mechanism can also be important in the frontogenesis process. Interestingly, in the bottom level of the model, the shelf-break now appears as a convergence zone, with flow up the slope meeting the offshore Ekman transport on the shelf itself. Comparison of the streamfunction fields in Figs 6 and 7 shows that at the later time, the streamlines converge into the jet along the shelf-break. The equation governing the streamfunction is (15), which, since there is no wind forcing, can be written as

72~

- - ~ 72~

po(O2\H/

0(~\

H

])

In the northern part of the basin, both the bottom friction torque and J E B A R terms play a role in concentrating the streamlines into the jet (see Fig. 8). However, in the southern part, these terms cancel each other out. This means that the J E B A R term suppresses the tendency, due to bottom friction, for the jet to "diffuse" in the direction of long, topographic wave propagation (this effect has been described by CSANADY, 1978). Figure 7

A long-time-scale, density-stratified shelf circulation model

Streamfunction

I~

Surface

I'[t~'''/~;'~;'~'~;'~-' / i/////'." !! Il J///;//I/

131

Density p

#/ I/7//?"

II I | f lI ' # I / II I I I

l!!lli!f/

il!ll/llll t I i. IIAIII

..-

ii Ill | I , IIAIII I I ! illlll I

I #

iill II ~{itl{~. ! IIilll

I

I I I |

I I III t ! II!l!lll ~I i I~I,i,I I I I I I IllII, I I III I i ! I i!l!il i I Ii "1,1'1' I i II IIIIII Ill ...!..,.s.ttq~t

I I ! I I ......

~ ......

Interval 0.20 Sv Time 53 days

Interval 0.100 Time 53 days

Surface Velocity

Bottom Velocity

g~z.z, , "

I °

. .,. . . ., . . .,. . . .° . . . .. . . . . . . . . ..

.

.

-= 0.10 m / s Time 53 days

.

.

.

.

I! .

.'.2.'2.2.

~. . . . . . . . . . . . . . . .

== 0.10 m / s Time 53 days - -

Fig. 6. The volume transport streamfunction, perturbation density in the uppermost level of the model, and the horizontal velocity field in the top and bottom levels 53 days into the integration. The contour interval for perturbation density is 0.1 kg m -3.

shows some transport into the model domain through the southern boundary. This transport feeds into the shelf-break jet and is driven across the contours by both the JEBAR and bottom friction terms. Another interesting feature of the streamfunction field is the "tongue" which is spreading along the shelf-break towards the eastern boundary.

f/H

132

R . J . GREATBATCH a n d A . GOULDING

Streamfunction 1~

Surface Density p

7777-

IIll/~

fl!/// .//. "

.,,,, IHJ//,

\~3,1, l ,"~' lllq :U:H ",

..I-

.I'"

I

I III I

1111

I Ill I I II!

I IIiiii "\ ,!Ii~

" I",',IWI

,' ii~ ~ ! I '""~;

c.-/,,,,,,m,, L.// / // / //fllllllil

I I! I,~

I I!III~I . .It. j lhl~

l l .t. ~ . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Interval 0.20 Sv Time 388 days

Interval 0.100 Time 388 days

Surface Velocity

Bottom Velocity

..4..
•

~

,

7.-.222.".2 ................

--

Time

=

0.10

m/s

388 days Fig. 7.

--*

=

Time

0.10

m/s

388 days

As Fig. 6 but for 388 days.

This can be seen in the northern part of the basin. The offshore transport associated with this "tongue" is driven by the JEBAR term arising from the along-isobath density gradients in the region. By 840 days (Fig. 9), it has extended all the way to the eastern boundary. This feature can also be seen in the density fields shown in Figs 7 and 9. Its propagation, which is in the opposite direction to that of long, topographic wave

A long-time-scale, density-stratified shelf circulation model

Friction

133

Jebor

i...........................

II I/ /I I/ ..................

Intervol 0.500E-02 Time 388 doys Fig. 8.

i

. . . . . . . . . . . . . .

Interval 0.500E-02 Time 388 doys

The bottom friction torque and J E B A R terms in equation (40) 388 days into the integration. The contour interval is 5 x 10 -3 m 2 s-1.

propagation, is due to advection of the density field by the velocities in the near bottom region (see Fig. 9). Comparison of the streamfunction fields shown in Figs 7 and 9 shows that at 840 days, the inflow through the southern boundary now extends further north. The total transport associated with this inflow is roughly 0.75 Sv (compared with the 2 Sv inflow through the northern boundary) and has increased from roughly 0.55 Sv at 388 days. At 840 days we can also see that part of the shelf-break jet is being forced offshore just north of the southern boundary. This forcing comes from the J E B A R term and is associated with the tendency for the light, inflow water to extend a little further offshore near the southern boundary than further north. This effect propagates northward, again in the opposite direction to long, topographic wave propagation, and ultimately has a dramatic effect on the model solution. By 1357 days (Fig. 10), it is already leading to the establishment of a "new" jet further offshore and in deeper water than the orginal one. The offshore transport we can see in Fig. 10 connecting the "old" and "new" jets is forced entirely by the J E B A R term. By 2001 days (Fig. 11), we clearly have two distinct jets with half the inflow transport through the northern boundary now circulating around the "tongue" extending to the eastern boundary and feeding into the offshore jet. Indeed, the presence of the original jet is more clearly seen in the surface velocity field (Fig. 11) than in the streamfunction. Interestingly, beneath the offshore jet we now have reversed flow (see Fig. 11). In addition, the bottom flow also has an upslope component associated with Ekman transport in the bottom level. At later times (not shown), a slow evolution

134

R.J.

GREATBATCI-I a n d A . GOULDING

Streamfunction

~

Surface

Density p

~,I7# '/-L;:>~> I16,//.-"~>'~GTill l:lil/:"

i~"l/~IW"'~"mJz :: flllll/;'/llll//L___<---~

"l' '""',,'""'\\';i ':',',,,i,, "'',,, ( ~' ",,'lib/ '' ' , ,,[-,;""r -, . . .-."

tl'""'g/~ ' "'""~ ",,,,," '"" ,' "'",,,,~, "~'~" ,,

l l l . 1 / . / / . /

~

| \ ~\llhl P\ f

~ \fill!It" \

l

I

t

I

•

t

~, i I illl Iiiii IIII! t" I! ,,IIIJll]l . . . . r' \ 1,1,111! lllilll I!l,lll!l

I i i

i

I

f

#liitii, ;tl ' "ll'l",l".~ i ....

'

'

t

I~ \

,tI ,,11)11 I II

I

,,.'

\ I I I I I

f t

ill

iI~ ,I II f ........I I,.:,~,.i,~,,,,

Interval

I

I i in

illl | I /lllllllllllllll .......

i<.:./"~Jllill ..............

0.20 Sv

Interval

84.0 days

Time

Surface Velocity

~<~

....

.

!

'~I ,l]' ' I I I! | ''1 ,Ij 'l/ 'i/ 'Jl . .' . .J. L

ii[

-= Time

J

I . I . I '

.

~

. .

.

840

0.100 days

Bottom Velocity

I / .....

" I f i

il

'

II~ll I I I I

Time

:'~

. .

.

. .

.

. .

.

.

.

.

.

"

........

~_~

;~',~;

. .

.

.

.

.

"-'"

|,

....

1 . . . I . . . I' I' I' .I . .' . . . . . . . . . . .

. .

. .

. .

. .

•

'

"'' " ' '

~''

. . " . . . . . . . . . . . . ........ .'. :. ?. :. . . ' ..' . L '. . . . . . . . . .

l

0.10 m/s 84.0 days Fig. 9.

-: Time

0.10 m/s 84.0 days

A s Fig. 6 b u t f o r 840 d a y s .

continues in which the light water at the surface gradually spreads further offshore, but the basic two jet structure remains. The dramatic splitting of the jet described above started at the southern boundary and extended northwards. It seems likely that the initial perturbation to the density field responsible for this effect is associated with the open boundary condition applied at the

A long-time-scale, density-stratified shelf circulation model

135

Streamfunction

/-/;,-;> I I t.'l I//I.1 I I ! / ////.."

! I l ! / /I; /"g/

r, " I

~tlll/t/ "-:----. ll~\ll ~t

IIII I , ~\111,q ; \ ~ ~.t11[q

\

\

\l'I',! t

I lhllI~-g:~A • I p I 11 X I I I iI!I~.,~ 1 Ill

\

-

\

I I I 1 i t

' i,!P ~u i i I I01!A ~1111 11 II 11!i " . I 11!! I t I ~'I I

I III ~ I I Illll iI )111 I I II VI,,~ 4 I ~,11Illl I [ Ill,,/ll~ lllll I I I ''I I I ' I II I'1111] I, .... i..L ~ hlJ IIIL..t l . I . .

Interval 0,20 Sv Time 1357 days Fig. 10.

The volume transport streamfunction at 1357 days.

southern boundary. However, its subsequent propagation northwards (in the opposite direction to that of long, topographic wave propagation) suggests that it is a "real" feature of our model. Indeed, it seems that any suitable perturbation to the density field could generate the effect. It should be noted, however, that the splitting of the jet is sensitive to the parameterization used for the bottom stress. In an experiment we have done in which bottom stress is parameterized in terms of bottom velocity, the splitting of the jet does not OCCUr.

We shall not describe this experiment further here. In a later paper we shall describe in detail a set of experiments designed to illustrate the development of a shelf-break jet.

5. S U M M A R Y A N D C O N C L U S I O N S

The work we have described was motivated originally by the need to find a computationally efficient method of modelling the circulation and the associated long-time-scale (seasonal and longer) variability on the Newfoundland and Labrador shelf. This shelf contains some of the richest fishing grounds in the world and the seasonal and interannual variability of the physical characteristics of the water on the shelf are believed to influence the seasonal and interannual variability of fish stocks. The model achieves its efficiency by integrating a reduced set of equations in which the horizontal momentum equations are reduced to a balance between the Coriolis, pressure gradient and friction terms. This filters out coastal trapped waves and allows a long time-step (of the order of i day) to be

136

R . J . GREATBATCH and A. GOULDING

Streamfunction

1~

[ 7f>~7>'~"7

"I

l

.

Surface .

.

.

.

.

i I #/1.#///~

i

: '

: "'" ',"," " I ~, \ ~,~ ^ll!!l!!l"\ \ l f

\

\

i

:

x l i I ILI!I!!I x I ,,ti ! ,i I,llllii|!llliii# I'"

'

::

, iii ,l~(tlli r • , :

,i iii Iii 'llliX I I ,, II I!!ll.i~tll!l t I I : " li'll/~ I'1 ' ' ' : | ~,lliJl/ilh Ill : i , I • III =,=.~I . . . . iI IiIi IllIlill/I II1~ ill ,, II Ixlil ,' I I I ,jill , I !l,.-l~ Jill ', I i t w i i l ~ i . III I ', t l I~\1111.,... : f ....... I.,..',,..,'J'~lu.,.,.,

.

.

.

III IIIII I I

Ii I il'i+l'i 1' il,ixxx\\kiliill/ f ! I,iXx\\\%,illlllll / t ,iI xx?~x,i/lllllll II tV..'~,ii

t t

I

I ,,, \tlllllllllll f ~ ~iiiiillllllllll

1 1

|t ,i illllllllllllll I lllllillllllllll

t1

,illllll

I i I i

i

IIIIIIIIIIIIIIIII I IIIIIIIIIIlllll I IIIIIIIIIIIIIIII I Illllllllllllllll

1 1 1 i

IIIIIIlllllllllll Illilllllllllllll i lillllllllllllll I IIIIIllllllllll I I IIIIlilIlUlIII I I IIIIIIlUlIIIIII

1 1 1 1 t t

iUlliillllllllll

I I I I I t

l

I ' . . . I . l J , l J l IIIU/IJ.I. t . t . . . . . . . . .

Interval 0 . 2 0 Sv Time 2001 days

Interval 0.100 Time 2001 days

Surface

Bottom

IIl|l ill|l, i

Velocity

I

''lJrfi i

s

|

'' I li' ' ''lI'

,

I ~' l''

t |

I

I

'

~ tl ~ tl

i ' ='

] ] ]

. : . : . ' . k . k L . Lt.!.: ..... .! "''ll]l' --- 0 . 1 0 Time 2001

m/s days Fig. 11.

p

I lllllll llllllll

II #1 l/lit_ ..... : I I .I I ! i I ~ t ~ . . . ~ ' ~ - i : I I ! I I, i /.~" --- _)~,~,1 :

i t t XX' +'+]"7 ~'~_[_A !\xxY\ ! liiU/t x\\%7 !1!11./+~.. I i . \ . IIIIll/ "

Density

. .

. .

. .

. .

'i

Velocity

. .

. .

. .

......... - - - -

~

. .

. . '

'

'

'

. . . . . . . . . . . . . . . . . . .

;.;.;..- ~.~.~.', '.'.:... -= 0.10 m/s Time 2001 days As Fig. 6 but for 2001 days.

used when integrating the temperature and salinity equations. These equations are maintained in their full generality. We have described the solution procedure in detail. In particular, at each time step a new density field is calculated from which the associated three-dimensional velocity field is diagnosed. This velocity field is then used in the temperature and salinity equations when calculating the density field at the next time step.

A long-time-scale, density-stratified shelf circulation model

137

The most important limitation of the model is the restriction to using a linear parameterization for the bottom stress. This can be in terms of either bottom velocity, vertically averaged velocity or bottom geostrophic velocity. However, on the long timescales for which the model is valid, the restriction to a linear parameterization is probably not a serious one (CsANADY, 1982). The model is also only valid in "wide" shelf regions; that is shelves for which the characteristic width is large compared to the internal Rossby radius of deformation. This is because information about baroclinic Kelvin waves is completely absent from the model equations. Fortunately, many shelf regions satisfy this requirement (notably, the Newfoundland and Labrador shelf, the Middle Atlantic Bight, the West Florida shelf, the South Australian shelf and the Weddell Sea shelf). We have illustrated the behaviour of the model by describing an experiment in which relatively light water is flushed into the model domain through the northern boundary and allowed to exit through the southern boundary. We described the boundary conditions used for this experiment in detail. A new method of handling closed boundaries (in our case, the eastern and western boundaries) was described. This successfully eliminates the boundary instability noted by KILLWORTH(1985) but has the possible disadvantage of forcing there to be no flow across sigma surfaces on the boundary (in a flat-bottomed model, this means that the vertical velocity is forced to be zero on the boundary). The "open" boundaries (in our case, the northern and southern boundaries) are treated in a manner similar to that used by BLUMBER~and MELLOR(1987). It was found to be important to include density mixing on the boundary (in particular, vertical mixing and mixing along the boundary) in order to prevent the development of spurious density gradients near the boundary. These can drive spurious flows through the J E B A R term. The offshore branch of the Labrador Current flows along the shelf-break and separates fresh, cold water on the shelf from relatively warm and saline water offshore (see Fig. 1). The model experiment we have described illustrates the development of a shelf-break jet given completely "unbiased" inflow conditions on the northern boundary. Indeed, the inflow water enters the model domain with uniform vertically averaged velocity and has a vertically uniform density profile which decreases linearly from the interior value at the eastern boundary. Furthermore the model ocean depth is a uniform 50 m along the inflow boundary. Within 1 year into the integration, a baroclinic shelf-break jet has developed, primarily in response to the offshore Ekman transport associated with the frictional stress in the bottom level of the model. Subsequently, however, both the J E B A R and bottom friction torque terms are responsible for concentrating the transport into the jet. As time progresses, the J E B A R term plays a dramatic role in splitting the jet and allowing the formation of a new jet offshore. The splitting of the jet is sensitive, however, to the parameterization used for the bottom stress (in this experiment, bottom stress is parallel and proportional to bottom geostrophic velocity). The importance of frictional stress associated with bottom geostrophic velocities in causing the offshore migration of an idealized density front has previously been noted by WRmnx (1989). Here, we have seen how this can also be important in the frontogenesis process itself. The model we have described is a generalization of that of HENDERSHOTrand RIZZOLI (1976) to allow for a three-dimensional density field. It is essentially a modification of those of HASSELMANN(1982) and KILLWORTH(1985) for application in a shelf/slope region. Acknowledgements--Thiswork has been funded by a grant from the Natural Sciences and Engineering Research Council of Canada, a grant from the NSERC/DFO Science Subvention Programme, a Department of Supply and

138

R . J . GREATBATCHand A. GOULDING

Services Contract and the Ocean Production Enhancement Network of Centre of Excellence. We are grateful to Moto Ikeda at the Bedford Institute of Oceanography for his encouragement, to Xiao Luo for running a number of experiments with the model and to Joy Simmons and Liz Crocker for typing the manuscript.

REFERENCES ASSELIN R. (1972) Frequency filter for time integrations. Monthly Weather Review, 100,487-490. BATTEENM. L. (1988) On the use of sigma coordinates in large-scale ocean circulation models. Ocean Modelling, 77, unpublished manuscript. BLUMBERGA. F. and G. L. MELLOR(1987) A description of a three-dimensional coastal ocean circulation model. In: Three-dimensional coastal ocean models, N. S. HEAPS, editor, American Geophysical Union, Washington, D.C., pp. 1-16. BRYAN K. (1969) A numerical method for the study of the circulation of the world ocean. Journal of Computational Physics, 4, 347-376. BRYAN K. and M. D. Cox (1972) An approximate equation of state for numerical models of ocean circulation. Journal of Physical Oceanography, 2,510-514. CLARKE A. J. and K. H. BRINK (1985) The response of stratified, frictional flow of shelf and slope waters to fluctuating large-scale, low frequency wind forcing. Journal of Physical Oceanography, 15,439-453. CSANADYG. T. (1978) The arrested topographic wave. Journal of Physical Oceanography, 8, 47-62. CSANADYG. T. (1982) Circulation in the coastal ocean, D. Reidel, Hingham, MA, U.S.A., 279 pp. FRIEDRICH H. and S. LEV1TUS(1972) An approximation to the equation of state for sea water, suitable for numerical ocean models. Journal of Physical Oceanography, 2, 514-517. GREENBERG D. A. and B. D. PETRIE (1988) The mean barotropic circulation on the Newfoundland shelf and slope. Journal of Geophysical Research, 93 (C12), 15541-15550. HASSELMANK. (1982) An ocean model for climate variability studies. Progress in Oceanography, 11, 69-92. HENDERSI.IO'IffM. C. and P. RIZZOLI (1976) The winter circulation of the Adriatic Sea. Deep-Sea Research, 23, 353-370. HUKUDAH., R. J. GREATBATCHand A. E. HAY (1989) A simple three-dimensional model of the circulation off Newfoundland. Journal of Geophysical Research, 94 (C9), 12607-12618. IKEDAM. (1986) Density-driven general circulation in a closed basin using a two-level model. Journal of Physical Oceanography, 16, 902-918. IKEDA M. (1987) Wind effects on the buoyancy-driven general circulation in a closed basin using a two-level model. Journal of Physical Oceanography, 17, 1707-1723. IKEDA M. (1988) A model study of wind- and buoyancy-driven ocean circulation. Journal of Geophysical Research, 93 (C5), 5078-5092. JAMES I. D. (1986) A front-resolving sigma coordinate sea model with a hybrid advection scheme. Applied Mathematical Modelling, 10, 87-92. JAMESI. D. (1987) A general three-dimensional eddy-resolving model for stratified seas. In: Three-dimensional models of marine and estuarine dynamics, J. C. J. NIHOUL and B. M. JAMART,editors, Elsevier Science Publishers, B.V., Amsterdam. KILLWORTH P. D. (1983) Some thoughts on the thermocline equations. Ocean Modelling, 48, unpublished manuscript. KILLWORTH P. D. (1985) A two-level wind and buoyancy driven thermocline model. Journal of Physical Oceanography, 15, 1414-1432. LYNCH D. R., F. E. WERNER,D. A. GREENBERGand J. W. LODER(1992) Diagnostic model for baroclinic, winddriven and tidal circulation in shallow seas. Continental Shelf Research, 12, 37-64. MELLOR G. L. and A. F. BLUMBERG(1985) Modelling vertical and horizontal diffusivities with the sigma coordinate system. Monthly Weather Review, 113, 1379-1383. MESINGERF. and A. ARAKAWA(1976) Numerical methods used in atmospheric models. GARP publication series No. 17, World Meteorological Organization. MIDDLETONJ. H., T. D. FOSTERand A. FOLDVIK(1987) Diurnal shelf waves in the southern Weddell Sea. Journal of Physical Oceanography, 17,784-791. M1DDLETON J. F. and D. G. WRIGHT (1990) Coastally trapped waves in a stratified ocean. Journal of Physical Oceanography, 20 (9), 1521-1527. PETRIE B. and A. ISENOR(1985) The near-surface circulation and exchange in the Newfoundland Grand Banks region. Atmosphere-Ocean, 23,209-227.

A long-time-scale, density-stratified shelf circulation model

139

PETP,aE B., S. AKEMaEAD, J. LAZIER and J. LODER (1988) The cold intermediate layer on the Labrador and Northeast Newfoundland shelves, 1978-1986. NAFA Science Council Studies, 12, 57-69. SARKISVAr~A. S. and V. F. IVANOV(1971) Joint effect of baroclinicity and bottom relief as an important factor in the dynamics of sea currents. Izvestiya of the Academy of the Sciences of the U.S.S.R. Atmospheric and Oceanic Physics, 7,173--188. SHAW P.-T. and G. T. CSANADY (1983) Self-advection of density perturbations on a sloping continental shelf. Journal of Physical Oceanography, 13,769-782. WRIGHT D. G. (1989) On the alongshelf evolution of an idealized density front. Journal of Physical Oceanography, 19,532-541. WRIGHT D. G., J. R. N. LAZIER and W. ARMSTRON~ (1988) Moored current and pressure data from the Labrador/ Newfoundland shelf, June 1985-July 1987. Canadian Data Report of Hydrography and Ocean Sciences No. 62, Bedford Institute of Oceanography, Dartmouth, Nova Scotia, Canada.

APPENDIX

1

Definition o f s y m b o l s used in the text )~, q~ z

longitude and latitude, respectively vertical coordinate measured upwards with the ocean surface at z = 0 and the ocean bottom at z = - H

a

z/H

u, v, w f~ H p p P0 P' = P - P0

eastward, northward and vertical velocity, respectively see equation (5) depth of the ocean pressure density representative density for sea water perturbation density buoyancy radius of the earth Coriolis parameter vertical eddy viscosity horizontal eddy diffusivity vertical eddy diffusivity surface wind stress bottom stress bottom friction coefficient volume transport streamfunction bottom velocity bottom geostrophic velocity vertically averaged velocity

b = gP'/Po a f v AH KH (r~, ~ ) (v~, ~ ) r ~p (ub, Vb) (Uby,/)bg) (~, ~)

APPENDIX

2

Parameterizing b o t t o m stress in terms o f b o t t o m velocity In this case, r b = por(ub vb) where (ub, Vb) = (U, v) evaluated at a = - 1 . In order to solve equation (15) for ~, we must express (Ub, Vb) (and hence rb) in terms of ~p and the buoyancy b. We begin by using equation (34) to write the m o m e n t u m equations (1) and (2) in the form

ifq = R + G + ~ ( v q o ) o /1-

(41)

140

R . J . GREATBATCHand A. GOULDING

where

R = - {ac@s~ (HA~ - H~B} + i ( H A , - H~B}}

(42)

- 1-~I1---~OPs+ iOPsl G= i= ~

and

q = u + iv. The

ao0 [cos tp O~.

(43)

0~J

surface and bottom boundary conditions are

lvqo=rs/po=(~+ir~)lpo

at

a=0]

and

[. 1

vqo=rq

at

(44)

o=-1

We now follow a procedure similar to that used by LYNCH et al. (1991) and define q 1 and

qz by

ifql=R +H~-f(vqla)° [ with

~l v q l o = r s / P 0 1 ~Vql a = rql

at at

a=0

[

•

(45)

a = -1j

and

ifq2 = with

1 + H1---g(vq2o)o

1 ~vq2~=O 1

vq2~,=rq2

at

o=0

at

(46)

cr=-I

Since R can be calculated from b, we can solve for ql at each time step. q2, on the other hand, is independent of the evolving density field and need only be solved for once, at the beginning of the integration. It follows that both ql and q2 are known at each time step. Now, since G is independent of the vertical coordinate, or,

q = ql + Gq2.

(47)

Vertically averaging (47), and using (12) to express g in terms of % gives i

_

G=[a~(-~O~+~a)-ql}/q2

_

(48)

and, hence,

q:ql+qZ{a~(-W~+~W~)-~l}/~12.

(49)

Evaluating (49) at a = - 1 now gives an expression for bottom velocity, qb = (Ub, VU), and hence bottom stress, ru = rqb, in terms of ~0and known quantities involving b. Substituting the expression for rb into (15) then gives the following equation, which can be solved for ~0.

141

A long-time-scale, density-stratified shelf circulation model a {ru 2 ~Pa \

0 [ru2~P,~cosdp\]

[ a fry 2a~p~

= PolOS, k " ]

a (rv 2aWl l

~(~)}

a

+ JEBAR- a{~-(~)-

f

~-~\--------~)j

(50)

where ul, vl, u2, v2 are all real, Ul + iv1 -----(ql -- q2ql/q2)b U2 + iv2 = q2b/q2 and subscript b denotes evaluation at 0 = - 1 . Once ~p is known, the three-dimensional (u, v) field can be obtained from (49) and F~ from (4). Finally, it should also be noted that -g[

~

,

.

where ~/is sea level [this follows from (43)]. Since G can be calculated from (48), it is possible to construct the sea level field predicted by the model, for comparison with observations.