Aquatic Ecology 36: 143–152, 2002. © 2002 Kluwer Academic Publishers. Printed in the Netherlands.

143

Mathematical models and computer programmes for the investigation of hydrophysical processes in Lake Shira V.M. Belolipetskii, S.N. Genova, L.V. Gavrilova and L.A. Kompaniets Institute of Computational Modelling, Siberian Branch of the Russian Academy of Sciences, Academgorodok, Krasnoyarsk, 660036, Russia (Fax: (3912) 432756; E-mail: [email protected]) Accepted 22 September 2001

Key words: salinity, simulation, temperature stratification, three-dimensional model, two-dimensional in vertical section, wind current

Abstract Physical processes determine to great extent the habitat of hydrobionts, as well as the transfer and sedimentation of substances, the intensity of pollution and rate of natural purification of water bodies. Mathematical models of different levels of complexity, developed to investigate hydrophysical processes in lakes, are discussed in this paper. The numerical algorithms and computer programmes described can be used to determine the influence of the morphometric characteristics and weather on the temperature regime of Lake Shira (Khakasia, Siberia) and the pattern of the wind currents. Examples of calculating the temperature regime in the context of a one- dimensional model and of calculating the parameters of wind currents in the water body of a simple geometrical form are given. The study suggests that the pattern of wind currents in Lake Shira is significantly affected by the density stratification, which depends not only on temperature but also on salinity. In order to construct a realistic pattern of currents a 3-D computer model of Lake Shira must be developed and used to estimate the validity of the two-dimensional and one-dimensional models. The present research can be further developed by extending the obtained algorithms to three-dimensional problems, taking into account the heat exchange, salinity and the geometry (bathymetric contours) of the water body. Calculations can be made for Lake Shira and the obtained data used in biophysical models.

Introduction Because of pollution and other anthropogenic impacts, the quality of water in inland lakes becomes much worse. Several factors and processes determine the ecological state of a water body, such as meteorological factors, hydrophysical processes, hydrochemical and hydrobiological factors, and external anthropogenic inputs. The physical parameters, e.g. temperature, govern the chemical and hydrobiological processes. The self- purification of a natural water body is largely determined by the temperature regime. Thermally stratified water currents have been the subjects of numerous studies. In some papers from the estuarine habitats, the water currents were explored for temperature, in addition to salinity (Khublarian, 1988). In other studies, temperature-stratified flows

were studied in reservoirs (Bocharov, 1994; Belolipetskii et al., 1994) including cooling-water reservoirs (Makarov, 1986) and lakes (Astrakhantsev et al., 1986; Belolipetskii, 2001, Dobrovolskaya et al., 1981; Dobrovolskaya & Simonov, 1986; Kompaniets & Gavrilova, 1999; Wang & Hutter, 2000). In all these studies, equations for the hydrodynamics and heat – mass transfer of different simplifications and different methods and numerical solutions were used for mathematical modelling of stratified flows. Earlier researchers investigated the mathematical modelling of hydro-thermodynamics processes in Lake Shira. This present study deals with the mathematical simulation of the hydrophysical processes in Lake Shira, including the lake’s thermal structure, wind and wind induced currents, salinity, etc. (see Belolipetskii, 2001). Mathematical models and com-

144 puter programmes of different levels of complexity can be used to investigate the hydrophysical processes in stratified water bodies. The one- dimensional model is used to determine the vertical distribution of temperature. The hydrostatic and non-hydrostatic models help to simulate the influence of wind and stratification on the pattern of water currents. The aim of the present investigation is to study a possibility of adaptation of the concrete onedimensional and two-dimensional, in vertical and horizontal planes, mathematical models for different problems relating to Lake Shira, and also to determine potential patterns of current in Lake Shira.

and stratification of a water body. The stratified fluid is a fluid the physical properties of which, when it is in a steady state, can change only in specific directions. Stratification can be due to various physical reasons. The most frequent is gravity, which causes the heterogeneity of liquid, termed as density stratification. The density stratification of the lake is associated with the non-uniform distribution of water temperature and salinity. The temperature and salinity fields are related to currents and turbulent mixing. The distributions of temperature and of salinity and water currents are interrelated. Therefore, when these processes are described mathematically, the motion equations and the heat-transport and salt-transport equations must be solved together.

