The numerical modelling of circulations in shallow lakes is a relevant tool for all environmental applications in which flow advection processes are of interest, e.g. for studies on nutrients, microorganisms, pollutants and sediment dynamics. While three-dimensional (3D) models are needed to properly describe the flow fields of basins with the main circulations in the vertical plane, two-dimensional (2D) models are commonly deemed to yield adequate results for lakes with prevailing horizontal circulations. However, the depth-averaged approximation is more limiting for wind-driven flows than for gravity-driven ones, such as rivers, as the driving force is a surface rather than a volume one, distributed along the depth through turbulence. In this work, the effects of such inaccuracy on the reproduction of circulation layouts are evaluated through compared simulations between a 2D Shallow Water solver and a 3D Reynolds-Averaged Navier-Stokes one. The models are first applied to a simple enclosed elliptical test basin and then to the real case of the Superior Lake of Mantua, a shallow fluvial lake in Northern Italy, thereby also investigating the influences of the interaction of wind with a riverine current and of a complex bathymetry on the compared results.

INTRODUCTION

Shallow lakes are inland basins in which the wind energy advects the whole water volume, turbulence being able to completely mix the water column. Therefore, they do not display the stratification typical of deep lakes. In these water bodies, the wind stress reaches the bottom, both through the produced turbulence and through wind-generated waves, due to transitional (or even shallow) water wave conditions often being achieved. This triggers an intense resuspension of nutrients, sediments and pollutants during storm events.

Shallow basins with regular shorelines and bathymetry display main circulations in the vertical plane, with windward currents in the surface layer and return flows in the bottom region. Horizontal circulations are instead dominant in lakes with prominent circulation vorticity (curl) sources (Józsa 2006, 2014; Toffolon & Rizzi 2009), such as irregular shorelines, complex bathymetry, unhomogeneous wind stress distribution and extended vegetation patches. These elements concur locally altering the equilibrium between the pressure gradient and the wind force, which generates perfect circulations in the vertical plane in a parallelepiped basin with the wind aligned with the shores. Areas in which one of the forces prevails locally are thus determined, forming gyres with concordant currents throughout the flow depth.

Three-dimensional (3D) models are needed to suitably describe the flow field in lakes with the main circulations in the vertical plane, as two-dimensional (2D) models would result in strong opposing currents being integrated as weak depth-averaged circulations. In addition to impairing studies in which the focus is on the advection of organisms or substances which concentrate near the surface or the bottom, depth-averaging provokes an underestimation of the bottom stress ‒ which is essential for investigations on resuspension dynamics ‒ and an overestimation of the surface setup (Krámer 2006). Hearn & Hunter (1990), conversely, state that ordinary 2D models result in smaller setups than 3D ones.

Depth-averaged models would be naturally applicable to shallow lakes with prevailing circulations in the horizontal plane, due to their low flow depth, well-mixed conditions along the vertical and resemblance to riverine currents to which they are efficiently applied (e.g. Costabile & Macchione 2015), as it is still performed in many recent studies (e.g. Schimmelpfennig et al. 2012; Fabian & Budinski 2013; Józsa 2014). This similarity is even stronger for fluvial lakes, i.e. basins in which the through-flowing current from the inlet to the outlet sections interacts with the wind force in shaping the circulation layout. Nevertheless, while gravity acts evenly on the entire water volume so that its effect is correctly accounted for in a 2D model, wind stress is applied on the free surface and is passed along the water depth by turbulence in the vertical plane. The depth-averaged schematisation, in which wind is considered as an even volume force, dividing the stress by the depth, is then a substantial approximation of the real physics. Some authors (Davies 1988; Hearn & Hunter 1988) proposed alternative 2D formulations in which the velocity distribution along the vertical is analytically approximated with convolution methods from the surface stress and slope, allowing an evaluation of bottom stress from an estimate of the bottom velocity rather than from the depth-averaged one, resulting also in more accurate calculations of the latter. This, however, relies on the vertical eddy viscosity distribution along the depth to be fixed throughout the domain, thus being suited to continental shelf seas, which are characterised by very large spatial scales, regular topographies and simple circulation patterns in the vertical plane. Problems arise when these models are applied to cases with complex bathymetries and flow fields or to transients with variations in the wind forcing in which the eddy viscosity distribution is significantly space- and time-dependent.

Hutter (1984) introduced a ratio to estimate whether a flow can be reasonably approximated with depth-averaged values, which expresses the degree to which irregularities in the velocity profile along the vertical are damped by the vertical eddy viscosity: 
formula
1
in which Th = H2/4νT,vert is the characteristic timescale of momentum transport by turbulence (νT,vert being the vertical eddy viscosity along the characteristic depth H), which damps the irregularities in the velocity profile in the vertical dimension, while Tadv = L/U is the characteristic timescale of the flow transport (U being the characteristic horizontal velocity along the characteristic length scale L), which instead advects the irregularities. When Δ ≪ 1, irregularities are more rapidly damped than they are advected, resulting in an almost uniform velocity profile along the vertical, except near highly-sheared boundaries (Józsa 2006), such as the logarithmic profile of turbulent water flow. However, there are shallow lakes with low Δ values, such as Lake Balaton (Hungary) (Józsa 2006), which are indeed characterised by main circulations in the vertical plane (e.g. Ciraolo et al. 2004). In addition, flow velocities in shallow lakes can increase by up to an order of magnitude during storm events, so that the same basin may be characterised by different Δ regimes. Furthermore, the approximation of evenly distributed wind stress along the water column, in addition to being increasingly inaccurate for growing flow depths, is also prone to fail for lower wind speeds, hence for inferior flow intensities, whereas Equation (1) suggests that depth-averaging works better for lower water velocities. Therefore, the Hutter criterion should not be used alone to state the convenience of the depth-averaged approximation for wind-driven flows.

This work focuses on the effects of the wind surface stress as volume force approximation ‒ adopted in 2D depth-averaged models ‒ on the simulation of flows in shallow lakes with the main circulations in the horizontal plane, such error being less understood and obvious than for basins with the prevailing flow cells in the vertical dimension. The results of a 2D and a 3D model, which rationally holds a more accurate representation of reality, are compared. Simulations are first performed for an enclosed elliptical fictitious basin case, to assess the differences under ideal topographical conditions. The models are then applied to a real-world shallow fluvial lake case (Superior Lake of Mantua, Northern Italy), investigating the effects of the complex bathymetry and of the interaction of wind forcing with a through-flowing current on the error induced by the depth-averaged approximation. This numerical study aims at shedding light on a category of basins to which 2D models are almost straightforwardly applied, highlighting instead relevant differences with the results of 3D models under a herein identified range of conditions.

