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.

*T*=

_{h}*H*

^{2/}4

*ν*is the characteristic timescale of momentum transport by turbulence (

_{T,vert}*ν*being the vertical eddy viscosity along the characteristic depth

_{T,vert}*H*), which damps the irregularities in the velocity profile in the vertical dimension, while

*T*=

_{adv}*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

*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: with: 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,

*S*the bed slope,

_{0}*S*the friction slope, computed from Manning's formula,

_{f}*τ*the wind surface stress,

_{s}*τ*the turbulent stress,

_{T}*ρ*the density and

*ν*the horizontal eddy viscosity. In the numerical scheme, the bed slope terms are upwinded (Bermúdez

_{T,hor}*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: in which

*ρ*= 1.225 kg/m

_{a}^{3}is the air density and

*c*the wind drag coefficient related to the wind intensity at 10 m height

_{10}*W*.

_{10}#### 3D model

*ν*the kinematic viscosity,

*ν*the global eddy viscosity and

_{T}*B*the summation of the external body forces. In the present simulations, Equation (6) comprises: the vertical gravity force

*f*, where

_{g}*γ*is the specific weight of water and ∀

*the volume of a grid cell; the horizontal wind force*

_{cell}*f*, which acts only over the surface cells,

_{s}*z*being the vertical mesh resolution; and the vegetation flow resistance force

_{g}*f*, 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,

_{p}*μ*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

*p*calculated by the 3D code as

_{s}*Δη*=

*p*/

_{s}*γ*, 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

*et al.*(2007), which has bottom elevation

*z*defined as function of the

_{b}*x*and

*y*coordinates: in which lengths are in metres and shores are placed where

*z*= −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

_{b}*η*= 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

*z*= −1.5 m was employed in a particular simulation, to determine the relevance of the variable bathymetry on results.

_{b}*c*and

_{10}*W*terms in Equation (4), hence

_{10}*τ*, to increase with fetch. A wind intensity over land

_{s}*W*= 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

_{10,l}*z*= 0.15 m (tall grass, see Józsa

_{0,l}*et al.*2007). Details about the wind stress distribution model and its implementation are present in Fenocchi (2015). A uniform wind stress distribution for

*W*= 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

_{10}*c*= 1.45 × 10

_{10}^{−3}having been estimated with the formula by Wu (1980), yielding

*τ*= 0.225 N/m

_{s}^{2}.

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/m^{1/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 *x _{g}* = 4 m and

*z*= 0.2 m, respectively (i.e. with aspect ratio

_{g}*r*=

_{g}*x*/

_{g}*z*= 20) and the free surface boundary being the water-at-rest elevation. The roughness height on the bottom was set to

_{g}*k*= 0.075 m, matching the original roughness in the 2D model according to the relation

_{s}*n*=

*k*

_{s}^{1/6/}26 (e.g. Christensen 1992).

#### Superior Lake of Mantua

*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.

Two wind intensity conditions were considered, based on the analysis of the local climate (Fenocchi 2015): 1) standard wind *W _{10,l}* = 1.94 m/s; and 2) ordinary storm wind

*W*= 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

_{10,l}*Q*= 20 m

^{3}/s, mean annual value for the basin, but simulations with reduced flow rate

*Q*= 5 m

^{3}/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.

*Shore Protection Manual*(CERC 1984). The aerodynamic roughness of land was still set to

*z*= 0.15 m, while that of the lotus flower cover to

_{0,l}*z*= 0.08 m. The resulting wind stress distribution for ordinary storm conditions is shown in Figure 4.

_{0,veg}*n*= 0.025 s/m

