THE NUMERICAL MODELLING OF TIDES IN A SHALLOW SEMI- ENCLOSED BASIN BY A MODIFIED ELLIPTIC METHOD by MULIA

A numerical model which renders possible the numerical solution of the nonlinear tidal equations by obtaining a solution for each individual tidal component is developed. For this purpose, a set of time independent non-linear equations for each tidal constituent is constructed. Each of these sets of equations is interrelated through the non-linear frictional terms, the approximation of which is accomplished by an iterative scheme. The method is tested for several models before it is applied to the real basin (Bight of Abaco). In order to evaluate the model and to construct the boundary conditions along the opening, a series of tidal observations were undertaken. The viability of the method is indicated by the fact that the results of computations using a coefficient of friction r = 0.0034 give good agreement with observations for all components and over all stations. STATEMENT OF PROBLEM The quantitative modeling of shallow water tides has long been of interest to oceanographers. A major consideration in this modeling has been 1) the representation of the frictional terms in the tidal equations and 2) the method of solution of the resulting equations. For physical reasons one expects the appropriate representation of the frictional terms to be non-linear. Frequently, a linear representation arising from a linearization of a quasi-quadratic non-linear representation has been assumed. This linearization can successfully be applied if the frictional force is small as in the case of a deep ocean. In shallow water, however, where the frictional force plays a significant role in the momentum balance (PR OU DM AN 1952; BOW DEN 1953), the linear theory has been inadequate. In particular it cannot describe several important aspects of shallow water tides, such as the interaction among tidal components and the associated generation of so-called shallow water tidal components. Non-linear theory gives deeper 2 Mar. Res. Indonesia Vol.21, 1978: 1-47


STATEMENT OF PROBLEM
The quantitative modeling of shallow water tides has long been of interest to oceanographers.A major consideration in this modeling has been 1) the representation of the frictional terms in the tidal equations and 2) the method of solution of the resulting equations.For physical reasons one expects the appropriate representation of the frictional terms to be non-linear.Frequently, a linear representation arising from a linearization of a quasi-quadratic non-linear representation has been assumed.This linearization can successfully be applied if the frictional force is small as in the case of a deep ocean.In shallow water, however, where the frictional force plays a significant role in the momentum balance (PR OU DM AN 1952;BOW DEN 1953), the linear theory has been inadequate.In particular it cannot describe several important aspects of shallow water tides, such as the interaction among tidal components and the associated generation of so-called shallow water tidal components.Non-linear theory gives deeper M. M. SIDJABAT finite differences.A stepwise integration is performed starting from initial conditions until a steady state solution is obtained.By a special arrangement of the grid system used in modeling the area under consideration, Hansen's method can treat all the non-linearities in the equation without many problems.Hansen's method has been recently employed extensively to calculate tidal oscillations in various types of basins by a number of authors (GRIJALVA 1962;BRETTSCHNEIDER 1967), and their results agree quite well with observations.The disadvantage of Hansen's method is that it does not readily give information on the propagation of each individual component of the tide or on the nature of the interaction among the tidal constituents.
The purpose of the present study is to develop a model which renders possible the numerical solution of the non-linear tidal equations by obtaining a solution for each individual tidal component.This is accomplished in the following way.It is generally recognized that the tides can always be resolved in terms of the principal tidal components and their higher harmonics regardless of the complexity of the structure of the area or. the non-linearity of the motions.This resolution allows one to construct a set of time independent non-linear equations for each of the tidal constituents.However, each of these sets of equations is not totally independent from the others.They are inter-related through the non-linear frictional terms which are here taken to be quasi-quadratic in the velocity.An iterative method is employed to produce an approximation of the interaction terms.At each stage of the calculation, an elliptic type equation is obtained for each tidal component, which can be solved by a numerical method similar to that employed by Defant and Sgibneva.The first step in investigating this approach is to design a basic elliptic method for a simple model basin, one of uniform depth and constant cross-sectional area, closed at one end with the other end open.One may then compare the numerical solution with the analytical solution obtained from simplified forms of the governing equations.The method is then tested on several idealized models before it is applied to a real basin.The convergence of the iterative scheme in approximating the interaction terms is investigated.
The method is then applied to a real semi-enclosed basin (Bight of Abaco).The tidal characteristics along the opening of the basin which are used as the boundary conditions for the numerical modeling, are determined from the experimental data by ordinary harmonic analysis.In order to check the validity of the proposed numerical model, these results are compared with the experimentally determined tidal topography of the area under study.This comparison allows the determination of an optimum value for the frictional coefficient.

