Introduction
 

The Gulf of Riga is a well-defined semi enclosed water basin (see fig. 1) connected with the Baltic Proper by two straits. The Irbe strait is rather wide (30 km) with the sill depth about 20 m, while the Muhu strait is rather narrow (5-10 km) and shallow (minimum depth 10-12 m). Both straits have a complex topology and depth distribution [1]. There are seven islands in the Gulf of Riga. Most of them are near the northern cost, except for Ruhnu, that is located in the central part. The depth in the central and southern part is 40-50 m, but in the northern part 15-20 m, average depth of the Gulf is 25.5 m [2]. The coastal line is rather smooth at the Latvian coast and Estonian coast near the Gulf of Pärnu, but it has the strongly varying curvature along the coast of Saaremaa. The quite complex bathymetry and boundary, complicity of the circulation and external influence characterizes Gulf as an interesting and difficult object for the simulation of the physical processes.

Figure 1. Gulf of Riga. Finite element mesh.

Adequate modelling of the physical processes is unsolved problem for the Gulf for a while. The physical structure of the Gulf is rather complicated and usually shows a three-dimensional behaviour. Development of the fully 3D operational hydrodynamical model requires proper understanding of the physics of the Gulf what can be achieved by means of development of the pilot models for the particular physical effects. This is especially important, because the verification of the 3D hydrodynamical models for large natural basins would be a problem itself.

The model of the vertical temperature and salinity structure of the Gulf of Riga carried out in [3] can be considered as the first step to the development of the 3D model and understanding of the vertical structure of the physical parameters. Next step must provide investigation of the horizontal structure of the Gulf and the development of an appropriate model for it. The main objective of this paper is to investigate the possibility to employ the shallow water equations for the calculation of the circulation in the Gulf of Riga. Calculations of the typical flow patterns as well as real-time calculations for selected time periods are also included in the paper.
 

Governing equations and numerical scheme

The shallow water equations (SWE) in conservative form [4] in the Boussinesq approximation stand as follows:
 

	(1)
	(2)

 where is the unit-width discharge, x is the elevation over a reference plane, h is the depth of the water column, m is the turbulent momentum exchange coefficient, g is the acceleration due to gravity, is the angular velocity of the Earth, C is the Chezy coefficient, is the wind velocity, is the drag coefficient for the air flow over water (C_d = 1.1× 10^-3 in the present paper), is the density of air, is the reference density of the water (r ₀ = 1000 kg× m^-3), is the deviation from reference density, p_a is the air pressure. These are standard SWE, except for the last term in the right hand side of (1), which represents additional pressure gradient due to the density gradient [4]. The equation of the state of the sea water [5] is used for determination of as a function of temperature (T) and salinity (S). The transport equations are used for T and S

(3)

 where c stands either for T or S, is velocity and D_c is the respective turbulent exchange coefficient.

The numerical scheme for solving (1)+(2) proposed in [6] is adapted here. The main advantages [6] of it are the unconditional stability and the absence of the spurious oscillations in the elevation field. An equation splitting is used at every time step to decouple the physical contributions. Convective transport is solved by means of using lagrangian scheme [7]. Spatial discretization of (1)+(2)+(3) is implemeneted using the standard Galerkin finite element method [8].

Density term in (1)+(2) causes interesting consequence. As density is the function of temperature and salinity, it is possible to write following equations:
 

	(4)
	(5)

Assuming non-diffusive case, summing up (3) for T and S, and using (2), it is possible to derive following transport equation for density:

(6)

Important conclusion from (6) is that at the steady state the streamlines must be parallel to the isolines of density as well as to the isolines of temperature and salinity.