Materials and methods Choice of models The lake The lake is located in Siberia (54◦ 29 N and 90◦14 E) and is an inland lake, without islands. The lake is saline and exhibits unique medicinal properties. The Son River flows into the lake from the south. As the inflow is very little, the river only influences the part of the lake near its mouth. The bathymetry of the lake is presented in Figure 1. The lake has an oval shape and is 9.4 km long and about 5 km wide. The surface area of the lake is 34.7 km2 and the maximal depth is ca. 24 m. The underground water input constitutes about 9% of the total water input of the lake. The temperature regime of the lake depends on weather conditions. The air temperatures at the shore and at the nearest weather station clearly differ (Figure 2), and this must be taken into account in the calculations. The water temperatures are much higher in the epilimnion (15–20 ◦ C) than in the hypolimnion (2–3 ◦ C). The salinity is higher in the hypolimnion. The summer stratification is very stable because of the temperature stratification and the high mineral content of the lake. The distinct density stratification is defined by both the decrease of the water temperature, and increase of salinity, with the lake’s depth. Because the patterns of the wind currents depend greatly on the stability of the stratification, they are expected to differ greatly in periods of weak and strong stratification. The factors affecting stratification In closed-basin water bodies like Lake Shira the dynamic processes are largely determined by the wind stress. The pattern of wind currents depends on the the direction and force of wind, geometry and the depth

To choose an adequate mathematical model, models of different levels of complexity have been compared. The one-dimensional model was used to determine the vertical distribution of temperature. The calculations show that the plane shallow-water model is inapplicable to Lake Shira and the density stratification must be taken into account. Hydrostatic and non-hydrostatic models have been used in the calculation of two-dimensional currents in a vertical section (the cross-section of Lake Shira). The parameterisation of the coefficient of vertical turbulent exchange was applied in accordance with the Prandtl-Obukhov formula. The influence of the wind force, stratification and geometry of the water body on the pattern of the current was evaluated. The test calculations of three-dimensional currents in water bodies of simple geometry (in the context of the hydrostatic model) have shown that in strongly stratified water bodies, in contrast to weakly stratified ones, the following three layers can be differentiated: the upper quasi-homogeneous layer, the lower layer, and the intermediate layer. The patterns of currents in these layers differ significantly. Method 1. The simplified one-dimensional model of the temperature regime in a closed-basin water body The temperature regime in a closed-basin stratified water body is formed by wind currents and heat exchange with the atmosphere. For the modelling of the vertical temperature stratification we consider the simplified one-dimensional statement. The distribution of the temperature over depth in the water body

145

Figure 1. The bathymetric map of Lake Shira (depth indicated in meters).

is determined based on the solution of the problem (Belolipetskii, 2001):
 ∂ ∂T ∂T F1 e−βz = , Kz +β ∂t ∂z ∂t cp ρ0 Kz

Fn ∂T =− ∂z cp ρ0

∂T =0 ∂z

for z = 0, (1) for z = H,

where t is time, z is the vertical coordinate directed downwards, T is the temperature of the water, Kz (z) is the vertical turbulent exchange coefficient, Fn is the entire heat flow through the free surface, FI is the short-wave solar radiation, β is the radiation absorption coefficient, cp is the specific heat capacity of water, ρ0 is the characteristic water density, H is the depth. Parameterisation of the coefficient of vertical turbulent exchange, heat flow and wind stress Turbulence considerably influences the heat and mass transfer. The semi-empirical models of turbulence are used to solve a wide range of problems. The most widespread semi-empirical model is the k-ε model of turbulence (Rodi, 1984). The problem becomes more complicated when the hydrobiological and hydrophysical processes are investigated concurrently. Thus, simpler parameterisations of the coefficient of vertical turbulent exchange should be used. Here we consider some of them. The