METHODS

Numerical models

2D model

Two-dimensional simulations were performed with the ORSA2D code (Petaccia et al. 2010), which integrates the Shallow Water Equations (SWE), written in conservative form. Wind stress was added to the momentum equations as a volume force term over the water depth. Effects of turbulence in the horizontal plane can be taken into account with a constant eddy viscosity formulation. The following system of continuity and momentum equations is solved with the Finite Volume (FV) method over unstructured grids, adopting the Roe (1981) scheme and either a first- or second-order upwind approach in space and time: 
formula
2
with: 
formula
3
in which t is the time, x the generic spatial coordinate, h the water depth, q = Q/x the unit discharge (Q being the total one), g the gravity acceleration, S0 the bed slope, Sf the friction slope, computed from Manning's formula, τs the wind surface stress, τT the turbulent stress, ρ the density and νT,hor the horizontal eddy viscosity. In the numerical scheme, the bed slope terms are upwinded (Bermúdez et al. 1998), the friction slope terms are treated in a semi-implicit way (Costabile et al. 2015) and the turbulent stresses are solved explicitly, while the local values of wind stress are given as steady constant external values. The second-order accuracy is obtained by means of a MUSCL procedure (Osher 1985), together with a Minmod limiter and a predictor-corrector approach (Hirsch 1988). The wind stress is ordinarily defined throughout this work as: 
formula
4
in which ρa = 1.225 kg/m3 is the air density and c10 the wind drag coefficient related to the wind intensity at 10 m height W10.

3D model

Three-dimensional simulations were carried out with the widely validated STAR-CCM+ v9.02 CFD (Computational Fluid Dynamics) package (CD-adapco 2014). Well mixed conditions, i.e. constant temperature and density, were adopted. The model integrates the conservative Reynolds-Averaged Navier-Stokes equations through a FV first-order or second-order upwind approach. The following continuity and momentum equations, written here with the application of the Boussinesq hypothesis for the Reynolds stress tensor term, are solved in sequence through a projection method, based on the SIMPLE algorithm: 
formula
5
with: 
formula
6
and: 
formula
7
in which is the mean velocity, the mean pressure, ν the kinematic viscosity, νT the global eddy viscosity and B the summation of the external body forces. In the present simulations, Equation (6) comprises: the vertical gravity force fg, where γ is the specific weight of water and ∀cell the volume of a grid cell; the horizontal wind force fs, which acts only over the surface cells, zg being the vertical mesh resolution; and the vegetation flow resistance force fp, which is applied only to the cells within the vegetated flow region in the Superior Lake of Mantua case, and is expressed with a porous Forchheimer-type closure, μ being the dynamic viscosity of water and K and β the intrinsic permeability and inertial parameter principal tensors, respectively.

The Realisable k-ɛ model (Shih et al. 1994) was used to parameterise νT. The free surface was represented as a solid horizontal boundary with a slip wall condition. In the case studies discussed below, the surface setups, obtained from the surface pressures ps calculated by the 3D code as Δη = ps/γ, were in fact much smaller than the vertical mesh resolution, so that moving boundary techniques such as Volumes of Fluids (VOF) did not prove necessary.

Case studies

This paper discusses steady-state conditions obtained by iterating until stable solutions were achieved from the given initial conditions.

Test elliptical lake

Preliminary tests were performed on the test elliptical basin introduced by Józsa et al. (2007), which has bottom elevation zb defined as function of the x and y coordinates: 
formula
8
in which lengths are in metres and shores are placed where zb = −1 m. The resulting bathymetry of the 1.4 km × 2.4 km basin is shown in Figure 1. Initial water-at-rest elevation was set to η = 0 m unless where specified, but η = ±0.7 m conditions were simulated as well (in analogy to Józsa et al. 2007), investigating the influence of the water depth h on the error of the depth-averaged model. A modified version of the basin with flat bottom and zb = −1.5 m was employed in a particular simulation, to determine the relevance of the variable bathymetry on results.
Figure 1

Bathymetry of the test elliptical lake introduced by Józsa et al. (2007).

Figure 1

Bathymetry of the test elliptical lake introduced by Józsa et al. (2007).

Most simulations employed a spatially heterogeneous wind stress distribution according to the Internal Boundary Layer (IBL) theory (e.g. Józsa 2014), which considers the variation of the wind turbulent profile downstream of an aerodynamic roughness change. In the present case, the transition is from land to water, causing both c10 and W10 terms in Equation (4), hence τs, to increase with fetch. A wind intensity over land W10,l = 10 m/s from NW (315 °) was considered, resulting in the wind stress field in Figure 2. The aerodynamic roughness of land for IBL wind stress calculations was set to z0,l = 0.15 m (tall grass, see Józsa et al. 2007). Details about the wind stress distribution model and its implementation are present in Fenocchi (2015). A uniform wind stress distribution for W10 = 11 m/s from NW, i.e. for approximately the average wind intensity resulting from the IBL model, was also tested, the wind drag coefficient c10 = 1.45 × 10−3 having been estimated with the formula by Wu (1980), yielding τs = 0.225 N/m2.
Figure 2

IBL-derived wind stress field for W10,l = 10 m/s from NW for the test elliptical lake.

Figure 2

IBL-derived wind stress field for W10,l = 10 m/s from NW for the test elliptical lake.

In 2D simulations, a refined grid mimicking as close as possible that of the 3D model in the planimetric dimension was first adopted, with 165,438 quadrilateral and triangular elements with vertices spaced ∼4 m one another. A coarser grid with 60,649 triangular elements, with ∼10 m vertices distance, was however used in later simulations. The bottom roughness was originally set to Manning's n = 0.025 s/m1/3 (Józsa et al. 2007). The effect of the algebraic turbulence model on results was also evaluated, enabling it in specific simulations and testing different νT,hor values.

The 3D model was run on a grid with 1,425,561 parallelepiped elements trimmed at the boundaries, the horizontal and vertical resolution being equal to xg = 4 m and zg = 0.2 m, respectively (i.e. with aspect ratio rg = xg/zg = 20) and the free surface boundary being the water-at-rest elevation. The roughness height on the bottom was set to ks = 0.075 m, matching the original roughness in the 2D model according to the relation n = ks1/6/26 (e.g. Christensen 1992).