Notation and Basic Equations
A Cartesian coordinate system is employed throughout.This system is oriented so that the origin and the x,y plane lie on the undisturbed water surface; z is positive upwards.Two dimensional horizontal vector quantities are capped with an arrow; two such quantities are the horizontal position vector: The symbols that will be used are as follows: The basic equations which will be used in the analysis in this study result from the application of the physical principles of conservation of momentum and mass to the tidal process.This application results in one vector and one scalar partial differential equation with kinematic and dynamic boundary conditions at the surface and bottom (LAMB 1932, KINSMAN 1965).Neglect of several terms in the vertical component of the momentum equation namely the vertical acceleration of water particles and the vertical component of the Coriolis force, leads to the elimination of the pressure and its replacement by the hydrostatic The Numerical Modelling Of Tides In A Shallow … (Mulia M. Sidjabat) approximation in the remaining equations.These equations are then vertically averaged to yield: In the first of the above equations the notation f x u (x,t) denotes the horizontal vector quantity (-fv, fu).In the long wave approximation, the inertial terms are frequently neglected since their relative magnitude is small compared to the leading terms in ( la ) (although these terms may be significant in a shallow bay).The area under consideration is a small basin (wavelength is large compared to the length of the basin), so that direct tide generation by moon and sun is negligible.The tide is regarded as generated in the open ocean and propagated into the bay.The specification of a unique solution requires the application of proper boundary conditions (RICHTMYER & MORTON 1967).The above equations are subject to two types of boundary conditions: 1. Fixed boundary (along the coast).The motion in the coastal region is fairly complicated and difficult to treat mathematically.
The boundary condition here is determined by the nature of the bottom topography along the coast, which will be too complex if real topography is considered.An attempt will be made to assume a hypothetical vertical boundary and apply the condition where n is a unit horizontal vector perpendicular to the boundary.This assumption is permitted so long as the transport through the hypothetical boundary is insignificant (MIYASAKI 1964); 2. Free boundary (along the opening).The characteristic of the tides along the opening are given (i.e., determined from tidal observations), where g(x,t) is a well behaved function.
Some consideration of the physical characteristics of the tides in shallow water encompassed by the above equations will be helpful.The tides of a basin are considered to be part of a co-oscillating system (DEFANT 1961) in which the period is determined by the tides in the outer sea.The character of the motion in the basin is influenced qualitatively by the geometry of the basin in several ways: 1.If the estuary is convergent, the amplitude tends to increase.2. Energy dissipation by bottom friction tends to reduce the amplitude of the incident wave.3.If the basin is small, the effect of the Coriolis force may be negligible.
Since the region under study is a shallow basin, the resistance force has an important effect, so that special care must be taken in evaluating the influence of this force.The time lag between the tide in the bay and that along the opening is a prominent manifestation of this property.BOW DEN (1953) pointed out that in a shallow bay, the frictional effect is more important than the Coriolis force to -the tidal motion.

Discussion of the Frictional Terms
The physical nature of bottom stress due to friction is difficult to ascertain.It has been customary to assume that this stress τ B , is approximated by the formula where u is the mean current (PROUDMAN 1952, BOWDEN 1953).This "quasi-quadratic" formulation makes the friction due to the bottom proportional to the square of speed and in a direction opposite to that of the flow.
Due to the intractability of the non-linear formulation, various attempts have been made to approximate the problem.Frequently, a linearized form of the friction is employed.Several alternate forms have been used, but.they are still essentially linearizations.BOWDEN (1953) andIPPEN (1966), for example, approximated the frictional coefficient from the mean velocity throughout the bay, a method which can only be applied to a basin whose boundary can be represented by a simple function.Some difficulties arise from this linearization.If one replaces r| u| u by the expression Ru, then the value of R depends upon the velocity u.The value of R here will be different from the value of r and it has a dimension (time) -1 .The other difficulty is that the interaction among the tidal components is absent.The linearized anlysis of tidal motion can lead to results which are reasonably accurate only when the frictional force is small, such as in the case of deep water.
When a tide propagates into a shallow basin, its motion is greatly distorted due to the bottom friction and shallow water constituents are generated.These constituents can only be generated by a non-linear frictional term and by other non-linearities in Eqs. ( la ) and ( lb ).The non-linear terms determine the mutual interaction of the tidal constituents and the propagation and shape of tides in shallow water.
The friction law requires the use of a friction coefficient, r, the numerical value of which depends upon the physical state of the medium as well as upon the nature of the bottom of the sea.A value of r = 0.003 has been found by actual observations of several workers (BOWDEN & FAIRBAIRN 1952, ROSSITER 1963) to be of the right order of magnitude for tidal currents.However, in practice, the value of r used for a particular area should be adjusted to compare with actual observations.

