Interaction between depth variation and turbulent diffusion in depth-averaged vorticity equations

Steady, depth-averaged, shallow water vorticity transport equations including advection, surface and bed shear stresses, and turbulent diffusion effects are written out in vorticity-velocity and stream function formalisms. The Boussinesq approximation is used to represent turbulent stresses in the effective stress tensor. We consider two different forms of the curl of the effective stress tensor: its complete form and the commonly used form neglecting the terms expressing interaction with variable water depth. After deriving the two equations in vorticity-velocity formalism, we recast the equations into stream function formalism, revealing all the internal effects associated with variable water depth. We examine the differences between the models through analytical solutions of the stream function equations for simple but realistic flows. The solutions are validated with CFD simulations.


Introduction
The shallow water equations are widely used to model and describe horizontal motions of environmental fluid flows.In the depth-integrated or depth-averaged form of the equations, the so-called effective stress tensor term includes the effects of molecular diffusion, turbulent diffusion and dispersion.Certain largescale ( 1 km) environmental fluid motions can be accurately described without considering the effective stress tensor [7,15,26,28,31].However, at smaller spatial scales, the contribution of the effective stress tensor increases among the different flow driving impacts.A possible example of a smaller spatial scale is a width of a river in a bend (100 m... 2 km), where the effect of secondary flow must be considered (with the dispersion part of the effective stress tensor) [3,12,13].Besides the secondary flow, the turbulent diffusion is another effect modelled with the effective stress tensor.The turbulent diffusion is commonly modelled employing the Reynolds stresses with the Boussinesq approximation that uses the eddy viscosity concept [4,5,12,13,19,33].
Interaction terms between the water depth variations-we call those depth-gradients (even though there are higher order derivatives of the depth function as well)-and the flow variables (depth-averaged velocities and vorticity) arise when the Boussinesq approximation is applied.There are only a few studies where all these joint effects have been considered.As an example, the interaction effects have been applied in [6,32] to model transient dissipating vortex structures over underwater trenches and slopes.On the other hand, there are many studies in which these interaction terms are simply neglected [14,17,21,23,25,30,34,36,37].Kuipers and Vreugdenhil [19] discuss neglecting the interaction effects by reasoning that the magnitude of their joint effect is negligible for large-scale steady or slowly varying motions.
The goal of this study is to reveal the impact of these interactions from the viewpoint of the equations and their solutions.Considering the equations, we compare not only the vorticity equations when the full and the simplified Boussinesq approximation is applied, but we also write out the equations in vorticity-velocity, and stream function formalisms.The solutions of the vorticity equations (employing the full and the simplified Boussinesq approximations) are compared for two different, realistic flow configurations.The first example is an open channel flow, while the second one is a wind-induced flow.
In Sect.2, we introduce some notations, and in Sect.3, two vorticity equations are derived.The two equations differ in how they consider the interaction between the depth gradient and the effective stress tensor (in other words, by the fully applied or the simplified Boussinesq approximation).In Sect.4, the Boussinesq approximation is applied to the turbulent stresses of the effective stress tensor with a constant eddy viscosity coefficient.We write down all the terms arising from the Boussinesq approximation with the velocity, vorticity and depth function triplets.In Sect.5, we turn the velocity-vorticity-depth function form of the vorticity equations into stream function-depth function form.We match the terms of the vorticity equations between the two formalisms and discuss the algebraic and geometrical contents.
In Sects.6-8, we derive a simple flow pattern and apply our vorticity balance equations to two realistic flows: a gravity-driven flow in a channel and a wind-driven flow at a shoreline.We examine the differences in the turbulent diffusion models (between the fully applied and the simplified Boussinesq approximation) by comparing analytical cross-sectional velocity profiles given as the solutions of the vorticity equations for each flow situation.In Sect.9, we describe model validation against simulations by computational fluid dynamics (CFD) for the two flow cases.