Superior Lake of Mantua

Further comparative tests were performed on the geometry of the Superior Lake of Mantua in the Po River Plain, Northern Italy. This is the most upstream of the three fluvial Mantua Lakes, originated at the end of the 12th century by the regulation of the Mincio River around the city of Mantua (Figure 3(a)). The lake hosts meadows of lotus flower (Nelumbo nucifera), a macrophyte characterised by cylinder-like stems spanning the whole water depth and both floating and emergent leaves and flowers. The largest lotus flower island, implemented in the models, occupies ∼10% of the lake surface (divided by trimming operations into eight blocks with interposed canals, see Figure 3(b)) and largely influences circulations (Fenocchi 2015; Pinardi et al. 2015; Fenocchi & Sibilla 2016). Due to the lack of topographical data, only the actual lacustrine part of the lake is modelled, omitting the upstream channelised transitional section from the Valli del Mincio wetlands, therefore introducing a fictitious inlet section. This affects a ∼ 500 m long upstream reach, in which the flow is however mostly river-like, except for drought discharges, so that the impairment of results is reduced. The bathymetry of the modelled area, relative to the control elevation of the basin η = 17.50 m above sea level, is shown in Figure 3(b). In all tests, the outflow is equally split between the two controlled outlets, the Vasarone and the Vasarina sluice gates, as in ordinary regulations. The influence of surface waves on circulations, owing to the increase of bottom friction, was proved to be negligible (Fenocchi 2015; Fenocchi & Sibilla 2016) and was therefore not included.
Figure 3

(a) Location of the Superior Lake of Mantua in the Po River basin and in the Mantua Lakes system; (b) bathymetry of the modelled area of the lake (adapted from Pinardi et al. (2015); ortophotos courtesy of the Lombardy Region).

Figure 3

(a) Location of the Superior Lake of Mantua in the Po River basin and in the Mantua Lakes system; (b) bathymetry of the modelled area of the lake (adapted from Pinardi et al. (2015); ortophotos courtesy of the Lombardy Region).

Two wind intensity conditions were considered, based on the analysis of the local climate (Fenocchi 2015): 1) standard wind W10,l = 1.94 m/s; and 2) ordinary storm wind W10,l = 10.10 m/s. A wind direction of 80 ° was adopted, typical of the summer season: such direction opposes the through-flowing current (Figure 3(b)), producing the most significant results. The reference discharge was set to Q = 20 m3/s, mean annual value for the basin, but simulations with reduced flow rate Q = 5 m3/s, typical drought discharge for the Superior Lake of Mantua, and null through-flow (enclosed lake conditions) were also performed in order to analyse the effect of through-flow on the differences between the models, i.e. to highlight the peculiarity of fluvial lakes.

In this case, multiple embedded IBL profiles (Krámer 2006) were considered in the calculation of the wind stress distribution, taking into account the multiple aerodynamic transitions between land, water and lotus flower covers. Wind sheltering by nearshore tall vegetation at the windward and leeward shores was also taken into account, according to the Backward- and Forward-Facing Step analogies (Ottesen Hansen 1979). Null wind stress was set over the lotus flower island, consistently with the simulation of summer conditions, as the leaf cover almost prevents wind from reaching the water surface (e.g. Józsa 2006). To account for correct fetch distances for westerly winds, the wind stress model considers the whole lake, extending upstream of the hydrodynamically modelled area, i.e. to the perimeter outlined in Figure 4 (the two upstream closed lines are vegetation islands). To take into account wind field irregularities, multiple radials were considered for wind stress field calculations, as suggested by the Shore Protection Manual (CERC 1984). The aerodynamic roughness of land was still set to z0,l = 0.15 m, while that of the lotus flower cover to z0,veg = 0.08 m. The resulting wind stress distribution for ordinary storm conditions is shown in Figure 4.
Figure 4

IBL-derived wind stress field for easterly ordinary storm wind for the Superior Lake of Mantua.

Figure 4

IBL-derived wind stress field for easterly ordinary storm wind for the Superior Lake of Mantua.

Two-dimensional simulations were performed on a grid with 45,336 triangular elements, with vertices spaced ∼12 m apart. The bottom roughness was still set to n = 0.025 s/m1/3 for bare bed areas. The flow resistance of the main lotus flower island was taken into account by an increased roughness coefficient, calculated with the Petryk & Bosmajian (1975) formula for emergent rigid vegetation: 
formula
9
in which n is the bare bed roughness, Cd = 1.1 the drag coefficient (very sparse canopy, see Wang & Wang (2011)), a = 0.12 m−1 the mean surveyed vegetation density for the lotus flower in the Superior Lake of Mantua (Fenocchi 2015) and hm,veg = 2.25 m the average flow depth in the vegetated region. A uniform normal velocity condition was set at the inlet boundary, while stage-discharge relations between the control elevation of the lake and the required flow rates were set at the outlets. Such water elevation and the inlet velocity were also adopted as initial conditions throughout the domain. The algebraic turbulence model was disabled, having been proved to be unsuitable in the test elliptical lake case, as discussed in the Results and discussion section.

In the 3D case, simulations were performed on a grid with 847,843 trimmed parallelepiped elements having xg = 6 m horizontal resolution and zg = 0.3 m vertical resolution (rg = 20), the free surface boundary being the control elevation. The bottom roughness height was still set to ks = 0.075 m. For the lotus flower island, the horizontal terms of the porous resistance tensors were calculated from the mean surveyed stems diameter d = 0.03 m and distance ls = 0.5 m (Fenocchi 2015) with Zinke's (2012) approach, adopting Gebart's (1992) formula for the intrinsic permeability, Khor = 8.75 × 10−2 m2, and Zinke's (2012) own equation for the inertial parameter, βhor = 3.85 × 10−2 m−1. Null resistance was set in the vertical dimension, given the reed-like nature of the stems. A uniform normal velocity and fixed discharge conditions were set at the upstream and downstream boundaries, respectively, the inlet velocity also being adopted as an initial condition. The 3D model of the Superior Lake of Mantua, detailed in Fenocchi (2015) and Fenocchi & Sibilla (2016), was indirectly validated against historical remote-sensed distribution maps of Chlorophyll-a by Pinardi et al. (2015).

RESULTS AND DISCUSSION