The lagrangian integration does not allow to implement the free-slip boundary conditions (BC) in the weak way just neglecting the contour integral [6] arising from the integration by parts of the convective terms in the weak form of (1)+(2). Therefore the condition =0 is used in [6] for the coastal boundaries of the computational domain. This boundary condition does not allow the convective transport of salinity and temperature along the boundary. Hence, changes in T and S at the boundary are due only to the diffusion in that case. There is no aim to model the thermal and concentration boundary layers, and no mesh refinement that would allow convective transport close enough to the boundary is anticipated. Therefore it is intended to use enhanced free-slip BC for the discharge at the closed boundary in the present paper. The system of linear algebraic equations is obtained after finite element discretization for each component of the discharge. Implementation of the free-slip BC has as a major consequence coupling of the both those systems. There are the two uncoupled equations before imposing BC at each node i of the computational mesh:

	(7)
	(8)

where q_i^x and q_i^y are the x and y components of the discharge at the node i, A_ij is the element of the stiffness matrix, b_i^x and b_i^y are the right hand sides of the discretized equations at the node i. The procedure of imposing the free-slip BC is developed as follows. Free-slip BC can be expressed as , where =(n_x, n_y)^T is the unit outward normal vector of the domain. Thus, applying BC at the boundary node i is enforced by replacing (7) and (8) by:

	(9)
	(10)

 where the first equation stands as BC itself and the second represents momentum equation for the tangential component of the discharge at point i. It is derived expressing q_i^x  and q_i^y from (7) and (8) and putting them into the expression for tangential component of discharge -n_yq_i^x+n_xq_i^y.The order of appearance of the new equations in the coupled system depends on the values of n_x and n_y. Diagonal predominance of the matrix is very useful for the rapid convergence of the iterative methods of solving the systems of the linear equations. Thus, it is convenient to rearrange equations according to the magnitudes of the diagonal elements.  q_i^x or q_i^y from (9)is posed on the main diagonal of the coupled stiffness matrix if, respectively, |n_x|<0.5 or |n_y|<0.5. One unsymmetric system of linear equations is obtained instead of two symmetric systems of linear algebraic equations imposing free-slip BC. A proper unsymmetric solver is therefore required to solve it. The biconjugate gradient squared method [8] has been seen to be very effective for this purpose, and the version of it from NSPCG package [9] has been used. The main disadvantage here is the increased time of calculations due to more computational effort required to solve unsymmetric system.

Obviously, in the absence of the separation the boundary contours are streamlines. The lagrangian integration of the convective transport along the boundary has been implemented creating a sorted list of the boundary elements, extending approach of [7] where integration at lagrangian step is implemented for inside of the domain.

 Two Dirichlet boundary conditions for the system (1)+(2) must be prescribed at the inflow but one at the outflow from the point of view of characteristic theory [10]. The simple imposition of discharge and elevation according to the measurements is quite standard practice for the river and tidal flow calculations [7]. However, when considering larger water bodies as gulfs, one needs to take into account also impact of Coriolis force, as well as wind stress. Semi-enclosed water bodies as Gulf of Riga are usually joined with the larger water basins by straits. Open boundary conditions therefore are closely related to the strait modelling i.e. development of the strait model or at least assumption of the structure of the motion in the strait is needed to impose open boundary conditions. There is available a large set of the measurements of the elevation [11], but measurements of the velocity are not enough to give the appropriate model forcing. It is the main reason why the Dirichlet boundary conditions for the elevation in the straits are proposed in the present work. Examining linearised 2D-equations of motion including Coriolis acceleration [12], one can determine that setting constant elevation across the strait leads to the zero velocities on this boundary at the steady state. In more general case it leads to the strong non-physical turning of the velocity vectors to the right from the direction of the strait at the boundary (see fig. 2a) and to the possible unstability and non-conservation of mass in the numerical solution. Hence, the measured elevation data can be imposed only at one point of the strait (coastal point). At the other points of the open boundary imposition of proper elevation distribution is required according to some hypothesis of the motion in the strait.

The surface long waves propagate along the strait (which is assumed as straight channel) as plane waves without influence of rotation of the Earth. There are several possibilities [12] of waves in the rotational case:

