Barotropic response in a lake to wind-forcing

. We report results gained with a three-dimensional, semi-implicit, semi-spectral model of the shallow water equations on the rotating Earth that allowed one to compute the wind-induced motion in lakes. The barotropic response to unidirectional, uniform winds, Heaviside in time, is de-termined in a rectangular basin with constant depth, and in Lake Constance, for different values and vertical distributions of the vertical eddy viscosities. It is computationally demonstrated that both the transitory oscillating, as well as the steady state current distribution, depends strongly upon the absolute value and vertical shape of the vertical eddy viscosity. In particular, the excitation and attenuation in time of the inertial waves, the structure of the Ekman spiral, the thickness of the Ekman layer, and the exact distribution and magnitude of the upwelling and downwelling zones are all signiﬁcantly affected by the eddy viscosities. Observations indicate that the eddy viscosities must be sufﬁciently small so that the oscillatory behaviour can be adequately modelled. Comparison of the measured current-time series at depth in one position of Lake Constance with those computed on the basis of the measured wind demonstrates fair agreement, including the rotation-induced inertial oscillation.


Introduction
Knowledge of water movements is a prerequisite for the study of a multitude of water quality problems of natural and artificial lakes. The computation of the current distribution in homogeneous lakes is commonly performed with the shallow water equations on the rotating Earth. In the present study, density variations are ignored (that is the case for most of the Alpine lakes in winter), so that the wind-induced motion is Correspondence to: Y. Wang (wang@mechanik.tu-darmstadt.de) restricted to the barotropic component for which the internal hydrostatic pressure gradient is exclusively due to the variation of the free surface over its equilibrium value. It is common knowledge that numerical codes integrating these spatially three-dimensional shallow water equations, explicitly in time, are conditionally stable; due to the explicit treatment of the friction (or diffusion) terms, the time step is limited by the smallest mesh size. For most problems of ocean dynamics, since space scales are usually large, the maximum usable time steps can equally be relatively large, and thus, time integration over physically relevant intervals is often possible. Since lakes are small, sufficient spatial resolution requires small mesh sizes limiting the economically justifiable integration time to values below the physically relevant intervals. Thus, implicit or semi-implicit integration routines are needed. For the latter, time steps are still limited by the value of the eddy viscosity and the space steps in the horizontal direction; but this restriction is not severe because the horizontal mesh sizes are usually fairly large. That is why, often, only semi-implicit schemes are applied instead of fully implicit or ADI schemes. There are also many ocean models which include implicit schemes for good reasons, especially in coastal ocean models or when a time-and space-dependent eddy viscosity is computed. Various implicit models have been developed and used by many researchers (e.g. Backhaus, 1983Backhaus, , 1985Backhaus, , 1987Blumberg and Mellor, 1987;Davies, 1987;Bleck et al., 1992;Davies and Lawrence, 1994;Ip and Lynch, 1994;Song and Haidvogel, 1994).
On the other hand, with regard to non-linear advection terms, if the second-order centered finite difference scheme is used (or other traditional high-order schemes, e.g. the spectral method used here), numerical oscillations always occur. The smaller the eddy viscosity is, the smaller the mesh size and time step for stable integration must be. Thus, less numerical oscillations are provoked, so that the small eddy viscosity can still assure stable numerical simulations. The attenuation with time of any oscillating motion increases with growing eddy viscosity. It follows that amplitude and persistence of oscillations are directly tied to the eddy viscosities. 368 Y. Wang et al.: Barotropic response in a lake to wind-forcing Values needed for stable integration are often larger than is physically justified, i.e. computed physical oscillations are faster attenuated than in reality. A study how the vertical eddy viscosity affects the time dependent and steady state wind-induced velocities is, therefore, urgent.
Wind-induced barotropic responses in lakes or the ocean have been studied by many researchers (e.g. Mortimer, 1974Mortimer, , 1980Simons, 1980;Hutter, 1982;Heaps, 1984;Sündermann, 1984). However, for most of these works, if the threedimensional field equations including the nonlinear advection terms are numerically solved, the inertial waves seem to be hardly observable in the numerical results, or they remain only for a fairly short time and then are rapidly damped out. This fact is due to the much larger eddy viscosities used in the models in order to restrain the numerical (not physical) oscillations and thus, assure numerical stability.
We employ the three-dimensional, hydrodynamic, semiimplicit, semi-spectral model SPEM, as originally developed by Haidvogel et al. (1991) and extended for semi-implicit integration in time by Wang and Hutter (1998). This software is able to cope with the above questions in a reasonable way. We propose three functional relations for the vertical distribution of the vertical eddy viscosity: (i) constant, (ii) sinusoidal half wave in an upper turbulently active layer extended to a constant value below and (iii) a linear increase with depth in the upper layer to a maximum continued by a constant value at this maximum below it. All these distributions are based on well justified arguments; here we use these distributions and vary their intensity.
We show that the offset of the surface current to the right of the wind, the Ekman spiral and the thickness of the Ekman layer depend strongly upon both the absolute value and the distribution of the vertical eddy viscosity, as do the windinduced inertial oscillations, the stored kinetic energy and the details in the up-and downwelling distributions of the motion. It is further demonstrated that the numerical values of the eddy viscosities need to be sufficiently small, i.e., at levels of their physical counterparts if the observed dynamics is to be reproduced by the computations. This demonstrates that robust software is needed if the barotropic circulation dynamics is to be adequately modelled. This paper, together with the previous work relating to the baroclinic circulation dynamics in lakes , makes up a complete pattern of forced motion response in enclosed lakes.
The paper is structured as follows: In Sect. 2, the numerical method is briefly introduced. Results are discussed in Sect. 3 for a rectangular basin with a constant depth, and in Sect. 4, for Lake Constance. Comparison of a series of measured results with the numerical simulation of that case is performed in Sect. 5. We conclude with a summary and some remarks in Sect. 6.