Comparisons between the results of the 2D and 3D models were performed considering the depth-integrated horizontal velocities from the latter. In the simulated cases, vertical velocities are negligible compared to the horizontal, being at most in the order of 10−5 m/s.

Test elliptical lake

A preliminary analysis on the effects of the order of accuracy of the numerical schemes and of the refinement of the computational mesh on results was first performed. The IBL wind stress distribution (Figure 2) case, which is considered as reference in the following discussion, was simulated with the 2D and 3D models employing both first- and second-order approaches. The first-order 2D code was run both on the finer grid, mimicking the 3D one, and on the coarser mesh. For each model, the differences between the results obtained by either of the solution schemes proved to be negligible, compared to the differences between the models themselves. The same holds for the different mesh resolutions in the 2D model. The flow field obtained with ORSA2D is also concordant to the one given by the 2D model by Józsa et al. (2007), enforcing the validity of the adopted SWE code. Figure 5 displays the depth-averaged velocity magnitude V profiles at Section A in Figure 1 for all these simulations.
Figure 5

Depth-averaged velocity magnitude profiles at Section A in Figure 1 for the reference IBL W10,l = 10 m/s from NW conditions obtained with the 2D and 3D models for different solution schemes and grids.

Figure 5

Depth-averaged velocity magnitude profiles at Section A in Figure 1 for the reference IBL W10,l = 10 m/s from NW conditions obtained with the 2D and 3D models for different solution schemes and grids.

Due to the negligible differences, it was decided to adopt the second-order scheme for STAR-CCM + , as it converges faster, and the first-order approach, together with the coarser grid, for ORSA2D, being computationally cheaper. Furthermore, second- and first-order schemes are the standard practice for 3D CFD and 2D SWE models, respectively, so that such choice is also relevant from a practical perspective.

The depth-averaged flow fields of the 2D and 3D simulations with the selected parameters for the reference IBL case are shown in Figure 6. A rough estimation of the Hutter ratio (Equation (1)) results in Δ ≈ 7 × 10−3, i.e. the flow characters should be theoretically well approximated by depth-averaged values. Both flow fields display a single, main clockwise cell resulting from the fetch-dependent wind stress, with a secondary counter-clockwise gyre penned to the southern end, as opposed to the behaviour observed for the uniform wind stress, characterised by two main counter-rotating cells (Podsetchine & Schernewski 1999; Józsa 2006). Substantial agreement between the models is found in the resulting depth-averaged flow directions, i.e. the circulation layouts are almost coincident. However, the flow intensities are different, as the 2D model simulates higher velocities in the pelagic region, whereas similar values are attained in nearshore areas. The same behaviour holds for the uniform wind stress, ruling out the influence of the variable stress on the discrepancy. The relative error of the depth-averaged velocity magnitude of the 2D model with respect to the 3D one, ΔV2D3D,rel = (V2DV3D)/V3D, is shown for the IBL reference case in Figure 7. The largest differences are due to the slight misalignment of the centres of the gyres between the models, while the described behaviour in the open-water and littoral regions is elsewhere clear.
Figure 6

Depth-averaged velocity fields in the test elliptical lake for the reference conditions obtained with the (a) 2D model and (b) 3D model.

Figure 6

Depth-averaged velocity fields in the test elliptical lake for the reference conditions obtained with the (a) 2D model and (b) 3D model.

Figure 7

Relative error ΔV2D3D,rel field between the depth-averaged velocity magnitudes from the 2D and the 3D models for the reference conditions on the test elliptical lake.

Figure 7

Relative error ΔV2D3D,rel field between the depth-averaged velocity magnitudes from the 2D and the 3D models for the reference conditions on the test elliptical lake.

Increasing the bottom roughness coefficient in the 2D model was first attempted to improve the agreement with the depth-averaged 3D results. This is standard practice, as 2D models should theoretically employ higher resistance coefficients than 3D ones as they lump turbulent energy dispersion in the vertical plane, directly reproduced by the latter, into the friction term (e.g. Morvan et al. 2008). Hearn & Hunter (1990) stated that in flows in which the Coriolis term is negligible, as for the elliptical basin and the Superior Lake of Mantua as well, the velocities obtained with 2D models are multiple of the depth-averaged ones from 3D models by a constant, so that equivalent flow fields could be obtained by means of a calibrated increased roughness coefficient.

With the test elliptical basin, however, the use of a single increased roughness coefficient throughout the bottom is not satisfactory. For the reference conditions, a closer agreement in the pelagic area (obtained for n2D ≈ 2 · n3D = 0.05 s/m1/3) is obtained at the cost of underestimating the velocity nearshore. Therefore, only compromise solutions could be pursued, such as the one in Figure 8(a) for n2D = 1.5 · n3D = 0.0375 s/m1/3. This has to be ascribed to the characteristics of the basin, in which the boundaries have a strong influence on circulations, as is typical of real-world small basins. A variable bed roughness increase within the 2D numerical domain is feasible in particular applications, but it would be meaningless from the point of view of evaluating the prediction capabilities of depth-averaged models, as it requires previous knowledge of the actual flow field for the specific conditions, possibly obtained from a 3D model.
Figure 8

Depth-averaged velocity fields simulated with the 2D model in the test elliptical lake for the reference conditions and (a) n2D = 1.5 · n3D = 0.0375 s/m1/3, (b) νT,hor = 105 · ν = 0.089 m2/s.

Figure 8

Depth-averaged velocity fields simulated with the 2D model in the test elliptical lake for the reference conditions and (a) n2D = 1.5 · n3D = 0.0375 s/m1/3, (b) νT,hor = 105 · ν = 0.089 m2/s.

Further attempts to improve the 2D results were performed by enabling the constant νT,hor turbulence model within ORSA2D, which allows a velocity damping which is not spatially uniform, but proportional to the gradient of the horizontal velocity. Only the effects of 2D turbulence due to shear layers and gyres in the horizontal plane are directly simulated, while those of turbulence due to shear layers along the vertical direction are still lumped into the friction coefficient. However, this approach did not lead to adequate results. Large values of the horizontal eddy viscosity, e.g. νT,hor = 105 · ν = 0.089 m2/s (Figure 8(b)), are needed to damp the flow velocities. Although coherent values to the ones of the 3D model are thus obtained both in the central and in the littoral areas of the basin, the circulation layout is altered, distortion growing together with νT,hor. This suggests that horizontal turbulent processes, the only ones that can be directly accounted for in depth-averaged models, are secondary in shaping the flow field. Little improvement of 2D simulations of wind-driven flows in shallow lakes can then be achieved by modelling them, even if more accurate and complicated turbulence models, such as k-ɛ, are applied.

