Modelling flows in shallow (fluvial) lakes with prevailing circulations in the horizontal plane: limits of 2D compared to 3D models

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.

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.

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_s calculated by the 3D code as Δη = 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.

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

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_g = 0.2 m, respectively (i.e. with aspect ratio r_g = x_g/z_g = 20) and the free surface boundary being the water-at-rest elevation. The roughness height on the bottom was set to k_s = 0.075 m, matching the original roughness in the 2D model according to the relation n = k_s^1/6/26 (e.g. Christensen 1992).

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,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 m³/s, mean annual value for the basin, but simulations with reduced flow rate Q = 5 m³/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 the 3D case, simulations were performed on a grid with 847,843 trimmed parallelepiped elements having x_g = 6 m horizontal resolution and z_g = 0.3 m vertical resolution (r_g = 20), the free surface boundary being the control elevation. The bottom roughness height was still set to k_s = 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 l_s = 0.5 m (Fenocchi 2015) with Zinke's (2012) approach, adopting Gebart's (1992) formula for the intrinsic permeability, K_hor = 8.75 × 10⁻² m², and Zinke's (2012) own equation for the inertial parameter, β_hor = 3.85 × 10⁻² m⁻¹. 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).

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⁻⁵ m/s.

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.

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.

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 = 10⁵ · ν = 0.089 m²/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 ΔV_2D−3D,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, |Δdir_s-b| = |dir(v_s)−dir(v_b)|, 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 W₁₀ = 11 m/s wind, IBL W_10,l = 10 m/s wind, IBL wind with η = +0.7 m, IBL wind with η = −0.7 m.

Table 1

Correlation coefficients between ΔV_2D−3D,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 .	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
|Δdirs-b|	0.322	0.373	0.308	0.429

τs distribution .	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
|Δdirs-b|	0.322	0.373	0.308	0.429

View Large

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.

A preliminary comparison was performed under no-wind and reference discharge Q = 20 m³/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⁻³ ÷ 5 × 10⁻² for increasing winds and discharges, so that the depth-averaged schematisation would still be acceptable.

The correlation coefficients of the relative velocity magnitude error ΔV_2D−3D,rel to the depth-averaged turbulence-related quantities and to |Δdir_s-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 ΔV_2D−3D,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 .	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
|Δ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 .	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
|Δdirs-b|	0.221	0.404	0.217	0.230	0.186	0.194	0.090

View Large

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−3D,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 |Δ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_3D)|, as listed in Table 3. The highest correlation coefficients are obtained to |Δdir_s-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)_2D−3D| 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 .	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
|Δ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 .	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
|Δdirs-b|	0.461	0.576	0.694	0.545	0.706	0.526	0.689

View Large

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.

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.

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.