• Poincare waves have nonlinear dispersion relation and a velocity component perpendicular to the boundary of the strait.
• Kelvin waves have the same dispersion relation and propagation speed as plane waves. There are no perpendicular component of the velocity for Kelvin waves and the distribution of the elevation across the direction of their propagation satisfies geostrophical relations [12] between elevation gradient and discharge.
• Rossby waves are present in straits with non-uniform bathymetry in their cross section.

(11)

a	b
Figure. 2. Dependence of flow on boundary conditions in the straits. Discharge vectors and  constant elevation lines.

Analysis of the application of SWE for Gulf of Riga

The dynamics of the water mass circulation in the Gulf of Riga are known to be determined by the following main factors:

• Wind is the main source of currents in the central part of the Gulf [1]. For the coastal regions it acts mainly as a source of the changes in the elevation.
• The circulation is significantly affected by the Coriolis force. The measure of the relative importance of the convective part and the Coriolis force in the equations of motion [12] is the Rossby number Ro=U_*/fL_*, where U_*and L_* are typical velocity and length scales of motion L_*/U_*is advective time scale. Flow is determined mainly by Coriolis force (geostrophic flow) if Ro<<1. For the central part of the Gulf of Riga Ro=0.1 m× s^-1/(1.22× 10^-4s^-1× 10⁵m)» 0.01. One can suspect therefore that in the central part flow will be mainly geostrophic. Generally it would not hold true for the shallow regions due to the higher velocity and for the straits due to the smaller length scale of the motion. Besides, above estimation will be true mainly for the steady state or slowly changing external forcing. Otherwise, the Rossby number based on the typical time scale of changes in forcing or on the local time scale instead of advective time scale must be constructed Ro=1/f. Then, the constraint >>1/f » 8000s » 2h stands as a condition for quasi-geostrophic flow at all time moments. Kelvin waves of any frequency satisfy a geostrophic relation across the strait, therefore implementation of the strait BC (see previous section) does not require condition Ro<<1, based on the local time scale.
• Bathymetry influences the flow patterns in two ways. Firstly, the distribution of shallow and deep regions determine the regions of motion dissipation and production. Secondly, in the fully geostrophic flow, lines of constant depth are at the same time the lines of constant elevation and the streamlines of motion [12].
• The main source of the kinetic energy dissipation is the bottom friction, mainly in the shallow regions.
• Turbulent dissipation is the another possible sink of the kinetic energy. The estimations for the turbulent viscosity coefficients vary widely (m =1-10⁴m^2× s^-1). The correct values for particular case can be determined by series of numerical experiments compared with measurements. However, it is possible to evaluate importance of viscosity term at least in the parts of Gulf where mainly geostrophic flow is suspected. The horizontal Ekman number E=m/ fL_*² measures the importance of the turbulent diffusion term against the Coriolis force. The estimation for the central part of the Gulf of Riga gives E=100/(1.22× 10^-4×10¹⁰)» 10^-4, i.e. viscosity has minor importance here at least for the small and medium values of µ. This is not generally true for shallow regions and straits due to the another length scale of motion.
• The influence of the Baltic proper to the structure of the motion in the Gulf is due to two factors, the level difference on the open boundaries, and the salinity difference between both water bodies. The variations in the measured levels occurs mainly with periods greater than one day (the changes of elevation due to tides is supposed to have minor importance for the Gulf of Riga). Maximum and minimum observed levels differ by » 2.5 m [11]. Although the measured elevations in different places of the Gulf are in rather good correlation, usually there occurs a level difference between different points, with the extreme values of » 40-50 cm. There is possible presence of a transit flow from strait to strait due to difference in elevations in both straits, as well as motion due to the change in elevation to fulfill the equation of continuity, e.g. if the elevation changes by 0.5 m per 1 day then flux through the strait must be as large as » 90000m³/s (surface area of the Gulf is » 16000 km²). The time of the propagation of the long surface waves from the Irbe strait to the Daugava river mouth is »3-4 h and is usually smaller than typical time of changes of the elevation. One can conclude therefore that applied elevation at the boundary propagates rather fast into whole the Gulf.
• Another main external influencing factor is the difference in salinity of the Gulf of Riga and the Baltic Proper. It causes a difference in density and can induce a density driven currents. Typical values of these differences concerning Irbe strait are D S=1.5-2PSU and Dr =1-1.5 kg× m^-3 [1].
• Influence of the air pressure to the elevation patterns can be estimated by comparing typical length scale of the atmospheric gyres (cyclones and anticyclones) » 1000 km with the size of the Gulf » 200 km. One can suspect then that difference in the air pressure between different locations in the Gulf is of minor importance and that all influence from the atmosphere comes from the wind forcing and applied elevation boundary conditions. The last may be a result of the actual air pressure.
• Differences in the temperature and salinity in the Gulf for autumn periods can reach 3 to 5°C and 1 to 2 PSU, respectively. Resulting density differences can induce currents and influence elevation distribution, e.g., difference in salinity 2 PSU can cause difference in elevation as large as m.