The mismatch between the 2D and 3D models has then to be attributed to spatially heterogeneous turbulent processes in the vertical plane. The reference IBL simulation was repeated on the modified basin with a flat bottom, yielding the same discrepancies, thus ruling out the influence of the variable bathymetry on the error. In the simulations with altered initial free surface elevation, instead a closer agreement between the models was obtained for the shallower water case, where flow circulation and wind stress transmission along the water depth are closer to the 2D schematisation. Larger dissimilarities were observed in the deeper water case, due to the converse reasons.

To identify the relevant parameters of the 3D simulations which trigger the differences with the 2D results, the relative error between the depth-averaged velocity magnitudes ΔV2D3D,rel was correlated to depth-averaged turbulence-related quantities from the 3D simulations: the turbulence intensity I = (2/3 k)1/2/, where k is the turbulent kinetic energy; the eddy viscosity νT; the non-dimensional horizontal and vertical vorticity components = (ωhor h)/V and = (ωvert h)/V, which are related to the shear layers in the vertical and horizontal planes, respectively, and are evaluated from the depth-averaged vorticity components ωhor and ωvert. The relative velocity error was also correlated to the absolute horizontal angle between the velocities at the surface and at the bottom of each stack of cells in the 3D model, |Δdirs-b| = |dir(vs)dir(vb)|, which is an efficient index of the local three-dimensional character of the flow field. The Pearson correlation coefficients are listed in Table 1 for the four discussed flow cases in the concave elliptical basin: constant W10 = 11 m/s wind, IBL W10,l = 10 m/s wind, IBL wind with η = +0.7 m, IBL wind with η = −0.7 m.

Table 1

Correlation coefficients between ΔV2D3D,rel and the selected variables from the 3D simulations on the test elliptical lake (non-significant values with p > 0.01 are in italics)

τs distribution Constant IBL IBL IBL 
η [m] + 0.7 −0.7 
I 0.360 0.342 0.368 0.522 
νT 0.271 0.261 0.188 0.061 
 0.623 0.810 0.960 0.743 
 0.321 0.100 0.376 0.163 
|Δdirs-b| 0.322 0.373 0.308 0.429 
τs distribution Constant IBL IBL IBL 
η [m] + 0.7 −0.7 
I 0.360 0.342 0.368 0.522 
νT 0.271 0.261 0.188 0.061 
 0.623 0.810 0.960 0.743 
 0.321 0.100 0.376 0.163 
|Δdirs-b| 0.322 0.373 0.308 0.429 

The highest correlations by far are obtained with , as is also noticeable from the strong resemblance between its depth-averaged field in Figure 9(a) and the velocity error map in Figure 7. In addition, the computed depth-averaged non-dimensional horizontal vorticity component is in all cases one order of magnitude larger than the vertical one. This means that the relevant turbulent effects responsible for the error of the 2D model lie in the vertical plane and are mostly related to the inaccurate assumption of uniform flow along the depth: such hypothesis, which approximates the unresolved shear layer with the logarithmic velocity profile of a uniform 1D flow, is the strict justification for Manning's resistance formula. However, in the present case, the shear layers produced a) by the wind surface stress, and b) by the variable velocity direction along the water column induced by the heterogeneous influence of wind along the depth, depart from the simplified assumption upon which the 2D model is based. The correlation coefficient to is higher for the IBL η = +0.7 m case, due to 3D effects along the vertical being more relevant in deeper waters, and conversely lower in the η = −0.7 m situation.
Figure 9

(a) Depth-averaged non-dimensional horizontal vorticity and (b) absolute flow direction difference between surface and bottom velocities |Δdirs-b| fields resulting from the 3D model of the test elliptical lake under the reference conditions.

Figure 9

(a) Depth-averaged non-dimensional horizontal vorticity and (b) absolute flow direction difference between surface and bottom velocities |Δdirs-b| fields resulting from the 3D model of the test elliptical lake under the reference conditions.

Correlations with |Δdirs-b| (Figure 9(b)) are also significant, showing that velocity damping by vertical turbulence also has an effect on the main gyre, which, despite being 2D on average, has in fact a significantly tilted rotation axis. Its centre shifts by ∼500 m from SW to NE from the surface to the bottom; the southern secondary gyre is similarly tilted. From the correlation coefficients for the IBL cases, it appears that the missing reproduction by the 2D model of the shear layers triggered by the varying velocity direction along the water column is more relevant for shallower depths, due to higher gradients being present.

Rather significant correlations are also obtained with the depth-averaged turbulence intensity: the lower values compared to may be due to this parameter being strongly dependent on the velocity magnitude, whereas the discrepancy between the models is independent from that.

The vorticity components were found to be correlated one another in all cases (the coefficient ranges r = 0.33 ÷ 0.54), although the turbulent flow is far from being isotropic. Most likely, as it occurs at a micro-scale in near-wall turbulence, also at the lake macro-scale vorticity is first produced by shear layers in the vertical plane and then partially shifted to the horizontal plane, contributing to gyres in a spatially heterogeneous fashion. This may explain why some improvements in the 2D simulations were achieved by reproducing the effects of the horizontal turbulence, however at the cost of an overall detrimental result.

The estimates of the free surface setups obtained from the models in the four conditions are very close. The differences between the maximum and minimum water elevations predicted by the 2D model were 6.5% ÷ 7.3% lower than those from the 3D model, except for the IBL case with η = +0.7 m, in which a 4.8% higher value was attained. A definite overestimation (Krámer 2006) or underestimation (Hearn & Hunter 1990) is therefore not present, showing that the 2D model is able to reproduce the global balance between wind and pressure forces, responsible for the overall surface setup.

This test case shows that 2D models can satisfactorily reproduce the depth-averaged circulation layouts obtained from 3D models on a simple basin with prevailing circulations in the horizontal plane, due to the correct reproduction of the global equilibrium between the two driving forces. Problems arise because of the non-uniform energy dissipation in the vertical plane owing to turbulent shear layers, which cannot be reproduced by 2D models, resulting in local overestimations of the depth-averaged velocities. The usefulness of 2D models would be then restricted to problems in which only the depth-averaged circulation layout is of interest, disregarding precise advection time scales and flow paths at specific depths. In fact, even in basins where circulations mainly develop in the horizontal plane, velocity variations in direction (see Figure 9(b)) and magnitude along the vertical can be significant.