Discussion of Coriolis Term
It has been generally stated that for a relatively narrow basin, the Coriolis force can be neglected (PROUDMAN 1952, DEFANT 1961).However, DRONKERS (1962) suggested that in a shallow bay with a wide opening, it may be important to consider the effect of this force on the tidal motion.In order to evaluate this question, one should consider the relative magnitude of the Coriolis force compared to other terms in the governing equations.
For the region under consideration it can be easily shown, that the non-linear acceleration is much smaller than friction term and that the ratio of the Coriolis term to the pressure gradient is approximately 0.2.This shows that the magnitude of the Coriolis force is substantially smaller than the leading terms.In the actual situation the relative magnitudes may differ from the ratio obtained above.Since the inclusion of this term into a numerical model will not introduce a major problem, it is just as well that it be taken into account.In doing this, one will be able to determine the relative importance of the Coriolis term.

METHOD OF ANALYSIS
In order to transform the tidal equations into an elliptic type equation for the surface elevation, one has to transform a set of nonlinear time dependent equations into several coupled sets of nonlinear time independent equations.This can be done by representing the surface elevation and current as harmonic functions of time and substituting into the governing equations.However, the resulting time-independent equations are not independent from each other.Each set contains a nonlinear term representing self interaction and interaction with other tidal components.
The harmonic representation of the surface elevation and current are as follows where the summation extends over both positive and negative values of the index n; ζ n (x) and U n (x) are complex functions of position and satisfy the relation holds for each tidal frequency ω n , and the sum ranges over all principal tidal frequencies and linear combinations there of; f' and U' are the residual surface elevation and current.These are random terms induced by the surface stress τ S .
By substituting the expression ( 4) and ( 5) into the tidal equations ( la ) and ( 1b ), the following equations are obtained: 9 The Numerical Modelling Of Tides In A Shallow … (Mulia M. Sidjabat) For simplicity of the notation, the x dependence of the variables ζ n (x) and U n (x) has been suppressed.The non-linear acceleration and Coriolis terms, for the time being, are not included.In addition, the surface elevation ζ n is regarded as small compared to the depth, which enables one to linearize the continuity equations.Also, it is assumed that the tidal component of the wind stress is negligible.Therefore, this stress contributes only to random terms in the above equations.
In order to construct the set of equations for each tidal component, further analysis of the frictional term is necessary.The frictional term in Eq. ( 6a) can be expressed as where T O , T^ and T2 are the terms that contain the zero, first and second order term with respect to e of the Taylor Series, where e is given by One can clearly see that the zero order term will contribute to the set of equations for a principal tidal component.Most of the terms generated by τ 1 will not contribute because the sum of ω m and ω n will never be zero (note that m ≠ -n).However, certain combinations of ω m , ω n and ω 1 may produce a frequency equal to a principal frequency, e.g., In a basin where M 2 is the dominant tide (as in the Bight of Abaco) λ is determined primarily by the M 2 fields.It follows that T 1 can produce a significant contribution only for the M 2 set of equations.Investigation of this case suggests that the size of this contribution is less than 25% of the contribution from T O .AS the principal objective is to demonstrate the viability of a modified elliptic method to model tides in a shallow basin and not to investigate all the ramifications of the assumption of quasiquadratic resistance, T 1 has been neglected in the computations.It should be pointed out, however, that T 1 could be included in the computation within the framework of the present method of analysis.
As in the case of T 1 , most of the terms of T 2 will not generate principal tidal frequencies, except for p = -r and q = -s, e.g., the constant part of T 2 , can be written as