Numerical method
The balance law of mass and the balances of linear momentum form the hydrodynamic field equations for the consid-ered fluid system. We apply these hydrodynamic equations in the shallow water approximation, with the Coriolis term and the hydrostatic pressure equation implemented. Thus, the field equations read Here a Cartesian coordinate system (x, y, z) has been used; (x, y) are horizontal, and z is vertically upwards, against the direction of gravity. The field variables (u, v, w) are the components of the velocity vector v in the (x, y, z) directions, φ is the dynamic pressure φ = p/ρ with pressure p and density ρ, f is the Coriolis parameter and g is the gravity acceleration. Eddy viscosities are taken into account by differentiating between the horizontal and vertical directions, ν h and ν v , respectively; thus, the anisotropy effects of the turbulent intensity are considered accordingly.
Wind stress forcing at the water surface, free-slip lateral boundary conditions, and linear bottom friction are used. Along all boundaries the normal component of the current is set equal to zero, which, at the free surface, corresponds to the rigid-lid approximation. With this method, however, a Poisson equation must be solved at every time step to obtain a pressure field which ensures ∇ · (hv) = 0, where h indicates water depth andv is the depth-averaged velocity (see Haidvogel et al. (1991) and Wang (1996) for details).
A semi-spectral model was designed with semi-implicit integration in time to solve the system of differential equations numerically. The model is based on the semi-spectral model SPEM developed by Haidvogel et al. (1991), in which the variation of the field variables in the vertical direction is accounted for by a superposition of Chebyshev polynomials, whereas a centered finite difference discretization on a staggered Arakawa grid is employed in the horizontal direction; it was extended by Wang and Hutter (1998) to account for implicit temporal integration. Due to the small water depths of lakes, in comparison to the ocean, the original SPEM model had to be altered to permit economically justifiable time steps in the computation of the circulation of a lake. In Wang and Hutter (1998) several finite difference schemes, implicit in time, were introduced; that scheme which uses implicit integration in time for the viscous terms in the vertical direction was the most successful one.
For a realistic irregular lake, a vertical topography following the σ -coordinate, and a horizontal shore following the curvilinear coordinate system is employed. By using the so- FIG. 1: Rotation angle Φ between wind and surface current (a) and absolute fluid velocity V surf at the free surface (b). Here Φ and V surf are plotted against ν v (solid curves), A (broken curves), A 0 (dotted curves), respectively, corresponding to the cases (i), (ii) and (   Here and V surf are plotted against ν v (solid curves), A (broken curves), A 0 (dotted curves), respectively, corresponding to the cases (i), (ii) and (iii).
called σ -transformation, a lake domain with varying topography is transformed to a new domain with constant depth, and this region is once again transformed in the horizontal coordinates by using conformal mapping, which maps the shore, as far as possible, onto a rectangle. Uniformity in grid size distribution is intended because numerical oscillations (instabilities) occur preferably on the small scales; however, it is difficult to achieve it in complex geometries. In such cases, to attain an uniform grid as far as possible, a bounding line, which deviates in some segments from the actual lake boundaries, is used for the conformal mapping. In these segments the actual boundaries can only be approximated by a step function, and the land areas must be masked with a special technique (Wilkin et al., 1995). For the computation in a rectangular basin, such a conformal mapping in the horizontal direction is not needed; however, for the Lake Constance, it is required, but not unique. In this case, a conformal mapping is performed for a bordering line which deviates in some segments from the actual shore line, so that a more uniform grid system is obtained. Then, in the numerical computations, the land areas within the grid system must be excluded by the aforementioned masking technique. This numerical code has proved its suitability in several lake applications, including diffusion problems , substructuring procedures  and wave dynamics in baroclinic wind-driven circulation dynamics . Here, the application of this code to the barotropic motions is performed. It will be seen in the following numerical results that almost all typical features in a homogeneous lake, e.g. inertial waves, Ekman spirals, can be reproduced by this model.

Barotropic motions in a rectangular basin
A rectangular basin of 65×17 km 2 extent with a 100 m depth corresponds to the space scale of Lake Constance; we assume homogeneous water, initially at rest, and subject to external wind forcing. This wind blows in the long direction of the rectangle, uniformly in space, Heaviside in time, and with strength 4.5 m s −1 (10 m above the water surface), corresponding to a wind stress of approximately 0.0447 Pa at the water surface. This is the wind-forcing we apply throughout the paper, except for Sect. 5. Integration starts at rest until a steady state is reached. We shall implement the discretization with 10-30 Chebyshev polynomials; the exact number depends on the magnitude of the vertical eddy viscosities (the smaller the vertical eddy viscosities, the more Chebyshev polynomials are needed to achieve stable and convergent numerical computations). We choose x = y = 1 km and let the numerical value of the horizontal eddy viscosity of momentum be constant ν h = 1 m 2 s −1 , because it turned out that the numerical values of ν h are not very crucial (Wang, 1996), and this value is generally reasonable (Csanady, 1978;Hutter, 1984). The insensitivity of the horizontal eddy viscosity can be expected and easily recognized because the horizontal variation of the velocity field is relatively small, in comparison with its vertical variation, except in the vicinity of shores and when tracer diffusion is considered.
3.1 Surface steady current in relation to distributions and amplitudes of the vertical eddy viscosity Three chosen vertical eddy viscosities are prescribed as follows: -case (i): ν v ∈ [0.005, 1] m 2 s −1 , -case (ii): FIG. 1: Rotation angle Φ between wind and surface current (a) and absolute fluid velocity V surf at the free surface (b). Here Φ and V surf are plotted against ν v (solid curves), A (broken curves), A 0 (dotted curves), respectively, corresponding to the cases (i), (ii) and (iii)   -case (iii): where A 0 ∈ [0.01, 1] m 2 s −1 .
In case (i), the vertical eddy viscosities are spatially constant, but their values are varied. Case (ii) corresponds to a z-dependence of ν v , where ν v is large in the upper layer, growing first with depth until it reaches a maximum in 30 m and then rapidly decreases to the value 0.01 m 2 s −1 at 45 m depth and below. This form corresponds with the vertical eddy viscosity distribution obtained by Svensson employing an one-dimensional k − ε model (Svensson, 1979) and by Güting, using a three-dimensional k − ε closure condition (Güting, 1998). Case (iii) assumes a linear variation of ν v from the free surface to 10 m depth below, in which ν v is kept constant at the value ν v = 0.03 m 2 s −1 (Csanady, 1980;Madsen, 1977). The maximum values of ν v can be varied by varying the free parameters A and A 0 . Choices like these are considered realistic and have, to some extent, been analysed with a linear equation system, in which the nonlinear advection terms are not included in finite depth oceans (see Heaps, 1984).
One typical result of the Ekman problem is the angle between the direction of the wind and the surface current in steady state. In the original problem, an infinite ocean with constant vertical eddy viscosity exists, and the surface current v surf is 45 • to the right of the wind. Results of our nonlinear calculations are summarized in Fig. 1. In case (i), the maximum value of is 42.5 • for ν v = 0.005 m 2 s −1 , and decreases rapidly with increasing vertical eddy viscosity to a value as low as 4 • for ν v = 1 m 2 s −1 . It can never reach 45 • as it does (independent of the value of ν v ) in the infinite ocean, because of the finiteness of the rectangular basin that induces a geostrophic flow which has the tendency to reduce . This simply indicates how significantly the circulating     motion in a finite basin can depend on the absolute values of the eddy viscosities. Even more surprising are the results of case (ii). Here varies non-monotonously with A: For A ≤ 0.12 it decreases with growing A, while it grows with A when A ≥ 0.12. Case (iii) presumes, for most values of A 0 , a very active turbulent near-surface layer and illustrates that for such a stiff situation, the typical dependence of can even be reversed. In Fig. 1b the moduli of the surface velocities v surf are shown as functions of the eddy viscosities. As expected, they rapidly decrease with increasing eddy viscosities. These results demonstrate how critical the distribution and absolute values of the eddy viscosities are; reducing the numerical diffusion to levels below the physical values is, therefore, compelling, if results are not to be physically falsified. Figure 2 shows the time series of the vertical velocity component w in 30 m depth at the four nearshore midpoints, as sketched in panel (b) of Fig. 2, subject to three different magnitudes of the constant vertical eddy viscosity, respectively. Immediately after the west wind sets in, the motion starts rapidly with an upwelling at the western end (Fig. 2a) and a downwelling at the eastern end (Fig. 2b), due to the direct effect of wind stress at the surface. In the northern midboundary point, there occurs an upwelling (Fig. 2c), while in the southern shore, a downwelling emerges (Fig. 2d). This behaviour is, clearly, due to the Coriolis force that causes a horizontal velocity drift to the right and therefore, is responsible for the downwelling (upwelling) at the southern (northern) shore on the northern hemisphere. The motion is characterized by oscillation, with a period of approximately 16.3 hours, which can obviously be identified as the inertial period using a Coriolis parameter approach for Lake Constance, f = 1.07 × 10 −4 s −1 . The smaller the vertical eddy viscosity is, the stronger the superimposing oscillations will be, and the longer the oscillations will persist. For the larger value ν v = 0.02 m 2 s −1 , only one day after the wind starts, the inertial oscillations are damped out, whereas the oscillations for ν v = 0.001 m 2 s −1 can still be clearly identified after six days. For this small eddy viscosity, due to the emerging strong oscillation, even a short-time downwelling (upwelling) at the western (eastern) shore can be observed.

Time series of the velocity components
In Fig. 3, the time series of the horizontal velocity components u (a) and v (b) for a medium value of ν v = 0.005 m 2 s −1 in the center of the basin, at different depths, are displayed. Steady motion is reached approximately after five days, but initially, an oscillating motion can be seen at all water depths, however, with decreasing amplitude as the depth increases. It is also interesting that the horizontal velocity that oscillates above 30-40 m, with the opposite phase from below this depth, corresponds approximately to the Ekman depth. This can be explained by the use of a linear equation system. It is known that for a linearized equation system, the horizontal motion consists of a drift current directly due to windstress, and a geostrophic current due to the surface slope (the surface pressure gradient) that can be described as follows: Neglecting the horizontal friction terms, the linearized horizontal momentum equations (2) and (3) then take the form with ν v = const. For the homogeneous case the dynamic pressure is independent of z. The local solution u(z, t), v(z, t) of (5) (depending only parametrically on x, y) may be written as   where u 1 (t), v 1 (t) is the geostrophic current which is a solution of due to the surface pressure gradient, whilst u 2 (z, t), v 2 (z, t) is the depth-dependent or frictionally induced velocity components With the suddenly imposed windstress τ 0 at the surface in the positive x-direction, and with a surface pressure gradient (a surface slope) ∂φ/∂x in the positive x direction and bottom stresses neglected, (7) and (8) have the local solutions of the geostrophic current and of the drift motion, respectively. The solution (10) was first given by Fredholm (Ekman, 1905). It is obvious that the two motions have opposite phase. It is of interest to note that such a combination of the currents can give a satisfactory approximation for the vertical current profile in regions far away from the lateral boundaries, even in those cases where nonlinear effects cannot be neglected in solving the horizontal circulation problem (Nihoul and Ronday, 1976). Above the Ekman depth, the drift motion (u 2 , v 2 ), induced directly by windstress, is dominant, whilst below this depth, the Ekman drift motion dies out and only the geostrophic current (u 1 , v 1 ) remains, which is independent of depth under homogeneous, hydrostatic conditions. From the solutions of such a simplified linear, horizontal momentum equation system, one can easily see that these two contributions to the current oscillate with opposite phase; at 30 to 40 m depth, the two oscillations possess comparably large amplitudes, and hence, the superimposing oscillation is hardly visible at this depth (see the 40m-depth curves in Fig. 3). Such oscillations have been discussed by Krauss (1979) and Csanady (1984). The same time series of the velocity components u and v, as shown in Fig. 3, are displayed again in the form of hodographs in Fig. 4 for 0, 5, 10, 20 and 60 m depths. The small circles along the curves mark time intervals of 1/4 of the inertial oscillations, which is approximately 4 hours. The arrows represent the velocity vectors in the steady state. This graph is very similar to the result obtained by Krauss (1979) in an infinite channel of rectangular cross-section. In the upper layer (0, 5, 10, 20 m), the wind-induced drift current, due to sudden wind, approximates the steady state in the form of inertial oscillations. The spiral in 60 m depth reflects the adaptation of the current field to the geostrophic current. The reflection of the water surface, due to the sudden wind, produces a sudden change of the pressure field, in the entire water column. The adaptation of the current field, which was zero, to the new condition, occurs in the form of inertial waves. Similar to those in the surface layer, these waves start together with the wind because the pressure gradient is felt immediately. In the hodograph for 20 m depth, one can see, at first, (perhaps during the first two hours) the same northwestward current as in 60 m depth; at this initial time, only the geostrophic current exists in 20 m, but when the surface  drift current reaches this depth, an abrupt change in direction occurs. It can be explicitly detected that the wind-induced current in the upper layer (in 0, 5, 10 and 20 m depth) oscillates with opposite phase of the geostrophic current in the lower layer.

Ekman spirals
In Fig. 5    the vertical velocity is negligible, for example, in the middle of a lake, the motion under homogeneous conditions can be considered a linear combination of (i) the geostrophic current, constant in the vertical, (ii) a bottom current deviation indicating the bottom Ekman layer, decaying exponentially away from the bottom which, together with the geostrophic current, satisfies the bottom boundary condition, and (iii) a surface-drift current (the Ekman layer) decreasing rapidly with depth and satisfying the surface dynamic boundary condition. In an infinite ocean with infinite depth, there exists only the surface-drift current. As is evident from Fig. 5a, for a constant vertical eddy viscosity, the drift current at the surface, in steady state, is directed 45 • to the right of the windstress (on the northern hemisphere). In a bounded rectangular basin, due to the surface slope under windstress, a geostrophic current shares the wind-induced motion. Subject to a western wind, a surface pressure gradient towards the east is built up, and hence, the caused geostrophic current points towards the north, as displayed in Fig. 5b; the current below the Ekman depth is approximately 60 m for ν v = 0.02 m 2 s −1 . If the northward geostrophic current is eliminated from the motion in Fig. 5b, a velocity field is obtained which is almost identical with that of the infinitely large ocean (Fig. 5a). If bottom friction cannot be ignored, the bottom Ekman layer occurs (Fig. 5c). The surface drift current, as well as the bottom current deviation, decreases rapidly with distance from the boundaries and a decay rate depending on the magnitude of the vertical eddy viscosity, as can be seen from a comparison of Fig. 5c,e. We also note that the northward geostrophic current automatically reduces the off-set angle of the surface current relative to that of the Ekman drift, a fact visible when comparing Figs. 5a and 5b. Finally, Figs. 5d,f display the Ekman spirals for a position near the midpoint of the southern shore (the distance from shore is 500 m) that is computed for two values of the verti-cal eddy viscosities. Here too, the value of the vertical eddy viscosity considerably affects the distribution of the velocity.

Total kinetic energy in relation to the vertical eddy viscosity and the wind-fetch
In Fig. 6, the time series of the total kinetic energy in the homogeneous rectangular basin of constant depth, subject to constant western (a) and southern (b) wind, respectively, are displayed, as previously computed with three vertical eddy viscosities, as indicated. As has been seen in the time series for the horizontal velocity, the inertial motions of the transient energy persist much longer, when ν v = 0.005 m 2 s −1 than when ν v = 0.02 m 2 s −1 , and this difference is even more distinct between ν v = 0.005 m 2 s −1 and ν v = 0.001 m 2 s −1 . The total kinetic energy stored in the basin for ν v = 0.001 m 2 s −1 is much larger than for ν v = 0.005 m 2 s −1 or ν v = 0.02 m 2 s −1 , due to the much larger energy input from windstress because of the much larger water velocity at the free surface (or more precisely, its component in the wind direction), and the likely smaller dissipation for the smaller vertical eddy viscosity. Comparison of Figs. 6a,b shows that the kinetic energy subject to the longitudinal, western wind is only slightly larger than that which is subject to the transverse, southern wind. The dependency of the total kinetic energy on the wind direction is worked out in terms of a polar diagram, which plots the variation of the total kinetic energy stored in the stationary circulation of the rectangular basin as a function of the direction of the wind. We compute the water motion subject to a constant wind from different directions (in intervals of 10 • ) for a fixed vertical eddy viscosity ν v = 0.02 m 2 s −1 . The polar diagram of the total kinetic energy is displayed in Fig. 7a. It is seen that the kinetic energy in the rectangular basin with constant depth depends only marginally on the wind direc-    Fig. 7. The dashed circles show the maximum value of the total kinetic energy, which is 9.71 × 10 8 N m for the case with consideration of the Coriolis force (a), but 8.57 × 10 9 N m for the case without the Coriolis effect (b). Fig. 8. Same as Fig. 7, but for a rectangular basin with larger aspect ratio. Here, the basin dimension is 60 km × 5 km × 100 m instead of the dimension of 65 km × 17 km × 100 m in Fig. 7. The dashed circles show the maximum value of the total kinetic energy, which is 9.71 × 10 8 N m for the case with consideration of the Coriolis force (a), but 8.57 × 10 9 N m for the case without the Coriolis effect (b).
tion. The kinetic energy attains its maximum for longitudinal, western wind, while for a transverse southern wind, the kinetic energy is minimum, although the maximum and the minimum differ from each other only by approximately 15%. As expected, the dependence of the kinetic energy on the wind direction is anti-symmetric about the E-W or S-N axis. It is also seen that the kinetic energy subject to southwest or northeast wind is larger than that for northwest or southeast wind. This loss of symmetry is due to the effect of the Coriolis force, which can be demonstrated by means of repeating the computations, but under the neglect of the Coriolis force. The analogous polar diagram of the total kinetic energy, without consideration of the Coriolis force, is displayed in Fig. 7b. In this case, the dependence of the kinetic energy on the wind direction is symmetric about the E-W or S-N axes. When the Coriolis force is ignored, the kinetic  energy is almost an order of magnitude larger. One likely reason may be that the water motion decays exponentially with increasing depth when the Coriolis force is accounted for, while the velocity decreases only linearly with depth, if the Coriolis force is neglected. The other reason may be that in the absence of the Coriolis force, the water motion at the free surface is mainly along the wind direction, hence, there is more energy input from the windstress. When the Coriolis force is neglected, the maximum kinetic energy occurs for a wind blowing in the diagonal direction of the basin. In fact, the difference in the total kinetic energy depends primarily on the wind-fetch, i.e. the integrated distance over which windaction takes place, taken along the wind direction. For diagonal winds, the wind-fetch is at maximum. In the presence of rotation, the integral for the fetch would probably have to be taken over the Ekman layer. This may be why the maximum of the kinetic energy is not reached for diagonal winds.
The stronger dependence of the total kinetic energy on the wind direction can be expected for a basin with a larger aspect ratio (e.g. a larger variation of the wind-fetch). In order to demonstrate this point, we repeat the computations for a basin with a larger aspect ratio. A basin dimension of 60 km × 5 km × 100 m is used instead of 65 km × 17 km × 100 m. The results are displayed in Fig. 8. In this case, the total kinetic energy depends on the wind direction much more intensively than in Fig. 7. Obviously, the basin must be very long and narrow before a directional dependence of the total kinetic energy is appreciable.
In the next section we will see that the response of Lake Constance, subject to a wind, is much more conspicuously dependent on the wind direction with regard to the total kinetic energy, even though its geometry is close to the rectangular basin in Fig. 7. This means that such a dependence on the wind direction depends not only on the wind-fetch, but also on the bathymetry of the lake basin.

Barotropic circulation in Lake Constance
Lake Constance consists of three basins: Obersee,Überlinger See and Untersee, but the Untersee is separated from the other two basins by a 5 km long channel; we shall be concerned here with the ensemble Obersee+Überlinger See, for brevity, also referred to as Lake Constance. It is approximately 64 km long and 16 km wide with has a maximum depth of 253 m and an approximate mean depth of 100 m, as shown in Fig. 9.
A study analogous to that above was also performed for Lake Constance with the same number of grid points (65 × 17) as for the rectangle. Computations were also done for ν h = 1.0 m 2 s −1 , ν v = 0.02 m 2 s −1 and ν v = 0.005 m 2 s −1 , respectively. In the second case, the larger number of polynomials, namely N = 25 instead of N = 10, was needed to achieve stable and convergent numerical integration. Let a spatially uniform, temporally constant wind with a wind speed of 4.5 m s −1 blow from Northwest (305 • (NW) in the longitudinal direction of the lake) until steady state conditions are established.

Horizontal distribution of the steady currents
The horizontal variation of the currents at the surface and at depth levels 10, 20 and 40 m is depicted in Fig. 10. The horizontal motion is displayed by arrow-diagrams. In the central part of the lake, the surface current is deflected by about 50-60 • to the right of the downwind direction and amounts to 5 cm s −1 , on average. The boundary currents somewhat off the northern and southern shore increase to a magnitude of 8-10 cm s −1 and are directed parallel to the boundary.
As is typical for drift currents, the direction of the motion in the open lake continuously turns to the right on with increasing depth and decreasing velocities. At 10 m depth and far away from the shores, the current deflection to the right of the wind exceeds 90 • , and the average velocity decreases to 3-4 cm s −1 . At the 40 m depth, the motion in the middle of the lake is nearly reflected to the upwind direction. However, above the 40 m depth, the shore currents are still basically    parallel to the shores with the exception of the northeastern shore, where the motion below the 20 m depth is almost in the upwind direction. Near the northeastern corner, a cyclonic gyre is visible.
In order to show the general influence of the depth configuration on the wind-driven circulation, the vertically integrated volume transport is displayed in Fig. 11   ume transport streamlines. The circulation consists, in principle, of two cells which rotate in such a way that the transport is downwind along the shores, orientated roughly in the direction of the wind. The central part of the circulation shows a net transport with an upwind component which results from the geographic current (caused and bounded by the bottom topography); with the bottom current prevailing over the surface drift current. In contrast, the drift current dominates the other currents in the shallow, nearshore region. Due to the Coriolis force, subject to a western wind forcing, the sloping bottom causes an intensification of the northward flow component at the northwestern shore as well as at the southeastern shore (Serruya et al., 1984;Wang, 1996). If there would be no Coriolis effect, the circulation would consist of two gyres, which are located almost exactly to the north and south of the Talweg, rotating in the clockwise and anticlockwise directions, respectively.

Ekman spirals and time series of the horizontal velocity
The vertical variation of the current depends strongly on the eddy viscosity. We display in Fig. 12 two steady Ekman spirals at the midlake positions of theÜberlinger See (above) and the Obersee (below), as they form, for an impulsively applied uniform wind from 305 • (NW) (approximately in the long direction), and as obtained with the two different indicated eddy viscosities. Those Ekman spirals are considerably affected by the ν v -values. The turning of the arrows which make up the spirals, also indicates that the surface Ekman boundary layer is thinner for the smaller value of the eddy viscosities (right panel) than for the larger value (left panel).
Equally interesting is the comparison of the time series of the horizontal velocity components u and v for various depths at the midlake positions of theÜberlinger See and the Obersee, as displayed by the hodographs in Fig. 13 and obtained with ν v = 0.005 m 2 s −1 . At both positions, transient oscillations can be discerned with the inertial period of approximately 16 h. The oscillations can be seen at all water depths, however, with decreasing amplitude as the depth increases. Furthermore, they die out before two inertial periods in theÜberlinger See, but only after 5-6 inertial periods in the Obersee. The reason is the smaller size of thë Uberlinger See and, therefore, the enhanced frictional resistance due to the lake bottom and the side shores. In the Obersee, the transient velocities oscillate around a slowly  FIG. 13: Hodographs, i.e., time series of the horizontal velocities in the middles of theÜberlinger See (a) and the Obersee (b) in 0, 5, 10 and 60 m depths. The motion is set up from rest. The small circles mark time intervals of approximately four hours. In order not to confuse the spirals at 20 m and 60 m depth the 20 m hodographs are shown separately in the panels (c) and (d) for the two positions.
Total kinetic energy [×10 9 Nm]   FIG. 13: Hodographs, i.e., time series of the horizontal velocities in the middles of theÜberlinger See (a) and the Obersee (b) in 0, 5, 10 and 60 m depths. The motion is set up from rest. The small circles mark time intervals of approximately four hours. In order not to confuse the spirals at 20 m and 60 m depth the 20 m hodographs are shown separately in the panels (c) and (d) for the two positions.
Total kinetic energy [×10 9 Nm]    varying mid-velocity at each depth; this is different from the corresponding behaviour in the rectangular basin where the motion oscillates approximately around the steady state (compare with Fig. 4). This difference is due to the complex topography of Lake Constance. From the spirals at larger depth (e.g. 20 m), during an initial time interval after the wind started, a sudden change in the direction of the motion can be seen when the wind-induced surface drift flow reaches this depth. Before this time, only the geostrophic current exists at this depth which, due to the restriction of the bottom topography, is along the Talweg. This is different from the behaviour of the flow in a basin with constant depth, in which the geostrophic motion is basically perpendicular to the wind direction or surface pressure gradient (90 • to the left on the northern hemisphere).

Total kinetic energy in relation to the vertical eddy viscosity and the wind direction
The time series of the stored, total kinetic energy in Lake Constance for NW and SW wind, ν v = 0.02 m 2 s −1 and ν v = 0.005 m 2 s −1 (Fig. 14) show that the inertial oscillations persist longer, and the value of the total kinetic energy is much larger for the smaller ν v -value (as shown before for the rectangular basin). More interesting is the comparison of the kinetic energies for the case of longitudinal NW wind and transverse SW wind (Fig. 14a,b). For the transverse wind, the temporal inertial oscillations persist much longer than those associated with the longitudinal wind, whereas the magnitude of the kinetic energy subject to the longitudinal wind is much larger. The variable response of the lake to constant wind forcing from different directions at 10 • interval is represented in the polar diagram of the kinetic energy uptake in Fig. 15a, computed for ν v = 0.02 m 2 s −1 . The stored, total energy depends very strongly on the wind direction. The geographical direction of maximum exposure coincides approximately with the longitudinal direction of the lake. The kinetic energy under a transverse wind is only one fifth of thekinetic energy for a longitudinal wind. The strong dependence of the total kinetic energy on the direction of the wind should be mainly due to the topography, and less to the dimension of the lake, since Lake Constance has a comparable aspect ratio to the rectangular basin in Fig. 7, in which the dependence on the direction of the wind is much weaker. This dependence of the total kinetic energy on the direction of the wind is similar, to some extent, to that obtained by Serruya et al. (1984), where only the kinetic energy of the two-dimensional vertically integrated net transport was calculated. The dependence of the kinetic energy on the wind direction, without consideration of the Earth's rotation, is also displayed in Fig. 15, in panel b. In this case, the total kinetic energy is much larger than the kinetic energy with the Coriolis force, while the dependence of the kinetic energy on the direction of the wind is somewhat weaker; here, the kinetic energy under a transverse wind is nearly half as large as the kinetic energy for a longitudinal wind.

Vertical distributions of the steady currents
The isotachs of the three velocity components in a steady state on the cross-section through the center of theÜberlinger See are presented in Fig. 16. The wind is blowing in the direction of the longitudinal axis of the lake, perpendicular to the plane of the graphs upwards. The most distinguished feature of the motion is a nearshore coastal jet in the direction of wind (Fig. 16a). The change of the horizontal velocity with depth happens primarly in the upper layer and in the middle of the cross-section (Fig. 16a,b). At lower depths, where the geostrophic current plays an important part in the motion, the velocity components do not significantly change with depth. Due to the effect of the Coriolis force, for a west wind, there occurs a downwelling along the southern shore (left side in Fig. 16c) and an upwelling along the northern shore (right in Fig. 16c), however, the vertical velocity component is much smaller than the horizontal components. In the zone far away from shore, the vertical velocity component is practically negligible. Figure 18 displays the transverse variations of the longshore transports. These are the depth integrals of the positive and negative velocities in the x-direction, as functions of the y-coordinates in the cross-section, with the wind direction (positive, broken lines) and against the wind (negative, dotted lines), in cross-sections through the centers of thë Uberlinger See and the Obersee, respectively. From Fig. 18, the nearshore coastal jet in the direction of the wind can be seen more clearly, not only in theÜberlinger See, but also in the Obersee. Close to the shore, the transport is almost only in the direction of the wind; far away from shore, where the motion is dominantly in the upwind direction, the transport is against it. This feature is due to the approximate parabolic shape of the cross direction.
If the positive and the negative transports are added, the net volume fluxes per unit width (the solid lines in Fig. 18) are obtained. These lines would be obtained with a depth integrated model. Obviously, the horizontal integration of the net volume flux along the cross-section must vanish in steady state.
The vertical distributions of the horizontally integrated transport at the centers of theÜberlinger See and the Obersee can be extracted from Fig. 19. It can be seen that the transport occurs primarly in the top 80 m, especially in the top 40 m. The transport in the wind direction (positive, broken lines) assumes its maximum at the surface, while the maximum of the transport against the wind (negative, dotted lines) occurs at approximately 20-30 m depth. Below 100 m there exists only the transport against the wind with very small values. The sum of the two fluxes for a fixed depth yields the net volume flux per unit depth, which is also displayed as solid lines in Fig. 19. Obviously, its vertically integrated total flux through a cross-section must vanish in a steady state. It is also obvious that these results cannot be obtained with the use of a vertically integrated model.
In Figs. 17, 20 and 21, graphs are displayed for the isotachs of the three velocity components in a cross-section of theÜberlinger See (Fig. 17) as well as the horizontal and vertical distributions of the longshore transport in the two indicated cross-sections of theÜberlinger See (Fig. 20) and the Obersee (Fig. 21), respectively. This is also the case in Figs. 16, 18 and 20, but now with a smaller vertical eddy viscosity ν v = 0.005 m 2 s −1 . Principally, they are very similar to those subject to ν v = 0.02 m 2 s −1 . The most distinguishing differences is that for the smaller vertical eddy viscosity, the surface layers of positive x-velocity are much smaller (Figs. 17 and 21) than in Figs. 16 and 19, while the absolute values of the velocity are much larger, and the two-dimensional depth integrated model is less justified (compare Figs. 18 and 20).

Comparison of measured values and computational results in Lake Constance
Unfortunately, there are insufficient data available which would allow validation of the model. Nevertheless, a partial check of the reliability of the numerical method was possible with data from 16 February 1993 (    weak. On 19 February, a strong wind was initiated; its amplitude reached 8 ms −1 and it persisted for the remainder of the period. For the same period, the water velocity and its direction at the "Mainauschwelle" (near the island Mainau) at 80 m depth were also measured; they are displayed as solid curves in Fig. 23a,b. In the first 80 hours, the water velocity was less than 1.1 cm s −1 , which was below the threshold of the current meter. At the fourth day, strong currents started with superimposed oscillations of a period of approximately 16 hours, which can be interpreted as inertial oscillations. The flow direction superimposed by oscillations is approximately 300-340 • (NW).
Here we simulate, numerically, the corresponding water motion with the measured wind input and check if the measured water velocity can be reproduced by the computed velocity field. In winter (here, in February), the water density in Lake Constance can be considered to be homogenous (Bäuerle et al., 1998). For this simulation the measured wind at the station "Boje Mitte" is used and extrapolated to the entire lake. This uniformity is certainly unlikely to be realistic, but we apply it due to lack of better knowledge. The shear stress at the water surface can be calculated with the aid of the classical drag formulas  where V wind is the wind speed 10 m above the water surface, ρ a the air density (ρ a = 1.225 kg m −3 ) and c 0 a dimensionless friction coefficient. The specification of c 0 is not unique and varies from author to author. A typical value for a weak or medium wind strength (V wind < 10 m s −1 ) is c 0 = 1.8×10 −3 (Lehmann, 1992). The wind speed measured at 4.4 m above the lake surface (Fig. 22) must be converted to the value at 10 m. Assuming a logarithmic wind profile, yields where z 0 is a measure for the corrugation of the water sur-face (roughness length). In the computations, we choose z 0 = 1.0 × 10 −4 m. The vertical eddy viscosity is assumed constant with value ν v = 0.02 m 2 s −1 and the simulation is started from rest, although at the initial time, a small motion must exist in nature.
The computed water velocity and its direction in 80 meters depth, at the position "Mainauschwelle" where the measured time series are plotted in Fig. 23, are also displayed in Fig. 23 as dashed lines. During the first three days, when a weak wind prevailed, the measured and computed current speeds could not be compared, but their directions could (Fig. 23b), and, they deviated considerably from one another. This is most likely due to the fact that the computations started from a state of rest, but there was some (small) motion in nature (however, not caught by the current meter) which was not considered in the simulation. A relatively strong wind is needed to bring the computed and the measured velocities together. This strong wind occurs after 80 hours when both measured and computed current speed and direction coincide quite well with one another. Only at the seventh day do the current speeds of the computations differ from the respective measured values. The computed current orientation at this depth lies between 300-340 • (towards Northwest), which is in fair agreement with the measured values (compare the solid and dashed lines of Fig. 23b), which are basically around 320 • , superimposed by oscillations. Given the measuring technique and the bold extrapolation of the wind over the entire basin, a better agreement can hardly be expected.
In Fig. 24a,b we have plotted time series of the two horizontal velocity components at several depths in the midlake position in the Obersee, as indicated in the inset. Evidently, the strong wind commencing after 80 hours generates equally strong oscillations; with a period of almost exactly 16.3 hours, they are most likely, inertial oscillations. On the other hand, comparison of the water velocity at 80 m depth at the Mainauschwelle (Fig. 23) with the driving wind suggests that these oscillations might simply be the direct response to the wind forcing. In order to confirm that the oscillations exhibited by the time series of water velocity are indeed inertial oscillations due to the Coriolis force, the same  computation but without consideration of the Coriolis force was repeated. For this case, the computed water velocities at the Mainauschwelle, at 80 m depth, and at the center of the Obersee, at several depths, are shown in Figs. 25 and 26. Comparison of Figs. 23 and 25 shows that without the Coriolis term, the computed velocity at the Mainauschwelle does not exhibit the oscillations shown by the measurements. Moreover, the direction of the current is toward 260 • , which is substantially different when one accounts for the Coriolis effects. Similar differences are also seen when Figs. 24 and 26 are compared. The computed current oscillations at the midpoint position of the Obersee are not established without the Coriolis force. These facts should be sufficient demonstration that the measured oscillations are indeed rotational effects due to the motion of the Earth.

Concluding remarks
In this paper, wind induced barotropic circulation in lakes was studied, first, from a more fundamental point of view, using a rectangular basin with constant depth, but later, with Lake Constance, a medium-size Alpine lake.      viscosities were varied in three different ways, as suggested by other studies, and the steady Ekman problem was solved. The following results were obtained, each demonstrating the significance of the eddy viscosity: -The direction of the surface water velocity at the midpoint of a rectangular basin relative to that of the wind, the so-called wind set-off, depends strongly on the absolute value and vertical distribution of the vertical eddy viscosity. It cannot be concluded that the current set-off (to the right on the northern hemisphere) decreases with increasing eddy viscosity. Depending upon the vertical distribution of the eddy viscosity, opposite or even nonmonotonic behaviour may occur (Fig. 1a).
-The surface current speed, however, decreases monotonically with an increase in the eddy viscosity (Fig. 1b).
-The amount of kinetic energy stored in the water depends equally upon the absolute value of the eddy vis-cosity (but less on its vertical distribution). Its transient behaviour from a state of rest to a steady state is characterized by oscillations that are attenuated in time, but equally, the smaller the vertical eddy viscosity is, the larger the oscillations. This is true for the rectangle (Fig. 6), as well as Lake Constance (Fig. 14).
-The wind-directional dependence of the total kinetic energy depends strongly on the wind-fetch and the lake bathymetry. For a rectangle with constant depth and a width to length ratio of 0.25, this dependence is less than approximately 10% (Fig. 7); for Lake Constance, it is substantial (Fig. 15).
-The Coriolis force implies a significant reduction of the total kinetic energy that can be established in a basin due to wind forces, when compared to the case when rotational effects are neglected (Figs. 7, 8 and 15).   -The vertical distribution of the water velocity, i.e., the Ekman spiral, depends considerably upon the absolute values of the eddy viscosities. Not only the Ekman depth is affected by these values, but equally so, the orientational distribution of the horizontal current with depth (Figs. 5, 12).
-Inertial oscillations are more easily established in open water than bounded channels; this demonstrates the significance of the frictional effects due to boundaries (Fig. 13).
All these results support the conjecture that it is important for adequate reproduction or prediction of observed current fields in lakes to use the correct orders of magnitudes of the numerical eddy viscosities. On the other hand, the eddy viscosity should be chosen according to physical considerations, and simultaneously, ensure numerical stability of a simulation. On the other hand, if the eddy viscosity desired by numerical stability is much larger than the values suggested by measurements, the numerical code needs to be adjusted, especially with some modern numerical treatments of the nonlinear advection terms. Inappropriate numerical schemes of these advection terms are often the reason for the production of numerical oscillations, which require large eddy viscosities in order to prevent their development and hence, assure numerical stability. At present, we are making an effort in this aspect. It is exactly the problem in many existing threedimensional lake circulation models, that in order to reach numerical stability, the eddy viscosities must be chosen to be much larger than physically acceptable so that the physical oscillations are also rapidly damped away or are even indiscernible. It has been clearly demonstrated by this model that observed inertial oscillations due to the rotation of the Earth persist long and are slowly attenuated. That these inertial waves are indeed a dominant effect is shown in a restricted comparison of the measured current at 80 m depth for Lake Constance. Neglecting the Coriolis effects does not yield re-sults that can be compared with the measured ones. This wave dynamic, in turn, determines the advective properties and thus, the transports of nutrients and pollutants both horizontally and vertically.