Superior Lake of Mantua

A preliminary comparison was performed under no-wind and reference discharge Q = 20 m3/s conditions. The 2D and 3D models yielded almost coincident flow fields and depth-averaged velocities, also within the lotus flower island, showing a successful matching of the adopted flow resistances, despite the different approaches. This result is expected, as the flow is purely gravity-driven, the ordinary condition for the application of depth-averaged models.

In this case, the Hutter parameter (Equation (1)) approximately ranges Δ ≈ 5 × 10−3 ÷ 5 × 10−2 for increasing winds and discharges, so that the depth-averaged schematisation would still be acceptable.

Figure 10 displays the results of the simulations with easterly ordinary storm wind (the wind stress field is shown in Figure 4) and reference discharge. Both models reproduce a prevailing influence of wind on the circulation layouts, which are almost equivalent, but the 2D model yields higher flow velocities in large part of the domain. This agrees with the behaviour found in the test elliptical lake case.
Figure 10

Depth-averaged velocity fields in the Superior Lake of Mantua for ordinary storm wind and reference Q = 20 m3/s discharge conditions obtained with the (a) 2D model and (b) 3D model.

Figure 10

Depth-averaged velocity fields in the Superior Lake of Mantua for ordinary storm wind and reference Q = 20 m3/s discharge conditions obtained with the (a) 2D model and (b) 3D model.

Stronger differences are present under easterly standard wind and reference discharge (Figure 11). The upstream-directed currents and the gyres reproduced by the 3D model, which result from the interaction between wind and riverine current in shaping the circulation layout, are almost absent in the 2D simulation, which returns instead a more river-like flow field, far less affected by wind. Higher depth-averaged velocities are generally attained in the 2D simulation, mainly due to the absence of weak flow areas inside gyres. The minor influence of wind in the results of the 2D model in the standard wind case is likely due to the small wind stress ( ∼2 orders of magnitude lower than in the storm wind case), which, uniformly supplied over the water depth, is insufficient to drive the advection process, so that the riverine current largely prevails. In the 3D model, in which the wind stress is unevenly distributed along the water column, the upper layers receive a forcing higher than the average one, which locally overcomes the momentum transport of the riverine current, triggering the formation of gyres.
Figure 11

Depth-averaged velocity fields in the Superior Lake of Mantua for standard wind and reference discharge conditions obtained with the (a) 2D model and (b) 3D model.

Figure 11

Depth-averaged velocity fields in the Superior Lake of Mantua for standard wind and reference discharge conditions obtained with the (a) 2D model and (b) 3D model.

Simulations with reduced Q = 5 m3/s and with null discharges show for storm wind forcing the same differences observed for the reference discharge. Conversely, more similar layouts between the models are obtained under standard wind conditions, almost coinciding in the enclosed-lake case (Figure 12). However, here the 2D model simulates lower depth-averaged velocities in most of the domain. The |Δdirs-b| map in Figure 13 shows that the 3D model predicts strong circulations involving the whole flow depth in those regions. The reason may be found in the intrinsic limitation of depth-averaged models, which cannot represent bottom currents springing in a direction different from the mean flow, so that the resulting global displacement of water volume is underestimated. Eventually, the strong balancing currents, which arise in the lake and involve the whole water column, are weaker in the 2D model prediction. Conversely, the 2D code reproduces slightly stronger depth-averaged velocities where the 3D model predicts strong variations of the flow direction along the water depth (Figure 13), i.e. mostly inside gyres, due to the reasons mentioned for the reference discharge and standard wind case.
Figure 12

Depth-averaged velocity fields in the Superior Lake of Mantua for standard wind and null through-flow conditions obtained with the (a) 2D model and (b) 3D model.

Figure 12

Depth-averaged velocity fields in the Superior Lake of Mantua for standard wind and null through-flow conditions obtained with the (a) 2D model and (b) 3D model.

Figure 13

Absolute flow direction difference between surface and bottom velocities |Δdirs-b| field resulting from the 3D model of the Superior Lake of Mantua for standard wind and null through-flow conditions.

Figure 13

Absolute flow direction difference between surface and bottom velocities |Δdirs-b| field resulting from the 3D model of the Superior Lake of Mantua for standard wind and null through-flow conditions.

The correlation coefficients of the relative velocity magnitude error ΔV2D3D,rel to the depth-averaged turbulence-related quantities and to |Δdirs-b| for the seven discussed conditions (no-wind with reference discharge; storm and standard winds with reference, reduced or null discharges) are listed in Table 2.

Table 2

Correlation coefficients between ΔV2D3D,rel and the selected variables from the 3D simulations on the Superior Lake of Mantua (non-significant values with p > 0.01 are in italics)

Wind Null Standard Storm Standard Storm Standard Storm 
Q [m3/s] 20 20 20 
I 0.674 0.378 0.545 0.314 0.216 0.255 0.204 
νT −0.045 0.109 0.072 0.008 0.051 −0.138 0.018 
 0.320 0.712 0.702 0.646 0.652 0.498 0.888 
 0.253 0.642 0.749 0.451 0.850 0.477 0.954 
|Δdirs-b| 0.221 0.404 0.217 0.230 0.186 0.194 0.090 
Wind Null Standard Storm Standard Storm Standard Storm 
Q [m3/s] 20 20 20 
I 0.674 0.378 0.545 0.314 0.216 0.255 0.204 
νT −0.045 0.109 0.072 0.008 0.051 −0.138 0.018 
 0.320 0.712 0.702 0.646 0.652 0.498 0.888 
 0.253 0.642 0.749 0.451 0.850 0.477 0.954 
|Δdirs-b| 0.221 0.404 0.217 0.230 0.186 0.194 0.090 

The correlation coefficients to the two non-dimensional vorticity components are much closer to one another than in the test elliptical lake case, even if the simulated is still one order of magnitude larger than . The weakest correlation to vorticity was obtained for the no-wind case, where the small differences between the models are mostly correlated to the turbulence intensity I, and hence to local turbulence effects which are not simulated by the 2D model. In the standard wind cases, the correlation of ΔV2D3D,rel to the depth-averaged non-dimensional horizontal vorticity is moderately higher than to the vertical one, while the opposite happens in storm wind simulations. These occurrences are likely due to the higher interaction between circulation and shear processes in the Superior Lake of Mantua case, which especially develops for stronger winds and is triggered by the irregular bathymetry and the through-flowing current.