MODELLING OF TIDES IN BASIN
Note that the restriction for n ≠ -m for the sum in the numerator is cancelled, and a positive number, δ, is introduced.Applying the Schwarz inequality for complex functions, one has Where 0 mn is the angle (real number) between the real vector U n and U m .Therefore, it is seen that the numerator of ( 8) is less than the denominator, which implies that In a similar manner, the higher order terms of ( 7) can presumably he demonstrated to be smaller than the second term.Therefore, to a good approximation, the frictional term is represented by In practical applications, one of the terms of the tidal expansion very often predominates over the others, which enables one to represent the frictional force by a single term that has the frequency of the principal component, BOWDEN (1953) and ORONKERS (1962) have approximated the non-linear frictional terms for a pure harmonic tide by expanding these terms in a Fourier series in time.The present scheme is essentially a generalization of this expansion to a multi-component almost periodic tide.
Now, the set of equations for a particular frequency which is obtained from above is Because λ(x) depends upon all the U n (x), the above set is both nonlinear and coupled to the other harmonic sets.On the other hand, if λ(x) is treated as a known function of position the set is linear and uncoupled.In this case, the solution for ζ n (X) and U n (x) can be obtained with comparative ease.This situation suggests the following iterative procedure for finding an approximate solution to the problem.Let the trial solution ζ°(x), Un(x) satisfy the frictionless system of equations At each state of the calculation, U n can be eliminated to give an elliptic type equation of the form where A n , B n , C n , D n , E n and F n are coefficients which can be computed for each x.Applying the boundary conditions ( 2 ) with appropriate substitution from ( 6b ), the above equation becomes subject to the boundary condition which makes the problem well posed of Neumann type (FORYSTHE & WASOW 1960).
The approximation of a non-linear acceleration term can be performed in a similar way to that used in the approximation of the factional term.Consider the inertial term MODELLING OF TIDES IN BASIN and let the solution of the tidal equations for a component of frequency ω k be sought.Then the terms that would contribute from the above expression are the ones the sum of whose frequencies are ω k , e.g., This relation is satisfied when M. M. SIDJABAT is beyond the scope of this paper.However, the author is inclined to judge the performance of the iterative scheme used here in a purely pragmatic way.The method was tried out in a number of models, and it worked satisfactorily.

NUMERICAL MODEL
For the numerical solution of boundary value problems, it is common to replace the differential equation by an approximate difference equation, and the continuous region in which the solution is desired by a set of discrete points.For this purpose, the entire region under consideration is divided into a network of equal squares of side L (Fig. la).This permits one to reduce the differential equations to a system of algebraic equations which involves many unknowns.
As has been mentioned previously, the frictional term is approximated from the previous velocity field.Then for given tidal component Eq. ( 1) can be written as follows Substituting (11a) and ( l 1 b ) into (10c), a single second order elliptic type equation for the complex surface elevation is obtained.