The comparison of the calculated data with the measurements are intended to do in the following three ways. Only the first of them is implemented in this paper.

• Comparison of the elevation at the coastal measurements stations [11] with calculated ones. Main difficulties arise from the local factors that must be investigated at each location and, if possible, screened out before comparison.
• Comparison of the velocity patterns. It is the most comprehensive and accurate verification, however, the measurements of velocity are very occasional. Therefore at the present stage they can serve mainly as the qualitative comparison subject.
• The large data set of salinity measurements are available [2] over a time period of more than 20 years covering most parts of the Gulf. The salinity distribution can reveal main flow features, as well as response of the system to the time dependent forcing. These factors allow to use salinity patterns as the main comparison source. Main difficulties here are rather rear observations at November-December contrary to summer periods [2], that can limit the choice of the time period for which the computations can be done. The influence of the interpolation of the initial data on the results of the calculation must be accounted, too.

Selection of typical situations

Analysis of the steady-state flow patterns would be of interest for classifying the possible flow structures of the Gulf. However, one have to note, that there is no sense for the steady-state situations of the density field because they can be reached during time periods which far extend periods of more or less uniform forcing by the wind and air pressure. Taking into account that the differences in the air pressure over the Gulf can be neglected, as the ²typical² steady-state situations are assumed the ones with the uniform wind over the homogeneous Gulf and prescribed elevation difference between the Irbe and Muhu sounds. The calculation series for 8 different wind directions (D a =45° ), wind speed 7 m/s were performed. Summary of long-term (1977-95) mean water levels in the LHMA monitoring stations for the wind velocities 7± 1m/s, directions a ± 22.5° are summarized in table 1, and plotted for the most characteristic stations for Irbe sound (Ventspils), western (Roja), southern (Lielupe), and eastern (Salacgriva) parts of Gulf on Fig.3. The scheme of the locations of internal stations see Fig. 4.

Figure 3. Long term mean water levels at selected stations for W=6-8 m/s.	Figure 4. Bathymetry of Gulf of Riga. Water level gauges of LHMA.

Wind direction	Ventspils	Kolka	Roja	Mersrags	Lielupe	Daugavgriva	Skulte	Salacgriva
0°	-3.47	-1.94	-1.88	-3.31	10.74	10.80	4.84	-0.05
45°	-6.93	-6.42	-8.04	-8.68	3.06	3.55	-2.49	-6.39
90°	-11.80	-12.04	-16.89	-17.84	-8.62	-7.90	-15.01	-17.64
135°	-5.37	-5.36	-10.32	-11.63	-4.83	-2.78	-7.39	-9.42
180°	6.86	8.14	6.96	5.15	12.15	15.59	11.66	10.56
225°	11.55	14.31	15.11	12.72	23.06	26.46	22.76	20.21
270°	11.27	14.25	14.74	11.95	26.24	28.64	24.56	19.88
315°	3.04	5.56	5.40	4.52	18.90	19.62	14.64	9.24

 Table 1. Average water levels at different monitoring station for different wind directions.