^{1/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: in which

*n*is the bare bed roughness,

*C*= 1.1 the drag coefficient (very sparse canopy, see Wang & Wang (2011)),

_{d}*a*= 0.12 m

^{−1}the mean surveyed vegetation density for the lotus flower in the Superior Lake of Mantua (Fenocchi 2015) and

*h*= 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.

_{m,veg}In the 3D case, simulations were performed on a grid with 847,843 trimmed parallelepiped elements having *x _{g}* = 6 m horizontal resolution and

*z*= 0.3 m vertical resolution (

_{g}*r*= 20), the free surface boundary being the control elevation. The bottom roughness height was still set to

_{g}*k*= 0.075 m. For the lotus flower island, the horizontal terms of the porous resistance tensors were calculated from the mean surveyed stems diameter

_{s}*d*= 0.03 m and distance

*l*= 0.5 m (Fenocchi 2015) with Zinke's (2012) approach, adopting Gebart's (1992) formula for the intrinsic permeability,

_{s}*K*= 8.75 × 10

_{hor}^{−2}m

^{2}, and Zinke's (2012) own equation for the inertial parameter,

*β*= 3.85 × 10

_{hor}^{−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

*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.

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.

^{−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,

*ΔV*

_{2D}_{−}

*= (*

_{3D,rel}*V*–

_{2D}*V*)/

_{3D}*V*, 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.

_{3D}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.

*n*≈ 2 ·

_{2D}*n*= 0.05 s/m

_{3D}^{1/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

*n*= 1.5 ·

_{2D}*n*= 0.0375 s/m

_{3D}^{1/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.

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.

*ν*= 10

_{T,hor}^{5}·

*ν*= 0.089 m

^{2}/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

*ν*. 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

_{T,hor}*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 *ΔV _{2D}*

_{−}

*was correlated to depth-averaged turbulence-related quantities from the 3D simulations: the turbulence intensity*

_{3D,rel}*I*= (2/3

*k*)

^{1/2}/, where

*k*is the turbulent kinetic energy; the eddy viscosity

*ν*; the non-dimensional horizontal and vertical vorticity components = (

_{T}*ω*)/

_{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

*ω*and

_{hor}*ω*. 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, |

_{vert}*Δdir*| = |

_{s-b}*dir(v*−

_{s})*dir(v*|, 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

_{b})*W*= 11 m/s wind, IBL

_{10}*W*= 10 m/s wind, IBL wind with

_{10,l}*η*= +0.7 m, IBL wind with

*η*= −0.7 m.

τ distribution
. _{s} | Constant . | IBL . | IBL . | IBL . |
---|---|---|---|---|

η [m]
. | 0 . | 0 . | + 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 | |

| Δdir|_{s-b} | 0.322 | 0.373 | 0.308 | 0.429 |

τ distribution
. _{s} | Constant . | IBL . | IBL . | IBL . |
---|---|---|---|---|

η [m]
. | 0 . | 0 . | + 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 | |

| Δdir|_{s-b} | 0.322 | 0.373 | 0.308 | 0.429 |

*η*= +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.

Correlations with |*Δdir _{s-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 m^{3}/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.

*Q*= 5 m

^{3}/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 |

*Δdir*| 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.

_{s-b}The correlation coefficients of the relative velocity magnitude error *ΔV _{2D}*

_{−}

*to the depth-averaged turbulence-related quantities and to |*

_{3D,rel}*Δdir*| 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.

_{s-b}Wind . | Null . | Standard . | Storm . | Standard . | Storm . | Standard . | Storm . |
---|---|---|---|---|---|---|---|

Q [m^{3}/s]
. | 20 . | 20 . | 20 . | 5 . | 5 . | 0 . | 0 . |

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 | |

| Δdir|_{s-b} | 0.221 | 0.404 | 0.217 | 0.230 | 0.186 | 0.194 | 0.090 |

Wind . | Null . | Standard . | Storm . | Standard . | Storm . | Standard . | Storm . |
---|---|---|---|---|---|---|---|

Q [m^{3}/s]
. | 20 . | 20 . | 20 . | 5 . | 5 . | 0 . | 0 . |

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 | |

| Δdir|_{s-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 *ΔV _{2D}*

_{−}

*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.*

_{3D,rel}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 |*Δdir _{s-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) _{2D}*

_{−}

*| = |*

_{3D}*dir(V*−

_{2D})*dir(V*|, as listed in Table 3. The highest correlation coefficients are obtained to |

_{3D})*Δdir*|, 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.

_{s-b}Wind . | Null . | Standard . | Storm . | Standard . | Storm . | Standard . | Storm . |
---|---|---|---|---|---|---|---|

Q [m^{3}/s]
. | 20 . | 20 . | 20 . | 5 . | 5 . | 0 . | 0 . |

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 | |

| Δdir|_{s-b} | 0.461 | 0.576 | 0.694 | 0.545 | 0.706 | 0.526 | 0.689 |

Wind . | Null . | Standard . | Storm . | Standard . | Storm . | Standard . | Storm . |
---|---|---|---|---|---|---|---|

Q [m^{3}/s]
. | 20 . | 20 . | 20 . | 5 . | 5 . | 0 . | 0 . |

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 | |

| Δdir|_{s-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.