- Email: [email protected]
Mathematical model of Upper Milford Haven—A branching estuary
-
Estuarine, Coastal and Shelf Science (1986) 23,403-418
Mathematical Model Haven-A Branching
V. Nassehi
of Upper Estuary
Milford
and D. J. A. Williams
Department of Chemical Engineering, Park, Swansea SA2 8PP, Wales
University College of Swansea, Singleton
Received 25 March 1985 and in revisedform 18 November 1985
Keywords: estuaries; tidal heights; salt intrusion; Haven estuary.
numerical
model; Milford
A mathematical model for hydrodynamics and salt intrusion in a onedimensional branching estuary has been developed. The mode1 satisfactorily predicts wafer levels and discharges throughout Upper Milford Haven estuary for spring and neap tides. There is also a measure of agreement between computed and observed salinities.
Introduction and other economic advantageshave resulted in the selection of the sheltered waters of Milford Haven for the location of industrial activities such asrefining petroleum and power generation. This growing industrial and commercial demand on tidal waters of the Haven requires careful management of water quality to avoid catastrophic effects on their ecology (Nelson-Smith & Case, 1984). Mathematical models which predict estuarine flow variables (water velocity and water surface elevation) and provide estimates of dispersion coefficients are useful tools for water quality management in this estuary. Although in general such models should recognisethe three-dimensional nature of flow in natural water systems, it is often possible from a consideration of the geometrical boundaries of an estuary to reduce the complexity and usesimplified models in either two or one dimensions (Dronkers, 1964). The outer part of this estuary (from Carr Jetty to St. Ann’s Head, Figure 1) has been modelled in fwo dimensions by Gunn and Yenigun (1985). In the present work we describe the development and use of a one-dimensional model for upper Milford Haven estuary.
Accessibility
Milford
Haven estuary-A
physical description
The physical environment of Milford Haven hasbeen described in somedetail by NelsonSmith (1965). The estuary is located in West Wales (U.K.) and it is of the drowned valley type. It comprises a network of channels connecting the inland waters of the Eastern and Western Cleddau and several smaller rivers to the sea(Figure 1). The estuary encloses 110 km of varied coastline and is 25 km wide at its mouth. 403 0272-7714/86/090403+
16 $03.00/O
0 1986
Academic
Press
Inc.
(London)
Limited
404
V. Nassehi arid D. J. A. Wiliiums
(Eastern
(Daucleddau)
Figure 1. Locations: a-St. Ann’s Head; &Port of Milford Haven; c-Carr Jetty; d-Hobbs Point; e-Pembroke Dock; f-Burton; g-Neyland; h-Lawrenny; i-Port Lion; j-Picton Point; k-East Hook; I-Slebech; m-Minwear Wood; n-Canaston Weir; o-Haverfordwest; pCrowhil1 Weir (0 water level observation station-1979 survey; 0 water level observation-1977 survey).
Beyond Carr Jetty and above the road bridge between Pembroke Dock and Neyland (Figure l), the narrowing main channel of the estuary (from this point landwards called Daucleddau) turns northwards until its limit at the junction of Eastern and Western Cleddau rivers at Picton Point, about 27 km from the estuary mouth. Within this upper part of the estuary water flow is predominantly one-dimensional. The mouths of other smaller rivers often form inlets which initially are of a similar width to the main channel, but usually are much shallower and generally their bedsare almost completely exposed at low water (West, 1978). These shallow inlets are thus thought to have negligible effect on the longitudinal momentum transfer. The tidal limit of spring tides for the Western Cleddau river is about 11 km [in the region of Crowhill Weir with a range of 3-9 m above ordnance datum (A.O.D.)] and about 8 km for the Eastern Cleddau river (in the region of Canaston Weir with a range of 4.2 m
Mathematical
model
of Upper
Miijord
405
Haven
A.O.D.) from Picton Point, respectively. The mean tidal ranges at Milford Haven are given as 6.3 m and 2.7 m for spring and neaps respectively (Admiralty Chart no. 3275). For both spring and neap tides there is a delay (max. 36 min) in the occurrences of high and low water between the port of Milford Haven and locations near to the tidal limits of the Western and Eastern Cleddau rivers (Nelson-Smith, 1965). The bed level of the Daucleddau varies from 15 m below ordnance datum (B.O.D.) near Lawrenny to about 5 m B.O.D. near Picton Point. The width of the main channel is of the order of 500 m. The Eastern and Western Cleddau rivers gradually decrease in width to a few meters towards their tidal limits. Except for a few pools there is very little water in these rivers at low water springs. Details of sediment distribution in the Daucleddau is reported by Williams (1971) and within the Cleddau rivers by West (1979). In the Cleddau rivers water mainly flows over a muddy bottom while in the Daucleddau reaches the bottom surface consists of sand and rock. During the 1977 survey of the upper estuary (West, 1978) the fresh water discharges of the Cleddau rivers were of the order of 1 m3 s- ‘. Another survey in spring of 1979 (Williams & Nassehi, 1980b) gives the respective mean daily fresh water flows of Western and Eastern Cleddau rivers as 6 m3 s- ’ and 5.5 m3 s- i. These flows are insignificant in comparison with the tidal flow in this estuary which is of the order of lo4 m3 s- ’ at Hobbs Point at mid-tide (Gunn & Yenigun, 1985). Although some vertical stratification has been reported in the uppermost reaches of the Milford Haven estuary (e.g. at Haverfordwest in the Western Cleddau and Minwear Wood in the Eastern Cleddau, Figure I), (Nelson-Smith, 1965) it is expected that seaward of these locations where the estuary channel is wider and deeper, vertical mixing will be almost complete. Recorded vertical salinity changes within the outer Daucleddau are less than 1%0on spring tides and less than 5%~ on neap tides (Gunn & Yenigun, 1985). For the purpose of over-all simulations of salt intrusion in the upper estuary by a one.dimensional model, localised vertical stratification effects near to the tidal limits of the Cleddau rivers have been ignored. Mathematical
model
The mathematical model consists of the numerical integration of the following governing equations (McDowell & O’Connor, 1977): (a) The one-dimensional hydrodynamical equations of continuity and motion (1)
-+---+--[- 1+- +H'p=o, dh dX
(b)
the one-dimensional
1
gA
aQ at
1 2
gAdx
Q2
lQlQn2
2y ax
pff4:3
A
salt transport
equation
as QdS t+--z-AG
i a
AtE [
=O. 8X
1
The variables are defined in the section entitled Notation.
(2)
(3)
406
V. Nassehi
Using the following
and D. 7. A. Williams
approximate
relationship p=o.75
equation (3) can be written
ss
1000,
(4)
as:
ah 1 aQ -+--+’ ax gA at
gA
+ --a Q2 +IQlQn’ ~2~413
[ 1 dx A
0*75H
as
2000t1.5 SG =O
(5)
whence equations (l), (3) and (5) can be solved in a combined manner. In the present model for reasons of convenience and mathematical rigour the nonlinear hyperbolic partial differential equations representing the hydrodynamical conditions in the estuary [eqns (1) and (5)] are solved by a finite difference technique (Gray, 1980) while the linear parabolic partial differential equation of the salt transport [eqn (3)] is solved by a Galerkin finite element method. After division of the solution domain into different segments according to the four point implicit finite difference scheme used herein, a variable f and its temporal and spatial derivatives can be approximated as (Liggett and Cunge, 1975): (6)
1
A!f(XJ)C+-& fi:: ++, +Sj+l-fi ’ i
This finite difference scheme is dissipative (ASCE, Task Committee, 1982) thereby damping out numerical dispersion; it is independent of variability of segment length and allows an extremely simple representation of boundary conditions (Cunge et al., 1980). Conclusions concerning stability and accuracy of the linearized hydrodynamical model based on this scheme are given by Liggett and Cunge (1975) as: (a) the scheme is unconditionally stable when 0.5 < f3< 1, (b) some parasitic oscillations appear for values of 0 near to 0.5 which disappear if 0 is chosen to be greater than 0.66 or the frictional resistance is enhanced, (c) the scheme is of second order accuracy when 8 = 0.5 and of first order accuracy for 0.5<8< 1. These conclusions can in general be extended to the cases of non-linear hydrodynamical models (Weaver, 1976). Application of the four-point implicit scheme to the hydrodynamical equations (1) and (5) results in the following two working equations: Q{;;-Qi’
hj:;
1
-hi+
+G1(hj;:+hj+‘)-G2=0,
’ + G,(Qi;;
+ Qi”)
These algebraic equations relate the prime unknowns
- G,=O.
(10)
of h and Q at the boundaries of every
Mathematical
model
sf Upper
Milford
Haven
407
segment. Coefficients G,, G,, G, and G, in equations (9) and (10) are relations in terms of space and time intervals, time weighting factor, hydrographic data, bottom friction, salinity and the values of the prime unknowns at previous and present (due to nonlinearity) time levels (e.g. see Williams & Nassehi, 1980a, for the case of 0= 1). Based on the working equations (9) and (10) a set of algebraic equations describing the hydrodynamical conditions along each branch of the Upper Milford Haven estuary is obtained. Added to this system of equations are the following hydrodynamical relations at the junction of the three branches: (a) the continuity condition is expressed as (Elliott, 1976): algebraic sum of the discharges across the boundaries of the junction
=
accumulation within the junction.
or
(11)
which after time-stepping
gives
Ce=A,h:+l - hi ; At (b) the momentum condition at the junction of relatively of equal free surface elevation at its boundaries
h,=h,,=h,,,=h,.
(12)
small rivers is the assumption
(13)
For larger junctions this assumption should be modified (Zi-Cai et al., 1983). The third working equation in the present model is derived by the application of the Galerkin finite element technique (Zienkiewicz, 1977) to the salt transport equation [eqn, (3)]. This equation includes two physical mixing mechanisms: convective longitudinal transport and longitudinal diffusion. These mechanisms are represented by the second and third terms, respectively in equation (3). It is known that in cases where there is strong convection [i.e. when the second term in equation (3) is dominant] the numerical solution of this equation by finite difference or finite element methods may become oscillatory and useless (e.g. see Holly & Usseglio-Polatera, 1984). To avoid this problem specially modified finite element solution schemes applicable to convective diffusion equations are widely used (e.g. see Morton & Parrot, 1980). However, in the present work preliminary results obtained using a standard Galerkin finite element solution of equation (3) did not show any parasitic oscillations and the third working equation of this model is based on the standard technique. Using linear Lagrangian shape functions and isoparametric elements (Zienkiewicz,
408
V. Nassehi and D. J. A. Williams
1977) for every element along the single channel of each branch, the Galerkin finite element form of equation (3) is
1/
-
\
dS,
I-
i-1
1
si
AAX +
6 si+1 (1
iI I1
1 -1 AE
1which after time stepping gives: 8Q
--y+3at
RAE A LIX Ax
eQ
si
=o
AX
Ah
7
BAE-
(14)
si+l
S;+’
-+1-ax 6Lt
=
-
ALIX -----6nt
OQ
Anx
l-8
e...v++.---Q--AE
2
Ax
3Llt
2
Aax
l-8
-+-
6nt
BAE AAx
2
! OQ ( eAE
3nt
1-Q AX
l-f? Q+nxAF
2
ALlX
---
nx
s.f+’ 1+1
1
l-8
l-8
Q+nxAe
6nt
2
ALiX
l-f?
--3nt
-Q-.xAE 2
-
l-8
The accuracy and stability of this scheme in relation to numerical diffusion and dispersion is presented by Smith et al. (1973). For practically large enough At values, 8 > 05 should be used (herein 0 = 0.67). This time step may be appreciably increased if Q/A (velocities) values are small and the physical dispersion coefficients (E) are large. Overall accuracy of the scheme depends on the total number of nodes (degree of mesh refinement) and is not affected by the use of second or third order elements instead of linear ones used in the present work. Based on the working equation (15) elemental matrix equations for all of the elements along the estuary branches are derived. In addition to this system is the extra equation arising from the salt continuity condition at the junction of the three branches. The necessarysalt continuity relation at the junction is: total amount of salt entering the junction
=
total amount of salt leaving the junction
+
accumulation of salt within the junction,
--
Mathematical
model of Upper Miljord
409
Haven
or (16)
CQS= ; jSc d’V,. If we assumethat mixing at the junction is complete then:
EQS=A~~
(17)
S,h,,
wherein
s, = s,,= s,,,= 41 + Sill 2
(18)
.
Using relation (18) and after time stepping we have: Q, j+l
sj+l,
-
Q:I+l
S;,+‘-Q;,+’
;4 (Sl, + S&j+ ’ hi+At’ -(S,,
S,:,+l
+ &)j hJ;
=
(19)
Let us write 2nt
j+l
al=yQl
9 c
2At a,=rQ,,
j+l
+e+’
(20)
c 2At
.+l
a3= A
a,,,
+hj+‘,
‘2 and a4
= + h,’ (S,, + S,,,)j.
Then after somealgebraic operations and multiplying the sidesby a large scalarparameter Aequation (19) gives: &a,
S{+‘-a,
S&+’
-a3
S{i’
+a,)=O.
(21)
The matrix equations basedon (15) are assembledtogether with those arising from application of (21) at the three boundaries of the junctions. The parameter A is chosen to be large (109) ensuring that in the finally assembledmatrix the coefficients applicable to salinities at the junction boundaries (i.e. S{” ‘, S{,+’ and S{,: ‘) are totally dominated by those of (21). In this way the constraint imposed by equation (19) is satisfied.
410
I,‘. Nassrhi
and D.
Application
3. A. Williuwe
of the model to the Upper Milford
Haven estuary
The division of the Daucleddau region and the Western and Eastern Cleddau rivers into segments (or elements) of variable length is shown in Figure 1. The shaded triangle contained by segment boundaries (1 l), (12) and (23) at Picton Point (Figure 1) represents the confluence of the Cleddau rivers. The solution algorithm assumesfirstly that at any time level the salinities along the estuary are given, therefore after insertion of the boundary conditions the system of algebraic equations representing the hydrodynamical relations becomesdeterminate and is solved using an iterative procedure (Williams & Nassehi, 1980~). The discharges so obtained are used to solve the salt transport equation (subject to salinity boundary conditions) at the next time level. These stepsare then repeated for a whole tidal cycle. Since the density term in the hydrodynamical equation of motion is relatively small the present algorithm which computes h and Q one n tout of phasewith S produces little error in their values. To obtain optimal accuracy without any danger of instability the time weighting factor B is chosen to be 0.67. The time increment (a t) is 250 seconds.Herein smaller (A t) values increase computational time without improving accuracy. For a typical spring tide of 125 h duration the present combined model converges after two tidal cycles. This required 270 secondsC.P.U. time allocation on an ICL 2966 computer. This time requirement is short enough to justify the useof the non-linear version of the hydrodynamical equation of motion [eqn. (5)] for any possible enhancement of the accuracy in the numerical solution. The model requires the following input data. (a) Hydrogruphic data. These data are values of cross-sectional area, hydraulic depth and water surface breadth at each segment boundary (or node) along the estuary. Herein values of A, H and B are arranged in double entry tables as functions of water surface elevation and location. The first row in these tables contains hydrographic data for each location corresponding to low water and the successiverows are in 0.2 m intervals. All of the geometrical data in the present work are with respect to O.D. Newlyn. Water surface breadth and low-water cross-sectional area landwards of Lawrenny are calculated utilizing echo sounding data of the 1977 survey of the Upper Milford Haven estuary (West, 1978). Corresponding data for the reach between Carr Jetty and Lawrenny (Figure 1) were found using Admiralty chart no. 3275. (b) Manning’s friction coefficient. Values of the friction coefficient ‘n’ were initially selected with the help of literature (Graf, 1971). Then using a trial and error procedure they were adjusted in somesections, in order to minimize the discrepancy between model output and field data observed during the spring tide of 25 April 1979 for the locations indicated in Figure 1. The finally accepted values of ‘ n ’ in the present work are 0.018 along the Daucleddau region and 0.02 along the reaches in the Western and Eastern Cleddau rivers. Boundary conditions. These data are water surface elevations and salinities observed at estuary mouth and fresh water flows and zero salinities at estuary head. (c)
(d)
Initial conditions.
At time t = 0 any set of arbitrary values of h, Q and S can be used,
Mathematical
model of Upper Milford
Haven
411
asthe initial conditions to start the solution algorithm. In practice, approximate low water values for h and Q, together with zero salinities are used. coeficients. Inherent in the one-dimensional approach to the modelling of unsteady state distribution of salt in an estuary is the difficulty of predicting dispersion coefficients. Although some progress has been made in theoretical analyses of mixing mechanismsin estuaries the relationships derived are complex (Smith, 1980; Chatwin & Allen, 1985) and cannot be readily used in practical cases.The longitudinal dispersion coefficient is best determined by field observation in conjunction with the direct useof the salt transport equation [eqn (3)]; however, the quantity of field data required for such a calculation is considerable. In order to avoid this problem someinvestigators have chosen a constant value for this coefficient throughout flow channel during the tidal cycle (West & Lin, 1983). In the present model a method based on an empirical formula suggestedby Thatcher and Harleman (1972) is used. This empirical relation in metric system of units is written as: (e) Dispersion
&=64K1
n(Q/A)
H5”6 + K2
(2%
values of K,, K, and m in this equation depend upon various physical factors and must be determined separately for each estuary and flow condition. In equation (22) the first term with K, = 1 is an analytical expression which allows for dispersion by velocity gradients over the flow depths (O’Connor & Thompson, 1976). In real flow systems, transverse velocity gradients, junctions and bends together with wind action produce extra dispersion and this is reflected by an increase in the value of K,. For example, Thatcher and Harleman (1972) used K, = 3 on tidal rivers on the eastcoast of U.S.A. but other investigators (e.g. seeMcDowell & O’Connor, 1977) have suggestedan increased (sometimesas much asten times) values for this constant in narrow and tortuous conditions. The second term of equation (22) is an empirical term which provides a contribution to Edue to the presenceof longitudinal density gradients. Results and discussion Following the calibration of the model using the spring tide of 25 April 1979 asthe input tide, the model was successfully validated by the insertion of water surface elevations at Carr Jetty observed during the neap tide of 2 May 1979 whilst the values of ‘ n ’ remained unchanged. The model was alsoused to predict the spring tide of 2 August 1977 with the samevalues of ‘ n ’ throughout the estuary. The comparisons between the observed water levels and computed results for these tides are given in Figures 2(a)-(d). These comparisonsshow a very close agreement between the observed data and computed water levels. These curves indicate that the incoming tidal wave maintains an almost symmetric shape throughout the Daucleddau region. However, as it propagates into the branches it is increasingly distorted. The maximum time lags between high water at Carr Jetty and the tidal limits of Cleddau rivers is about 20 min. Computed dischargesfor the spring tide of 25 April 1979 for various locations together with the corresponding depths mean velocities are given in Figures 3(a)-(d). Whenever available these figures include instantaneous velocities recorded by West during 1977 survey of the estuary (West, 1978) and the measureof agreement between these values is shown.
412
V. Nassehi
and D. J. A. Williams
7
4
t
(0 ) Locotlon
I Port
4
Lion
(b)
Location
1
Plcton
Point 0,
4 t 30 30
// 0
I
0
5
‘iYI-u.Ll 0
I2
3456789101
0
(c ) Location
I I I I I23456789101112
I
I
I
I
I
I
I
I.
I23456789101112
h : Lawrenny
I
I
I
I
I
I
I
1:
I
0
Time
I
2
3
4
5
6
7
8
9
IO
II
I2
(h )
Figure 2. Comparison of computed and observed water levels observed during the spring tide of 25 April 1979, H the spring tide of 2. August 1977, 0 water levels observed 1979; the time of origin is 7.30 a.m. B.S.T.; model -- - - model output for the neap tide).
surface elevations (0 water water levels observed during during the neap tide of 2 May output for the spring tides,
Mathematical
model of Upper Miljord
Haven
413
I : 0 ) LocatIon
3000
I : Port
LlOrl
t
f
( c ) Locotlon
k
East
Hook I-O
0.5
0.0
1 -0/
-0.5
1 t I
I
I
I
I
I
I
1
I
I
I
I
I
4.0
f
Y) E La c :: ?
( b ) Locotlon
1 Plcton
pomt
t
( d 1 Locotlon
I.0
1 Slebech
0.5
o-o
\ -0.5
3000
t
,
,
0
I
23456
,
,
,
,
]
]
,
,
,
,
7
8
9
IO
II
12
bl 0 Time
Figure 3. Computed discharges, 25 April velocities, 2 August
discharges 1979; ---1977).
(h
I
I 2
I 3
I 4
I 5
I 6
I 7
I 8
I 9
I IO
I Ii
I 12
)
and velocities (time of origin is 7.30 a.m. B.S.T.: ~ velocities, 25 April 1979; 0 observed instantaneous
The constants (i.e. K,, K, and m) in equation (22) were adjusted to find appropriate 1: values. This involves a trial and error procedure in which model output is compared with the observed salinity data for upper Milford Haven reported for tides of similar amplitude by Williams and Jolly (1975) and West (1978). The finally accepted values for these constants under various conditions are: For the spring tide of 25 April 1979 and corresponding fresh water flows (6 m3 s- ‘) for Eastern and Western Cleddau rivers, K, = 10, K, = 30 and m = 2.
-1.0
414
I’. Nassehi
and D. J. A. William
I
‘**
-
150 100 50
-
/
Locot~on a ,-. \
o
Corr
c
Jetty /
‘\
/ \
\
/
\
/
\
-
-1
./
\
,’
\
\
I
/
;
\d/
\J/
IrllrlllllYI 12345670
0
-E
/
(a)
9
IO
II
12
I
0
0
1
123456789
I
2
IO
3
4
5
6 Time
Figure 4. Dispersion coefficients (----E b, etc. correspond to the time of maximum
7
8
9
IO
II
12
II
I;
f (h)
values, spring or minimum
tide of 25 April 1979). Points of computed Q/A values.
a,
For the neap tide of 2 May 1979 and corresponding fresh water flows (6 m3 s- ‘) for Eastern and Western Cleddau rivers K, = 10, K, = 26 and m = 2. It was noted that while the variation of fresh water flows within the ranges observed for 1977 and 1979 surveys (1m3s-’ - 6 m3 s- ‘) have virtually no effect on the hydrodynamics of the estuary, salt distributions in the uppermost parts of the Eastern and Western Cleddau rivers are affected by such changes in the input. In practice for both spring and neap tides the insertion of lower fresh water discharges into the model requires an increase in the selected values of K, for the upper reaches. The optimal values for K, are 35 and 30 for spring and neap tides of 1977, respectively. The computed dispersion coefficients for the spring tide of 25 April 1979 for three different locations are given in Figures 4(a), (b) and (c). At seaward locations [Figures 4(a)
Mathematical
model
of Upper
Miljord
40
Haven
415
40 (b)
30
30
20
20
IO
IO
zE 0 2 .(I 0
0 Distance
from Carr Jetty
(km)
(Daucleddau
cn 40
and West
40
Cleddau)
(d 1
30
30
20
20
IO
0
IO
5
IO
15 Distance
20 from Carr Jetty
0 (km)
(Daucleddau
and East Cleddau)
Figure 5. Simulated and observed salinities [m observed data, Williams and Jolly; A observed data, West; low tide, 1.5 h after low tide, - 3 h after low tide, - - - - high tide; (a) and (c), spring tide; (b), (d) neap tide].
and (b)] the variation in Eduring tidal cycles always correspond to the variation of (Q/A) at these locations. However, at Slebech [Figure 4(c)], this correspondence is not observed reflecting the influence of the considerable longitudinal density gradient which exists at this location during a tidal cycle. Computed salt concentrations throughout the Daucleddau, Western Cleddau and Eastern Cleddau regions at low, high and two intermediate levels for the spring and neap tides are shown in Figures 5(a)-(d). The comparison between computed salt concentrations and available data indicates that the model is capable of simulating salt intrusion in the Upper Milford Haven estuary. These simulated salt concentrations are in general in phase with the predicted water surface elevations at the seaward locations in the estuary
416
V. Nassehi
(a)
Location
& D. J. A. William
1 Picton
Pod
(b)
Locatlon
I Slebech
3
4tl 0 I
2
345
6
7691011
0
12 Time
Figure 6. Spring tide of 25 April simulated salinities).
I 2
3 4
5 6 7 8 9 1011
12
(h)
1979 (-
predicted
water
surface
elevations,
----
[for typical case see Figure 6(a)]. The phase difference between these variables at the tidal limits of Eastern and Western Cleddau rivers [for typical case see Figure 6(b)] is probably due on the flood, to the combined effects of the distortion of the tidal current and the ponding of the fresh water inflow and on the ebb, the delayed discharge of salt water remaining in various branches and pools into the Eastern and Western Cleddau’s main flow channels. Summary The present combined model provides a reasonable prediction of tidal dynamics and simulation of salt intrusion in Upper Milford Haven which is a branching one-dimensional tidal water system. The use of a finite difference scheme for the solution of the hydrodynamical equations in conjunction with the simultation of salt transport by a finite element technique which is readily extended to branching domains makes the model suitable for the practical cases of one-dimensional tidal water networks. Acknowledgements The authors gratefully acknowlege the assistance of members of the Department of Chemical Engineering, University College of Swansea, and the Welsh Water Authority, during the 1979 survey of the Upper Milford Haven estuary. Notation x 6
= distance from mouth of estuary (space variable), = time (time variable),
-.-
Ax At h
Q B A H n g S P 4 hc VC ii ij
I, II, III
Mathematic
model
of Upper
Milford
Haven
417
= spatial increment, = time increment, = water surface elevation with respect to a datum, = volumetric discharge, = breadth of water surface averaged over a segment, = cross-sectional area, = hydraulic depth, = Manning’s friction coefficient, = acceleration due to gravity, = depth mean longitudinal salinity (ppt), = depth mean water density, = surface area of a triangular junction, = water surface elevation at a triangular junction, = volume of junction with surface area A,, = one-dimensional dispersion coefficient, = time-stepping weight factor e.g. S{‘” = OS{’ i + (1 -@S{ = indices relative to points in the space-time plane = indices relative to the sidesof a triangular junction area.
References A.S.C.E. Task Committee 1982 Modelling hydraulic phenomena: a glossary of terms. A.S.C.E. Journal o.f Hydraulics Division, 108 (HY7), 845-853. Chatwin, P. C. &Allen, C. M. 1985 Mathematical models of dispersion in rivers and estuaries. Annual Rwiew of Fluid Mechanics 17,119-149. Cunge, J. A., Holly, F. M., Jr & Verwey, A. 1980 Practical aspects of Computational River Hydraulics. Pitman Publications, London. Dronkers, J. J. 1964 Tidal Computations in Rivers and Coastal Waters. North Holland Publishing Company, Amsterdam. Elliott, A. J. 1976 A numerical model of the internal circulation in a branching estuary. Special Report 54, Chesapeake Bay Institute, The John Hopkins University, Baltimore, MD. Graf, W. H. 1971 Hydraulics of sediment transport. McGraw Hill, New York. Gray, W. 1980 Do finite element models simulate surface flow? Finite Elements in Water Resources, Proceedings, Third International Conference, section 1, p. 122. S. Y. Wang et al., editors; Rose Printing Co. Inc., Tallahassee, U.S.A. Gunn, D. J. & Yenigun, 0.1985 Modelling oftidal motion in shoaling waters; the Estuary of Milford Haven, Estuarine Coastal and Shelf Science 21,337-356. Holly, F. M., Jr & Usseglio-Polatera, J. M. 1984 Dispersion simulation in two-dimensional tidal flow, A.S.C.E.,Journal of Hydraulic Engineering 110 (7), 905-927. Liggett, J. A. & Cunge, J. A. 1975 Numerical methods of solution of the unsteady flow equations. In Unsteady Flow in Open Channels (Mahmood, K. & Yevjevich, V., eds). Water Resources Publications, Fort Collins, Colorado, U.S.A. McDowell, D. M. 8s O’Connor, B. A. 1977 Hydraulic Behaviour of Estuaries, The McMillan Press Ltd. Morton, K. W. & Parrott, A. K. 1980 Generalised Galerkin methods for first order hyperbolic equations, Journal of Computational Physics 36, (2), 249-270. Nelson-Smith, A. 1965 Marine Biology of Milford Haven: The Physical Environment, Field Studies, 2 155-188. Nelson-Smith, A. & Case, R. M. 1984 Aspects of the shore and sublittoral ecology of the Daucleddau estuary (Milford Haven). ZoologicalJournal of Linnean Society 80, 177-190. O’Connor, B. A. &Thompson, G. 1976 A mathematical model of chloride levels in the Wear estuary (U.K ) International Symposium on Unsteady Flow in Open Channels, Paper Hl, BHRA Fluid Engineering, Cranfield, U.K. Smith, I. M., Faraday, R. V. &O’Connor, B. A. 1973 Rayleigh-Ritz and Galerkin finite elements for diffusion convection problems. Water Resources Research (9), 593406. Smith, R. 1980 Buoyancy effects upon longitudinal dispersion in wide well-mixed estuaries. Philosophical Transactions of zhe Royal Society of London, Section A, 296.
418
V. Nassehi
& D. 3. A. Williams
Thatcher, M. L. & Harleman, D. R. F. 1972 A mathematical model for the prediction of unsteady salinity intrusion in estuaries. Ref. 144, Department of Civil Engineering, M.I.T. Cambridge, MA, U.S.A. Weaver, T. J. 1976 Instability in tidal flow computational schemes. A.S.C.E.,J ournal of Hydraulics Division 102 (HY5), 569-580. West, J. R. 1978 A report on the hydrodynamics and water quality of the Upper Milford Haven, vols 1 and 2. Department of Civil Engineering, University of Birmingham. West, J. R. 1979 A report on the hydrodynamics and water quality of the Upper Milford Haven, vol. 3. Department of Civil Engineering, University of Birmingham. West, J. R. & Lin, C. H. 1983 An evaluation of a moving coordinate system model of salinity instrusion into the Mersey Estuary. Water Pollution Control, 402-417. Williams, B. R. H. 1971 The Milford Haven. Brixham Lab. Report, BL/C/1357, ICI Ltd. Williams, D. J. A. &Jolly, R. 1975 Survey of pollution and mathematical modelling of pollutant distribution and transport. Report by Department of Chemical Engineering, University College of Swansea. Williams, D. J. A. & Nassehi, V. 1980a Mathematical tidal model of the Tay Estuary, Proceedings of the Royal Society of Edinburgh, Section B 78,171-182. Williams, D. J. A. & Nassehi, V. 1980b A one-dimensional mathematical tidal model of Upper Milford Haven Estuary. Report to Welsh Water Authority, Department of Chemical Engineering, University College of Swansea. Zi-Cai, L., Ling-Jia, Z. & Hui-Li, W. 1983 Difference methods of flow in branch channel A.S.C.E.,Journal of Hydraulic Engineering, 109 (3), 424-447. Zienkiewicz, 0. C. 1977 Finite element method, 3rd ed. McGraw-Hill, New York.
Estuarine, Coastal and Shelf Science (1986) 23,403-418
Mathematical Model Haven-A Branching
V. Nassehi
of Upper Estuary
Milford
and D. J. A. Williams
Department of Chemical Engineering, Park, Swansea SA2 8PP, Wales
University College of Swansea, Singleton
Received 25 March 1985 and in revisedform 18 November 1985
Keywords: estuaries; tidal heights; salt intrusion; Haven estuary.
numerical
model; Milford
A mathematical model for hydrodynamics and salt intrusion in a onedimensional branching estuary has been developed. The mode1 satisfactorily predicts wafer levels and discharges throughout Upper Milford Haven estuary for spring and neap tides. There is also a measure of agreement between computed and observed salinities.
Introduction and other economic advantageshave resulted in the selection of the sheltered waters of Milford Haven for the location of industrial activities such asrefining petroleum and power generation. This growing industrial and commercial demand on tidal waters of the Haven requires careful management of water quality to avoid catastrophic effects on their ecology (Nelson-Smith & Case, 1984). Mathematical models which predict estuarine flow variables (water velocity and water surface elevation) and provide estimates of dispersion coefficients are useful tools for water quality management in this estuary. Although in general such models should recognisethe three-dimensional nature of flow in natural water systems, it is often possible from a consideration of the geometrical boundaries of an estuary to reduce the complexity and usesimplified models in either two or one dimensions (Dronkers, 1964). The outer part of this estuary (from Carr Jetty to St. Ann’s Head, Figure 1) has been modelled in fwo dimensions by Gunn and Yenigun (1985). In the present work we describe the development and use of a one-dimensional model for upper Milford Haven estuary.
Accessibility
Milford
Haven estuary-A
physical description
The physical environment of Milford Haven hasbeen described in somedetail by NelsonSmith (1965). The estuary is located in West Wales (U.K.) and it is of the drowned valley type. It comprises a network of channels connecting the inland waters of the Eastern and Western Cleddau and several smaller rivers to the sea(Figure 1). The estuary encloses 110 km of varied coastline and is 25 km wide at its mouth. 403 0272-7714/86/090403+
16 $03.00/O
0 1986
Academic
Press
Inc.
(London)
Limited
404
V. Nassehi arid D. J. A. Wiliiums
(Eastern
(Daucleddau)
Figure 1. Locations: a-St. Ann’s Head; &Port of Milford Haven; c-Carr Jetty; d-Hobbs Point; e-Pembroke Dock; f-Burton; g-Neyland; h-Lawrenny; i-Port Lion; j-Picton Point; k-East Hook; I-Slebech; m-Minwear Wood; n-Canaston Weir; o-Haverfordwest; pCrowhil1 Weir (0 water level observation station-1979 survey; 0 water level observation-1977 survey).
Beyond Carr Jetty and above the road bridge between Pembroke Dock and Neyland (Figure l), the narrowing main channel of the estuary (from this point landwards called Daucleddau) turns northwards until its limit at the junction of Eastern and Western Cleddau rivers at Picton Point, about 27 km from the estuary mouth. Within this upper part of the estuary water flow is predominantly one-dimensional. The mouths of other smaller rivers often form inlets which initially are of a similar width to the main channel, but usually are much shallower and generally their bedsare almost completely exposed at low water (West, 1978). These shallow inlets are thus thought to have negligible effect on the longitudinal momentum transfer. The tidal limit of spring tides for the Western Cleddau river is about 11 km [in the region of Crowhill Weir with a range of 3-9 m above ordnance datum (A.O.D.)] and about 8 km for the Eastern Cleddau river (in the region of Canaston Weir with a range of 4.2 m
Mathematical
model
of Upper
Miijord
405
Haven
A.O.D.) from Picton Point, respectively. The mean tidal ranges at Milford Haven are given as 6.3 m and 2.7 m for spring and neaps respectively (Admiralty Chart no. 3275). For both spring and neap tides there is a delay (max. 36 min) in the occurrences of high and low water between the port of Milford Haven and locations near to the tidal limits of the Western and Eastern Cleddau rivers (Nelson-Smith, 1965). The bed level of the Daucleddau varies from 15 m below ordnance datum (B.O.D.) near Lawrenny to about 5 m B.O.D. near Picton Point. The width of the main channel is of the order of 500 m. The Eastern and Western Cleddau rivers gradually decrease in width to a few meters towards their tidal limits. Except for a few pools there is very little water in these rivers at low water springs. Details of sediment distribution in the Daucleddau is reported by Williams (1971) and within the Cleddau rivers by West (1979). In the Cleddau rivers water mainly flows over a muddy bottom while in the Daucleddau reaches the bottom surface consists of sand and rock. During the 1977 survey of the upper estuary (West, 1978) the fresh water discharges of the Cleddau rivers were of the order of 1 m3 s- ‘. Another survey in spring of 1979 (Williams & Nassehi, 1980b) gives the respective mean daily fresh water flows of Western and Eastern Cleddau rivers as 6 m3 s- ’ and 5.5 m3 s- i. These flows are insignificant in comparison with the tidal flow in this estuary which is of the order of lo4 m3 s- ’ at Hobbs Point at mid-tide (Gunn & Yenigun, 1985). Although some vertical stratification has been reported in the uppermost reaches of the Milford Haven estuary (e.g. at Haverfordwest in the Western Cleddau and Minwear Wood in the Eastern Cleddau, Figure I), (Nelson-Smith, 1965) it is expected that seaward of these locations where the estuary channel is wider and deeper, vertical mixing will be almost complete. Recorded vertical salinity changes within the outer Daucleddau are less than 1%0on spring tides and less than 5%~ on neap tides (Gunn & Yenigun, 1985). For the purpose of over-all simulations of salt intrusion in the upper estuary by a one.dimensional model, localised vertical stratification effects near to the tidal limits of the Cleddau rivers have been ignored. Mathematical
model
The mathematical model consists of the numerical integration of the following governing equations (McDowell & O’Connor, 1977): (a) The one-dimensional hydrodynamical equations of continuity and motion (1)
-+---+--[- 1+- +H'p=o, dh dX
(b)
the one-dimensional
1
gA
aQ at
1 2
gAdx
Q2
lQlQn2
2y ax
pff4:3
A
salt transport
equation
as QdS t+--z-AG
i a
AtE [
=O. 8X
1
The variables are defined in the section entitled Notation.
(2)
(3)
406
V. Nassehi
Using the following
and D. 7. A. Williams
approximate
relationship p=o.75
equation (3) can be written
ss
1000,
(4)
as:
ah 1 aQ -+--+’ ax gA at
gA
+ --a Q2 +IQlQn’ ~2~413
[ 1 dx A
0*75H
as
2000t1.5 SG =O
(5)
whence equations (l), (3) and (5) can be solved in a combined manner. In the present model for reasons of convenience and mathematical rigour the nonlinear hyperbolic partial differential equations representing the hydrodynamical conditions in the estuary [eqns (1) and (5)] are solved by a finite difference technique (Gray, 1980) while the linear parabolic partial differential equation of the salt transport [eqn (3)] is solved by a Galerkin finite element method. After division of the solution domain into different segments according to the four point implicit finite difference scheme used herein, a variable f and its temporal and spatial derivatives can be approximated as (Liggett and Cunge, 1975): (6)
1
A!f(XJ)C+-& fi:: ++, +Sj+l-fi ’ i
This finite difference scheme is dissipative (ASCE, Task Committee, 1982) thereby damping out numerical dispersion; it is independent of variability of segment length and allows an extremely simple representation of boundary conditions (Cunge et al., 1980). Conclusions concerning stability and accuracy of the linearized hydrodynamical model based on this scheme are given by Liggett and Cunge (1975) as: (a) the scheme is unconditionally stable when 0.5 < f3< 1, (b) some parasitic oscillations appear for values of 0 near to 0.5 which disappear if 0 is chosen to be greater than 0.66 or the frictional resistance is enhanced, (c) the scheme is of second order accuracy when 8 = 0.5 and of first order accuracy for 0.5<8< 1. These conclusions can in general be extended to the cases of non-linear hydrodynamical models (Weaver, 1976). Application of the four-point implicit scheme to the hydrodynamical equations (1) and (5) results in the following two working equations: Q{;;-Qi’
hj:;
1
-hi+
+G1(hj;:+hj+‘)-G2=0,
’ + G,(Qi;;
+ Qi”)
These algebraic equations relate the prime unknowns
- G,=O.
(10)
of h and Q at the boundaries of every
Mathematical
model
sf Upper
Milford
Haven
407
segment. Coefficients G,, G,, G, and G, in equations (9) and (10) are relations in terms of space and time intervals, time weighting factor, hydrographic data, bottom friction, salinity and the values of the prime unknowns at previous and present (due to nonlinearity) time levels (e.g. see Williams & Nassehi, 1980a, for the case of 0= 1). Based on the working equations (9) and (10) a set of algebraic equations describing the hydrodynamical conditions along each branch of the Upper Milford Haven estuary is obtained. Added to this system of equations are the following hydrodynamical relations at the junction of the three branches: (a) the continuity condition is expressed as (Elliott, 1976): algebraic sum of the discharges across the boundaries of the junction
=
accumulation within the junction.
or
(11)
which after time-stepping
gives
Ce=A,h:+l - hi ; At (b) the momentum condition at the junction of relatively of equal free surface elevation at its boundaries
h,=h,,=h,,,=h,.
(12)
small rivers is the assumption
(13)
For larger junctions this assumption should be modified (Zi-Cai et al., 1983). The third working equation in the present model is derived by the application of the Galerkin finite element technique (Zienkiewicz, 1977) to the salt transport equation [eqn, (3)]. This equation includes two physical mixing mechanisms: convective longitudinal transport and longitudinal diffusion. These mechanisms are represented by the second and third terms, respectively in equation (3). It is known that in cases where there is strong convection [i.e. when the second term in equation (3) is dominant] the numerical solution of this equation by finite difference or finite element methods may become oscillatory and useless (e.g. see Holly & Usseglio-Polatera, 1984). To avoid this problem specially modified finite element solution schemes applicable to convective diffusion equations are widely used (e.g. see Morton & Parrot, 1980). However, in the present work preliminary results obtained using a standard Galerkin finite element solution of equation (3) did not show any parasitic oscillations and the third working equation of this model is based on the standard technique. Using linear Lagrangian shape functions and isoparametric elements (Zienkiewicz,
408
V. Nassehi and D. J. A. Williams
1977) for every element along the single channel of each branch, the Galerkin finite element form of equation (3) is
1/
-
\
dS,
I-
i-1
1
si
AAX +
6 si+1 (1
iI I1
1 -1 AE
1which after time stepping gives: 8Q
--y+3at
RAE A LIX Ax
eQ
si
=o
AX
Ah
7
BAE-
(14)
si+l
S;+’
-+1-ax 6Lt
=
-
ALIX -----6nt
OQ
Anx
l-8
e...v++.---Q--AE
2
Ax
3Llt
2
Aax
l-8
-+-
6nt
BAE AAx
2
! OQ ( eAE
3nt
1-Q AX
l-f? Q+nxAF
2
ALlX
---
nx
s.f+’ 1+1
1
l-8
l-8
Q+nxAe
6nt
2
ALiX
l-f?
--3nt
-Q-.xAE 2
-
l-8
The accuracy and stability of this scheme in relation to numerical diffusion and dispersion is presented by Smith et al. (1973). For practically large enough At values, 8 > 05 should be used (herein 0 = 0.67). This time step may be appreciably increased if Q/A (velocities) values are small and the physical dispersion coefficients (E) are large. Overall accuracy of the scheme depends on the total number of nodes (degree of mesh refinement) and is not affected by the use of second or third order elements instead of linear ones used in the present work. Based on the working equation (15) elemental matrix equations for all of the elements along the estuary branches are derived. In addition to this system is the extra equation arising from the salt continuity condition at the junction of the three branches. The necessarysalt continuity relation at the junction is: total amount of salt entering the junction
=
total amount of salt leaving the junction
+
accumulation of salt within the junction,
--
Mathematical
model of Upper Miljord
409
Haven
or (16)
CQS= ; jSc d’V,. If we assumethat mixing at the junction is complete then:
EQS=A~~
(17)
S,h,,
wherein
s, = s,,= s,,,= 41 + Sill 2
(18)
.
Using relation (18) and after time stepping we have: Q, j+l
sj+l,
-
Q:I+l
S;,+‘-Q;,+’
;4 (Sl, + S&j+ ’ hi+At’ -(S,,
S,:,+l
+ &)j hJ;
=
(19)
Let us write 2nt
j+l
al=yQl
9 c
2At a,=rQ,,
j+l
+e+’
(20)
c 2At
.+l
a3= A
a,,,
+hj+‘,
‘2 and a4
= + h,’ (S,, + S,,,)j.
Then after somealgebraic operations and multiplying the sidesby a large scalarparameter Aequation (19) gives: &a,
S{+‘-a,
S&+’
-a3
S{i’
+a,)=O.
(21)
The matrix equations basedon (15) are assembledtogether with those arising from application of (21) at the three boundaries of the junctions. The parameter A is chosen to be large (109) ensuring that in the finally assembledmatrix the coefficients applicable to salinities at the junction boundaries (i.e. S{” ‘, S{,+’ and S{,: ‘) are totally dominated by those of (21). In this way the constraint imposed by equation (19) is satisfied.
410
I,‘. Nassrhi
and D.
Application
3. A. Williuwe
of the model to the Upper Milford
Haven estuary
The division of the Daucleddau region and the Western and Eastern Cleddau rivers into segments (or elements) of variable length is shown in Figure 1. The shaded triangle contained by segment boundaries (1 l), (12) and (23) at Picton Point (Figure 1) represents the confluence of the Cleddau rivers. The solution algorithm assumesfirstly that at any time level the salinities along the estuary are given, therefore after insertion of the boundary conditions the system of algebraic equations representing the hydrodynamical relations becomesdeterminate and is solved using an iterative procedure (Williams & Nassehi, 1980~). The discharges so obtained are used to solve the salt transport equation (subject to salinity boundary conditions) at the next time level. These stepsare then repeated for a whole tidal cycle. Since the density term in the hydrodynamical equation of motion is relatively small the present algorithm which computes h and Q one n tout of phasewith S produces little error in their values. To obtain optimal accuracy without any danger of instability the time weighting factor B is chosen to be 0.67. The time increment (a t) is 250 seconds.Herein smaller (A t) values increase computational time without improving accuracy. For a typical spring tide of 125 h duration the present combined model converges after two tidal cycles. This required 270 secondsC.P.U. time allocation on an ICL 2966 computer. This time requirement is short enough to justify the useof the non-linear version of the hydrodynamical equation of motion [eqn. (5)] for any possible enhancement of the accuracy in the numerical solution. The model requires the following input data. (a) Hydrogruphic data. These data are values of cross-sectional area, hydraulic depth and water surface breadth at each segment boundary (or node) along the estuary. Herein values of A, H and B are arranged in double entry tables as functions of water surface elevation and location. The first row in these tables contains hydrographic data for each location corresponding to low water and the successiverows are in 0.2 m intervals. All of the geometrical data in the present work are with respect to O.D. Newlyn. Water surface breadth and low-water cross-sectional area landwards of Lawrenny are calculated utilizing echo sounding data of the 1977 survey of the Upper Milford Haven estuary (West, 1978). Corresponding data for the reach between Carr Jetty and Lawrenny (Figure 1) were found using Admiralty chart no. 3275. (b) Manning’s friction coefficient. Values of the friction coefficient ‘n’ were initially selected with the help of literature (Graf, 1971). Then using a trial and error procedure they were adjusted in somesections, in order to minimize the discrepancy between model output and field data observed during the spring tide of 25 April 1979 for the locations indicated in Figure 1. The finally accepted values of ‘ n ’ in the present work are 0.018 along the Daucleddau region and 0.02 along the reaches in the Western and Eastern Cleddau rivers. Boundary conditions. These data are water surface elevations and salinities observed at estuary mouth and fresh water flows and zero salinities at estuary head. (c)
(d)
Initial conditions.
At time t = 0 any set of arbitrary values of h, Q and S can be used,
Mathematical
model of Upper Milford
Haven
411
asthe initial conditions to start the solution algorithm. In practice, approximate low water values for h and Q, together with zero salinities are used. coeficients. Inherent in the one-dimensional approach to the modelling of unsteady state distribution of salt in an estuary is the difficulty of predicting dispersion coefficients. Although some progress has been made in theoretical analyses of mixing mechanismsin estuaries the relationships derived are complex (Smith, 1980; Chatwin & Allen, 1985) and cannot be readily used in practical cases.The longitudinal dispersion coefficient is best determined by field observation in conjunction with the direct useof the salt transport equation [eqn (3)]; however, the quantity of field data required for such a calculation is considerable. In order to avoid this problem someinvestigators have chosen a constant value for this coefficient throughout flow channel during the tidal cycle (West & Lin, 1983). In the present model a method based on an empirical formula suggestedby Thatcher and Harleman (1972) is used. This empirical relation in metric system of units is written as: (e) Dispersion
&=64K1
n(Q/A)
H5”6 + K2
(2%
values of K,, K, and m in this equation depend upon various physical factors and must be determined separately for each estuary and flow condition. In equation (22) the first term with K, = 1 is an analytical expression which allows for dispersion by velocity gradients over the flow depths (O’Connor & Thompson, 1976). In real flow systems, transverse velocity gradients, junctions and bends together with wind action produce extra dispersion and this is reflected by an increase in the value of K,. For example, Thatcher and Harleman (1972) used K, = 3 on tidal rivers on the eastcoast of U.S.A. but other investigators (e.g. seeMcDowell & O’Connor, 1977) have suggestedan increased (sometimesas much asten times) values for this constant in narrow and tortuous conditions. The second term of equation (22) is an empirical term which provides a contribution to Edue to the presenceof longitudinal density gradients. Results and discussion Following the calibration of the model using the spring tide of 25 April 1979 asthe input tide, the model was successfully validated by the insertion of water surface elevations at Carr Jetty observed during the neap tide of 2 May 1979 whilst the values of ‘ n ’ remained unchanged. The model was alsoused to predict the spring tide of 2 August 1977 with the samevalues of ‘ n ’ throughout the estuary. The comparisons between the observed water levels and computed results for these tides are given in Figures 2(a)-(d). These comparisonsshow a very close agreement between the observed data and computed water levels. These curves indicate that the incoming tidal wave maintains an almost symmetric shape throughout the Daucleddau region. However, as it propagates into the branches it is increasingly distorted. The maximum time lags between high water at Carr Jetty and the tidal limits of Cleddau rivers is about 20 min. Computed dischargesfor the spring tide of 25 April 1979 for various locations together with the corresponding depths mean velocities are given in Figures 3(a)-(d). Whenever available these figures include instantaneous velocities recorded by West during 1977 survey of the estuary (West, 1978) and the measureof agreement between these values is shown.
412
V. Nassehi
and D. J. A. Williams
7
4
t
(0 ) Locotlon
I Port
4
Lion
(b)
Location
1
Plcton
Point 0,
4 t 30 30
// 0
I
0
5
‘iYI-u.Ll 0
I2
3456789101
0
(c ) Location
I I I I I23456789101112
I
I
I
I
I
I
I
I.
I23456789101112
h : Lawrenny
I
I
I
I
I
I
I
1:
I
0
Time
I
2
3
4
5
6
7
8
9
IO
II
I2
(h )
Figure 2. Comparison of computed and observed water levels observed during the spring tide of 25 April 1979, H the spring tide of 2. August 1977, 0 water levels observed 1979; the time of origin is 7.30 a.m. B.S.T.; model -- - - model output for the neap tide).
surface elevations (0 water water levels observed during during the neap tide of 2 May output for the spring tides,
Mathematical
model of Upper Miljord
Haven
413
I : 0 ) LocatIon
3000
I : Port
LlOrl
t
f
( c ) Locotlon
k
East
Hook I-O
0.5
0.0
1 -0/
-0.5
1 t I
I
I
I
I
I
I
1
I
I
I
I
I
4.0
f
Y) E La c :: ?
( b ) Locotlon
1 Plcton
pomt
t
( d 1 Locotlon
I.0
1 Slebech
0.5
o-o
\ -0.5
3000
t
,
,
0
I
23456
,
,
,
,
]
]
,
,
,
,
7
8
9
IO
II
12
bl 0 Time
Figure 3. Computed discharges, 25 April velocities, 2 August
discharges 1979; ---1977).
(h
I
I 2
I 3
I 4
I 5
I 6
I 7
I 8
I 9
I IO
I Ii
I 12
)
and velocities (time of origin is 7.30 a.m. B.S.T.: ~ velocities, 25 April 1979; 0 observed instantaneous
The constants (i.e. K,, K, and m) in equation (22) were adjusted to find appropriate 1: values. This involves a trial and error procedure in which model output is compared with the observed salinity data for upper Milford Haven reported for tides of similar amplitude by Williams and Jolly (1975) and West (1978). The finally accepted values for these constants under various conditions are: For the spring tide of 25 April 1979 and corresponding fresh water flows (6 m3 s- ‘) for Eastern and Western Cleddau rivers, K, = 10, K, = 30 and m = 2.
-1.0
414
I’. Nassehi
and D. J. A. William
I
‘**
-
150 100 50
-
/
Locot~on a ,-. \
o
Corr
c
Jetty /
‘\
/ \
\
/
\
/
\
-
-1
./
\
,’
\
\
I
/
;
\d/
\J/
IrllrlllllYI 12345670
0
-E
/
(a)
9
IO
II
12
I
0
0
1
123456789
I
2
IO
3
4
5
6 Time
Figure 4. Dispersion coefficients (----E b, etc. correspond to the time of maximum
7
8
9
IO
II
12
II
I;
f (h)
values, spring or minimum
tide of 25 April 1979). Points of computed Q/A values.
a,
For the neap tide of 2 May 1979 and corresponding fresh water flows (6 m3 s- ‘) for Eastern and Western Cleddau rivers K, = 10, K, = 26 and m = 2. It was noted that while the variation of fresh water flows within the ranges observed for 1977 and 1979 surveys (1m3s-’ - 6 m3 s- ‘) have virtually no effect on the hydrodynamics of the estuary, salt distributions in the uppermost parts of the Eastern and Western Cleddau rivers are affected by such changes in the input. In practice for both spring and neap tides the insertion of lower fresh water discharges into the model requires an increase in the selected values of K, for the upper reaches. The optimal values for K, are 35 and 30 for spring and neap tides of 1977, respectively. The computed dispersion coefficients for the spring tide of 25 April 1979 for three different locations are given in Figures 4(a), (b) and (c). At seaward locations [Figures 4(a)
Mathematical
model
of Upper
Miljord
40
Haven
415
40 (b)
30
30
20
20
IO
IO
zE 0 2 .(I 0
0 Distance
from Carr Jetty
(km)
(Daucleddau
cn 40
and West
40
Cleddau)
(d 1
30
30
20
20
IO
0
IO
5
IO
15 Distance
20 from Carr Jetty
0 (km)
(Daucleddau
and East Cleddau)
Figure 5. Simulated and observed salinities [m observed data, Williams and Jolly; A observed data, West; low tide, 1.5 h after low tide, - 3 h after low tide, - - - - high tide; (a) and (c), spring tide; (b), (d) neap tide].
and (b)] the variation in Eduring tidal cycles always correspond to the variation of (Q/A) at these locations. However, at Slebech [Figure 4(c)], this correspondence is not observed reflecting the influence of the considerable longitudinal density gradient which exists at this location during a tidal cycle. Computed salt concentrations throughout the Daucleddau, Western Cleddau and Eastern Cleddau regions at low, high and two intermediate levels for the spring and neap tides are shown in Figures 5(a)-(d). The comparison between computed salt concentrations and available data indicates that the model is capable of simulating salt intrusion in the Upper Milford Haven estuary. These simulated salt concentrations are in general in phase with the predicted water surface elevations at the seaward locations in the estuary
416
V. Nassehi
(a)
Location
& D. J. A. William
1 Picton
Pod
(b)
Locatlon
I Slebech
3
4tl 0 I
2
345
6
7691011
0
12 Time
Figure 6. Spring tide of 25 April simulated salinities).
I 2
3 4
5 6 7 8 9 1011
12
(h)
1979 (-
predicted
water
surface
elevations,
----
[for typical case see Figure 6(a)]. The phase difference between these variables at the tidal limits of Eastern and Western Cleddau rivers [for typical case see Figure 6(b)] is probably due on the flood, to the combined effects of the distortion of the tidal current and the ponding of the fresh water inflow and on the ebb, the delayed discharge of salt water remaining in various branches and pools into the Eastern and Western Cleddau’s main flow channels. Summary The present combined model provides a reasonable prediction of tidal dynamics and simulation of salt intrusion in Upper Milford Haven which is a branching one-dimensional tidal water system. The use of a finite difference scheme for the solution of the hydrodynamical equations in conjunction with the simultation of salt transport by a finite element technique which is readily extended to branching domains makes the model suitable for the practical cases of one-dimensional tidal water networks. Acknowledgements The authors gratefully acknowlege the assistance of members of the Department of Chemical Engineering, University College of Swansea, and the Welsh Water Authority, during the 1979 survey of the Upper Milford Haven estuary. Notation x 6
= distance from mouth of estuary (space variable), = time (time variable),
-.-
Ax At h
Q B A H n g S P 4 hc VC ii ij
I, II, III
Mathematic
model
of Upper
Milford
Haven
417
= spatial increment, = time increment, = water surface elevation with respect to a datum, = volumetric discharge, = breadth of water surface averaged over a segment, = cross-sectional area, = hydraulic depth, = Manning’s friction coefficient, = acceleration due to gravity, = depth mean longitudinal salinity (ppt), = depth mean water density, = surface area of a triangular junction, = water surface elevation at a triangular junction, = volume of junction with surface area A,, = one-dimensional dispersion coefficient, = time-stepping weight factor e.g. S{‘” = OS{’ i + (1 -@S{ = indices relative to points in the space-time plane = indices relative to the sidesof a triangular junction area.
References A.S.C.E. Task Committee 1982 Modelling hydraulic phenomena: a glossary of terms. A.S.C.E. Journal o.f Hydraulics Division, 108 (HY7), 845-853. Chatwin, P. C. &Allen, C. M. 1985 Mathematical models of dispersion in rivers and estuaries. Annual Rwiew of Fluid Mechanics 17,119-149. Cunge, J. A., Holly, F. M., Jr & Verwey, A. 1980 Practical aspects of Computational River Hydraulics. Pitman Publications, London. Dronkers, J. J. 1964 Tidal Computations in Rivers and Coastal Waters. North Holland Publishing Company, Amsterdam. Elliott, A. J. 1976 A numerical model of the internal circulation in a branching estuary. Special Report 54, Chesapeake Bay Institute, The John Hopkins University, Baltimore, MD. Graf, W. H. 1971 Hydraulics of sediment transport. McGraw Hill, New York. Gray, W. 1980 Do finite element models simulate surface flow? Finite Elements in Water Resources, Proceedings, Third International Conference, section 1, p. 122. S. Y. Wang et al., editors; Rose Printing Co. Inc., Tallahassee, U.S.A. Gunn, D. J. & Yenigun, 0.1985 Modelling oftidal motion in shoaling waters; the Estuary of Milford Haven, Estuarine Coastal and Shelf Science 21,337-356. Holly, F. M., Jr & Usseglio-Polatera, J. M. 1984 Dispersion simulation in two-dimensional tidal flow, A.S.C.E.,Journal of Hydraulic Engineering 110 (7), 905-927. Liggett, J. A. & Cunge, J. A. 1975 Numerical methods of solution of the unsteady flow equations. In Unsteady Flow in Open Channels (Mahmood, K. & Yevjevich, V., eds). Water Resources Publications, Fort Collins, Colorado, U.S.A. McDowell, D. M. 8s O’Connor, B. A. 1977 Hydraulic Behaviour of Estuaries, The McMillan Press Ltd. Morton, K. W. & Parrott, A. K. 1980 Generalised Galerkin methods for first order hyperbolic equations, Journal of Computational Physics 36, (2), 249-270. Nelson-Smith, A. 1965 Marine Biology of Milford Haven: The Physical Environment, Field Studies, 2 155-188. Nelson-Smith, A. & Case, R. M. 1984 Aspects of the shore and sublittoral ecology of the Daucleddau estuary (Milford Haven). ZoologicalJournal of Linnean Society 80, 177-190. O’Connor, B. A. &Thompson, G. 1976 A mathematical model of chloride levels in the Wear estuary (U.K ) International Symposium on Unsteady Flow in Open Channels, Paper Hl, BHRA Fluid Engineering, Cranfield, U.K. Smith, I. M., Faraday, R. V. &O’Connor, B. A. 1973 Rayleigh-Ritz and Galerkin finite elements for diffusion convection problems. Water Resources Research (9), 593406. Smith, R. 1980 Buoyancy effects upon longitudinal dispersion in wide well-mixed estuaries. Philosophical Transactions of zhe Royal Society of London, Section A, 296.
418
V. Nassehi
& D. 3. A. Williams
Thatcher, M. L. & Harleman, D. R. F. 1972 A mathematical model for the prediction of unsteady salinity intrusion in estuaries. Ref. 144, Department of Civil Engineering, M.I.T. Cambridge, MA, U.S.A. Weaver, T. J. 1976 Instability in tidal flow computational schemes. A.S.C.E.,J ournal of Hydraulics Division 102 (HY5), 569-580. West, J. R. 1978 A report on the hydrodynamics and water quality of the Upper Milford Haven, vols 1 and 2. Department of Civil Engineering, University of Birmingham. West, J. R. 1979 A report on the hydrodynamics and water quality of the Upper Milford Haven, vol. 3. Department of Civil Engineering, University of Birmingham. West, J. R. & Lin, C. H. 1983 An evaluation of a moving coordinate system model of salinity instrusion into the Mersey Estuary. Water Pollution Control, 402-417. Williams, B. R. H. 1971 The Milford Haven. Brixham Lab. Report, BL/C/1357, ICI Ltd. Williams, D. J. A. &Jolly, R. 1975 Survey of pollution and mathematical modelling of pollutant distribution and transport. Report by Department of Chemical Engineering, University College of Swansea. Williams, D. J. A. & Nassehi, V. 1980a Mathematical tidal model of the Tay Estuary, Proceedings of the Royal Society of Edinburgh, Section B 78,171-182. Williams, D. J. A. & Nassehi, V. 1980b A one-dimensional mathematical tidal model of Upper Milford Haven Estuary. Report to Welsh Water Authority, Department of Chemical Engineering, University College of Swansea. Zi-Cai, L., Ling-Jia, Z. & Hui-Li, W. 1983 Difference methods of flow in branch channel A.S.C.E.,Journal of Hydraulic Engineering, 109 (3), 424-447. Zienkiewicz, 0. C. 1977 Finite element method, 3rd ed. McGraw-Hill, New York.