Prandtl-Obukhov formula (Belolipetskii & Genova, 1997) is: 
 
2  √ ∂u 2 ∂v   (0.05h)2 + 1 − Ri ∂z ∂z Kz =   Kmin for Kz > Kmin , (2) for Kz ≤ Kmin , 

2  g ∂ρ ∂v ∂u 2 Ri = , + ρ ∂z ∂z ∂z where u and v are the horizontal components of the wind current velocity, ρ is water density, Ri is the Richardson number, Kmin is the background value of the coefficient of the vertical turbulent exchange (determined from the field data), h is the thickness of the quasi-homogenous layer, determined by the first calculation point from the surface where h = zk .


2
2
  ∂v g ∂ρ 2 ∂u (0.05zk ) + − ∂z |zk ∂z |zk ρ0 ∂z |zk ≤ Kmin . To evaluate the coefficient Kz we define the vertical gradients of velocity from the approximate Eckman solution for wind currents

146

Figure 2. Air temperature in 1999.

uek =

v ek

e−αx [(τX + τY ) cos(αz)+ 2ρ0 αK0

+(τy − τx ) sin(αz) , (3)

e−αx  (τy − τx ) cos(αz)+ = 2ρ0 αK0

+(τx + τy ) sin(αz) ,

where τx , τy are the components of the shear friction exerted by the wind on the water surface, K0 is the characteristic value of √ the coefficient of vertical turbulent exchange, α = f/K0 . From (3) we find 2
ek 2
ek 2
∂u ∂v τ + = e−2αz , (4) ∂z ∂z ρ0 K0  τx2 + τy2 . The depth of the upper where τ = quasi-homogenous layer can be evaluated from the following condition: the thickness of the layer where K0 ((∂u/∂z)2 + (∂v/∂z)2 ) decreases eπ -fold is equal to the thickness of the upper quasi-homogenous layer h1 . Then from (4) it is clear that  K0 h1 = π . 2f Supposing K0 = Kz (0) we find K0 =



(0.05π)2τ

ρ0 4f 2 + (0.05π)4gρz|z=0

.

For Kz the following formula is obtained 
2 g ∂ρ τ 2 . (5) e−2αz − Kz = (0.05h1) ρ0 K0 ρ0 ∂z Munk & Anderson (1948) have proposed a correction that takes into account the stratification of the form Khl Kz = √ , (6) 1 + 3.33 Ri where Khl is the coefficient of the vertical turbulent exchange for a homogeneous fluid. For the definition of Khl , formula (2) (for Ri=0) is used, so it is clear from formula (6) that:    ∂u  (7) Kz = (0.05h)2   (1 + 3.33 Ri)−1/2 . ∂z In the methods described, the velocity gradient and the stratification define the intensity of the vertical turbulent exchange. At the first stage the simplest variant of the equation for the state of the salt water is used in the form: ρ = ρ0 · (1 − 1.5 × 10−4 T ).

(8)

In the complete problem statement it is necessary to calculate the water salinity S and use the equation of state ρ = ρ(S, T ) (Gill, 1982):

147 

ρ(S, T ) = ρw + S 0.824493 − 4.0899 × 10−3 T + +7.6438 × 10−5 T

− 8.2467 × 10−7 T 3 +

+5.3875 × 10−9 T 4 + S 3/2 (−5.72466×

×10−3 + 1.0227 × 10−4T − 1.6546×  (9) ×10−6 T 2 + 4.8314 × 10−4 S 2 , where ρw is defined by the simplified formula ρw = [1 − 0.68 × 10−5 (T − 4)2 ] × 103 . An important parameter influencing the temperature regime of the water body is heat flow. The total heat flow through the free surface is calculated from the well-known relations (Belolipetskii et al., 1994; Belolipetskii & Genova, 1997). The wind stress is determined by Davtyan’s formula:     τ = ρa (0.9 + 0.17 W 2 ) × 10−3 W 2  W 2 , where ρa is the air density, W 2 = (wx , wy ) is the vector of the wind velocity, and the index 2 corresponds to the height at which the wind velocity is measured −2 m. Method 2. The numerical algorithm for the calculation of wind currents in stratified water bodies. The two-dimensional non-hydrostatic model in a vertical section Here we propose the numerical algorithm for the calculation of wind currents in closed-basin water bodies. The basic assumptions are as follows: slow currents; rigid lid approximation; the approximation of the boundary layer; and Boussinesq’s approximation. Taking into account these assumptions, we have 1 ∂p ∂ ∂u ∂u =− + Kz , ∂t ρ0 ∂x ∂z ∂z 1 ∂p ∂ ρ ∂w ∂w =− + Kz + g, ∂x ρ0 ∂z ∂z ∂z ρ0

(10)

∂u ∂w + = 0, ∂x ∂z ∂T ∂T ∂T ∂ ∂T +u +w = Kz , ∂t ∂x ∂z ∂z ∂z ∂Ci ∂Ci ∂Ci ∂Ci ∂ (11) + + (w + wgi ) Kz + fc , ∂t ∂x ∂z ∂z ∂z ρ = ρ(T , Ci ), where u, w are the components of the velocity of water current in x, z directions respectively, wgi is the

fall velocity of the pollutant, p is pressure; Ci is the concentration of the ith pollutant, fc comprises the stirring-up and sedimentation, g is gravitational acceleration. The initial-boundary problem is solved. The initial conditions as well as the boundary conditions are given: on the free surface (z = 0) τx ∂u =− , w = 0, Kz ∂z ρ0 K

∂T Fn =− , ∂z cp ρ0

Kz

∂Ci = −qC0 i , ∂z

on the rigid wall (z = H ): FH ∂T = u = 0, w = 0, Kz , ∂z cp ρ0

Kz

∂Ci = qCHi . ∂z

Here FH is the heat exchange with the bottom of the water body; qC0 i and qCHi are pollutant flows on the free surface and on the bottom respectively. The conditions of the non-passage and the absence of heat transfer are set on the lateral surface. Here we consider the algorithm for the solution of the stated problem for the variable coefficient of vertical turbulent exchange using the stream-function and vortex as the new variables. The numerical algorithm for the solution of the problem involves three stages (Belolipetskii, 2001). At the first stage, the vortex of velocity + is determined, and the semi-implicit method is used for the time sampling:
 ∂ ∂+n+1 g ∂ρ +n+1 − +n = + Kz + ,t ∂z ∂z ρ0 ∂x
 ∂Kzn ∂y n ∂ ∂Kzn ∂wn + · − · , ∂z ∂x ∂z ∂z ∂z where t is the time step, and the time layers are denoted by the indexes n, n + 1, ωn+1 = +(tn+1 , x, z), Kzn = Kz (tn , x, z), un = = u(tn , x, z, ), wn = w(tn , x, z), tn+1 = = tn + ,t. At the second stage, the stream-function (ψ) and the velocity components (u, w) are obtained using the derived vortex: ∂ψ ∂ 2ψ ∂ψ ∂ 2ψ , w=− . + = −+n+1 , u = 2 ∂x ∂z2 ∂z ∂x At the third stage, the processes of heat and pollutant transfer in the system are simulated (9). The numerical realisation of this problem involves the method of splitting into physical processes, the

148 sweep method, and the method of successive overrelaxation. Method 3. The numerical algorithms for the investigation of wind currents in stratified water bodies. Three-dimensional currents. Hydrostatic model The three-dimensional mathematical model of the hydro-thermodynamics of a water body in Boussinesq’s approximation of slow currents and in the approximations of hydrostatics and boundary layer has been used for Lake Shira (Astrakhantsev et al., 1986). The velocities calculated taking the Coriolis force into account are approximately equal to those calculated ignoring it (Kompaniets & Gavrilova, 1999). That is why the basic equations can be written disregarding the Coriolis force as
 z ∂ g ∂u ∂u ∂ρ ∂η = − dz, Kz +g ∂t ∂z ∂z ∂x ρ0 ∂x 0


 z ∂v ∂ ∂η g ∂v ∂ρ (12) Kz +g = − dz, ∂t ∂z ∂z ∂y ρ0 ∂y 0

∂v ∂w ∂y + + = 0, ∂x ∂y ∂z where u, v, w are the components of the vector of the current velocity, η is the elevation of the free surface. The boundary conditions for the equation of motion are set on the undisturbed surface (z = 0): τy τx ∂v ∂u ∂η = − , Kz =− , = w, Kz = ∂z ρ0 ∂z ρ0 ∂t and on the bottom of the lake (z = H (x, y)):     ∂u ∂v Kz u V  , Kz = −k1 = −k1 v V  , ∂z ∂z w = −u

∂H ∂H −v , ∂x ∂y

H u(x, y, z, t) dz, 0

 v (x, y, t) =

1 H (x, y)

∂T cos(n, z) = 0, ∂z on the water-body surface (z = 0) Fn ∂T =− . Kz ∂z cn ρ0 Kz

Here Ax , Ay , are the coefficients of horizontal turbulent diffusion. At the initial moment of time the velocity distribution, the elevation of the free surface and the temperature distribution are set. Integrating Equations (10) from z = 0 to z = H (x, y) we obtain the system of shallow water equations for the determination of the depth averaged velocities and of the elevation of the free surface  ∂η τx u|V | ∂u −g = −k1 + + F1 , ∂t ∂x H ρH τy ∂η  v |V | ∂v −g = −k1 + + F2 , ∂t ∂y H ρ0 H ∂η ∂ ∂ + ( uH ) + ( v H ) = 0, ∂t ∂x ∂y F1 = −

g ρ0 H0

H

H v(x, y, z, t) dz. 0

The condition of sliding is set on the lateral surface. The equation of heat transfer has the form

z dz

0

F2 = −

g ρ0 H0

(14)

∂ρ dz, ∂x

0

H

z dz

0

where k1 is the non-dimensional coefficient and V = ( u, v ) is the vector of depth averaged velocities 1  u(x, y, t) = H (x, y)

∂ 2T ∂T ∂T ∂T ∂T +u +v +w = Ax 2 + ∂t ∂x ∂y ∂z ∂x (13)
 2 ∂ T ∂T ∂ Kz . +Ay 2 + ∂z ∂z ∂y It has the following boundary conditions: at the rigid boundaries ∂T ∂T Ax cos(n, x) + Ay cos(n, y)+ ∂x ∂y

∂ρ dz. ∂y

0

For the homogeneous fluid F1 = F2 = 0. The condition of non-passage is set on the lateral surface. At the initial moment of time the initial distribution of depth averaged velocities, the elevation of the free surface and the initial temperature distribution are set. The numerical algorithm for the solution of the three-dimensional problem is based on the three-stage method. At the first stage, the horizontal constituents of the depth-averaged components of velocity and the elevation of the free surface are calculated based on the set of Equations (14), taking into account the boundary conditions. At the second stage, u, v, w are

149 calculated from the obtained values of the average velocity and the level of the surface of Equations (12). Finally, at the third stage, the new fields of temperature and density are found from the equation of heat transfer and the equation of state. The numerical solutions are found by the method of finite differences. The explicit conditionally stable McCormack difference scheme and the implicit unconditionally stable difference scheme are used to solve the equation for inhomogeneous shallow water (14).

Discussion The definition of the hydrophysical regime of a water body is an important multi-parametric problem that can be solved by different methods. Lake Shira is a significantly density-stratified water body. This feature must be taken into account in choosing an adequate mathematical model. Several approaches to modelling the water flow in a water body have been reported in the literature. In some well-known papers, models of the stratified currents involve the temperature stratification only. To illustrate, the two- dimensional, in a vertical section (x, z), model with a two-layer temperature stratification has been used for elongated water bodies (Dobrovolskaya et al., 1981; Dobrovolskaya & Simonov, 1986). The layer interface is in the thermocline, which is taken as infinitesimal, with the density of every layer as taken equal. The twodimensional XZ-models involving temperature stratification are also often used to model the stratified currents in water bodies (Lavelle et al., 1991; Edinger & Buchak, 1983; Wang & Kravitz, 1980; Wang, 1985; Dunbar & Burling, 1987). The so-called equations of shallow water for the uniform fluid, which are twodimensional in a horizontal section, are used when the length of the water body significantly exceeds its depth (Lake Onega, Caspian Sea, Chudskoe Lake). The depth-averaged velocities are used to determine the horizontal circulation of water (Baklanovskaya et al., 1979; Dobrovolskaya et al., 1981; Andreev & Vorobjeva, 1986; Perminov & Chechel, 1986). However, the three-dimensional models of stratified currents must be used to determine real wind currents. Dobrovolskaya et al. (1981); Volkova et al. (1981); Astrakhancev et al. (1986), Dobrovolskaya & Simonov (1986), Wang & Hutter (2000) use simplified mathematical 3-D models of wind currents in Boussinesq’s assumption, and in the assumptions of hydrostatics. Of considerable interest is the work by

Wang & Hutter (2000), in which the three-dimensional model is used to study the wind-induced transfer of admixtures in temperature-stratified lakes. The examples of wind currents are given for stratified Lake Constance (length – 64 km, width – 16 km, maximal depth – 250 m). The coefficients of turbulent exchange used in the calculation are in the form of partly constant distribution over depth. The zones of water circulation depending on temperature stratification of the water body are determined for the fixed wind, but the influence of the wind force on the pattern of the current has not been studied. Having analysed the above works, we have chosen several models to simulate the wind currents in Lake Shira: the one-dimensional model, the model for shallow water, and two models for the modelling of wind currents in the vertical section: the hydrostatic model (Kompaniets & Gavrilova, 1999) and the non-hydrostatic model (Belolipetskii, 2001). The parameters and variables used in this paper are listed in Table 1. Based on the detailed weather data, the calculations of the thermocline in Lake Shira were conducted using the one-dimensional model (1). The comparison of the calculation data with the observation evidence (Figure 3) suggests that in the deep-water area the temperature profile can be determined using the one-dimensional approximation. The influence of the wind is taken into account in determining the coefficient of vertical turbulent exchange according to formula (5). The basic calculations were made using the twodimensional models in a vertical section. They involved a simplest schematisation of Lake Shira as a rectangle, 5000 m long and 21 m deep, for the temperature stratification (Figures 4b and 5b). Coefficient Kz was defined by formula (7). Figures 4a and 5a show the isolines of the stream-function for the wind 5 m s−1 . For the summer temperature profile (Figure 4b) ,ρ/ρ0 ≈ 2 × 10−3 , according to formula (8) and ,ρ/ρ0 ≈ 4.3 × 10−3 , according to formula (9) (salinity varies from 16 to 18 g l−1 , ,ρ = ρmax − ρmin ). The two-dimensional, in the vertical section, non-hydrostatic model (10–11) was used to make the calculations (Figure 4a) for the summer temperature profile. Two circulation areas were formed for the characteristic summer stratification (Figure 4a), whether or not salinity was taken into account. In spring and autumn, when variations in the temperature of the water over depth are insignificant, the variations in the density over depth in Lake Shira depend on the vertical distribution of salinity. When the

150 Table 1. Basic variables and parameters Ax , Ay (m2 s−1 ) cp (kcal kg−1 ◦ C−1 ) g = 9.81 (m s−2 ) f = 2ω sin ϕ (s−1 ) Fn (kcal m−2 s−1 ) H (m) Kmin (m2 s−1 ) k1 = g C−2 Kz (m2 s−1 ) P (kg m−1 s−2 ) T (◦ C) u, v, w (m s−1 ) wgi (m s−1 ) wx , wy (m s−1 ) η (m) ρ (kg m−3 ) ρ0 = 103 (kg m−3 ) ρa = 1.3 (kg m−3 ) τx , τy (kg m−1 s−2 ) ψ (m2 s−1 ) + (s−1 )

Coefficients of horizontal turbulent diffusion Water heat capacity Gravity constant Coriolis parameter, where ω is the angular velocity of the earth and ϕ is the geographical latitude Total heat flow through the free surface Depth of water body Background value of the coefficient of the vertical turbulent exchange (≈ 2 × 10−6 –2 × 10−4 ) Bottom drag coefficient, where C = H 1/6 /km , km = 0.02 is Manning’s coefficient Coefficient of vertical turbulent exchange Air pressure Water temperature Components of the water current in the x-, y-, z-directions Fall velocity Horizontal components of the wind 2 m above water surface Free surface elevation Density of water Density of water at 4 ◦ C Density of air Horizontal components of the shear frictions exerted by the wind on the water surface Stream-function Vortex

Figure 3. The temperature distribution over depth. Observed and calculated profiles.

temperature distribution follows the pattern shown in Figure 5b, ,ρ/ρ0 ≈ 0.6 × 10−3 and ,ρ/ρ0 ≈ 1.9 × 10−3 , according to formula (9) (salinity varies from 16 to 18 g l−1 ). For the weak stratification (Figure 5a) one circulation area is formed if density is calculated according to (8) and two circulation areas are formed if salinity is taken into account, according to (9) (in calculations the salinity of the water varies over depth from 16 g l−1 to 18 g l−1 ). Hence, even small relative variations in density (of the order of 10−3 ) lead to a change in the pattern of the current. Similar calculations were carried out using the two-dimensional analogue of the three-dimensional hydrostatic model (12)–(14) (with formula (8)). The results of calculations for the temperature stratification (Figure 4b) show that if the wind velocity is |W | < 8 m s−1 , then these results agree well with those obtained using the non-hydrostatic model. But if |W | > 8 m s−1 , then we have current velocities differing over depth, and the transition from the twocirculation current to the one-circulation current is realised if |W | is 9 m s−1 for the hydrostatic model and if |W | is 10.5 m s−1 for the non-hydrostatic model. If

151

Figure 4. The diagram of stream-lines (a) for the initial temperature (b).

Figure 5. The diagram of stream lines (a) for initial temperature (b).

salinity is taken into account, the wind force for the realisation of the one-circulation current is 1–1.5 m s−1 greater. The bathymetry also influences the formation of the wind-induced current. If the profile of the bottom is not rectangular, one circulation area is formed at a stronger wind than in the case with a rectangular water body (calculated using the non- hydrostatic model). Our calculations show that the shallow-water model cannot be applied to Lake Shira over the entire depth. However, when a two-circulation zone is formed in the lake, the shallow-water equation may possibly be used to determine the current in the layer above the layer where the density abruptly increases. The test calculations of three-dimensional currents (in the context of the hydrostatic model) have shown that in stratified water bodies the following three layers can be differentiated: the upper quasi-homogeneous layer, theintermediate layer and the lower layer. The patterns of currents in these layers differ significantly. That is why investigations into the spatial-temporal dynamics of the components of an ecosystem must take into account the qualitative difference in the currents in the study areas. In the 3-D problem the effects of the wind force, stratification and the profile of the lake

bottom may become even more significant. For example, calculations carried out using the 3-D model show the emergence of several sub-circulation water zones within the main circulation. In the work of Astrakhantsev et al. (1986) the model was even successfully applied to a weakly stratified Lake Ladozhskoye. Thus, the pattern of wind currents in Lake Shira is significantly affected by the density stratification, which depends not only on temperature but also on salinity. In order to construct a realistic pattern of currents, a 3-D computer model of Lake Shira must be developed and used to estimate the validity of the two-dimensional and one-dimensional models. This research can be further developed by extending the obtained algorithms to three-dimensional problems, taking into account heat exchange, salinity and geometry of the water body. Calculations can be made for Lake Shira and the data obtained used in biophysical models.

Acknowledgements The work was done under the support of the Ministry of Education of Russia and CRDF (grant No. REC-

152 002), RFFR (No. 99-05-64695), INTAS (97-O-519), FSP ‘Integration’ (project No. 73).

References Andreev OA and Vorobjeva LV (1986) Wind and drain circulation of the Ladoga Lake (numerical experiments). In: Modelling and Experimental Investigation of Hydrological Processes in Lakes. Nauka, Leningrad, pp. 17–21 Astrakhantsev GP, Oganesyan LA and Ruhovec LA (1986) Threedimensional mathematical model of hydrothermic of the closed reservoir. In: Modelling and Experimental Investigation of Hydrological Processes in Lakes. Nauka, Leningrad, pp. 13–16 Baklanovskaya VF, Palcev BV and Chechel II (1979) On boundary problems for Saint-Venant system of equations on a plane. JCM&MPh 19: 708–725 Belolipetskii VM (2001) Numerical modelling of wind currents in stratified reservoirs. Water Res 2: 133–137 Belolipetskii VM and Genova SN (1997) Numerical modelling of hydrothermal processes in well-drained basins. Russ J Numer Anal Math Modelling 12: 319–334 Belolipetskii VM, Genova SN, Tugovikov VB and Shokin YuI (1994) Numerical Modelling of Problems of hydro-Ice-Thermics of Currents. SB RAS, Novosibirsk Bocharov OB, Vasiliev OF and Ovchinnikova TE (1994) Twodimensional vertical model of hydrothermical processes in temperature-stratified long length reservoir. Proc Acad Sc 339: 327–330 Dobrovolskaya ZN, Epichov GP, Korjavov PP and Moiseev NN (1981) Mathematical models for calculation of dynamics and quality of complex water systems. Water Res 3: 32–51 Dobrovolskaya ZN and Simonov AI (1986) Mathematical modelling of the currents in stratified reservoir. In: Modelling and Experimental Investigations of Hydrological Processes in Lakes. Nauka, Leningrad, pp. 6–10 Dunbar DS and Burling RW (1987) A numerical model of stratified circulation in Indian Arm. British Columbia. J Geophys Res 92: 13075–13105

Edinger JE and Buchak EM (1983) Developments in LARM2: A longitudinal-vertical, time varying hydrodynamic Reservoir Model. In: Technical Report E-83-I, US Army Engineer Waterways Experimental Station. Vicsburg, Miss Gill AE (1982) Atmosphere-Ocean Dynamics. Vol. 2. Academic Press, New York, pp. 256–278. Khublarian MG and Frolov AP (1988) Modeling of Intrusion Processes in Estuary and Underground Water-Bearing Horizonts. Nauka, Moscow Kompaniets LA and Gavrilova LV (1999) Numerical algorithm for calculation of the wind currents of the stratified liquid. In: Dynamics and Thermics of Rivers, Reservoirs and Littoral Region. Moscow, pp. 51–54 Lavelle JM, Cokelet ED and Cannon GA (1991) A model study of density intrusion into and circulation within a deep, silled estuary: Puget Sound. J Geogr Res 96: 16779–16800 Makarov II, Sokolov AS and Shulman SG (1986) Modeling of Hydrothermics Processes of Water Reservoir-Cooler of Heat and Atomic Power Plants. Energoatomizdat, Moscow Munk WH and Anderson ER (1948) Notes on the theory of the thermocline. J Marine Res 7: 25–29 Perminov SM and Chechel II (1986) The investigation of flows, horizontal water circulation and water exchange between the eastern and western parts of Northern Caspian. In: Hydrophysics of North Caspian. Nauka, Moscow, pp. 6–10 Rodi W (1984) Models of turbulence of environment. In: Prediction Methods for Turbulent Flows. Mir, Moscow Volkova GB, Kvon VI and Filatova TN (1981) Numerical modelling of wind flows in Chudskoe Lake. Water Res 3: 91–99 Wang DP (1985) Numerical study of gravity currents in a channel. J Phys Oceanogr 15: 299–305 Wang DP and Kravitz DW (1980) A semi-implant two-dimensional model of estuarine circulation. J Phys Oceanogr 10: 441–454 Wang Y and Hutter K (2000) Methods of substructuring in lake circulation dynamics. Adv Water Res 23: 399–425