Finite element simulation of water circulation in the North Sea

C. A. Brebbia and P. W. Partridge*

Department of Civil Engineering, University of Southampton, Southampton S09 5NH, UK (Received February 1976)

The modelling of tidal effects, storm surges and currents in large bodies of water is considered. The solution is attempted using the evolutionary shallow water equations with velocities and wave heights as unknowns. Two finite element simulation models are described based on six noded triangular elements. Special consideration has been given to the adequacy of the models which w e r e applied to the North Sea only after extensive tests in channels. Results for velocities and wave heights are compared and discussed. A set of conclusions on the applicability and scope of the models is presented.

Introduction
This paper is concerned with the modelling of tidal effects, storm surges and current patterns in large bodies of water. The solution is attempted using the shallow water equations, which are evolutionary equations with velocities and wave heights as unknowns. They require the initial conditions as well as the boundary conditions to be known. The solution of these equations is usually found by applying a numerical technique. The method used is of fundamental importance. In a finite element or finite difference approach the grid size will determine the type of phenomenon which can be investigated. In addition grid size relates to stability criterion and accuracy in evolutionary problems. The refinement of a model, though desirable in principle, may demand a large number of parameters which require more experimental data. These data can be difficult to obtain and produce a new type of error affecting the confidence one can have in the results. The analyst usually has to compromise between having a sophisticated model or a practical one, giving reliable results for the variables under consideration. In addition, large models are expensive to run. We describe here two finite element models. Both models have been developed using six noded triangular elements, but one is based on an implicit integration scheme, the other in an explicit one (the former allows for elements with curved sides). Special consideration has been given to the adequacy of the models and only
*Present address: Federal University, Porto Alegre, Brazil.