Notations
In this paper, we use scalar-and vector fields of the 2D (horizontal) plane, in Cartesian coordinates.We clarify a few notations below.By the gradient of a scalar field F(x, y), we mean the following column vector: . ( The Poisson-bracket of the scalar fields F(x, y) and G(x, y) is the skew-symmetric, bilinear operator.We use the notations of cross product and curl in a not formally correct way, meaning scalar-valued quantities by them.Let a and b be vector fields with components a x (x, y), a y (x, y) and b x (x, y), b y (x, y), respectively.We identify the cross product and the curl with the following determinants: The symbol • is used for both vector and matrix operations: The first two operations above (from the left to the right) are the dot product of two vector fields and the divergence of a vector field, respectively.The third operation is the right multiplication of an M four-dimensional square matrix by a (column) vector, resulting in a column vector.The fourth operation is the left multiplication of the previous result by a row vector, resulting in a scalar quantity.The fifth example is the divergence of the tensor N.
3 The vorticity equations

The depth-averaged shallow water momentum equations
Let h(x, y) denote the water depth as shown in Fig. 1 and defined as the sum of the vertical distance between the still water level and the water bed, h 0 (x, y), and the local vertical deviation of the (steady) water surface η(x, y): h(x, y) = h 0 (x, y) + η(x, y). (5) Thus, the bathymetry-the shape of the water bed-is represented by h 0 (x, y).Let denote the Cartesian horizontal components of the 3D velocity field, where ū(x, y, z) and v(x, y, z) are the local time-averaged velocity components, u (x, y, z) and v (x, y, z) are the fluctuations around the means.The depth-averaged velocity components, U (x, y) and V (x, y), are defined as [4,9,12,13,19,22,33,35] The steady, depth-averaged shallow-water momentum equations and continuity equation are [4,12,13,19,33] where f , g and ρ are the Coriolis coefficient, the gravitational acceleration and the density of the fluid, respectively; τ wx (x, y) and τ wy (x, y) are the components of the wind stress field (surface stresses), τ bx (x, y) and τ by (x, y) are the components of the bed stress field (bottom stresses), T x x (x, y), T xy (x, y), T xy (x, y) and T yy (x, y) are elements of the effective stress tensor.Equations ( 8)-( 9) in vectorial form read as where

The vorticity form of the momentum equations
Taking the curl of the momentum equation ( 10) yields the scalar vorticity equation (by definition (3)): where we used the ∇ ×(∇η) = 0 identity.We introduce the depth-averaged vorticity ζ(x, y), which represents the rotation of a vertical fluid filament of depth h(x, y) [4] and interpreted as the following scalar by definition (3) The left-hand side of the vorticity equation ( 13) can now be written as [17,37] Here, V • ∇ζ is the advection of the depth-averaged vorticity by the flow, a vorticity transport term.The term ζ ∇ • V is a 2D counterpart of the 3D vortex-stretching, a vorticity source proportional to the divergence of the depth-averaged velocity field.This divergence is caused by depth variations as seen by expanding the continuity equation (11) as Rewriting the vorticity equation (13) with Eqs. ( 15) and ( 16) results Equation ( 17) describes a steady vorticity balance between the advection of the depth-averaged vorticity, the vortex stretching, the curl of the Coriolis force, the curl of the surface and the bed stresses, the curl of the effective stress tensor, and their interactions with the depth-gradients.Depth gradient can occur by changes in both the water surface and the bed topography.In this study, only the gradually changing depth functions are considered.In other words, we require from h that to be differentiable at any location.Therefore, we do not discuss any steeply varying steady flows (to which the conservative form of the shallow water equations would be necessary).

Basic equations of the study
The vorticity definition (14), the vorticity equation (17) and the continuity equation (11) form one system of equations for our study: We will refer to system (18) as FDGI : full depth-gradient interactions.
With the name above we mean the full depth-gradient interactions in the curl of the effective stress tensor term (the last term on the right-hand side of the vorticity equation of system (18)).We consider another system of equations where the interaction terms between the depth-gradient and the effective stress tensor are neglected, i.e., from which System (18) can be replaced by We will refer to system (20) as SDGI : simplified depth-gradient interactions.
System (20) provides our other base model.We will use the names FDGI and SDGI not only for the reformulations but also for the solutions of systems ( 18) and (20), respectively.We will also write out the difference FDGI − SDGI in order to emphasize the distinction between our two models, and so the physical background and mathematical formulation of the interaction effects between the effective stress tensor and the depth gradients.For the description of this difference, the curl of the effective stress tensor term is resolved as The first term on the right-hand side of Eq. ( 21) is found in SDGI.The second and third terms contain the depth-gradient interaction effects.The right multiplication of the effective stress tensor with the depth gradient, T • ∇h appears in both terms.This operation can be seen as a scale and rotation of the depth gradient vector by the effective stress tensor.In the second and third terms on the right-hand side of Eq. ( 21), we find the curl and the cross product of this transformed depth gradient, respectively.The difference between our models is characterized by these last two terms:

The Boussinesq approximation and its application
In this section, we discuss briefly the components of the effective stress tensor T. We apply the Boussinesq approximation to the Reynolds stresses and rewrite the curl of the effective stress tensor terms (the last terms on the right-hand sides) of the vorticity equations of FDGI and SDGI.We are aware that the validity and efficiency of the Boussinesq approximation are debatable [19,24,27].However, the Boussinesq approximation is widely used for various flow cases [1,4,5,8,10,12,13,[17][18][19]30,31,33].Our goal is not to investigate the validity or the accuracy of the Boussinesq approximation but to clarify the depth gradient interaction effects arising from its application.

The effective stress tensor
The elements of the effective stress tensor T are [4,12,13,19,33] where ν is the kinematic viscosity.The effective stresses consist of three different parts.The terms express the effects of the laminar viscous friction according to the widely used Newtonian fluid model.The terms are referred to as turbulent friction or Reynolds-stresses.The terms are the differential advection terms (also called as dispersion stresses) that occur due to the depth-averaging process.These terms express the effects of deviations (fluctuations) between the local mean 3D velocities and the depth-averaged velocities in a vertical sense.This effect results in horizontal (streamwise and lateral) momentum transport.The dispersion terms are considered in some general studies [5,12,13] and specific cases when the vertical variability of the velocities is important, such as in models of river bends [3,11].

Application of the Boussinesq approximation
A widely used expression of the Reynolds-stresses is the so-called Boussinesq eddy viscosity concept (Boussinesq approximation) [1,4,5,8,10,12,13,[17][18][19]30,31,33], which assumes that the Reynolds-stresses are in the same form as the Newtonian viscosity law for the laminar viscous friction terms (24).We apply the Boussinesq approximation for the Reynolds stresses and neglect the differential advection terms in the effective stress tensor components (23) to get Here, ν e is the eddy viscosity coefficient.The more significant effect of turbulent diffusion than molecular diffusion is achieved by the eddy viscosity coefficient.The value of the eddy viscosity coefficient can be larger by a few orders of magnitudes than the kinematical viscosity coefficient.We consider ν e to be constant, to clarify see the impacts of the depth-gradient related interactions.The Boussinesq eddy viscosity concept is essential in environmental hydrodynamics, even though its validity and efficiency are debatable [19,24,27].By introducing the following velocity gradient tensor: we obtain the following theorem.

Theorem 1
The curl of the effective stress tensor term (21), with the Boussinesq eddy viscosity approximation (27), yields the following form with the depth-averaged vorticity ζ and the depth-averaged velocity V: Proof The proof is provided in Appendix.
The single term in the right-hand side of Eq. ( 29) is called the diffusion of the vorticity.The diffusion of the vorticity is the only term that represents the curl of the effective stress tensor term in many different applications, e.g., [14,17,21,23,25,30,34,36,37].This simplification can be motivated by computational reasons and justified by the negligible effect of the depthgradient interactions with the Reynolds stresses [19].However, as we will show, there are depth-gradient interactions in Eqs.(29)(30)(31), due to hidden depth-gradient influences of the variables.We can rewrite now Eq. ( 22) which expresses the difference between our models, and so the interaction effects between the depth gradient and the effective stress tensor: Among the terms of the right-hand side of Eq. ( 32), we find the curl and the cross product of the transformed depth gradient vector by the velocity gradient tensor ∇V • ∇h, the dot product of the depth gradient and the vorticity gradient, and the scale of the vorticity by the Laplacian of the depth function and the magnitude of the depth gradient.

Stream function formalism of the vorticity equations
The depth-averaged velocity and vorticity, V and ζ carry the effects of the variable water depth by their definitions ( 7) and ( 14) [12,13,19,33].Consequently, the velocity, vorticity, and their derivatives carry implicit depth-gradient effects.For the same reason, the advection of the vorticity V • ∇ζ and the diffusion of the vorticity ζ terms include implicit depth-gradient effects as well.To reveal all the hidden influences of the depth variations in Systems FDGI and SDGI-( 18) and ( 20), respectively-we turn them to stream function formalism.

The stream function
We define the depth-integrated velocity field q (also known as specific flow rate or specific discharge) as and its stream function ψ(x, y) [4,[14][15][16] as The depth-averaged velocity field becomes One advantage of the stream function formalism is that the continuity equation is automatically satisfied.We can express the divergence of the depth-averaged velocity field with the stream function (35) from the continuity equation ( 16) as The Poisson-bracket {•, •} are defined by Eq. ( 2).The value of the Poisson-bracket at a certain location is determined by the intersection angle of the contour lines of its arguments.The value of (36) depends on the intersection angle of the ψ = const.streamlines and the h = const.contour lines of the depth function, the so-called isobath lines [17].The expression is zero where contour lines coincide, which also implies zero divergence for the depth-averaged velocity.Substituting the stream function ( 35) into the depth-averaged vorticity formula ( 14) results in the following expression for the depth-averaged vorticity in terms of the stream function and the depth function [14,36]: The influence of the depth gradient is contained in the second term of the right-hand side of (37).Its sign is determined by the relative orientation of the depth gradient and the gradient of the stream function.If the streamlines are perpendicular to the isobath lines, then the second term in the vorticity formula (37) is zero.

The stream function form of the vorticity equations
Theorem 2 Systems FDGI and SDGI-( 18) and (20), respectively-can be rewritten by applying the effective stress formulae (29)- (31), the stream function (35) and the vorticity formula (37) in the following forms: and respectively.
Proof The proof is provided in Appendix.On the right-hand side of Eq. ( 39), the term tr refers to the trace of the matrix product of the Hessians of the depth function and the stream function.
The difference of the models expressed by Eq. ( 32) can be rewritten as We illustrate the transition from the vorticity-velocity form to the stream function form on FDGI.Figures 2  and 3 are transition diagrams on which we can see the hidden depth-gradient effects revealed.The advection of vorticity V • ∇ζ -in which the depth-gradient effects are hidden-contributes to all the terms in the stream function form.All these terms are associated with Poisson brackets.The same terms are expressed and discussed in [4].Similarly, the diffusion of vorticity ζ contributes to all the terms on the right-hand side.
The stream function formalism provides a geometric point of view of the vorticity equation.All the bilinear operators can be related to intersection angles between various contour lines of scalar fields, as seen in Fig. 4. In Fig. 4 The geometric interpretation of the bilinear operators of FDGI the following sections, we restrict our attention to a very simple flow configuration, where all of these contour lines are coincide.We apply FDGI and SDGI to this idealized flow configuration and solve the corresponding simplified equations.We illustrate the differences between FDGI and SDGI by those simple but realistic solutions.

The concept of uniform, isobathimetric, straight flows and reduction in the stream function equations
We call the flow uniform if the magnitude of the depth-integrated velocity (33) |q| = |∇ψ| is constant along a streamline, i.e., the contours of |∇ψ| coincide with the streamlines: The flow is isobathimetric if the streamlines coincide with the isobath lines: The flow is straight, if the streamline curvature κ is zero: Consequently, all the terms related to the advection and the Coriolis force in FDGI and SDGI are zero, and all the contour lines are drawn in Fig. 4 coincide (ϑ = 0 for all cases).It follows from that the scalar products have the form In this flow configuration, there is no vorticity transport in the streamwise direction.On the other hand, the lateral vorticity transport is maximal and it occurred only by the turbulent diffusion including the depth gradient interactions.
We prescribe uniform, isobathimetric, straight flow in the following way: the streamlines and isobaths are both parallel to the y coordinate axis.The left panel of Fig. 5 provides a plan view of the streamlines and the isobaths.Based on this flow configuration and the definition of the stream function (34), the following ansatz is introduced for the stream function and depth function: Ansatz (52) transforms FDGI and SDGI [Eqs.(38) and (39)] into the ordinary differential equations , respectively, where = d/dx.The unknown function is the depth-integrated velocity profile q(x), and the depth function h(x) (cross-sectional profile) plays the role of a variable coefficient.A cross-sectional profile and a velocity profile can be seen in the right panel of Fig. 5.We turn the stream function equations into dimensionless form.Let L denote the characteristic width of the flow domain, where x * is the dimensionless coordinate.The depth function and the depth-integrated velocity are given as and respectively, where H is the characteristic depth and V is the characteristic velocity of the flow.The characteristic velocity can be the mean velocity, averaged to the cross-sectional area, where Q denote the volume flow rate.The stress fields have the form where τ w and τ b are characteristic values of the stresses.
The dimensionless form of FDGI and SDGI-by dropping the star notation-is respectively, where the dimensionless numbers are introduced.
The difference of the models expressed by Eq. ( 41) can be rewritten as The net interaction between the depth gradients and the turbulent diffusion results in lateral vorticity transport.That is linear in the depth-integrated velocity and its first and second derivatives (the first and the second velocity gradients), and nonlinear in the cross-sectional profile and depth gradients.
In the following two sections, we solve Eq. (60) for FDGI and (61) for SDGI for two simple flow cases: a gravity-driven channel flow and a wind-driven shoreline flow.

Solution for open channel flow with simple geometry
A steady state uniform flow in a straight channel with congruent cross sections in the streamwise (y) direction is considered.Our goal is to determine and compare the depth-integrated velocity profiles q(x) for some given cross-sectional profiles h(x) by solving FDGI and SDGI.The uniform, isobathimetric, straight flow concept can be applied, and FDGI and SDGI describe a lateral vorticity balance between the turbulent diffusion and the bed shear stresses of the channels.
First, we consider a rectangular channel (see left panel, Fig. 6) where there is no depth gradient, and the equations therefore simplify to the same form and give the same solutions.Then, we consider channels with depth variations (see middle and right panels, Fig. 6) expressed by prescribed cross-sectional profiles.
The bed stress field, τ b , is assumed to be a linear function of the depth-integrated velocity q [28,30], as where C b is a linear friction coefficient with dimension [s −1 ].Although the bed stresses are most likely related to the depth-averaged velocities V, Eq. (64) also applies [26].Here, we apply the formula (64) because the channels have no significant depth variations and they are bounded by vertical walls.
Fig. 6 The cross section profiles of the channels The relation between the characteristic units used for the dimensionless form of τ b and q-formulae (57) and (59)-is given as from which the second dimensionless number in (62) simplifies to The dimensionless bed stress curl term has the form where formulae (64), ( 33), ( 14), ( 37) and (52) were used, respectively.We do not consider wind stresses for the channel flow problem, i.e., Applying Eqs. ( 67) and (68) to Eq. (60) for FDGI and Eq.(61) for SDGI, and rearranging leads to and which are third-order linear homogeneous ordinary differential equations for the velocity profile q(x).Equation (70) for SDGI has a compact structure with respect to the vorticity: where ζ = q h − h q h 2 .Because of this structure, SDGI is integrable for arbitrary cross-sectional profiles h(x).The solution is where c 1 − c 3 are integration constants.
For the constant depth function (the constant is chosen to the unity), and both FDGI and SDGI simplify to The solution is Fig. 7 Cross section h(x) of the rectangular channel and velocity profiles q(x) with q 1 = −1 velocity gradient at x = 1 for different κ dimensionless numbers A rectangular channel is described with the unity depth function as We prescribe the following conditions at the middle of the channel and its right boundary: where q 1 is the velocity gradient at the right endpoint (x = 1).Applying the foregoing conditions, solution (75) becomes Figure 7 shows several velocity profiles obtained for different κ dimensionless numbers.The velocity gradient at x = 1 is chosen to be q 1 = −1.The solution curves reflected about the y axis-a similar boundary value problem can be solved on −1 ≤ x ≤ 0. Larger κ leads to smaller flow velocity, which is as expected because larger κ implies a wider channel or higher bed friction coefficient [see Eq. ( 62)] from which energy dissipation and smaller velocities follow.
In order to include depth variations in the channel cross section, we apply the following Gaussian function for the depth function, where δ is a shape parameter.The general solution of FDGI and SDGI with the depth function (79) is and respectively, where c 1 − c 6 are integration constants, M(•, •, •) and U(•, •, •) are Kummer's functions [2].The channel is defined as Fig. 8 The cross section of the channel and the velocity profiles with q 1 = −1 velocity gradient at x = 1 for different κ dimensionless numbers.Solutions FDGI and SDGI-( 80) and ( 81), respectively-are denoted by dashed and solid curves, respectively.On the left panel, the channel has smaller depth-gradients (δ = 0.1 in depth function ( 82)) than on the right panel, where δ = 0.5 The boundary conditions are the ones for the rectangular channel: Solution (81) of SDGI is denoted by solid curves in Fig. 8. Solution (80) of FDGI is represented by the dashed curves.On the left panel, there are smaller depth gradients, while on the right panel the depth gradients are higher, and the shape parameter of the depth function ( 82) is δ = 0.1 and δ = 0.5, respectively.In Fig. 8, it appears that the simplified diffusion effects-SDGI-overestimate the velocities for the open channel flow compared to the FDGI solution as we can see in the dashed and solid curves-FDGI and SDGI solutions, respectively-in Fig. 8.Given that the difference of the turbulent diffusion expressions lies in the absence (SDGI) or presence (FDGI) of the depth-gradient interactions with the turbulent stresses, the difference between the solutions increases with the depth-gradients (see right panels of Fig. 8).The difference between FDGI and SDGI lessens with the increase in the dimensionless number κ, i.e., with a wider channel or higher bed friction.
This result agrees with the physical expectations, namely the difference between solutions is higher when turbulent diffusion has greater effect on vorticity balance.For κ = 1, there is about 20% difference between FDGI and SDGI solutions when δ = 0.5 (right panel of Fig. 8).

Solution over sloping shoreline with wind forcing
We assume an idealized sloping shoreline with congruent cross sections in the streamwise (y) direction representing a coastline section of a lake or sea.The cross section can be seen in Fig. 9. Our goal is to compare the velocity profiles, q(x), given as solutions of FDGI [Eq.( 60)] and SDGI [Eq.( 61)] for a wind-driven recirculating flow.The flow is assumed to be isobathimetric along the straight shoreline.This flow represents the half-section of a large-scale eddy, from the shoreline to the middle of the eddy.Such large-scale circulating flows can be found in [17,18,26,30,31].
Let the bed stress field has the following linear dependence on depth-averaged velocity: where c b is a linear friction coefficient with dimension [ms −1 ].Usually, a quadratic formula is used to relate bed stress to depth-averaged velocity in numerical models.Even so, the linear relationship expressed by Eq. ( 84) has also been utilized in certain studies, see [28,30] for examples.The relation between the characteristic units used for the dimensionless form of τ b and V is given as from which the second dimensionless number in (62) reads as Applying formulae ( 84), ( 14), ( 37), ( 33) and ( 52), the dimensionless bed stress curl term reads as The wind stress τ w is prescribed to be constant and parallel to the shoreline: The wind stress curl term has the form The dimensionless number related to the wind stress term in Eq. ( 62) may be expressed as Substituting Eqs. ( 87) and (89) into Eqs.(60) for FDGI and (61) for SDGI yields and respectively.We obtain third-order linear nonhomogeneous ordinary differential equations for the velocity profile q(x).The inhomogeneity on the right-hand side corresponds to the wind forcing.The sloping shoreline is described by the following linear depth function: where δ is the slope.Applying depth function (93), the general solution of Eqs.(91) for FDGI and (92) for SDGI is and We describe a circulating domain of water in between the shore and the open water zone.The center of the eddy is placed at x = 1.Solutions (94) and (95) give zero velocities at the shore, where x = 0, and naturally satisfy a "no-slip" condition.At the center of the eddy, we prescribe zero velocity with arbitrary velocity gradients: where q 1 is the velocity gradient and q 11 is the derivative of the velocity gradient.
The SDGI solutions (solid curves) underestimate the peak velocities compared to the FDGI solutions (dashed curves) as seen in Fig. 10.This tendency is opposite to the channel flow solutions and can be caused by the wind-forced, nonhomogeneous nature of the shoreline problem.The increasing difference between FDGI and SDGI with greater depth gradients (δ = 0.5) is also observable but less than at the channel flow problems, especially for the curves corresponding to κ = 0.1.However, the difference between the solutions is significant at this parameter value.For the solutions with zero velocity gradients at x = 1 (see middle panels of Fig. 10), the equal values of the dimensionless numbers result in approximately the same solutions.These values imply equal weights for the wind and bed shear stresses which are greater than the turbulent diffusion effects by one order of magnitude.Nonzero boundary velocity gradients (see upper panels of Fig. 10) result in adjacent solutions for κ = 1.

Validation with CFD models
We now solve the channel and shoreline flow problems with CFD models.Our goal is not only to compare the analytical solutions with the numerical ones but also to alter the parameters to give realistic values for the dimensionless numbers.With such values of κ and λ, we can see the difference between our models of different turbulent diffusion expressions for realistic flow situations.
We chose the Mike21 software by DHI [8] to solve numerically the dynamic shallow water equations.Previously, Mike21 has been successfully verified and extensively applied to many riverine and coastal environments [18,20,26,29].Therefore, we give only a brief description of the solver here, and for further details, we refer to the model documentation [8].The numerical solution is obtained by an element-centered, Godunovtype finite-volume scheme [8,10] with second-order accurate time integration.The maximum time step was determined by the Courant-Friedrich-Levy (CFL) stability condition ensuring that CFL is not larger than 0.8 [1].The spatial domain was discretized with an unstructured triangular mesh that could capture complex bathymetry and curved shorelines.The spatial resolution was determined by a mesh refinement test to obtain mesh-independent results.
The numerical model solves the incompressible Reynolds-averaged and depth-integrated shallow water equations-consists of two momentum equations and the continuity equation, respectively: where g, ρ and ν e are the gravitational constant, the density of the fluid, and the eddy viscosity coefficient, respectively.Compared to the earlier steady shallow flow Eqs. ( 8) and ( 9), the unsteady shallow flow Eqs. ( 8) and ( 9) are in depth-integrated form, and the key variables depend on both space and time.The water depth h and the free surface elevation η are also having time dependence, as the stress fields τ b and τ w .In Mike21, the bed stress field is given by the following quadratic function of the depth-integrated velocity components: where n is the Manning-Strickler roughness coefficient.The eddy viscosity is determined using the Smagorinsky scheme [1] with a coefficient of 0.28.Given that Smagorinsky turbulence model is similar to our SDGI model, we are therefore able to relate the numerical solutions to the SDGI solutions.
The flow is steady with a constant volume flow rate, which is the cross-sectional integral of the depthintegrated velocity profile.We use a compatibility condition between the numerical and analytical velocity profiles, which requires their integral (the total volume flow rate, Q) must be the same for comparison purposes.The volume flow rate is determined by the numerical model, and the integral condition results in an expression for the velocity gradient q 1 at the boundary point.

Open channel flow
For the rectangular channel, we apply the integral condition on the solution (78) as from which the velocity gradient at the boundary is Condition (101) serves as a compatibility condition between the analytical velocity profile (78) of the flow in the rectangular channel and the numerical velocity profile for the same channel.The volume flow rate used by the numerical simulation is Q = 1 m 3 s −1 , which corresponds to a bed slope of 0.0007.The width of 3 .On reaching steady state, the cross-sectional velocity profile is determined and then, non-dimensionalized according to the mean velocity.We fit the numerical and analytical model by the dimensionless number κ and the dimensionless velocity gradient q 1 .For the dimensionless number κ, we use its general definition (62), in which the characteristic value of the bed stress, τ b appears.In this way, we can use characteristic bed stress directly from the numerical model.The bed stress and the eddy viscosity are varying-cross section wiselythrough the computational grid, and we use their mean values for calculating κ.We obtain κ = 3.23.The integral condition (101) results in q 1 = −3.04(dimensionless).The solutions can be seen in Fig. 11.
The same procedure has been used to validate the FDGI solution (80) and the SDGI solution (81) corresponding to the Gaussian cross-sectional profile (82).The compatibility condition with the CFD model, has been applied to the SDGI solution.The dimensionless number has the value κ = 4.08, and the (dimensionless) velocity gradient at the endpoint is q 1 = −2.15.The shape parameter for the Gaussian cross-sectional profile (82) is δ = 0.5.It can be seen in Fig. 12 that there is a good correspondence between the numerical solution (dots) and the SDGI solution (81).The FDGI solution (80) is represented by the dashed curve.
We can say that both the numerical and analytical models with the simplified Boussinesq approximation overestimate the maximum depth-integrated velocities in the middle of the channel.The complete application of the turbulent diffusion formulas, including interaction effects with the depth gradients, results in such solutions which predict lower velocity maxima.

Wind-driven flow
In the case of the wind-driven shoreline flow problem, we compare our FDGI solution (94) and SDGI solution (95) with simulations by a CFD model.The numerical model domain is rectangular horizontally with a width of 2 [m] and a length of 100 [m].The transverse slope of the bed is 0.5, and the Manning-Strickler roughness coefficient is n = 0.02 sm − 1 3 .The wind forcing is homogeneous in space, being parallel with the longer axis, and having a speed of 2 ms −1 (which is converted to wind stress using the empirical formula proposed by Wu [8]).This setup results in a simple circulating flow pattern consisting of a single eddy.The flow is parallel with the shoreline in the middle section of the longer axis.On the shallow side, the flow is oriented in the same direction as the wind; on the deeper side, the flow is oriented in the opposite direction.The center of Fig. 12 The cross section of the Gaussian channel and the velocity profile of the numerical solution (dots), SDGI solution (81) (solid curve), and FDGI solution (80) (dashed curve) Fig. 13 The cross section of the shoreline and the velocity profile of the numerical solution (dots), SDGI solution (95) (solid curve), and FDGI solution (94) (dashed curve) the eddy is located at 1.15 [m] from the shoreline.We calculate the velocity profiles at this region and compare them with the analytical models in Fig. 10.
The numerical data correspond to κ = 0.30 and λ = 16.44.By applying the compatibility condition to the SDGI solution, the following linear formula is obtained relating the velocity gradients: A possible realization of the solution is seen in Fig. 13, where the boundary velocity gradients have been chosen from the function q 1 (q 11 ) above to reach an optimal correspondence between the shape of the numerical and the SDGI profiles (dots and solid curve, respectively).For the wind-driven shoreline problem, both SDGI and FDGI solutions are not as close as in the open channel flow case.According to the numerical simulation, the velocity profile along the cross section is nearly symmetrical, while our solutions are skewed, resulting in larger velocities at the shallower side.The skewness of the SDGI solution is lower, and thus, it is closer to the numerical solution, as expected.The difference between the numerical model and analytical solutions arises from simplifications in the analytical solutions; namely, application of a linear formula for bottom stress and constant eddy viscosity.In the numerical model, the bottom shear stress and the eddy viscosity are both significantly higher in the shallow zone compared to the deeper one.As a result of these two counterbalancing effects, the velocity can increase by less intensity crosswise compared to the analytical solutions.Nevertheless, the SDGI solution gives a sufficiently good approximation for our purposes, namely, to reveal the impact of the depth-gradient related interactions.

Conclusions
The vorticity balance is commonly used to analyze the various mechanisms responsible for large-scale circulation patterns.For this purpose, the vorticity-velocity-depth-function form of the shallow water equation is used; however, this formulation cannot consider all the impacts of the bathymetry, more precisely, the depth-gradient related interactions.
Two important advances are presented in this paper.The first is the deduction of the stream function-depth function form [Eq. ( 38)] of the vorticity equations beside the more known vorticity-velocity-depth-function form.We showed that the stream function-depth function formalism is able to reveal all the depth-gradient effects that are hidden in the case of the vorticity-velocity-depth-function formalism.To our knowledge, these depth gradient-related terms have been presented for the first time.A further advantage of the stream function formalism is that it provides a clear geometrical point of view for the vorticity equation based on the relative position of streamlines, isobath lines, and further isolines.We derived and compared two mathematical models applying the Boussinesq approximation.The SDGI model represents the classical vorticity formalism by neglecting some depth gradient effects.In contrast, the FDGI model is considered a full model by incorporating all the depth gradient interaction terms.
The second advance concerns the application of the two models to gravity-driven open channel flow and wind-driven shoreline flow.The channel flow involved a homogeneous solution of the vorticity equations with a single dimensionless number κ, whereas the shoreline flow involved an inhomogeneous solution with two dimensionless numbers κ and λ.Solutions obtained for these realistic flow problems using the two models were found to differ significantly in terms of the velocity profiles.Interactions related to the depth gradient could not necessarily be neglected without incurring significant error.
To validate our solutions, the two flow situations were also simulated using CFD.Similar results were obtained between the analytical models and the CFD simulations, even though the analytical models were based on constant eddy viscosity and linear bed-friction formulae assumptions, while the CFD models used variable eddy viscosity coefficients and quadratic bed-friction formulae.The CFD simulations provided a closer match to SDGI solutions because the depth-gradient effects were neglected in the both models.
The newly introduced dimensionless numbers in the present study have practical implications.Our sensitivity analyses revealed that the depth-gradient effects could vanish depending on the flow and geometry characteristics.On calculating, the bottom stress and eddy viscosity coefficient, the dimensionless numbers κ and λ can be used to assess whether the full depth-gradient interactions in the Boussinesq approximation are required in order to gain a more accurate solution than the shallow water equations can provide.by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Availability of data and materials
The numerical datasets generated during the current study are available from the corresponding author on request.

Conflict of interest
The authors have no competing interests to declare that are relevant to the content of this article.
Ethical approval This declaration is not applicable.

Proof of Theorem 1
The resolved form of the curl of the effective stresses (21) in detail: We substitute the Boussinesq eddy viscosity concept (27) into Eq.(105); and after that we insert the depthaveraged vorticity (14 The result is known as the diffusion of the depth-averaged vorticity. With the Boussinesq eddy viscosity concept (27), the interaction terms (106) and (107) turn into and respectively.We introduce the velocity gradient tensor (28), ∇V, and the following quantities:

Proof of Theorem 2
We apply the stream function (35) on the resolved formula for the advection and stretching of the depth-averaged vorticity (15), from which The resolution of the diffusion of vorticity (29)- -by the depth-averaged vorticity expression (37) yields where Now, we resolve the interaction terms in (30)- -by applying the stream function (35) and the depth-averaged vorticity expression (37).The first term has the form The second interaction term can be rewritten as where with which The third interaction term transforms as The resolution of the interaction terms (31)-

Fig. 1
Fig. 1 Definition sketch of the water depth h(x, y)

Fig. 2 Fig. 3
Fig. 2 The transition diagram between the advection related terms of the vorticity-velocity-depth function and the stream functiondepth function formalisms.The advection of vorticity (only implicit depth-gradient effects) is placed on the left, and the vortex stretching (both explicit and implicit depth-gradient effects) is placed on the right

Fig. 5
Fig. 5 On the left panel, the plan view of the streamlines and isobaths (straight black and grey lines) can be seen.The right panel shows the cross section of the water body represented by the depth function h(x), and the distribution of the depth-integrated velocities (the velocity profile)

Fig. 9
Fig.9The cross section profile of the shoreline

Fig. 11
Fig.11The cross section of the rectangular channel and the velocity profile of the numerical solution (dots) and solution (78) (solid curve) Open access funding provided by Budapest University of Technology and Economics.The research reported in this paper is part of Project No. BME-NVA-02, implemented with the support provided by the Ministry of Innovation and Technology of Hungary from the National Research, Development and Innovation Fund, financed under the TKP2021 funding scheme.This work has been supported by the Hungarian National Research, Development and Innovation Fund under contract NKFI K 137726.The research presented in the article was carried out within the framework of the Széchenyi Plan Plus program with the support of the RRF 2.3.1 21 2022 00008 project.Furthermore, Project No. TKP-6-6/PALY-2021 has been implemented with the support provided by the Ministry of Culture and Innovation of Hungary from the National Research, Development and Innovation Fund, financed under the TKP2021-NVA funding scheme.