Water level in between Kolka and Ventspils means was set as boundary condition for the Irbe strait. Absence of representative time-series of water level for Muhu sound (Virtsu station) caused necessity to produce more artificial ² characteristic² boundary condition there. Assumption of plane mean water surface of the northern part (neglecting Lielupe and Daugavgriva stations with possible large fresh-water influence to water level) lead to a least square solution for the water level difference between Irbe and Muhu in table 2.

Wind direction	0°	45°	90°	135°	180°	225°	270°	315°
D h (Irbe-Muhu), cm	1.59	2.73	4.08	2.41	-2.37	-4.19	-2.44	0.20

Table 2. Applied elevation difference between Irbe and Muhu straits for different wind directions.
The numerical parameters C=100 m^1/2/s, m=50 m²/s are used for all calculation series.
 

Discussion

The discharge and elevation for selected three wind directions (S, NE, and W) are plotted on figs. 5-7.
 

Figure 5.Water elevation and discharge fields. South wind 7 m/s. Elevation difference Irbe- Muhu -2,37 cm.	Figure 6.Water elevation and discharge fields. North-East wind 7 m/s. Elevation difference  Irbe-Muhu -2,73 cm.	Figure 7. Water elevation and discharge fields. West wind 7 m/s. Elevation difference Irbe- Muhu -2,44 cm.

The main features of the computed circulation patterns are as follows:

• Rather obviously, the elevation is lower at the upwind coast and higher at the downwind coast (e.g. for the west wind fig. 7 the lower elevations are at the western coast, higher at the eastern coast). The elevation isolines nearly follow the bathymetry lines (see also fig. 4) and streamlines in the central part; that confirm the qualitative concept of the geostrophic flow. Furthermore, due to employed hypothesis of the motion in the straits the flow between them are close to geostrophic, too. The difference between maximum and minimum elevation is » 8-10 cm.
• There are at least one cyclonic and one anti-cyclonic vortex present in each of the eight computations. One gyre of size » 60-70 km is located in the southern part of the Gulf. Its direction depends on wind, cyclonic (counter-clockwise) circulation is found for N, NW, W (see fig. 7) winds, anti-cyclonic (clockwise) for S (see fig. 5), SE, E winds. The calculations indicate presence of two vortexes of the opposite direction in the southern part of the Gulf for the NE (see fig. 6) and SW winds.
• The series of smaller vortexes with the same direction as of the southern one extend northward along the western coast of the Gulf, especially for S (fig. 5), E, SE winds. They may be responsible for large level fluctuations in Mersrags.
• Another main vortex is located in between Kihnu and Ruhnu islands. Its rotation is opposite to the southern vortex. The anticyclonic circulation (figs. 6-7) is more intensive there in comparison with cyclonic one.
• The flow patterns are quite complicated in the Irbe strait itself. Following the complex bathymetry lines there is possible existence of vortexes which requires a careful data analysis when trying to calculate water exchange values from velocity field. The extension of the computational domain further into the Baltic Proper may be required.
• The comparison of the velocity and elevation distribution for the calculations with the opposite wind directions shows that flow patterns are nearly opposite to each other, except for the strait zone and zone near the cape of Kolka. Hence, the elevation distribution, the directions of vortexes, and, even the water flows through the straits are nearly opposite for such two calculations.
• Increasing the coefficient of turbulent viscosity causes the flow deviations from the geostrophic and tends to the smoothing of the vortexes. The decreasing of magnitudes of the velocities takes place, too. The influence of the change of the coefficient of viscosity to the distribution of elevation is negligible at the coastal zone and small in the center of the Gulf. The similar evaluation was done for the Chezy coefficient, its decreasing lowers the magnitude of the velocity. It means that precise value of the viscosity coefficient can be determined by comparing results of the simulation with the measured patterns of velocity paying more attention to the flow structure, whilst the calibration of the Chezy coefficient requires mainly the comparing of the magnitudes of the velocities.
• There are both in - and outflowing currents at the steady-state flow patterns through the Irbe strait, while the unidirectional flow is always present through the Muhu strait. The distinct transit flow Muhu-Irbe is obtained for N, NE (see fig.6) winds, but it is less developed in reverse direction for W, SW, S (see figs. 5, 7) winds. Summary of the discharges through the both straits is given in fig. 8 and can be explained by the superposition of the applied wind stress and elevation difference. One should note that the maximum in - and outflows through the Irbe strait are computed for N, NE winds with inflow along Latvian but outflow along Estonian coast (see fig. 6) while minimum water exchange would be expected for W, NW (see fig. 7) winds, which are most favourable for the development of the vortexes in the Irbe strait.

  Figure 8. Calculated steady-state discharges through the straits and applied water level difference for different wind directions.