Such interaction is highlighted in the 3D simulations by the higher correlation coefficients between the two depth-averaged non-dimensional vorticity components compared to the test elliptical lake, whose range here is r = 0.48 ÷ 0.87. While these coefficients are almost constant with discharge for standard wind, increasing values are obtained for decreasing discharges under storm wind conditions. This shows how the riverine current and the weak winds hinder the regular transmission of wind stress along the water depth, which on the other hand is more easily attained under strong wind conditions and leads to a more isotropic turbulent flow.

The correlation coefficient to |Δdirs-b| increases with discharge and is higher under standard wind conditions (Table 2). This highlights that the depth-averaged assumption of constant velocity direction along the vertical is less adequate when the interaction between wind and riverine current in shaping circulations is stronger.

However, the analysis of the discrepancies between the two models cannot be limited here to the differences between the velocity magnitudes, as different flow directions are obtained in some conditions. The same parameters of the previous analysis were hence also correlated to the absolute values of the difference between the depth-averaged velocity directions obtained from the 2D and the 3D models, |Δdir(V)2D3D| = |dir(V2D)dir(V3D)|, as listed in Table 3. The highest correlation coefficients are obtained to |Δdirs-b|, proving that the error in the depth-averaged velocity estimation by the 2D model is significant where the differences in the velocity direction along the flow depth are relevant. In addition, correlation coefficients to , despite being much lower than those for the velocity magnitude differences, are still higher than those to (except for the no-wind case). This suggests once again that the inaccuracy of the 2D model is mostly linked to turbulent processes in the vertical plane, in particular to the shear produced by varying velocity direction along the water column.

Table 3

Correlation coefficients between |Δdir(V)2D3D| and the selected variables from the 3D simulations on the Superior Lake of Mantua (all values are significant with p < 0.01)

Wind Null Standard Storm Standard Storm Standard Storm 
Q [m3/s] 20 20 20 
I 0.040 0.122 0.359 0.144 0.177 0.239 0.374 
νT −0.157 0.297 0.120 0.365 0.134 0.150 0.081 
 0.136 0.215 0.412 0.121 0.379 0.295 0.232 
 0.162 0.130 0.283 0.115 0.232 0.253 0.153 
|Δdirs-b| 0.461 0.576 0.694 0.545 0.706 0.526 0.689 
Wind Null Standard Storm Standard Storm Standard Storm 
Q [m3/s] 20 20 20 
I 0.040 0.122 0.359 0.144 0.177 0.239 0.374 
νT −0.157 0.297 0.120 0.365 0.134 0.150 0.081 
 0.136 0.215 0.412 0.121 0.379 0.295 0.232 
 0.162 0.130 0.283 0.115 0.232 0.253 0.153 
|Δdirs-b| 0.461 0.576 0.694 0.545 0.706 0.526 0.689 

Similar estimates of the free surface setups were still obtained from the two models in all conditions, confirming the considerations drawn for the elliptical lake case.

Even if the complicated flow fields resulting from the real-world geometry of the Superior Lake of Mantua make the interpretation of results more difficult than for the test elliptical basin, some conclusions can be drawn from the analysis of this fluvial lake case. Two-dimensional hydrodynamic simulations of shallow fluvial lakes with prevailing circulations in the horizontal plane and complex bathymetry may face the following limitations. (1) When strong interactions between the wind and the riverine current in shaping circulations are present, i.e. for low wind intensities and significant discharges, the assumed uniform distribution of the wind stress along the water depth can cause an underestimation of the influence of wind on circulations, resulting in a flow structure with much less developed circulations. (2) When wind is the dominant driving force, i.e. for strong intensities or low discharges, the error appears to be limited to different velocity magnitudes. Under strong wind conditions, higher velocities can be simulated, due to the lack of velocity damping by the spatially heterogeneous vertical turbulence, as was observed also in the elliptical lake test case. For low winds and discharges, lower velocities can be conversely simulated, due to the impossibility of reproducing any flow direction variation over the water column, which in turn causes an underprediction of the displaced water volume and of the resulting balancing currents over the whole flow depth. (3) When the riverine current is dominant, valid results are produced, as the mainly gravity-driven flow is effectively simulated by depth-averaged models.

CONCLUSIONS

The application of 2D models to wind-induced flows in shallow lakes with prevailing circulations in the horizontal plane, though intuitively valid, can result in various levels of inaccuracy, due to the treatment of wind as a volume force over the flow depth, neglecting dispersion by vertical turbulence, and to the synergetic, depth-averaged motion of the water column. In these basins, the significant three-dimensional characters of shear, unrestricted by stratification, condition circulation dynamics and impair depth-averaged simulations. The discrepancies with 3D simulations are higher for complex bathymetries, low wind intensities and deeper waters, where conditions depart strongly from the ideal 2D ones. These errors cannot be reduced by any expedient, as the spatial heterogeneity of turbulence and shear in the vertical plane strongly depends on the flow field, which is not previously known.

In ordinary shallow lakes, the error of 2D models is possibly limited to an inaccurate velocity magnitude prediction, while the circulation layout is generally correctly reproduced. In fluvial lakes, when strong interactions between the wind and the through-flowing current in determining circulations are present, 2D models result instead in a strong underestimation of the effect of wind, and produce flow fields which are far more similar to riverine ones. In addition to the inaccurate reproduction of the depth-averaged flow fields, this may result in an underprediction of the water residence times in fluvial lakes, as the development of water-trapping gyres is underestimated. As a practical example, the export of nutrients would be overrated, resulting in an underpredicted eutrophication potential. To an extreme consequence, the error of depth-averaged models may lead to the misinterpretation of fluvial lake environments for enlarged river ones, ignoring the influence of wind, thus leading to significant mistakes.

Therefore, use of 3D models is especially supported for shallow fluvial lakes with main circulations in the horizontal plane and no prevailing influence of wind or through-flowing current, i.e. typically for low wind intensities, which are characteristic of many basins in lowland regions. This conclusion is rather unexpected, as the main limits to the application of depth-averaged models to fluvial lakes appear to be the same features which make them similar to rivers, for which these models were first developed.

ACKNOWLEDGEMENTS