MODELLING OF TIDES IN BASIN
For the points which are located on the boundaries, a central difference approximation cannot be employed, unless so-called fictitious boundary points are introduced (Fig. 1c).By so doing, one is able to set the relation between the fictitious points and the interior points by using the boundary conditions ( 15 The matrix appearing in connection with the finite difference equations of a simple net-system, has a simple structure.It is a symmetric matrix and has many zero elements.This enables one to find the solution of the matrix by a triangularization method ( FORSVTHE & WASOW 1960).However, if the real basin is taken into consideration, the irregularities of the boundaries and the bottom structure upset the symmetry, making the application of the triangularization method complicated.The other alternative is to employ a simple Gauss elimination method.However, elimination is ordinarily coded by actually storing the coefficient matrix for which the storage capacity of the computer becomes a problem.To overcome the storage problem, the scheme of inverting a matrix representing the equations under consideration has to be modified.A brief illustration is described as follows.Let A be a n x n matrix with coefficients represented by where each C k is an m-by-m matrix of sept-diagonal form (the non-zero elements are all on its main diagonal and the 3 adjacent diagonals on each side), and S k , N k are tri-diagonal matrices (VARGA, 1962).This matrix resembles the matrix associated with the problem considered in this paper.The matrix A then is reorganized and reduced into A', in such a way that the diagonal elements of A are on the central column of A' By doing this, one is able to reduce the storage requirements in computer for inverting the matrix.The geometrical structure of the area under study and the grid size of the net system representing it are the main factors in determining the degree of reduction in the size of the matrix.When the matrix has been reduced, a new algorithm is used to invert the new matrix A'.
Although elimination methods are rarely used for solving elliptic difference equations, the author has examined a modified Gauss elimination which can be employed to solve the problem considered in this paper.An iterative scheme, such as an over-relaxation method (FOX 1950. YOUNG 1962) was also tried to solve the simultaneous equations.However, it was observed that computations converge very slowly ( 7 0 iterations) and still indicate significant discrepancies from the true solution (for the case of the rectangular basin) compared to the Gauss method.
From the computed ζ r e , and ζ i m , the velocity can be obtained for each grid point.The real and imaginary parts of the velocity, U and V, are given by MODELLING OF TIDES IN BASIN At each stage of the calculation for each component, the computed velocity components are stored on disk or tape.When the computations for all components have been completed, these velocity components can be retrieved to evaluate the parameter λ(x) which is then used for the next iteration calculations.The computations are repeated to obtain the solution for ζ n and U n for each constituent, until the solutions converge to within a prescribed tolerance.

COMPUTATIONS
The numerical method described in the previous chapter is written in the Fortran IV language for IBM 360/365.Grid points representing the area under consideration are coded in such a way that the program can be applied to any form of basin.In order to verify this program, it is, first tested for a rectangular basin of constant depth.Using a linearized form of the friction term, the numerically computed solution for this rectangular basin can be compared to the analytical solution.The convergence of the computation is also investigated.Following this, computations for the simple model of the Bight of Abaco are performed in order to gain insight into the effect of the depth distribution and coastline configuration to the tidal motion; and also for further investigation of the convergence of the computation.Finally, computations with the real boundary conditions along the opening are conducted for the Bight of Abaco.

Analytical Solution
The analytical solution, which serves as the standard for comparison with the numerical method, is for the linear form of the hydrodynamical equations.For the sake of simplicity, the Coriolis force is not considered.The solution is developed for a one dimensional model, e.g., a channel with closed end and of uniform depth.The governing equations which arise from the linearization of ( la ) and ( lb ) for one dimensional problems are

Numerical Solution
In order to verify the numerical method to be employed to compute tides in the Bight of Abaco and to test the iterative scheme described in the previous chapter, a rectangular basin of uniform depth was chosen as a model.The basin was open at one end and closed at the other.The size and depth of the model were taken as comparable with the real basin under study.The size of the basin was 46 x 46 km.The tidal constituent M 2 was used as a forcing function along the opening.The numerical parameters used in computation are as follows: The rectangular basin is represented by 70 grid points.Ordinarily, the size of the matrix representing the equations considered here will be 140 x 140, however, by the suggested scheme, this size can be reduced to 140 x 40.
The amplitude at each grid point along the entrance was set equal to 17 cm, which is taken as the average amplitude of M 2 component along the entrance of the Bight of Abaco, and the phase is set equal to zero.By setting the Coriolis parameter f equal to zero, and allowing for the linear form of the frictional term, [The value of R (in 22a) used in this computation is an arbitrary value (0.0003) which is irrelevant to the model] the result of the computations can be compared with the analytical solution which was obtained by the method discussed previously.The method of numerical computation was confirmed because at every grid point the computed tidal amplitude corresponds to the analytical solution.Between the corresponding points, the deviation in amplitude is less than 0.01 cm and in phase is less than 0.05 degree.
The computation with nonlinear friction was performed; it took 4 iterations to achieve the convergence of the solutions (after 4 iterations the difference between the solution of two consecutive computations is very small (< 0.1%)), however, there is no way to justify whether the solution is the true one.The numerically computed solution for both computations with linear and non-linear friction were plotted in Fig. 2. In this figure, the computation with non-linear friction was first conducted.Then, from the resulting velocities, the average along the channel was used to evaluate the coefficient of friction for the computation with linear friction.As one may observe, there is a slight difference between the amplitude distribution of the two computations, but a more significant difference is observed in the phase distributions.If the depth of the basin is set larger than 5 m, this difference becomes less significant.
Computations with and without the Coriolis force were also made; the results indicate that the effect of this force is very small (1 to 3%) compared to the total tidal amplitude.

Simplified Model
The grid system for the Bight of Abaco can be seen from Fig. 3.In setting this grid system, several assumptions were made.The Bight was assumed to be closed at the northern end, since the channels between the islands there are narrow and shallow.In other words, the inflow and outflow of water between these islands due to the tides was assumed negligible.This assumption is supported by the tide observation at station 20 where the amplitudes and phases differ considerably from those observed at the northern end of the Bight.For practicality, this opening should also be treated as closed, for the size of it is less than the mesh size' o f the grid system (4572 m) used in the numerical modeling of the region.
The above assumption was also adopted for the southeast boundary of the Bight, where the marls are located.Then, the area under consideration has only one opening, that is the western entrance.A boundary line was drawn along this entrance, in such a way that all the complicated land distribution along the entrance lay outside of the net system.The numerical modeling of the area would be too complicated if this line had to be drawn outside of the land areas, although the setting of the boundary values might be easier.The numerical treatment of this type of problem has apparently not been investigated.
The net-system has 155 grid points, 16 of which are located along the free boundary.The numerical values used here are practically the same as in the previous model, except for the mesh size.The boundary values for the grid points along the boundary are made uniform; the amplitude is 17 cm, and the phase is set equal to zero.The depth is taken to be uniform.
For this model, it is observed that the number of iterations required for convergence of the solutions is larger than that for the rectangular basin.This may be mainly due to the irregularity of the boundaries.In order to speed the convergence of the solution, a new algorithm was introduced in approximating the non-linear form of the frictional term.In place of Eq. (3.7), the equation was used.The parameters p 1 and p 2 are known as weight parameters, which satisfy the relation The effect of using the best combination of p 1 and p 2 reduces the number of iterations necessary for convergence of the solution by a factor 4. Fig. 4 illustrates the number of iterations for 3 different combinations of p 1 and p 2 .It can be seen that p 1 = 0.75 and p 2 = 0.25 is the best combination of the three.There values were used for the rest of the computations.
A computation was carried out using the frictional coefficient r equal to 0.003.From the computed tidal amplitude and phase, the cophase and corange lines were constructed for the model.It was found that the tidal amplitude increases as the tide approaches the northern M.M. SIDJABAT boundary of the model, and at the same time the rate of change in phase decreases.This phenomenon is mainly a result of reflection of the tidal wave by the corresponding boundary.Near the boundary, the size of the amplitude is twice as large as the size of the incoming wave.As the reflected wave propagates away from the reflecting boundary, the amplitude decreases again due to the damping effect of the bottom on the model.If a frictionless model is considered, a pure standing wave is obtained with zero phase everywhere in the model.A computation was also performed for the model where the real depth was considered.The boundary conditions are the same as before.The use of real depth in the computation results in significant changes in the amplitude and phase distributions.One can see that, due to the small depth along the entrance, the tidal amplitude drops drastically, but as the tide enters the deeper part of the model its amplitude increases.In this model the smallest amplitude is observed right after the entrance.This can be explained by similar arguments to those considered for the model of uniform depth, but in this case the reflected waves apparently have small effects at the regions near this entrance.The difference in phase distribution in the 2 models is also noticeable.Due to the shallower depth near the coast for the model with real depth, the phases are retarded compared to the model with uniform depth.Both models demonstrate that the coast and bottom configuration have very significant effects on the propagations of tides.

Real Model
The grid system used for this model is essentially the same as the one for the simple model.The real depths are taken into consideration; the depth distribution which was obtained from the Hydrographic Chart (H.O.5990) -published by U.S. Naval Oceanographic Office in 19H3, is illustrated in Fig. 3. Again, the numerical parameters used here are also identical to those in the previous model except for the number of tidal components considered in the computations.The boundary values along the opening were obtained from the observed tides from the 5 tidal stations that were installed along this entrance.The proper setting of the tidal constants along this opening is the most important factor in obtaining an accurate solution.These 5 stations were not exactly located along the line that is used as a free boundary of the model.Since the bottom configuration along this entrance is so complicated and because of the limited number of observations in the neighborhood of the opening, it would be difficult to obtain the boundary values based on these observations alone.

MODELLING OF TIDES IN BASIN
The following procedure was therefore adopted to carry out the determination of the boundary values for the numerical computations.First, corresponding to each of these 5 stations, a standard point was determined along the free boundaries of the model in such a way that the point was the projection of the station to the opening line of the grid system.Then, the approximated tidal constants at these 5 points were determined by visual inspection, from which the boundary values for each grid point along the opening can be obtained by cubic spline interpolation (through the aid of curves).This interpolation was performed for each tidal component both in amplitude and phase.Using these initial guesses of the boundary values, the computations were performed as prescribed previously.These computations give the solution of each component for all grid points in the model.From these results., the computed tidal amplitudes and phases for each boundary station were obtained by polynomial interpolation.By doing this, the numerically computed tidal constants can be compared to the observed ones.It is recognized that for this stage of calculation, the difference between the computed and observed tides are significant.The source of these large deviations is due mainly to the improper setting of the boundary condition along the opening.These boundary conditions can be improved by readjusting the tidal constants at those 5 standard points in accordance with the numerical results.Following this, the computations were performed again, and it was observed that the results of computations were improved.The above procedure can be repeated several times until the observations and computations at the stations along the entrance agree within a prescribed tolerance.In the computations considered here, the third readjustment of the boundary value was found to be satisfactory.This final setting of the boundary conditions is shown in Fig. 5.The difference between the computed and observed values for both amplitude and phase for the stations along the opening is less than 19f.For the adjustment of the boundary condition the value of r = 0.003 was used on the computations.
When the proper boundary condition along the opening had been determined, a number of computations using different values of r (from 0.002 to 0.006) were conducted.It appears that this small) change of the value of r does not require the readjustment of the boundary conditions.
It was observed that the best value for the frictional coefficient r is 0.0034.Using this value, the tidal amplitude and phase are computed for the 5 tidal components considered in this study.From the results, the corange and cophase for each constituent were constructed (Figs. 6 -7).One can observe that the appearance of the amplitude and phase distribution depends upon the period of the tidal component.The semi- diurnal components have the same qualitative distribution and the same thing holds true for the diurnal components.However, there is a noticeable difference between corange of K 1 and O 1 components.This difference is mainly due to the different amplitude distribution of K 1 and O 1 along the opening (Fig. 5), which in turn may result from errors in the harmonic analysis for these constituents due to the shortness of the records and the meteorologically induced noise.