• The calculated water levels are less than observed at almost all observation stations. Average D h between calculated and observed values for all observation stations are summarized in table 3. There, and further the water level is referred as water level above the one at the Irbe boundary. The water level difference at the Kolka station can be explained by the sensitivity of the water level with respect to location near the cape of Kolka. Water level differences at another stations most probably are due to the deviation of the water density from the mean density of the Gulf at these coastal stations, i.e. higher density at the Roja station and lower (according to dominant spreading of riverine water along the eastern coast) at all others. Besides, water gauges may be located at the places with permanent two-layer stratification (riverine water above saline). The plots of the observed water levels (shifted to the respective values from table 3) together with calculated ones are given in fig. 9 for characteristic western (Mersrags), southern (Daugavgriva) and eastern (Skulte) stations.

  Figure 9a. Calculated and measured (shifted) elevations at Mersrags station.	Figure 9b. Calculated and measured (shifted) elevations at Daugavgriva station.	Figure 9c. Calculated and measured (shifted) elevations at Skulte station.

•  Comparison indicates reasonable qualitative and ever quantitative agreement, with exception for E, SE winds. This disagreement most probably is due to non-satisfied assumption about the plane water surface of the northern part of the Gulf, which has caused the overestimate of the applied Virtsu water level for these wind directions.

  Station	Kolka	Roja	Mersrags	Lielupe	Daugavgriva	Skulte	Salacgriva
  D h (cm)	0.40	-1.33	0.64	9.33	11.06	5.76	3.23

   Table 3. Difference between observed and calculated mean water levels at different stations.

Real-time calculations

Real-time calculations are performed for two one-and-half months long time periods in autumns for which the water-level time-series at Virtsu was available, from 5-Oct-93 until 16-Nov-93, and from 16-Nov-94 until 31-Dec-94. The calculations are performed assuming homogeneous Gulf, applying uniform wind field from the measurements at Daugavgriva and water elevations measured at Ventspils and Virtsu. All the forcing parameters were available two to four times daily, linear interpolation was done in between the measurements. The histogram of the water level at the straits for both periods are shown on fig. 10 but the wind direction distribution is depicted on fig. 11.

Figure 10. Histogram of applied water levels for calculation periods.	Figure 11. Wind direction probability for calculation periods.

The time periods are quite different from the point of the forcing conditions (the summary see table 4). The water level of the Gulf is almost half a meter lower in the first period with the higher atmosphere pressure and prevailing S winds. The observed water level fluctuations are higher at the Muhu strait, the average water level is slightly higher at Irbe for autumn-93, but vice-versa for autumn-94.

	X-XI/93	XI-XII/94