after extensive tests on channels 1 were they applied to North Sea studies. The North Sea is an important and busy seaway, especially since the discovery of gas and oil. From the numerical point of view, the area is wellconditioned, North Sea topography being regular and changes in depth gradual. Nevertheless, the models can and have been applied to different regions (e.g., the Solent in England 2. Other models of the North Sea exist: an explicit finite difference scheme by Heaps a, an implicit three-node finite element one by Grotkop 4 and a quartic quadrilateral finite element model by Davis and Taylor 5. The present model is based on the shallow water equations which are vertically averaged versions of Navier-Stokes equations, and take into consideration tides, bottom friction, advective forces, coriolis, wind tangential stresses and atmospheric pressure gradients. Results for velocities and wave heights in the North Sea are compared and discussed. A set of conclusions on the applicability and scope of the models is presented, indicating areas where further work is required.

Shallow water equations

The evolutionary equations used in marine and certain types of estuarial modelling are called the shallow water equations. They are a vertically integrated version of Navier-Stokes momentum equations and the continuity equation which acts as a constraint condition. In addition initial and boundary conditions have to be fulfilled. The different assumptions involved

Appl. Math. Modelling, 1976, Vol 1, September 101

Simulation of water ckculation in the North Sea. C. A. Brebbia and P. W. Partridge

,X 3 v 3

-q
D

X I vI

.....

Figure I Geometrical notation for the shallow water equations

are treated in detail elsewhere 2'6. The two shallow water momentum equations are:

c~- + VI-~x1 + Vz~x 2 = B1

(1)
(~V2 (~V2 V (~V2 = B 2

position of the seabed, the shape and variation in shape of the coastline, friction between the seabed and the water, hence the material of the seabed, the meteorological conditions, including wind, etc. Although the circulation of the earth and the astronomical forces of the sun and the moon act on the water as body forces, the main cause of tidal water movements in areas such as the North Sea is the driving force caused by tidal motion of the water on the boundaries of the area under consideration. The shape of the land surface containing the body of water is usually very complex, in some cases not even static, though the effects of erosion generally occur over too large a period of time to be important. Bottom friction is introduced in the model via Chezy coefficients. The inadequacy of using constant Chezy coefficients for all the model is evident. The different materials making up the seabed have different frictional resistances as the water depth and the velocities change. It must be pointed out that bottom friction and wind are of great importance in the movement of shallow water. The main causes of inaccuracies in tidal predictions are the wind forces and atmospheric pressure variations, which are important for large areas such as the North Sea.

where Bl
= ~'~V2 - - g ~ X 1 rl

Boundary and initial conditions

(2)
272 s

B 2 ----- - - ~ ' ) V 1 - - g O x 2

--

Vii are the averaged velocities:

=-

v,

fi

'L

l) i d x 3

(3)

H is the total depth, H = r / + h, where q is the wave height above a certain datum plane and h is the depth from the datum to the bottom of the sea. x3 is the coordinate in the vertical direction (Figure 1), f2 = 2e~ sin is the coriolis coefficient, q~ is the latitude and co the angular rotation of the earth, g is gravity, p the water density and p, the atmospheric pressure. The surface and bottom stresses are written as:

The solution of equations (1) and (5) require the knowledge of the corresponding boundary and initial conditions. The boundary conditions of the model are of two types: (a) fixed or land boundaries such as those given by the coastlines, where the normal velocities are zero, and the tangent velocity can be set free; (b) open boundaries where the elevation of the sea level (or the normal component of velocity) is prescribed. The determination of the initial conditions requires the knowledge of the free surface position at t = 0. Usually this knowledge is not possible and the models have to be started with zero elevation and zero velocity conditions. This is called a 'cold start'.

Finite element model

In order to build finite element models the two momentum equations (1) and continuity (5), including influx-type boundary conditions have to be written in the following weighted residual way:

zi

= ~ ~-(W

7 W~

2 1 p

Jr-

W2)1/2

i = 1;2

(4)

= -

(V2~+ V~) i/2

i = 1,2

ff,ov, (
ff{

~11 + V2~x 2 - B 1 & + Vl~x

~V1}

bV, d A = O
(6)

c is the Chezy coefficient, W~ are the wind speed components and ? is a parameter related to atmospheric density p, (usually given as a constant multiplied by p,). In addition equations (1) have to satisfy the vertically integrated continuity equation, i.e.,

f;f ~V2+ V,,~V2 ~Y2 ex~+Y2~ x2-B2 } bVzdA=O

+ ? I(HVi) +

6H dA = f(HV. H~)3H dS

~t +

(HV,) +

(HV2) = 0

(5)

The continuity equation is usually integrated by parts to render a simpler expression. This integration gives:

The systems of equations (1) and (5) describe the movement of large bodies of shallow water. The factors affecting the movement are many: the morphology and

+ HV2 ~x2

&

j dA
(7)

= fHV.6HdS

102

Appl. Math. Modelling, 1976, Vol 1, September

Simulation of water circulation in the North Sea. C. A. Brebbia and P. W. Partridge

F2

(10)

F.
v

b
Six node elements. (a) Straight sides: (b) curved sides

or more simply,

Figure 2

MQ + KQ

= r

(11)

The above weighted residual statements [equations (6) and (7)3 are the starting point for the finite element models. Assume that over an element the same interpolation function applies for the 1/1, V2 and H unknowns, i.e. Vl
=

Formula (11) is valid for each unconnected element. The next stage is to assemble all the element equations into a global system and impose boundary conditions in H and V~. To eliminate proliferation of notation the global system will be defined with the same notation as equation (11).

Time integration
Two time integration schemes were used, one an implicit and the other an explicit scheme. The implicit integration procedure is the trapezoidal rule. Starting with:
MQ + KQ = F

~)~1,

V2 = ~bt~2,

H = 4)/t"

(8)

is the interpolation function and I~'~,/-P are nodal values of Vi, H. In what follows six noded triangular finite elements with curved boundaries were used in order to define the boundaries better (Figure 2). These elements are called isoparametric and can be formulated by a simple coordinate transformation, the details of which have been given by Connor and Brebbia 7. Curved elements have the important feature that they tend to eliminate the spurious forces that may be generated on the boundaries by straight side elements joining at an angle. Substituting equations (8) into (6) and (7), one obtains: M~ + K ~ - f~MW2 + G1 n n + F I = 0

(12)
Q _ Q, + Q o 2
(13)

one assumes:

Q=Qt-Qo At
F-

Fo + F, - 2

Hence equation (12) becomes:

(9)
or

(2M+K)Q,=(Fo+

2 M -K)Qo F,)+ (~

(14)

MI~2 + f~MW1 + KI~2 + G2H n + F2 = 0 and

X*Qt M/-P - CI ~ - C2 V~2+ FH = 0

= F*

(15)
(16)

The recurrence relationship is then: Q, = (K*)- 1F*

where

dpdA

G,= g fr~r,,dA
F, = 4, T ~ P~

M= fCp~dA
+ J~ --if(W, + W2) '/2 dA
i=1,2

,, dA

c, = f4,SH4,dA,
and ( ),i = ~x,'

Fn =

fHv. TdA

The K* matrix to be inverted generally is a large nonsymmetric banded matrix of size approximately three times the number of nodes by six times the element band width (i.e., the maximum difference between element nodal point numbers plus one). The computer program has been optimized by taking boundary conditions into account in such a way that the corresponding rows and columns are eliminated from the element matrices before assembling. This significantly reduces the maximum size of the global matrix. It was also advantageous to store the matrix in a one-dimensional form such that only one and not two addresses need to be evaluated each time an element of the array is accessed. The explicit time integration scheme used was the well-known fourth order R u n g e - K u t t a scheme.

( ) = ~t

Equations (9) can be written as

North Sea model

The above finite element formulation has been applied to model the North Sea. This is a shallow sea varying in depth from under 50 m in the south to 400 m in a trench off the coast of Norway. Depths were obtained from Admiralty charts Sections were drawn at different angles across the whole region to determine

1v~ (B o

I" il[
tiM
C

V[ H "J

--

--

C2

Appl. Math. Modelling, 1976, Vol 1, September

103

Simulation of water circulation in the North Sea: C. A. Brebbia and P. W. Partridge

I-ERWICK BERGEN

ALTIC SEA

On the land boundaries, the no-slip boundary condition 1'11= V2 = 0 is specified. This assumption simplifies the necessary computing and is reasonable as the North Sea is a regularly shaped region. Because of the imposition of this condition, curved boundary elements, which are more expensive to run, were not necessary. By contrast when modelling the Solent, curved sided boundary elements were used allowing the tangential velocity to remain free 2.

Stability and accuracy

The smallest stability limit for the North Sea as given by the Friedricks-Lewy-Courant condition, occurs for an element off the Norwegian coast. It gives: At ~ 450 sec
I r i 3

(18)

Km

250

Figure 3

Finite element mesh for the North Sea

the best locations for nodes and the more accurate way of representing the bottom topography. Elements were carefully positioned in order to obtain the best possible representation of the topography using a predetermined number of nodes. The final mesh comprises 228 nodes and 97 six noded elements as shown in Figure 3. The tidal characteristics of the N o r t h Sea are complicated, tidal amplitudes vary from zero to six metres and high water times change throughout the cycle around each of the three amphidromic points. Tidal heights for boundary conditions were taken from the charts of co-tidal lines. The waveheight forcing functions were specified in the form (Figure 4)

~ = ao/sln[- ~- +e

(17)

on each of the extremes of the tidal boundaries, with the intermediate heights being linearly interpolated. This formulation approximates the ff~ost important tidal component for the North Sea. These curves have been taken from the Admiralty Tide Tables. They are assumed to be referred to the same datum since other information is not available. Because of this, the results presented in this paper may not be quantitatively correct but the comparison between models are still valid. The charts of co-tidal lines give the approximate tidal range and high water times for the interior points on the grid. The circulation pattern and velocity magnitudes for parts of the Sea may be seen in the tidal stream atlases 8,9. Wind and storm surges were not modelled as this series of tests was carried out to investigate the general performance of the model. After several tests it was decided to take a continuous tidal boundary in the Northern part from node 41 to 1 and 1 to 11 (Figure 4), otherwise instability originated from the Shetland Islands element, showing that one element is inadequate to represent a discontinuous tidal boundary properly. Tidal conditions were also specified at the Dover Strait (points 226 to 228) and after a number of trials also for the Baltic Sea boundary (nodes 86-87-109). It was found that inclusion of the Baltic Sea improved the waveheight results.

The worst case on a tidal boundary gives At ~<650 sec. This is important as instabilities always start at these boundaries. The average value is around 900 sec and for the shallow southern North Sea the criterion suggests a limiting time step of less than 2000 sec. The explicit programme which uses a fourth-order R u n g e - K u t t a procedure was run with a time step of 600sec. For the implicit programme instead a time step of 30min was used. To obtain stable results with both models requires the application of special techniques. For this work three different techniques were used. The first and simplest of them is to work always with a constant value of friction over all the region, starting with a low value (c = 10ml/Z/sec) and increasing it by 10 over two to four cycles. This technique did not give good results and the solution tends to become unstable for large values of c, i.e., small values of friction. In addition it is unrealistic o assume that the Chezy coefficient will be the same over all the domain. The second technique was to prescribe a higher order of friction for the elements on tidal boundaries and a smaller value for internal elements. This is because there is a general tendency for the tidal boundary to generate disturbances. These perturbations may be due to a number of factors and cause the transmission of short waves through the system. The specification of higher values of friction for the tidal boundary elements reduce the propagation of errors. It was decided to apply a Chezy coefficient c = 20 ml/2/sec at the boundaries, which produces stable results.
i

6'

E) 2 O

Hours: 8"5 Before Bergen;

65

415

2.5 O-51 1"5 High water Dover , Kirkwall;

2[5 After , Dover

Figure 4 Tidal boundary conditions for the North Sea. . . . . . . ,

, Lerwick; . . . . .

104

Appl. Math. Modelling, 1976, Vol 1, September

Simulation of water circulation in the North Sea. C. A. Brebbia and P. W. Partridge

The third different stabilizing technique is to start with a re~tlistic value of friction from the beginning and try to remove the short waves by 'numerical smoothing'. This smoothing can also be applied every time the level of friction is decreased. The operation needs to be carried out for a number of steps and after the perturbations have been removed the solution does not need to be smoothed any longer. This numerical smoothing consists of taking for the next time step not the actual nodal values just obtained but weighted averages. These averages are calculated by weighting every nodal value by a constant and adding to it the weighted values of the neighbouring points. For the programmes described here the node under consideration has half the weight and the other half is distributed among six neighbouring points proportionally to their area of influence. (For boundary nodes only 3 neighbouring points are taken into consideration.) In general, the coarser the mesh the more relative weight the central node will have. The actual weight used does not seem to be too important provided that the coefficient for the node under consideration is reasonably large by comparison with the coefficients for the neighbouring nodes. A simpler way of weighting may be for instance, to multiply the solution vector by the mass matrix. Smoothing has been successfully applied by the authors in small estuarial areas, such as the Solent in England but it is less necessary for the North Sea as the system is more stable.

080 t

~
I I

-O80
,:

O80

8t
L I I I

-080 O 0-80 Velocity, VI (m/see)

Figure 5 Velocity ellipses for (a) implicit and (b) explicit models. A, Node 51; B, node 57; C, node 105; D, node 148

c = 20 (ml/Z/sec) and values of internal friction of 40 and 60. The solution was initiated from the implicit model results for c = 20 throughout instead of direct cold start. (This was done simply to save computer costs.)

Results
Velocity ellipses are useful to find out if steady state has been reached, if there are any disturbances present in the system which are likely to cause instability, the magnitude of drift velocities, the changes in velocity magnitudes due to changes in Chezy's coefficients, time step, etc. Two comparable sets of results obtained using the same tidal boundary Chezy coefficient (20 mX/2/sec) and the same internal friction (c = 40) are shown in Figure 5. It can be seen "that the ellipses for the explicit scheme tend to be smaller (this is also true for other points) than those obtained with implicit integration. The drift is also less in the explicit method. This may be due to the large time step used in the implicit solution. As a smoother circulation pattern and smaller drifts are obtained one could conclude that the explicit solution is more accurate in this case. There are also slight differences in the shape of the ellipses. The large drift (in the order of 10 to 20cm/sec for some points) may also be explained by the coarseness of the grid. A typical plot of velocities over all the Sea is also shown in Figure 6 (results are from the implicit programme, the Chezy coefficient is 20 on tidal boundaries and 40 inside), where they are compared against results published in the tidal stream atlases available in Britain. The general trend of the velocities compares well. We should be aware, however, that the tidal stream atlases velocities are smoothed out, the observations are made in only the top layer of the water. The programme instead yields depth averaged velocities given by equation (3). Hence the currents worked out by the programme are not exactly as the atlases currents. The programme results are being affected by local fluctuations in depth to a greater degree than the figures in the atlases. Graphs showing the waveheight solutions for the same tests and at the same points as the velocities are shown in Figure 7. There are comparatively high variations in tidal range and high water times. The flatness of the waveheight curves at nodes near an amphidromic point is also noticeable.

Tests
Many tests were run with different time integration schemes, friction coefficients and smoothing schemes but only two of them will be presented for brevity. The first, Test 1, uses implicit integration (At = 30 min) and the second, Test 2, explicit fourth-order R u n g e - K u t t a ( A t - 10min). Test 1 was started with a Chezy coefficient of 10 ml/2/sec, after two tidal cycles this was increased to 15 and after another two to 20. Then the friction 20 was left for the elements on the tidal boundary but the value of the internal friction was decreased to = 40 over four tidal cycles. Finally, the internal Chezy coefficient was taken as a variable given by: c = 151og~(0.9H) [in mU2/sec] (19)

This formula gives low friction in the interior of the North Sea (e.g. for H = 55 m, c = 60). The same friction was applied during 6 more tidal cycles to obtain repetition of results. In addition the results were numerically smoothed over 3 h (i.e., 6 steps) after change of Chezy's coefficient, then the smoothing was stopped. For the level of friction given by formula (19) the velocity ellipses tend to increase in size and their drift is accentuated as the level of friction is reduced. As the friction is variable, the results are not being obtained under constant conditions of damping and the ellipses do not quite close even after three or four tidal cycles. It is surprising that good results were reported by Davis and Taylor 5 after only three tidal cycles from cold start, using this variable friction formula. Test 2 was run with tidal boundary elements friction

Appl. Math. Modelling, 1976, Vol 1, September

105

Simulation of water circulation in the North Sea. C. A. Brebbia and P. W. Partridge

/_..

Conclusions
. . . . \

/.C-",, ' ~ ' ,

(,

"

"

'

'

'

"

..

\, ,,,
\','~ ,

.
",

,
.

.
,

::.) \-:-/

For a constant level of friction throughout the grid the results became unstable for a value of Chezy's coefficient of 60 ml/2/sec unless special procedures are used. Sudden reductions of the level of friction (i.e., increases in c) cause a small shock to be transmitted through the system. Short waves are also generated on the tidal boundary by the discrete changes in the imposed tidal height at each time step. A way of damping out those effects is by decreasing the c coefficients in elements on the tidal boundary and this technique has been used for all the results shown in this paper. Another damping procedure is the numerical smoothing, for which the model could be started with a realistic level of friction and the results numerically smoothed until the disturbance caused by

4.0
/-.' . .

20 A

Figure 6a Velocity vectors obtained with the implicit programme

2 hours after high water at Dover '-

Iv
. 0

,
60 12.O
b

'oL
40

>~

o 2.

O
ABERD

6"0 12-O Time after low water (Bergen)

FixTure 7 Wavelength for (a) implicit and (b) explicit models. A,

Point 51 B, point 57: C, point 105; D, point 148

GREAT
BRITAIN ~

"t

2,am

~IBANTWERP

Figure 6b Velocity vectors from the Admiralty charts 2 hours

after high water
2,Om I.Om

Finally the co-range lines chart for the North Sea has been computed from results obtained using the implicit programme with friction c = 20 ml@sec on tidal boundary and for interior elements c as given by formula (19) and shown in Figure 8. If the co-range lines are compared against the results shown in Admiralty chart 5058 the agreement is reasonable. Similar curves were also reported by G r o t k o p 4 and Nihoul 6.

Figure 8 Co-range tidal lines

106

Appl.

Math.

Modelling,

1976,

Vol

1, S e p t e m b e r

Simulation of water circulation in the North Sea: C. A. Brebbia and P. W. Partridge the cold start has been reduced. Since it is not necessary to alter the friction parameter after this the introduction of further short waves may in some cases be avoided. For the results reported in this paper the numerical smoothing technique was only used once for a short time, to stabilize the implicit model results when the Chezy coefficient changed from 40 to the value given by formula (19). Both time integration schemes, implicit and explicit, give similar results; however, the fourth-order R u n g e - K u t t a is more accurate using At = 10 min which is the highest allowable time step in this case. But it should be also pointed out that the implicit programme (At = 30 min) needs less than half the computer time required for the explicit programme. For evolutionary processes of the type here described computer time can be very expensive and the programme should be further optimized before undertaking production runs. It seems, however, that for a problem with the dimensions of the North Sea and a finite element grid similar to the one used here, the explicit programme may be more convenient to use. This is not the case when the domain is smaller (for instance for the Solent) or the grid very fine. In other words implicit schemes are more expensive per timestep than explicit ones but allow for larger timesteps. This can be of interest in problems where the time step may be increased considerably over a simpler scheme, without significantly affecting the accuracy of the results. The viability of the six nodes finite element circulation model for the North Sea has been established. Second order elements of this type are especially suitable to represent accurately the topography of the region (i.e., variable depth a n d curved boundaries).

References
1 2 3 4 5 6 7 8 9 Partridge, P. W. and Brebbia, C. A. Proc. Am. Soc. Mech. Eng. (J. Hydraulics Div.) 1976, in press Brebbia, C. A. and Partridge, P, W., in 'Mathematical Models for Environmental Problems', (Ed. C. A. Brebbia), Pentech Press, London, 1976, p 141 Heaps, N. S. Mere. Soc. R. Sci., Liege, 1972, 6, 143 Grotkop, G. Comput. Methods Appl. Mech. Eng., 1973, 2 Davis, J. M. and Taylor, C. Proc. Finite Element Methods in Flow Problems Conf., Swansea, 1973 Nihoul, J. C. J. (Ed.), 'Modelling of Marine Systems', Elsevier, Amsterdam, 1975 Connor, J. J. and Brebbia, C. A., 'Finite Elements for Fluid Flow', Bunerworths, London, 1976 Hydrographer for the Navy, 'Tidal Streams Atlas North Sea-Southern Portion', HMSO, London, 1962 Hydrographer of the Navy, 'Tidal Streams Atlas North Sea-Flamborough Head to Portland Forth', HMSO, London, 1962

Appl. Math. Modelling, 1976, Vol 1, September

107