Tidal Observations
The tidal observations described in this paper were made in the Bight of Abaco, Bahamas, by R. L. Snyder of Nova University and J. H. Filloux of Gulf General Atomic Corporation, San Diego.A chart of the Bight is shown in Fig. 8. Tide gauges designed and built by J. H. Filloux (1969) were used.These were shallow water versions of a small, high-resolution, self-contained and bottom mounted instrument designed to record oceanic tidal fluctuations anywhere in the ocean.
The observations were conducted from February 5 to March 24, J969; and 15 tidal gauges were installed throughout the region.A second series of observations was run from July 10 to August 17, 1969.For accurate positioning of the tidal gauges, Decca Hi-Fix was used.Five gauges were placed on the inside entrance of the Bight, the information from which was used to set the boundary condition along the opening for the numerical computation.Inside the Bight, 10 gauges were installed.The observed tidal amplitudes and phases obtained from these gauges were used as comparisons to the computed ones.In order to gain a qualitative understanding of the tidal propagation through the northern and western openings, 4 gauges were planted outside the corresponding entrances.In addition to the tide gauges, four weather stations were installed.The observations were made from R.V. GULFSTREAM and R.V. L.F.R. BEL'LOWS, research vessels of the Physical Oceanographic Laboratory of Nova University.

Harmonic Analysis
As mentioned previously, the surface elevations were inferred from the pressure changes measured at the bottom of the sea.Hence, the atmospheric pressure fluctuations were also incorporated in the tidal records, e.g., The Numerical Modelling Of Tides In A Shallow … (Mulia M. Sidjabat) where p is the pressure recorded by a tidal gauge and p a is the atmospheric pressure which was also recorded simultaneously at the surface of the sea.However, since the only significant period in the atmospheric fluctuations that is also studied in the tidal fluctuations is S 2 (period of 12 hours), the above correction has been applied directly to the analyzed constants for this constituent.It was assumed that the atmospheric pressure is uniform everywhere in the Bight of Abaco.Therefore, information from only one record atmospheric pressure was needed.In order to evaluate the hydrodynamical model considered in this study, it is necessary to numerical computations with observations.There are several ways to perform this comparison, and one of these is to calculate the variance between the computed tides and the observed tides.This variance will be defined as the mean square value of the difference between the computed tides and the observed ones.This variance can be represented in four different ways.