Average wind speed (m/s)	5.5	5.7
Prevailing wind direction	S,SW	SW,W
Average air pressure (mbar)	1020	1011
Average water level at Irbe (cm)	-25	21
Water level fluctuations at Irbe (cm)	-50/0	-5/51
Average water level at Virtsu (cm)	-29	22
Water level fluctuation at Virtsu (cm)	-64/14	-10/61
Irbe-Virtsu level difference (cm)	-39/38	-30/21

 Table 4. Summary of forcing data for calculation periods.

Appropriate initial conditions are one of the problems for the real-time calculations. The measurements of the velocity covering all the Gulf for the fixed time moment are impossible, therefore the proper initial distribution of the velocity must be obtained in some different way. The influence of the initial conditions on the further development of the system must be evaluated. The time of the relaxation of the total kinetic energy after rapid changes in the applied external forcing can serve as the criterion of such an influence. The fast decreasing of the kinetic energy takes place after ending of the wind forcing (fig. 13). The time of relaxation is approximately 1 to 1.5 days and is small comparing to the typical calculation time period of 20 to 50 days. The extreme case of the rapid change of the wind direction by 180 degrees has been investigated, too (fig. 13).

Figure 13.Changes of the total kinetic energy of the Gulf after switch-off of the wind  forcing (dashed) and after change of the wind direction by 180° (solid).

The very fast decrease of the kinetic energy determined by wind takes place at the first 4 to 5 hours. Following 3 to 6 days are necessary to reach a nearly steady state determined by the new wind direction. The total calculated time of relaxation is 5 to 7 days, that is not rather small. However, so rapid change in forcing is rarely observable and for slower change this time is smaller. Above considerations lead to conclusion that influence of the initial conditions is not large and they can be obtained from the steady state calculation done with the constant values of the wind and applied elevation in the straits equal to the initial ones of the real-time calculations.

The comparison of the computed and measured water levels at the two stations for which the observations were available (Daugavgriva and Skulte) during calculation period are shown on fig. 12 and summarized in table 5. The disagreement between average water elevations is less than for steady-state calculations, mainly due the lower average wind velocities for both periods (5.6 m/s instead of 7 m/s for steady state calculations).

Figure 12a. Water level time-development at Daugavgriva and Salacgriva stations. X-XI/93.

Figure 12b. Water level time-development at Daugavgriva and Salacgriva stations. XI-XII/94.

 However, general agreement indicated in fig. 12 (calculated results are shifted to the values from table 3) is mainly due to the rather fast propagation of the applied boundary elevations into the whole Gulf with 3-4 hours delay. Occasional events of extreme water levels are not fetched by the model. Generally, the comparison indicates the necessity to include density field in the calculation and, possible, enhancement of the detalization of the calculation mesh near the coast. There are also the indications of the necessity to account for atmosphere pressure, for instance, high atmosphere pressure (~1030 mbar), and, probably, reasonable pressure gradients during 11-16-XI-93 in combination with NW winds caused an overestimation of the water levels at both stations.

	X-XI/93	XI-XII/94
Average water level at Daugavgriva (cm)	-17	32
The same, calculated (not shifted)	-25	21
Water level fluctuations at Daugavgriva (cm)	-60/31	-7/94
The same, calculated	-51/9	-4/49
Average water level at Salacgriva (cm)	-25	26
The same, calculated (not shifted)	-25	21
Water level fluctuations at Salacgriva (cm)	-56/24	-12/68
The same, calculated	-51/12	-4/50

 Table 5. Summary of water level comparison for Daugavgriva and Salacgriva stations.

The flow patterns during the real-time calculation indicate all the spectrum of features described in the previous section including various interswitching between the described characteristic flow patterns. However, the structure of vortexes lead to interesting conclusions regarding similar spreading of water masses during time periods with different forcing as the two under consideration. The illustration of below features see fig. 14 with flow-paths of passive tracer from sources at southern part of the Gulf, Irbe strait, Northern part of the Gulf, and Pärnu bay.

Figure 14. Spreading of passive tracer from selected locations for real-time calculation  periods X-XI/93 (a) and XI-XII/94 (b).