The authors would like to thank J. Józsa and T. Krámer for providing some of their own results on the test elliptical lake case. The contributions of the anonymous reviewers are also gratefully acknowledged.

REFERENCES

REFERENCES
Bermúdez
A.
Dervieux
A.
Desideri
J. A.
Vázquez
M. E.
1998
Upwind schemes for the two-dimensional shallow water equations with variable depth using unstructured meshes
.
Computer Methods in Applied Mechanics and Engineering
155
,
49
72
.
CD-adapco
2014
STAR-CCM+ v9.02 User Guide
.
Melville
,
USA
.
CERC (Coastal Engineering Research Center)
1984
Shore Protection Manual
.
4th edn
.
Waterways Experiment Station, US Army Corps of Engineers
,
Vicksburg
,
USA
.
Christensen
B. A.
1992
Replacing hydraulic radius in Manning's formula in open channels
. In:
Channel Flow Resistance: Centennial of Manning's Formula
.
Yen
B. C.
(ed.).
Water Resources Publications
,
Highlands Ranch
,
USA
, pp.
271
287
.
Ciraolo
G.
Lipari
G.
Napoli
E.
Józsa
J.
Krámer
T.
2004
Three-dimensional numerical analysis of turbulent wind-induced flows in the Lake Balaton (Hungary)
. In:
Shallow Flows – Research Presented at the International Symposium on Shallow Flows, Delft, Netherlands, 2003
.
Uijttewaal
W. S. J.
Jirka
G. H.
(eds).
CRC Press/Balkema
,
Leiden
,
The Netherlands
, pp.
661
669
.
Costabile
P.
Macchione
F.
2015
Enhancing river model set-up for 2-D dynamic flood modelling
.
Environmental Modelling & Software
67
,
89
107
.
Costabile
P.
Macchione
F.
Natale
L.
Petaccia
G.
2015
Flood mapping using LIDAR DEM. Limitations of the 1-D modeling highlighted by the 2-D approach
.
Natural Hazard
77
,
181
204
.
Fabian
J.
Budinski
L.
2013
Horizontal mixing in the shallow Palic Lake caused by steady and unsteady winds
.
Environmental Modeling & Assessment
18
,
427
438
.
Fenocchi
A.
2015
Circulation Dynamics in a Shallow Fluvial Lake – The Case of the Superior Lake of Mantua
.
PhD Thesis
,
Department of Civil Engineering and Architecture, University of Pavia
,
Pavia
,
Italy
.
Gebart
B. R.
1992
Permeability of unidirectional reinforcements for RTM
.
Journal of Composite Materials
26
,
1100
1133
.
Hirsch
C.
1988
Numerical Computation of Internal & External Flows, Volume 2: Computational Methods for Inviscid and Viscous Flows
.
J. Wiley & Sons
,
Chichester
,
UK
.
Hutter
K.
1984
Fundamental equations and approximations
. In:
Hydrodynamics of Lakes
.
Hutter
K.
(ed.).
CISM Lectures No. 286, Springer
,
New York
,
USA
, pp.
1
38
.
Józsa
J.
2006
Shallow Lake Hydrodynamics – Theory, Measurement and Numerical Model Applications
.
Mundus-Euroaquae lecture notes, Department of Hydraulic and Water Resources Engineering, Budapest University of Technology and Economics
,
Budapest, Hungary
.
Józsa
J.
2014
On the internal boundary layer related wind stress curl and its role in generating shallow lake circulations
.
Journal of Hydrology and Hydromechanics
62
,
16
23
.
DOI: 10.2478/johh-2014-0004.
Józsa
J.
Krámer
T.
Napoli
E.
2007
The impact of terrain roughness and water level changes on wind-induced shallow lake circulation patterns
. In:
Proc. Fifth International Symposium on Environmental Hydraulics
,
Tempe
,
USA
,
4–7 December 2007
.
Krámer
T.
2006
Solution-adaptive 2D Modelling of Wind-induced Lake Circulation
.
PhD Thesis
,
Department of Hydraulic and Water Resources Engineering, Budapest University of Technology and Economics
,
Budapest, Hungary
.
Morvan
H.
Knight
D.
Wright
N.
Tang
X.
Crossley
A.
2008
The concept of roughness in fluvial hydraulics and its formulation in 1D, 2D and 3D numerical simulation models
.
Journal of Hydraulic Research
46
,
191
208
.
Osher
S.
1985
Convergence of generalized MUSCL schemes
.
SIAM Journal on Numerical Analysis
22
,
947
961
.
Ottesen Hansen
N. E.
1979
Effects of boundary layers on mixing in small lakes
. In:
Hydrodynamics of Lakes
.
Graf
W. H.
Mortimer
C. H.
(eds).
Developments in Water Science 11, Elsevier
,
Amsterdam
,
The Netherlands
, pp.
341
356
.
Petaccia
G.
Soares-Frazão
S.
Savi
F.
Natale
L.
Zech
Y.
2010
Simplified versus detailed Two-dimensional approaches to transient flow modeling in Urban Areas
.
Journal of Hydraulic Engineering
136
,
262
266
.
Petryk
S.
Bosmajian
G.
1975
Analysis of flow through vegetation
.
ASCE Journal of the Hydraulics Division
101
,
871
884
.
Roe
P. L.
1981
Approximate Riemann solvers, parameter vectors and difference schemes
.
Journal of Computational Physics
43
,
357
372
.
Schimmelpfennig
S.
Kirillin
G.
Engelhardt
C.
Nützmann
G.
2012
Effect of wind-driven circulation on river intrusion in Lake Tegel: modeling study with projection on transport of pollutants
.
Environmental Fluid Mechanics
12
,
321
339
.
Shih
T. H.
Liou
W. W.
Shabbir
A.
Yang
Z.
Zhu
J.
1994
A New k-ε Eddy Viscosity Model for High Reynolds Number Turbulent Flows-Model Development and Validation
.
NASA Technical Memorandum 106721, Institute for Computational Mechanics in Propulsion and Center for Modeling of Turbulence and Transition, NASA Lewis Research Center
,
Cleveland, USA
.
Zinke
P.
2012
Application of a porous media approach for vegetation flow resistance
. In:
Proc. River Flow 2012
.
Murillo Muñoz
R.
(ed.).
San José, Costa Rica, 5–7 September 2012
,
CRC Press/Balkema
,
Leiden
,
The Netherlands
,
Vol. 1
, pp.
301
308
.