4.
Average mean square difference due to all components over all stations for different value of r, e.g., The value of the variances discussed above are given in Tables 1 to IV.The tabulated variances in Tables I and II indicate that the computed tides are in good agreement with the observed ones both in amplitude and phase for the two given components.It is observed that the variance varies from station to station as well as from component to component.A significant deviation is noted for station 8.Although the computed amplitudes seem to agree very well with the observed amplitudes, the differences in phase are very noticeable.These discrepancies are probably the failure of the numerical model to adequately represent the complex bottom topography just south of the Grand Bahama Island.On the other hand, observations and data analysis may contribute errors to a certain degree.From Table I I I , one can see that the relative magnitude of the variance for S 2 and N 2 are larger than for the other 3 components.This fact could be due to several factors other than numerical modeling.The amplitude of S 2 and N 2 are small compared to the amplitude of M 2 .K 1 and O 1 and they are also small compared to the noise level M. M. S SIDJABAT accompanying the observations.A small error in the observed values for these components may produce a large proportional deviation from the computed values.It is also relevant to point out that in determining the amplitude and phase of the tidal records, a conventional harmonic analysis was employed without making any further correction for the interference of other tidal constituents that have frequencies in the neighborhood of the principal components.'These neighbor components were not taken into account in the harmonic analysis because of the noise level and record length.Although the analyzed constants for K 1 really represent a combined curve for K 1 and P 1 (and similarly S 2 for S 2 and K 2 ), this does not affect the principal aspect of this study if the observations are simultaneous for the whole region.However, the data for stations 8 and 26 are about three months later than the others and therefore the data for these stations are not as reliable in the calculations.Fortunately, station 26 is located at a relatively non-sensitive position of the area.Since the contribution of station 26 to the overall solution in the Bight is therefore reduced, this also suppresses the effect of the error on the solution.Also, M 2 , the principal constituent which mainly determines the frictional force, is not affected by the change in period.Error estimates for the observed tidal amplitudes and phases were not made because of time limitations.
In the following, several important factors that bear on the accuracy of the hydrodynamical and numerical modeling are discussed.
1.The hydrodynamical equations that are used in the model have been simplified to some extent in order to simplify the numerical modeling.2. The assumption that along the hypothetical boundaries does not correspond to the actual conditions.3. The inaccuracy of the depth distribution.4. Since the numerical modeling entails the representation of a continuous system by a discrete one, an error can be introduced.This error can be particularly significant if the scale of the change in the variables of the problem is locally small compared to the mesh size. 5.The computations are subject to round-off error.Since the iterative scheme is used, each computed value will influence the values which are yet to be determined.In reality, this distribution of errors is often more or less a random process, and the effect of the errors will partially cancel out.One of the advantages of the method considered in this study over Hansen's method is that this method can readily give information on the tidal characteristics of a tidal component.In Hansen's method, the tides for each grid point are computed as a time series, from which one has to perform a harmonic analysis in order to obtain the tidal components.On the otherhand, a disadvantage of this method is its limitation to a system that can be represented by harmonic functions of time.
The numerical method was written in the Fortran IV language for IBM 360/365-G compiler.Storage requirements were 150 K bytes and machine time approximately 20 minutes for a 5 component calculation with four different values of friction coefficient.
P i g u r e 1 .Grid points for an a r b i t r a r y region ( a ) and coding systems used in the n u m e r i c a l m o d e l ( b a n d e ) .T h e s y m b o l s u s e d i n t h i s f i g u r e a r e d e f i n e d a s f o l l o w s : where the notation of E, VV, N, S and C stand for East, West, North, South and Central points, respectively (Fig. lb).
a,b) or (16 a,b) at the boundary points.To obtain more insight into this, consider a point on the east boundary.The boundary conditions in a finite difference form are where an overbar denotes a fictitious point.Hence, the expression for the fictitious points are give by Then, the set of difference equations for the east boundary point, by using Eqs.( 19a ) and ( 19b) can be written as follows MODELLING OF TIDES IN BASIN Thus, in obtaining solutions for the surface displacement, ζ,. in a region under consideration, Eqs. ( 17 ) and ( 18 ) must hold at every interior point, and Eqs. ( 15 a,b) and ( 16 a,b) must be satisfied at any appropriate fixed boundary point.

Figure 2 .
Figure 2. Phase and amplitude distribution of the computed solutions for M2 corresponding to a rectangular basin with uniform depth.The lines with triangles are for the computation with linear friction and the lines with the circles are for the nonlinear friction.
Figure 3. Depth distribution of the Bight of Abaco.

Figure 4 .
Figure 4. Rate of convergence of the solution for the amplitude and phase for a model basin at a particular grid point (at the middle of the basin) for different combinations of pi and p 2 :

Figure 5 .
Figure 5.The final setting of the amplitude and phase distribution of the 5 constituents along the opening based on the observed and computed tides at the stations that are located along the entrance.The symbols used in this figure are defined as follows:

Figure 8 .
Figure 8. Geographical position of the Bight of Abaco.

Figure 9 .
Figure 9. Distribution of average mean .squaredeviation due to all components over all stations with respect to the frictional coefficient r (the solid line for the consideration of Coriolis force and the dashed line computations without Coriolis force).(Note: the scale of r should be multiplied by √% and the curve with Coriolis force should be disregarded -see Foreward).

Table I .
Variance due to M 2 for each station