• Due to closed circulation pathlines in the southern part of the Gulf the mixing of it with other water masses is low and water spreading northwards occurs mainly in coastal zone (with preference to the eastern coast).
• The water entering the Gulf through the Irbe strait almost never can be found in the northern part of the Gulf. Flowing around the cape of Kolka it is captured by the series of vortexes along the western coast of Gulf, between Kurzeme and Greta grund, the main sill of the Gulf. The convergence of the vortexes near Mersrags do not allow these water masses to spread in the southern part of the Gulf.
• Absence of permanent longshore currents in the Pärnu bay leads to weak water exchange of this water body with the rest of the Gulf.
• The most dynamic are the water masses of the northern part of the Gulf. Due it is more shallow, the velocities would reach higher values there. They can also left the Gulf not only through Muhu strait but also along the northern coast of Irbe strait.

Conclusions

The formulation of the shallow water model for the Gulf of Riga has been given together with the physical situation analysis in the present paper. The attempt to find out the conditions of applicability of the model without accounting of the density field is performed, too.

The shallow water model can be applied for:

• Qualitative investigations of the reaction of the Gulf to applied meteorological forcing.
• Obtaining the typical flow patterns.
• Real-time calculations for the autumn-winter periods.

Acknowledgements

This work has been in part supported by grant No. 93.311 of Latvian Science Council and the Finnish Institute of Marine Research within the framework of the ² Gulf of Riga Project² funded by the Nordic Council of Ministers. Authors owe Dr. Davide Ambrosi (CRS4, Italy) and Viesturs Berzins (Latvian Fisheries Research Institute) for many fruitful discussions and useful advices concerning, respectively, numerical realization and forcing data interpretation. We acknowledge Drs. Tarmo Kõuts (Estonian Hydrometeorological Agency) and Jevgeniy Zaharèenko (LHMA) for the access to the respective water level measurements.
 

References

1. V. Berzins. Hydrology. In: Ecosystem of the Gulf of Riga between 1920s and 1990s. Ed. Ojaveer. Tallinn, 1995.
2. V. Berzins, U. Bethers, J. Sennikovs. Gulf of Riga: bathymetric, hydrological and meteorological databases, and calculation of the water exchange. Proc. of the Latvian Academy of Sciences. (7/8):107, 1994.
3. J. Sennikovs, U. Bethers. Modelling of the vertical temperature and salinity structure of the Gulf of Riga. Latvian Journal of Physics and Technical Sciences, (1):19-41, 1995.
4. T. S. Murty, Z. Kowalik. Numerical Modelling of Ocean Dynamics. World Scientific Publishing, Singapore, 1993.
5. UNESCO 10th report of the joint panel on oceanographic tables and standards. UNESCO Tech. Papers in Mar. Sci., 36, 1981.
6. D. Ambrosi, S. Corti, V. Pennati and F. Saleri. Numerical simulation of unsteady flow at Po river delta. Journal of Hydraulic Engineering of the American Society of Civil Engineers (ASCE), vol. 122,  pp.735-743, 1996
7. J. P. Benque, J. A. Cunge, J. Fenilet, A. Hauguel, F. M. Holly. New method for tidal current computation. Journal of Waterway, Port, Coastal and Ocean Division ASCE, 108:396-417, 1982.
8. A. Quarteroni and A.Valli. Numerical approximation of partial differential equations. Springer-Verlag, 1994.
9. W. D. Joubert, D. R. Kincaid, T. C. Oppe. NSPCG User’s guide, Version 1.0. Center for Numerical Analysis The University of Texas at Austin, April 1988.
10. G. B. Witham. Linear and non-linear waves. Wiley, New-York, 1974.
11. Meteorological bulletins of the Latvian Hydrometeorological Agency. Tables TGM-1, TGM-1M. Technical report, LHMA, 1974-1995.
12. J. Pedlosky. Geophysical Fluid Dynamics. Springer-Verlag, 1982.