Hydraulic Reference Manual Iber

Twodimensionalmodellingoffreesurfaceshallow
waterflow
HydraulicReferenceManual
Iberv1.0
19.07.2014
Two-dimensional modelling of free surface shallow water flow
REFERENCE MANUAL
1. PRESENTATION ...................................................................................... 5
2. HYDRODYNAMIC MODULE ...................................................................... 7
2.1. Introduction ........................................................................................7
2.2. Hydrodynamic Equations .......................................................................7
2.3. Bed Friction .........................................................................................8
2.4. Superficial Friction due to Wind ..............................................................9
2.5. Effective Stresses.................................................................................9
2.6. Hydrodynamic boundary Conditions ...................................................... 11
2.6.1. Closed Boundaries ........................................................................ 11
2.6.2. Open Boundaries .......................................................................... 13
2.7. Internal boundary conditions ............................................................... 14
2.7.1. Gates .......................................................................................... 15
2.7.2. Weirs .......................................................................................... 15
2.7.3. Gate-Weir Combination.................................................................. 16
2.7.4. Local Losses ................................................................................. 16
2.8. Infiltration ........................................................................................ 17
2.8.1. Green-Ampt ................................................................................. 17
2.8.2. Horton ........................................................................................ 18
2.8.3. Linear ......................................................................................... 19
2.9. Initial Abstraction............................................................................... 19
3. TURBULENCE MODULE .......................................................................... 21
3.1. Introduction ...................................................................................... 21
3.2. Scales of Turbulence in Shallow Waters ................................................. 22
3.3. Constant Turbulent Viscosity................................................................ 23

Iber. Turbulent Flow in Shallow Waters
3.4. Parabolic Profile of turbulent Viscosity ................................................... 23
3.5. Mixing Length Model ........................................................................... 24
3.6. The Rastogi and Rodi K-E model (Rastogi & Rodi, 1978) .......................... 24
3.7. Dimensional Analysis of the turbulent terms in the shallow water equations 25
4. NON-STATIONARY SEDIMENT TRANSPORT MODULE ............................ 28
4.1. Sediment Conservation Equations......................................................... 28
4.2. Bed Transport.................................................................................... 29
4.2.1. Stress Partition ............................................................................. 29
4.2.2. Bedload Discharge ........................................................................ 29
4.2.3. Bed slope correction ...................................................................... 31
4.2.4. Avalanche model .......................................................................... 32
4.2.5. Considerations on a non-erodible layer ............................................ 32
4.3. 2D Turbulent Suspended Load Module................................................... 32
4.3.1. Turbulent Suspended Load Equation ................................................ 32
4.3.2. Entrainment/Deposition (E-D) term Calculation ................................. 33
4.3.3. Sedimentation Velocity .................................................................. 35
5. NUMERICAL SCHEMES .......................................................................... 37
5.1. Calculation Mesh ................................................................................ 37
5.2. Discretization in finite volumes for the 2D-SWE Equations ........................ 38
5.2.1. Discretization of the convective flux terms ........................................ 39
5.2.2. Discretization of the bed slope source term ...................................... 42
5.3. Discretization of the Transport Equations in the K-E Turbulence Model, and in
the Suspended load Transport Model ........................................................... 43
5.3.1. Depth-averaged Transport Equation ................................................ 43
5.4. Discretization of the Exner Conservation of Sediment Equation ................. 45
5.4.1. Non-erodible bed layer .................................................................. 46
5.5. Managing Wet-Dry Fronts .................................................................... 47
Page 2 of 56
6. REFERENCES ......................................................................................... 50
7. NOMENCLATURE ................................................................................... 52
Page 3 of 56
IBER
TurbulentFlowinShallowWaters
PRESENTATION
Page 4 of 56
1. PRESENTATION
Iber is a numerical model for simulating turbulent free surface unsteady flow and environmental
processes in river hydraulics. The ranges of application of Iber cover river hydrodynamics, dambreak
simulation,floodzonesevaluation,sedimenttransportcalculationandwaveflowinestuaries.
AtthismomentIberhas3maincomputationalmodules:ahydrodynamicmodule,aturbulencemodule
andasedimenttransportmodule.Allofthemworkinafinitevolumenonstructuredmeshmadeupof
triangleorquadrilateralelements.Inthehydrodynamicmodule,whichconstitutesthebaseofIber,the
depthaveragedtwodimensionalshallowwaterequationsaresolved(2DSaintVenantEquations).The
turbulent module allows including turbulent stressess in the hydrodynamic calculation, in this way
allowingitsusefordifferentturbulentmodelsforshallowwatersatdifferentdegreesofcomplexity.The
mostrecentversionofIberalsoincludesaparabolicmodel,amixinglengthmodelandakmodel.The
sedimenttransportmodulesolvesthebedloadandtheturbulentsuspendedloadtransportequations,
basedonthesedimentmassbalanceevolutionatthebed.
This manual describes in detail the equations and models used in the different Iber computational
modules,includingtheusednumericalschemes.
InordertocitetheIbermodelintechnicaldocumentsandpublicationsthefollowingreferenceshould
beused:
Blad,E.,Cea,L.,Corestein,G.,Escolano,E.,Puertas,J.,VzquezCendn,M.E.,Dolz,J.,Coll,A.(2014).
"Iber: herramienta de simulacin numrica del flujo en ros". Revista Internacional de Mtodos
NumricosparaClculoyDiseoenIngeniera,Vol.30(1)pp.110
Page 5 of 56
IBER
HYDRODYNAMICMODULE
Page 6 of 56
2. HYDRODYNAMICMODULE
2.1. Introduction
ThehydrodynamicmodulesolvesthedepthaveragedShallowWaterEquations(2DSWE),alsoknownas
thetwodimensionalSt.VenantEquations.Theseequationsconsiderahydrostaticpressuredistribution
and a relatively uniform distribution of the depth velocity. The hydrostatic pressure hypothesis
reasonablycompliesinriverflowsandestuarywavecurrents.Likewise,thehypothesisofuniformdepth
velocitydistributionusuallycompliesinriversandestuaries,eventhoughtheremaybezonesinwhich
thisrequirementisnotfulfilledduetolocalthreedimensionalflowsorsalinewedges.Inthesecasesitis
necessarytostudytheextensionofthesezonesanditspossiblerepercussioninthemodelresults.At
thismoment,thenumericalmodelsbasedonthe2DSWEarewidelyusedincoastalandriverdynamics
studies,floodzonesevaluationandsedimentandcontaminationtransportcalculations.
2.2. HydrodynamicEquations
The hydrodynamic module solves the conservation of mass and momentum equations in the two
horizontaldirections.
h hU x hU y
MS
t x y
hU x hU 2x hU x U y h exx h xy
e
Zs s,x b, x g h 2
gh 2 sin U y MX
t x y x 2 x x y
hU y hU x U y hU 2y Zs s,y b, y g h 2 h exy h eyy
gh 2 sin U x MY
t x y y 2 y x y

Wherehisdepth,Ux,Uyarethedepthaveragedhorizontalvelocities,gisthegravityacceleration,Zis
thefreelayerelevation,sisthefreesurfacefrictionduetowindinducedfriction,bisthebedfriction,
isthewaterdensity,istheearthsrotationangularvelocity,isthelatitudeofthestudiedpoint,exx,
exy,eyyaretheeffectivehorizontaltangentialstressess,andMs,Mx,Myarerespectivelythetermsof:
masssource/drainandmomentum,whichareusedtomodelprecipitation,infiltrationanddrainage.
Thefollowingtermsareincludedinthehydrodynamicequations:
Hydrostaticpressure
BedSlope
Page 7 of 56
ViscousandTurbulenttangentialstresses
BedFriction
Superficialfrictionduetowind
Rainfall(Precipitation)
Infiltration
Likewise,wet/dryfronts,bothstationaryandnonstationary,thatmayappearintheareaofstudyare
modelled. These fronts are fundamental in the fluvial and coastal flood zones modelling (rivers and
estuaries).Inthiswaythepossibilitytoevaluatetheextensionoffloodzonesinriversisintroduced,asis
themovementofwavefrontsinestuariesandcoastalzones.
2.3. BedFriction
Thebedwieldsafrictionforceoverthefluidthatisequivalenttothefrictionwithawall,eventhough,in
generalthebedsroughnesshasagreatervalueinthecaseofriversandestuaries.
Bed friction has a double effect on the flow equations. On one side it produces a friction force that
opposes the mean velocity, and on the other side, it produces turbulence. Both effects can be
characterized by the friction velocity uf, which is a form of expressing the beds tangential stress in
velocityunits:
b
uf

Wherebisthebedfrictionforcemoduleandisthewaterdensity.
In depth averaged models its not possible to calculate friction velocity by means of standard wall
functions, as it is done in wall type boundaries, because these equations are solved in a vertical
direction. This means that its necessary to relate friction velocity uf with depth averaged velocity by
meansofafrictioncoefficient.Thebedstresscanbeexpressedas:
2
b u f2 Cf U
WhereCfisthebedfrictioncoefficient.Therearedifferentexpressionstoapproximatethiscoefficient.
Mostofthemassumeuniformchannelflowinalogarithmicprofileofdepthvelocities.
Page 8 of 56
Unlike1Dmodels,in2Dmodelsthehydraulicradiusisnotdefinedasthewetsectionareaoverthewet
perimeter,sincein2Dmodelsitmakesnosensetodefineatransversalsection.Takingacolumnoffluid
oflengthxanddepthh,thehydraulicradiuswouldbecalculatedas:
A h x
Rh h
Pm x
Meaningthatin2Dmodelshydraulicradiusanddepthareequivalent.
The bed friction is evaluated using the Manning formula, which uses the Manning n coefficient as
parameter.TheManningformulausesthefollowingroughnesscoefficient:
n2
C f g 1/3
h
2.4. SuperficialFrictionduetoWind
Thefrictionforceduetowindoverthefreesurfacecanbecalculatedconsideringthewindvelocityat10
metersofheightandadraggingcoefficient,usingtheVanDornequation(VanDorn,1953).
s Cvd V10 2
Whereisthewaterdensity,V10isthewindvelocityat10metersofheightandCvdisthesurfacedrag
coefficient.BydefaultadragcoefficientofCvd=2.5106isconsidered.
2.5. EffectiveStresses
Thehorizontaleffectivestressesthatappearinthehydrodynamicequationsincludetheeffectsofthe
viscousstresses,oftheturbulentstressesandthedispersiontermsduetothenonhomogeneityofthe
depthvelocityprofile.
ije ijv u'i u' j D ij
v
Where ij aretheviscousstresses, u' i u' j aretheturbulentstresses(alsocalledtheReynoldsstresses),
andDijarethelateraldispersionterms:
1 Zs
U i u i U j u j dz
h Z b
D ij
Page 9 of 56
Thedispersiontermsarenotconsideredinthe2DSWEequations(rememberthehypothesisofuniform
depthvelocityprofile),becauseitisimpossibletocalculatetheminageneralwaywithadepthaveraged
model.Theirimportancewillbegreateraslessuniformisthedepthvelocityprofile.Atypicalsituation
in which these terms can gain importance is in channels with elbows or small curvature radii, like in
confluencechannels(Figure1)
Q1
Q
3
Q
2
Figure1.Secondaryflows(left)anddepthvelocityprofile(right).Maincausesofthedispersionterms.
Theviscousstressesarecalculatedfromthefluidskinematicviscosity, as:
U U j
ijv i
x x
j i
Ingeneral,theorderofmagnitudeoftheviscousstressesismuchlowerthantherestofthetermsthat
appearinthehydrodynamicequations,exceptnearthewallsandinlaminarflow.
Theturbulentstressesaremuchlargerinmagnitudethantheviscousstresses,especiallyinrecirculation
zones,whereturbulenceproductioniselevated.Inthecaseoftwodimensionalshallowequationsthe
turbulentstressesconstitutethreenewunknownstocalculate,addedtothedepthandvelocitiesUx,Uy
they mean a total of 6 unknowns. This is what is known as turbulence closure problem, because it is
necessarytosolve3equationswith6unknowns.Duetothis,itisnecessarytouseaturbulencemodel
that allows calculating these turbulent stresses. Most of the models calculate the turbulent diffusive
termsfromthefollowingexpression:
u'i u' j U i
t
x j x j x j
Where t istheturbulentviscositycalculatedintheturbulencemodule.Theproblemliesinthenon
existence of a universal turbulence model, that may allow calculating in a precise form the turbulent
stresses. Through time some models have been developed with different degrees of complexity. The
BoussinesqformulationisusedbyallturbulencemodulesinIber.
Page 10 of 56
2.6. HydrodynamicboundaryConditions
Inatwodimensionalproblemitisnecessarytodistinguishtwotypesofboundaries:openandclosed.
The closed boundaries, also called walltype boundaries, are waterresistant, meaning that they dont
allowwaterthroughthem.
2.6.1. ClosedBoundaries
Thepresenceofawalltypeboundarygeneratesalateralfrictionforceinthefluidinasimilarfashion
thanthefrictionexertedbythebedsroughness.Ataclosedboundarythefollowingconditionscanbe
used:
Slipcondition(nulltangentialstress)
WallfrictionCondition(wallfunctions)
Theslipconditionisequivalenttounderminingthefrictionstressgeneratedbythewalltypeboundaries
over the fluid. In hydraulic engineering, and especially in river engineering, the surface contact with
lateralboundariesismuchlowerthanthesurfacecontactwiththebed,duetotheseparationbetween
horizontal and vertical scales, so the friction forces in the wall boundaries can be undermined. In this
casetheslipconditionwouldhold.
In problems where the horizontal and vertical dimension are similar (narrow section channels) this
roughness force can have some importance in the flow development, even though in general the
influenceissmall.Ifyouwanttoconsiderthelateralfriction,awalltypefrictioncanbeintroduced;this
consists in imposing a tangential force in the opposite direction to the flow. In this case, Iber
distinguishes between smooth turbulent regime and rough turbulent regime in function of the walls
roughnessandtheflowvelocityclosetothewall.
Thewallfrictionvelocity(u*)isdefinedinfunctionofthewallfriction( w )as:
w
u*

Thetangentialvelocityatthewallcanbeexpressedasafunctionofthefrictionvelocity,theroughnesss
depthandthewalldistanceas:
u*
u Ln E y y
y u*

Page 11 of 56
Whereyistheperpendiculardistancetothewall,andEisaparameterwhichvaluedependsontheflow
characteristics. To calculate E, Iber considers smooth turbulent flow conditions, rough turbulent flow
conditionsandatransitionbetweensmoothandrough(Table1).
KSu* u*
Regimentype KS u Ln E y

SmoothTurbulent K S 5 E 9.0
30
RoughTurbulent 5 < KS 70 E=
KS
1
SmoothRoughTransition KS 70 E=
0.11 + 0.033 K S
Table1.WallFriction.
Theturbulenceisdefinedassmoothwhenthisrelationshipcomplies:
KSu*
KS 5

Where Ks is the walls roughness height, a measure of the walls roughness, with length units. In this
conditionthetangentialvelocityatthewallcanbeexpressedasafunctionoffrictionvelocityandthe
kinematicviscosity,as:
u* yu
u Ln 9.0 *

Roughturbulentregimeisappliedwhenthisrelationshipcomplies:
K Su*
K S 70

Inthiscondition,thetangentialvelocityatthewallcanbeexpressedasafunctionofthefrictionvelocity
andthebedsroughnesssheight,as
u* y
u Ln 30
KS
Inthetransitionbetweenregimes,thetangentialvelocityatthewallcanbeexpressedasafunctionof
frictionsvelocity,kinematicviscosityandroughnesssheight,as:
Page 12 of 56

u y
u * Ln
0.11 + 0.033 K
u*
S

2.6.2. OpenBoundaries
In open bondaries, different types of boundary conditions can be introduced. So that the 2DSWE
equationsarecorrectlypresentedfromamathematicalpointofview,thenumberofconditionstouse
intheopenboundarieswilldependiftheboundaryisaninletoranoutlet,itwillalsodependonthe
typeofregimeattheboundary(fast/slow).Inaninletboundaryitisnecessarytointroduce3boundary
conditionsiftheregimeissupercritical(oneforeachofthe3St.Venantequations),whileiftheregime
is subcritical, it is enough to present 2 conditions. In an outlet boundary, it is enough to introduce a
singleconditioniftheregimeissubcritical;meanwhile,itisnotnecessarytoestablishanyconditionif
the regime is supercritical. If the user introduces fewer conditions than those necessary, from a
mathematicalpointofviewtheequationswillbeundeterminedandnocorrectsolutionwillbefound.
Theconditionswhichcanbeenteredbytheuserareheight,velocitycomponents,oracombinationof
both.InIberdifferentoptionsareconsideredtointroduceboundaryconditions,asexplainedinTable2.
Itiscommoninriverhydraulicsfortheflowtobeinsubcriticalregimeinthecontoursofthemodeled
section. In this case, what is often done is to introduce free surface depth or water level in the
downstreamboundary.Intheupstreamboundarythemostcommonistoenteratotalinflow(incubic
meters per second cms) in the direction of the flow, which in general, lacking more precise data, is
assumedtobeperpendiculartotheinletboundary.Itisalsopossibletointroduceupstreamvelocities
(inm/s)orspecificflow(m2/s),thoughthisoccurslessoften.Inthecaseoftotalflowintheinlet,aunit
flowdistribution(m2/s)isentered,accordingtothefollowingexpression:
h 5/3
qn Q
h dy
5/3
Whereqnisthenormalspecificdischarge(m2/s)ateachpointoftheinlet,andQisthetotalinletflowat
mentionedboundary.Theintegralinthedenominatorisextendedalongthelengthoftheboundary.
Besides height, other possibilities can be considered at the outlet, such as weirs or rating curves. The
weirconditionestablishesarelationshipbetweenoutflowandheightateachpointofthecontour:
q = Cd (ZS ZW )1.5
Page 13 of 56
WhereCdis theweirsdischarge coefficient,Zsisthewaterlevelupstream oftheweirandZwisthe

waterlevelontopoftheweir.TheusermustintroduceCdandZw.
The rating curve contour condition establishes a general relationship between the outflow and the
water level en each contour point. This relationship is introduced by the user using a table defining
valuesofspecificdischargeandwaterlevel.
ThesetofconditionsimplementedinIberinopencontoursareshowninTable2.
Boundary Regime Boundaryconditions
Total discharge in normal direction to the

Subcritical/Critical
boundary
Discharge
Total discharge in normal direction to the
Supercritical
boundaryandaveragevelocity
Inlet Specific discharge in normal direction to

Subcritical/Critical
theboundary

a) Specificdischargeinnormaldirectionto
Unitdischarge theboundaryandwaterdepth
Supercritical
b) Specificdischargeinnormaldirectionto
the boundary and water surface
elevation
a)Waterdepth
b)Watersurfaceelevation
Subcritical
Outlet c)Weir(elevationanddischargecoefficient)
d)Ratingcurve
Supercritical/Critical Notnecessarytointroduceanycondition
Table2.ContourConditionsimplementedinopencontours.
2.7. Internalboundaryconditions
Theinternalconditionsareusedtomodelhydraulicstructureslikegates,weirsorbridgeswhichchange
theconditionsofthesystem.
Page 14 of 56
The internal boundary conditions implemented in Iber can be used to model the following flow
conditions:
Flowbelowagate
Flowoverafreefallingweir
WeirGatecombination
Localloss
2.7.1. Gates
Theequationusedisthegatedrainage,whichworksonfreeorsubmergedflow.Thedatatoinsertare
thedischargecoefficient,thegatesbedlevel,thegateopeninganditswidth.Bydefault,thedischarge
coefficientwillhaveavalueofCd=0.6.
ZU
ZD
h
ZB
Figure2.Schemeandequationsofthegateinternalcontourcondition.
(ZD ZB ) / (ZU ZB ) DischargeEquations
FreeFlow 0.000.67 Q = Cd B h 2g (Z U ZB )
Transition 0.670.80 Q = Cd B h 6g (ZU ZD )
SubmergedFlow 0.801.00 Q = Cd B h 2g (ZU ZD )
2.7.2. Weirs
Therectangularweirdischargeequationisused,whichcanworkforfreeflowandsubmergedflow.The
datatoinputaretheweirsupperlevel,thedischargecoefficientandtheweirslength.Bydefaultthe
dischargecoefficienthasavalueofCd=1.7.
Page 15 of 56
ZU
Zw ZD
ZB
Figure3.Schemeandequationsoftheweirinternalcontourcondition.
(ZD ZW ) / (ZU ZW ) DischargeEquation
FreeFlow <0.67 Q = Cd B (ZU ZW )1.5
SubmergedFlow >0.67 Q = 2.598 Cd B (ZD ZW ) (ZU ZW )0.5
2.7.3. GateWeirCombination
Inthiscasetheconditionincludesthetwopreviousones,sobothparametersmustbeused,thegates
andweirs.Thetotaldischargedflowistakenasthesumbelowthegateandtheflowabovetheweir.
ZU
Zw
ZD
Q Qcompuerta Qvertedero
h
ZB
Figure4.Schemeandequationsofthegateweirinternalcontourcondition.
2.7.4. LocalLosses
In this case, there is a local energy loss in the flow transference between two finite volumes of value
H= v2/2g. The St. Venant Equations are the mathematical expressions of the conservation of mass
andmomentum,awaytoconsiderthislossofenergyisoverthefrictionslope.Forthis,anadditional
termisaddedtothefrictionslope,Sf,throughtheedgeofthefinitevolume,thisaddedtermisequalto
H/V,whereVistheelementsvolume.Inthisway,theenergylossthroughthetwoelementswillbe
Page 16 of 56
equal to H+ SfL, where L is the distance between elements centre on both sides of the boundary
wherethelocallossisimplemented.
v2
H
2g
H
S 'f S f
V
Figure5.Schemeandequationsofthelocallossinternalboundarycondition.
2.8. Infiltration
Tosimulaterainfallprocesses,itcanbenecessarytoconsiderwaterinfiltrationinnonsaturatedterrain
for calculating runoff. Modeling infiltration of surface waters is especially important in rainfallrunoff
models.
Infiltrationisconsideredinthemodelusinganegativetermintheconservationofmassequation(loss
ofwater):
h hU x hU y
=-i
t x y
Whereiistherealinfiltrationrate,measuredastheminimumbetweenthepotentialinfiltrationratef
(soilsinfiltrationcapacityateveryinstant,whichdependsonsoilcharacteristicsandconditions),andthe
amountofreadilyavailablesurfacewaterforinfiltration.
h
i = min ( f , )
t
To calculate the potential infiltration 3 infiltration common models are used: The GreenAmpt model,
theHortonmodelandlinearmodel.
2.8.1. GreenAmpt
Theinfiltrationrate,expressedinm/s,iscalculatedineachcellusingtheGreenAmptformula(Chow,
Maidment, & Mays, 1988), in it; it is assumed that there exists a saturated front that separates the
saturated soil region, immediately below the terrain, and the nonsaturated region, where there is
suction.
Page 17 of 56
Asthewaterinfiltrates,thesaturatedfrontdescendsandthedepthofthesaturatedregionLbecomes
bigger.Theinfiltrationratefiscalculatedas:
h t F
f k s 1 F f dt L L0 + i
L0 F 0
Whereksthesaturatedpermeabilityofthesoil,histhedepth,isthesuctioninthenonsaturatedsoil
region,thehumiditycontentchangeofthesoilasthesaturatedfrontincreases,ithesoilsinitial
humiditycontent, thesoilstotalporosity,andListhedepthofthesaturatedregion.Theinfiltration
rateisequaltothepotentialinfiltrationrateaslongasthereisenoughsurfacewatertoinfiltrate.
Theparameterstointroducebytheuserforthismoduleare:
Thesoilssaturatedpermeability(ks)
Thenonsaturatedregionssuction()
Soilseffectiveporosity(drainable)(e)
Soilsinitialeffectivesaturation(Se),definedas:
i - r
Se =
e
Where r isthesoilsretentioncapacity(nondrainablehumidity)and i thesoilsinitialhumidity.The
soilsporosity isequaltothedrainableporosityplusthesoilsretentioncapacity( = e r ).
Starting with the effective porosity and the soils initial effective saturation, the change in humidity
contentcanbecalculatedasthesaturationfrontadvanceswiththeformula:
= - i = - r - e Se e 1 Se
AlltheparametersoftheGreenAmptequationcanbeintroducedasspatialvariables(differentforeach
elementofthemesh).
2.8.2. Horton
IntheHortonmodel,thepotentialinfiltrationrateiscalculatedas:
f f c f 0 - f c exp k t
Page 18 of 56
Where t is time, starting with the beginning of the precipitation. The user must introduce as model
parameters an initial infiltration rate ( ), the infiltration rate at infinite time ( ) and the constant k,
whichdefinesthetemporalvariationofthepotentialinfiltrationrate.AlltheparametersoftheHorton
infiltrationequationarespatiallyvariable(calculatedforeachelementofthemesh).
2.8.3. Linear
The linear model considers an initial abstraction P0 (volume per unit of area), and afterwards some
continuousconstantlosses(volumeperunitofareaandunitoftime).Thevalueoftheinitialabstraction
andthecontinuouslossescanvaryperelement.

Figure6.Temporalevolutionoftheinfiltrationrateaccordingtothelinearmodel
2.9. InitialAbstraction
IftheGreenAmptorLinearinfiltrationmodelsareusedtocalculatetheinfiltrationlosses,thereisalso
the possibility to consider an initial abstraction. This initial abstraction can represent other processes
likethecanopyretention(becauseofvegetation)orlanddepressions,ortheinfiltrationcapacityofdry
soilswithanelevatedporosity.
Theinitialabstractionisdefinedasavolumeperunitofarea,andasthathaslengthunit.Thisvalueis
subtractedfromthewaterthatreachesthesoil,beingfromrainfallorfromsurfacerunoff.Thatwayit
canbeusedeitherinareaswithorwithoutrainfall.
Page 19 of 56
IBER
TURBULENCEMODULE
Page 20 of 56
3. TURBULENCEMODULE
3.1. Introduction
Agreatnumberofstudiesinhydraulicengineeringimpliestheanalysisofopenchannelflows,manyof
whichcanbeconsideredshallowflows,understandingthisaslowdepthwithsmallverticalhorizontal
relationship.Practicallyallopenchannelflowsareturbulent.Whilelookingatariver,smallswirlsappear
and disappear with an almost chaotic movement, showing the complexity of the turbulent motion.
Theseswirlsarethemaincauseofmixingprocesses,sotheyplayanimportantroleinthediffusionof
solublesubstancesandsuspendedsolids,amongothers.
Despite the fact that all flows in river engineering are turbulent, in many cases the turbulence is not
largeenoughtohaveanotoriousinfluenceinthemeanvelocity.Suchisthecaseofrivers,estuariesand
in general in coastal zones with a smooth enough geometry so that no recirculation zones appear.
Nevertheless,eveninthistypeofsituationsitsimportanttoconsiderturbulencewhilemodellingthese
watersystems,sinceitplaysanimportantroleintheprocessesoftransportandmixingofcontaminants
andsediments.Thediffusionofheat,ofasolute,orsuspendedsedimentisproducedbasicallybecause
ofturbulence,exceptinlaminarflow,whichgenerallydoesnotoccurinhydraulicengineering,evenless
inriversandestuaries.Theturbulentdiffusioncoefficientisseveraltimeshigherinmagnitudethanthe
molecular diffusion coefficient. Because of this, it is necessary to evaluate previously the turbulent
kineticenergytocalculatethediffusiveflow.
OneofthemaincharacteristicsofIberistheinclusionofdiverseRANStypeturbulencemodels,which
aresolvedintheturbulencemodule.Thefollowingturbulencemodelsareincludedforshallowwaters,
enumeratedbyorderofcomplexity:
ConstantTurbulentViscosity
ParabolicModel
MixingLengthModel
k(kaepsilon)modelfromRastogiRodi(RastogiyRodi,1978)
Usingmodelswithdifferentlevelsofcomplexityallowsselectingthemostadequateforeachcasestudy,
considering the complexity of the flow and the model. In general the mixing length model provides
satisfactoryresultsinriversandestuaries,thiswaydiscouragingtheuseofotherturbulencemodelsin
these cases. In hydraulic structures such as open flow channels with marked elbows and recirculation
Page 21 of 56
zones,itmaybenecessarytouseatleastonemixinglengthmodelandonekmodeltotesttheresults.
The choice of the turbulence model that best fits each case is based on user experience, considering
that the more complex the model is, the longer it takes to calculate and more complex the equation
solutionis.
The objective of the turbulence models is to calculate the Reynolds stresses. Models based on the
Boussinesqhypothesis(theonesenumeratedbefore,usedinIber),theReynoldsstressesareevaluated
fromthefollowingexpression:
U i U j 2
uiu j t k ij
x
j x i 3
Theturbulencemodelcalculatestheturbulentviscositytobeusedinthepreviousexpression.
3.2. ScalesofTurbulenceinShallowWaters
Oneofthemaincharacteristicsofshallowflowsisthesplitupbetweenhorizontalandverticalscales,
this happens because the vertical extension of the fluid (limited by depth) is much smaller than the
horizontalextension.Thisscalesplitupisapplicableinthespatialdimension,velocities,andofcourse
the turbulence. Speaking of turbulence, its main effect is the fragmentation between tridimensional
(swirls)andbidimensionalturbulentstructures.The3Dturbulencespatialscaleislimitedbyitsdepth,
and thus its structures are much smaller than those associated to the 2D turbulence, which are only
limited by the horizontal scale. The 3D turbulence is originated mainly by the bed friction, while the
originsofthe2Dturbulencearethehorizontalvelocitygradients
Itsimportantfortheturbulencemodeltoincludenotonlytheeffectsofthe3Dturbulence,formedby
thebedfriction,butalsothe2Dturbulence,formedbythehorizontalvelocitygradients.Intheshallow
water models, the bidimensional characteristics of the flow are implicitly considered in the transport
equations when considering a homogenous depth velocity profile, while the tridimensional causes of
turbulenceareusuallyincludedusinganothertermwhichdependsonthebedtangentialstress.Inthe
same way, even when considering 3DSWE, the turbulence model should consider turbulence
anisotropyinthehorizontalandverticaldirections.
Iber uses the following turbulence models. All of them are shallow water depth averaged turbulence
models.
Page 22 of 56
3.3. ConstantTurbulentViscosity
Theorderofmagnitudeoftheturbulentviscositycanapproximateddependingontheconsideredflux.
Several published articles propose these turbulent viscosity approximated values in function of the
consideredflux.Thisfocusisverysimple,andcannotbeconsideredasanadequate,norrealisticinany
case,turbulencemodel,sinceitdoesnotconsiderspatialchangesofturbulentviscosity.Itsimportant
to highlight the importance of spatially variable turbulence viscosity to the mean velocities.
Furthermore,theexistingtablesonlyprovideapproximatedvalues.Becauseofthesereasons,theuseof
thismethodisdiscouraged,sinceitcanleadtoexcessivelylargeerrors,generallycomingfromtheuse
oflargevaluesofturbulentviscosityandfromnotconsideringitsspatialvariability.
Another important downfall of the method is produced while modeling nonstationary flow regimes,
sinceinthesecases;turbulencevariesthroughspaceandtime.
3.4. ParabolicProfileofturbulentViscosity
Thismodelassumesaparabolicdistributionofturbulentviscosityindepth,thiswayconsideringamean
viscosityindepthwhichcanbecalculatedfromthefollowingexpression:
t 0.068 u f h
Wherehisthetotaldepthandufisthebedfrictionvelocity,calculatedfromthebedtangentialstress,
fromtheexpression:
f
uf

IftheManningformulaisusedtocalculatethebedfriction,thefollowingformulaforturbulentviscosity
canbeused:
t 0.068 g n U h 5/6
Thismeansthattheturbulentviscositydependsonthelocaldepth,thedepthaveragedvelocityandthe
manning coefficient. For better results, a coefficient is introduced to adjust the turbulent viscosity,
whichcanbehandledbytheuser.
t C m 0.068 g n U h 5/6
Page 23 of 56
3.5. MixingLengthModel
In the shallow waters mixing length model, the turbulent viscosity is calculated based on the local
characteristicsoftheflowusingthefollowingexpression:
2
u
t min 0.267 h, d wall 2S ij S ij 2.34 f
2
Where =0.41 is the von Karman constant. Its a relatively simple algebraic model that produces
acceptable results in flows where the turbulence is locally generated mainly from bed friction. It
considers the turbulence caused by horizontal velocity gradients, but it does not consider convective
transportorturbulencedissipation.Usingthemodelinflowswithstrongrecirculationzones,theresults
havelargeerrors.
3.6. TheRastogiandRodiKEmodel(Rastogi&Rodi,1978)
In this model k represents the turbulent kinetic energy and the rate of dissipation of the turbulent
energy. The model considers turbulence caused by bed friction, velocity gradients, dissipation and
convectivetransport.Itsolvesthefollowingshallowwaterequationsforkmodel:
k U x k U y k
t k
3
2 t S ij S ij c k u f
t x y x j k x j h

U x U y
c 1 2 t S ij S ij c u f c 2
4

t
2

t x y x j
x j k h2 k
k2 b 1
t c c k c f1/2 c 3.6c 3/2 1/2
k c 2 c cf
U2
Withthefollowingconstants:
c 0.09 c 1 1.44 c 2 1.92 k 1.0 1.31
Where k is the kinetic turbulent energy, is the turbulence dissipation rate, Sij is the deformation
sensor.Thetermswhichincludethebedfrictionvelocityuf areresponsibleformodelingtheformingof
turbulenceduetobedfriction.
The k model is relatively sophisticated. The results for shallow water turbulent flows are relatively
good,beingoneofthemostusedmodels,wheneverturbulencelevelsareimportant.Nevertheless,its
Page 24 of 56
levelofcomplexitydoesnotguaranteecorrectresultsforalltypesofflows.Likeanyotherturbulence
model,theresultsfromthekmodelmustbecriticallyanalyzed,forthis,userexperienceinturbulent
flowisessential.
3.7. DimensionalAnalysisoftheturbulenttermsintheshallowwaterequations
If the shallow water equations are turned dimensionless, the following nondimensional numbers are
obtained:
U 1 H UL UL
F T Rl Rt
gh Cf L t
Which respectively refer to the relationship between water mass inertia (convective forces) and
pressure forces (F), bed friction force (T), laminar tangential stresses (Rl) and turbulent tangential
stresses (Rt). The relative importance of these associated processes to each dimensionless number is
inverselyproportionaltothevalueofeachnumber.Forexample,thelargeradimensionlessnumberis,
thelessimportantitwillbefortheprocessitrepresents.Thatway,foralargelaminarReynoldsvalue,
theflowisturbulentandthusthelaminarforceslooseimportanceintheflowdevelopment.Likewise,
the importance of turbulent stresses in the mean velocity will depend on the value of the turbulent
Reynoldsnumber,whichdependsontheturbulentviscosity.Aparabolicmodelcanbeusedtoestimate
theorderofmagnitudeoftheturbulentviscosity,withthefollowingexpression:
1 1
t u f h g n h 5/6 U 0.21 n h 5/6 U
6 6
TheManningformulahasbeenusedtoestimatethebedfrictionvelocityuf.Thisestimationwillbemore
preciseincaseswheretheturbulenceisgeneratedmainlybybedfriction,asisthecaseinrivers,and
will provide some errors in cases where the turbulence is originated from horizontal shear stress, for
example in cases of recirculation zones. In any case, using this approximation, the turbulent Reynolds
numbercanbeexpressedas:
UL 4.8 L
Rt
t n h 5/6
Thisexpressioncanbeusedasafirstevaluationoftheimportanceofturbulentstressinthevelocityand
depth. For example, if we model a river section with a depth of 10m, a width of 400m, a manning
coefficientof0.025andameanvelocityof0.5m/s,theapproximateturbulentviscosityis0.02m2/sand
Page 25 of 56
the turbulent Reynolds number is close to Rt ~ 11000, which is a pretty large value, so the turbulent
stressesareexpectedtohaveasmalleffectonthemeanflow.
According to the expression for Tb, the importance of bed friction grows in shallow water flows, and
loosesimportanceastherelationshipbetweendepthandthehorizontaldimensiongrows.
Page 26 of 56
IBER
SEDIMENTTRANSPORTMODULE
Page 27 of 56
4. NONSTATIONARYSEDIMENTTRANSPORTMODULE
Thesedimenttransportmodulesolvesthenoncohesivesedimentnonstationarytransportequations.
The equations include the bedload transport equations and the suspended sediment transport
equations, coupling the bedload and the suspended load through a sedimentationrise term. The
sedimenttransportmoduleusesthevelocity,depthandturbulencefieldsfromthehydrodynamicand
turbulencemodules.Thebedloadiscalculatedusinganempiricalformula,chosenfromtheMeyerPeter
Muller or the Van Rijn Methods. The suspended load transport is modeled from a depth averaged
turbulenttransportequation.
Hydrodynamic
+
Turbulence
Bed-Load Suspended Load
Sediment
Conservation
Bed variations

Figure7.NonStationarysedimenttransportscheme.
4.1. SedimentConservationEquations
ThebedlevelvariationiscalculatedwiththeExnersedimentconservationequation:
Z b q sb,x q sb,y
1 p D -E
t x y
Wherepistheporosityofthesedimentthatformsthebedlayer,Zbisthebedlevel,qsb,xandqsb,yarethe
twosolidflowcomponents.ThedifferenceDEisthebalancebetweenthebedloadandsuspendedload.
Page 28 of 56
4.2. BedTransport
4.2.1. StressPartition
Thetotalbedloadstressintheriverbedisgeneratedbythegrainsroughness(whichisproportionalto
thesedimentdiameter)andtheshapeofthebedlayer(ripples,dunes,antidunes).Onlythegrainstress
contributes the bedload movement. Thus, before calculating the bed sediment flow it is necessary to
estimatethisgrainstress.Todothis,IberusestheEinsteinformulabasedonthetotalstress:
1.5
n K1/6
s (m)
= s
*
bs
*
b ns K s 2 3 D s
n 25
WherenisthetotalManningcoefficient,nsistheequivalentManningcoefficientduetothegrain,Dsis
thesedimentsdiameter,Ksthegrainsroughnesssheight(measuredfromthegrainsdiameter), b is
* *
the bed total stress, bs is the bed stress due to the grain, b , bs are dimensionless total and grain
stress,calculatedwith:
b bs
*b = *bs =
s g Ds s g Ds
Where s isthesedimentsdensity,and isthewaterdensity.Iberusesavalueof Ks 2.5 Ds
4.2.2. BedloadDischarge
TheBedloaddischargeiscalculatedbasedonempiricalformulas.Inthelatestversionofthemodel,two
wellregardedandwellknownformulasareused:
MeyerPeterMller
VanRijn
MeyerPeterMller(MeyerPeter&Muller,1948)
The original MeyerPeter & Mller equation, deduced from gravel beds of diameters up to 30mm,
calculatestheBedloaddischargewiththefollowingexpression:
q*sb = 8 *bs - *c = 8 *bs - *c

3/2 3/2

Wherethedimensionlessdischargeiscalculatedas:
Page 29 of 56
q sb
q*sb
s
1 g Ds
3
In the case of flat beds, a dimensionless critical bed stress is considered of c

0.047 . On the
contrary,itisnecessarytomakeabedslopecorrection.Detailedinsection4.2.3
Othercorrectionshavebeenmadetothisformula,(Wong,M,2003),(Wong&Parker,2006),resulting
in:
q*sb = 3.97 *bs - *c

3/2

If the bed is flat, c 0.0495 , is considered. On the contrary, it is necessary to make a bed slope

adjustment,detailedinsection4.2.3.ThiscorrectedversionisincludedinIber.
VanRijn(vanRijn,L.C.,1984)
InthevanRijnFormulathebedloaddischargeiscalculatedfromthefollowingexpressions:
T 2.1
T 0.3 q 0.053 0.3
*
sb
D*

T1.5
T 0.3 q 0.100 0.3
*
sb
D*
Where T is a dimensionless parameter that measures the excess of bed friction above a critical value
whichdefinesthemovementthreshold.
*bs - *c
T
*c
Thedimensionlessdiameterisdefinedas:
1/3
gR s -
D* Ds 2 con R

Page 30 of 56
4.2.3. Bedslopecorrection
Wheneverthebedisnotflat,theequationsmustbecorrectedtoincludetheeffectsofgravity,eitherto
increase the bedload discharge in a positive slope, or to reduce it in an adverse slope. The correcting
factor is applied over the critical stress of movement threshold, as detailed by Apsley and Stansby
(Apsley&Stansby,2008),wheretheypresentandgeneralizetheworksofmanyotherauthors,citingat
least(Dey,2003)and(Wu,2004).
Toconsiderthebedslopebothinthemovementthresholdandinthebedloaddischarge,thesediments
weight, due to the beds slope, is combined in a vector form with the bed stress to get an effective
stress.Sibisaunitvectorinthedirectionofthemaximumslope;thedimensionlesseffectivestressis
definedas:
*bs,eff = *bs + Do sin b
Where is the maximum slope, with the horizon, and D0 is a parameter depending on the particle
shape.Sothatthemovementbeginsintheabsenceofflow,whenever isequaltotheinternalfriction
angleofthematerial( ),theparameterD0isdefinedas:
*c,0
Do
tan
*
Where c,0 isthedimensionlesscriticalstressfortheflatbed.Ontheotherside,theeffectivecritical
stressisreducedproportionallytothegravitycomponentnormaltothebedslope.
*eff,crit = *c,0 cos
Where *c,0 is the same as before. Starting at this point, the bedload discharge equations presented
above are used, substituting the (critical and bed) stresses by effective stresses, and getting the
sediment discharge, which is a function of the fluids stress and the bed slope in each of the x, y
directions.
Theformulasusedareentirelyinvectorformforthebedloaddischarge,andthusabletoconsiderany
flowdirectionwithrespecttothemaximumslope
Page 31 of 56
4.2.4. Avalanchemodel
Apsley and Stansby (Apsley & Stansby, 2008) also proposed including an avalanche sliding model to
avoidthoseslopeshigherthanthebedmaterialsfrictionangle.Forthat,ifanangle betweentwo
finite volumes is higher than the friction angle , then a unit bedload discharge appears from the
upstreamtothedownstreamelementequalto:
0.5 L2 (tan tan )

q aval =(1 p)
cos t
WhereListhemaximumhorizontaldimensionoftheadjacentfinitevolumes.
4.2.5. Considerationsonanonerodiblelayer
Whilecalculatingbedloaddraggingandchangesinthebedlayersection,thepossibilitytoincludearock
layer,ornonerodiblesurface,belowwhichthebederosioncannotgofurther.
4.3. 2DTurbulentSuspendedLoadModule
4.3.1. TurbulentSuspendedLoadEquation
Thesuspendedloadtransportequationusesthevelocity,depthandturbulencefieldscalculatedinthe
hydrodynamic and turbulence modules. The suspended load is modelled using a depthaveraged
equation.TheequationusedintheIberis:
hC hU x C hU y C C Dsx Dsy
t h E D
t x y x j Sc,t
x j x y
Where C is the depthaveraged concentration of suspended solids, Ux, Uy are the horizontal depth
averagedvelocitycomponents, t istheturbulentviscosity,isthemoleculardiffusioncoefficientfor
suspended solids, and Sc,t is the Schmidt number, which relates the moment turbulent diffusion
coefficientwiththesuspendedturbulentdiffusioncoefficient.
The terms Dsx, Dsy model the suspended sediment dispersion due to the nonhomogeneous vertical
velocity profile and sediment concentration. Normally its effect is undermined in 2D shallow water
models,despiteitsimportancemayberelevantwhentheconcentrationsandvelocitiesvaryalongthe
depthofanelement,asisthecaseinopenchannelswithslamcurvatureradii.
Page 32 of 56
ThetermsEandDrespectivelymodelthosebedloadgrainswhichbecomesuspended(entrainment)and
deposition of suspended sediments to the bed layer. The difference represents a balance, and thus a
couplingofthebedloadandsuspendedload.
4.3.2. Entrainment/Deposition(ED)termCalculation
Three formulas are implemented for calculating the entrainment/deposition term (ED):Van Rijn (van
Rijn,1987),Smith(Smith,J.D.,1977)andAriathuraiandArulanandan(Ariathurai&Arulanandan,1978).
Thefirsttwoarevalidonsandbedlayers,whiletheAriathuraiisvalidforcohesivebedlayers.Thethree
formulasareespeciallyrecommendedinthelatestASCEManualofSedimentTransport,amongthem,
themostwidelyusedistheVanRijnFormula.
VanRijn
IntheVanRijnformula(vanRijn,1987),thetermEDisevaluatedfromthefollowingexpression
E D Ws c*a ca Ws C* C
Where isacoefficientthatrelatesthemeansuspendedparticleconcentrationandtheriverbedload
concentration, which value comes from the Rouse profile for the depth sediment concentration
distribution,Wsisthesedimentgrainvelocity,Cisthedepthaveragedsuspendedloadconcentration,
C* is the depthaveraged suspended load concentration at equilibrium condition (suspended load

transport capacity), c a and c*a are respectively the instantaneous concentration and the equilibrium

concentrationataheightz=aabovetheriversbedlayer,whereaistheheightofthelayerinwhichthe
transport is produced (theoretical separation limit between the bedload and suspended load). This
heightmustbeapproximatedfromthegraindiameter.Thecoefficient iscalculatedfromthevertical
concentrationdistribution(Rouseprofile),fromthefollowingintegral:
h-a
= ws
a = 3 D50
h h-z a k u *
a

z h-a
dz
Where=0.41isthevonKarmanconstant.
TheequilibriumconcentrationnearthebedlayerproposedbyVanRijn(vanRijn,1987)is:
Page 33 of 56
D50 T1.5
c*a 0.015
a D*0.3
1/3
gR
a = ks ks 3 Ds D* D 2

Smith
ThisformulaissimilartotheVanRijn,itsonlydifferenceliesintheformulatocalculatetheequilibrium
concentration,forwhichSmith(Smith,J.D.,1977)proposesthefollowingformula:
1.56 103 T
c
*
a
1 + 2.4 103 T
a = 26.3 *s - *c Ds + k s k s 3 Ds
AriathuraiyArulanandan
For cohesive soils the formula proposed by Ariathurai and Arulanandan (Ariathurai & Arulanandan,
1978)isused,inittheerosiondependsonthedifferencebetweenthetangentialstressandthecritical
tangentialstresstostarttheerosion ce ,itsalsodependantofavalueMrepresentativeoftheerosion
rate(equivalenttoanerosionratewhen b = 2 ce ):

E = M b 1
ce
In cohesive soils a modification is introduced for calculating D, which considers the deposition critical
tangentialstress cd .Inthiscase:
D P Ws C
With:

P 1 b if b < cd and P 0 otherwise
cd
Page 34 of 56
4.3.3. SedimentationVelocity
Thesedimentationvelocityiscalculatedbasedonitsdiameter(vanRijn,1987):
R g D50
2
Ws = D50 < 10-4 m
18
Ws =
10
D50

1+0.01 D*3 -1 10-4 m < D50 <10-3m
Ws =1.1 R g D50 10-3m < D50
Page 35 of 56
IBER
Shallowwatersturbulentflow
NUMERICALSCHEMES
Page 36 of 56
5. NUMERICALSCHEMES
All the equations described so far: the hydrodynamic equations (twodimensional shallow water
equations), the turbulence equations and sediment transport equations are solved using the finite
volume method. The finite volume method is one of the most widely and commonly used in
computationalfluiddynamics.InthissectionthenumericalschemesusedinIberarebrieflydescribed.
InthereferencespresentedinSection6moredetaileddescriptionsoftheusednumericalschemesmay
befound.
ThecharacteristicsoftheusednumericalschemesinallIbermodulesare:
FiniteVolumeSchemes,consideredinanintegralandconservativeway.
NonstructuredMesh.Meshesmadeupofelementsof3or4sides.
Abletosolveunsteadyflow(insubcritical,supercriticalandvariableregimes)
Abletosolvegraduallyvariedflow(movinghydraulicjumps,nonstationaryshockwaves)
It solves the hydrodynamic equations using the upwind High Resolution extension of the Roe
schemes(higherorderandnonoscillatory)
Upwindtreatmentofthebedlayerslopeterm.
Centredtreatmentoftheothersourceterms
FirstandSecondOrderschemesinspace.
Timeexplicitschemes
Nonstationary wetdry fronts treatment using stable and conservative schemes (conservation of
mass).
5.1. CalculationMesh
Tosolveadifferentialequationusingthefinitevolumemethodfirstisnecessarytospatiallydiscretize
the study domain. To do that, the study domain is divided in relatively small cells (calculation mesh).
Iberworkswithnonstructuredmeshofelementswhichmayhave3or4sides.Irregularelementsof3
or4sidesmaybecombinedinsidethesamemesh.Themainadvantageofworkingwithnonstructured
meshesistheeaseofadaptationtoanygeometry,sinceitisnotnecessary forthemeshtohaveany
Page 37 of 56
type of internal structure or organization. This attribute makes them especially functional for river
hydraulics.

Figure8.Exampleofanonstructuredmeshmadeupoftriangularelements
5.2. Discretizationinfinitevolumesforthe2DSWEEquations
For the discretization using the finite volume method, Iber Works with the twodimensional shallow
waterequationswritteninaconservativevectorform,as:
w Fx Fy
Gk
t x y k
WheretheconservativevectorvariablewandthevectorflowtermsFx,Fyrepresentthefollowing:

qx qy
h 2
q x gh 2 qxqy
w qx Fx Fy
h 2 h
2
qy q q q y gh 2
x y

h h 2
AndtheGktermsrepresentthesourcetermsincludedinthehydrodynamicequations,aspresentedin
section2.2.
Tomakethisspatialdiscretizationofthemassandmomentumconservationequationswiththefinite
volume method, the integral of the differential equations is solved at each of the meshs cells. This
procedure is especially advantageous to solve the conservation equations, since the equations are
solved in an integral form, which allows formulating the conservative methods in a simple way. The
Page 38 of 56
temporalandspatialDiscretizationofthetwodimensionalshallowwaterequationsintheverticalform
isexpressedas:
w in 1 w in
A i Fx n x Fy n y dL G k,i A i
t Li
k
5.2.1. Discretizationoftheconvectivefluxterms
Todiscretizethefluxterms,theupwindGodunovtypeconservativeschemesareused.TheupwindRoe
schemeisimplementedforafirstandsecondorderprecisioninspace.Inthosezoneswithrecirculation
floworwithlargespatialgradientsvelocitiesitisnotadvisabletousethefirstorderschemetosolvethe
hydrodynamicequations,sinceitpresentshighlydiffusivevelocityfields.
A conservative formulation of the equations, solved with upwind Godunov type, provides a good
solutionfortransonicshocks,soitisarecommendedmethodtomodelhydraulicjumps,dambreaksor
shockwaves.
Inthehydrodynamicequationsdiscretization,thecontourintegralcorrespondingtotheconvectiveflux
termsiscalculatedbasedonanumericalfluxfunctionas:
F n
Li x x Fy n y dL LR (w L ,w R ,n ij )
jK i
Where ij is a numerical flux function defined for each edge LR, where L and R are the left and right
nodes of the corresponding edge. A detailed description of the formulation of the numerical flux for
upwindGodunovtypeschemescanbefoundinthereferencesofsection6.
FirstorderupwindRoescheme
IberusestheupwindRoescheme,whichsolvesthenumericalfluxwiththefollowingexpression:
ZL + ZR 1
LR J LR w R w L
2 2
Z
Z Fx n x Fy n y J= J = X D X -1 J = X D X -1
w
WhereZLyZRrepresentthefluxnormaltothecontouratbothsidesoftheLRedge.Thematrix|J|LRis
theabsolutevalueofthejacobianmatrixoftheflowZ,evaluatedinthemeanRoestate,definedas:
Page 39 of 56
hL + hR h L U x,L h R U x,R h L U y,L h R U y,R

h = hL hR c = g =
U x
=
U y
2 hL hR hL hR
ThevaluesandvectorsemofthejacobianmatrixJ,canbewrittenas:
1 2 c n 2x + n 2y 2 n x U
n U
x y

y 3 2 - c n 2x + n 2y
1 1 1

e1 U x c n x e 2 - c n y e 3 U x - c n x
U c n U
y c n y x y - c n y
For its implementation in the numerical flux calculation in Iber, the difference in states (wRwL) is
decomposedtotheleftandrightoftheedgeconsideringthevectorbaseem:
3
wR - wL =
m=1
m e m
Writingthenumericalfluxas:
Zi + Z j 1 3
LR
2
-
2
m=1
m m e m
Withthecoefficientscalculatedas:
hR - hL 1
1 =
2
+ U x,R h R - U x,L h L n x + U y,R h R - U y,L h L n y - U
2c
x x
n h - h
n + U
y y R L
1
c
R L
h - h n U h - U h - U
2 = U y,R h R - U y,L h L - U y x
x,R R x,L L
h - h n
x R L y
hR - hL 1
3 =
2
U x,R h R - U x,L h L n x + U y,R h R - U y,L h L n y - U
2c
x x
n h - h
n + U
y y R L
Thedescribedschemeisoffirstorderinspacescheme.Theupwindconvectivefluxisequivalentfroma
mathematicalpointofviewtoaddingadiffusionterm(whichisgenerallycallednumericalorartificial
diffusion) with a diffusive (numerical) coefficient proportional to the meshs size. Which is why it is
convenient to use fine meshes to reduce the error introduced by this numerical diffusion, it is also
convenienttousehigherorderschemes.
Page 40 of 56
NumericalFluxextensiontosecondorder
The scheme described is the first order scheme due to the introduced numerical diffusion at the
convective flux discretization. Despite it, this scheme is widely used as the default scheme in many
computationalfluiddynamicschemes,duetoitsnumericalstability.Whenmoreprecisionisrequiredin
aregularsizemesh,itisnecessarytoappealtohigherorderschemes.InIbershydrodynamicmodule,
theorderofprecisioncanbeincreasedthroughanextensionofthesecondordernumericalfluxandits
TVD(TotalVariationDiminishing)limitation:
Zi + Z j 1 3
LR
2
-
2

m=1
m m (1- m (1- m )) e m
Where m m t / d LR and d LR thedistancebetweentheLandRelements.
Thefollowingfluxlimitatorsareincluded:
Minmod(bydefault)
Superbee
VanLeer

Whichareafunctionofthe rm parameter,anindicatorofthevariablejumpbetweenupwindedgeand

calculationedge:
(1- )
m m m
rm i,j = upwind
(1- )

m m m
i,j
Theimplementedexpressionsare:
Minmod: (r)=max 0,min r,1
Superbee: (r)=max 0,min 2r,1 ,min r,2
0 si r 0

VanLeer: (r)= 2r
1+r si r>0
Page 41 of 56
5.2.2. Discretizationofthebedslopesourceterm
Iberusesacentreddiscretizationofallthesourcetermsexceptofthebedslopeterm.Themaincause
ofusinganupwinddiscretizationofthebedslopeinsteadofacentreddiscretizationisthatitcalculates
exactlythehydrostaticsolutioninanirregularbathymetry,thiswayavoidingspuriousoscillationsinthe
freesurfacewatersandinthevelocities.Theseoscillationsareingeneralsmall,buttheycanbecome
larger in magnitude in irregular bathymetry problems, such is the common case in river and coastal
hydraulics.
Whenever the upwind numerical schemes are used to discretize convective flow, the upwind
discretizationofthebedslopeprovideshasbetterpropertiesandprovidesmorepreciseresultsthanthe
classiccentredschemes.
TheuseddiscretizationofthebedslopesourceterminafinitevolumeCicanbeexpressedas:
Si = S dA =
Ci
S
jKi
ij
WhereSijisanupwinddiscretizationofthebedslopesourcetermineachedgeoftheconsideredfinite
volumeandiscalculatedas:
0
| nij | hL hR -1 -1
S ij = -g zb , R zb , L I X | D | D X n x ,ij
2 2 n
y ,ij
Thesameasintheconvectiveflux,toimplementthisdiscretizationinIber,itisdecomposedusingthe
vectorsemas:
3
Sij = mem
m 1
1
1C cij zij | n ij |
2
2C 0
1
3C cij zij | n ij |
2
Page 42 of 56
5.3. DiscretizationoftheTransportEquationsintheKETurbulenceModel,and
intheSuspendedloadTransportModel
The conservation equations used in the k turbulence model and in the suspended load transport
modelcanbeexpressedas:
Fj
S
t x j
Where = C h is the conserved variable, C is the nonconserved variable, corresponding to the
considered model (sediment concentration, kinetic turbulent energy, turbulent dissipation rate), S are
thesourcetermswhichmodelgenerationordestructionprocessesoftheconservedvariable,andFis
theconvectiveanddiffusiveflux,whichmaybeexpressedas:
C
Fj C h u j e h
x j
Whereeistheeffectivediffusivitycoefficient,whichincludesthemolecularandturbulentdiffusion.In
general, turbulent diffusion is many times larger than the molecular diffusion, meaning that this term
canbeomitted.
In the k turbulence model it is necessary to solve 2 transport equations; one due to the turbulent
kineticenergy(k)andtheotheronefortheturbulencedissipationrate().
5.3.1. DepthaveragedTransportEquation
Tospatiallydiscretizethetransportequationusingthefinitevolumemethodthedifferentialequationis
integrated in each cell of the computational mesh. This process is especially advantageous for solving
the conservation equations, since they are being solved in an integral form, which allows the
formulation of natural conservative methods. Integrating the convectiondiffusion equation in a two
dimensionalcellweobtain:
C h i C h i
n 1 n
A i C h u dA e h C dA E - D dA
t Ai Ai Ai
Where =C.h is the average value at the cell and S=ED is the source term due to
generation/destruction processes of the considered variable. Applying the divergence theorem to the
secondandthirdtermweobtain:
Page 43 of 56
C h i C h i
n 1 n
A i C h u n dL e h C n dL E - D dA
t Li Li Ai
Where the integrals around the cells area are extended to the cells edges. This equation is the
conservation equation expressed in an integral form. The integrals that appear in the integral form
conservationequationaresolveddiscretely,thiswaygetting:
C h i C h i
n 1 n
t
Ai C h u n L
j K i
ij ij h C n L E - D A
j K i
ij ij i i
Wherethesummationsareextendedtoalltheedgesthatformtheboundaryofthecell.Thesubindex
ijidentifiesthecommonfaceofthecellsiandj.Eachtermofthesummationsrepresentafluxofthe
considered variable that enters the cell through the corresponding face, that way the sum of all the
fluxesthroughallthefacesthatformthecontourofthecessareequaltothebalanceofwhatgoesout
minuswhatgoesin,forexample;thenetfluxoutsidethecell.Infunctionoftheusedinterpolationto
calculatethefluxthroughthecellscontours,especiallytheconvectiveflux,differentnumericalschemes
canbefound.
The diffusive flux between cells is calculated using a centred discretization of second order, without
numericalstabilityproblems.Theconvectivefluxdiscretizationismoreproblematicfromthenumerical
stabilitypointofview.Inthenextparagraphssomeimplementednumericalschemesarepresentedin
thecodetodiscretizetheseterms.
FirstOrderUpwindscheme
Intheupwindschemesthedirectionoftheinformationflowisconsideredtomakethediscretizationof
theconvectiveflux.Thediffusivefluxisdiscretizedusingacentredschemeofsecondorder.Inthecase
oftheconvectiveflux,theinformationtravelsinthesamedirectionasthewatervelocity,forexample
upstream to downstream. The simplest upwind method is the first order upwind scheme, calculated
using:
C h u n ij u n ij C h ij un,ij C h ij
The velocity is discretized in a centred manner, while the value of C.h is discretized in an upwind
manner,consideringtheupstreamnodevalue:
Page 44 of 56
u n ,ij u n ,i (1 ) u n ,j
C h ij C h i si u n ,ij 0
C h ij C h j si u n ,ij 0
Where isalinearinterpolationcoefficient.Foranequidistantdiscretizationthevalueof is0.5.
Thefactofupwindingtheconvectiveterminthefirstorderschemeisequivalentfromthemathematical
pointofviewasaddingadiffusiveterm(usuallycallednumericaldiffusion)withadiffusivity(numerical)
coefficient n proportional to the mesh size. The numerical diffusion is a consequence of the used
techniquefornumericalstabilization.Thedesiredsituationinanumericalschemeisforittobestable
withtheintroductionofaslittleaspossiblenumericaldiffusion.Thisiswhyitssoconvenienttousefine
meshestoreducetheerrorintroducedbynumericaldiffusion
SecondOrderUpwindscheme
The previously mentioned scheme is of first order because the numerical diffusion introduced in the
discretizationoftheconvectiveflux.DespitethisitiswidelyusedinComputationalFluidDynamic(CFD)
codesasthedefaultscheme,duetoitsnumericalstability.Whenmoreprecisionisrequiredinaregular
sizemesh,itisnecessarytoappealtohigherorderschemes.Todothisasecondorderschemeoftype
MUSCL(MonotonicUpstreamSchemeforConservativeLaws)hasbeenimplementedIber,whichmakes
alinearreconstructionofeachmeshelementofnonconservedvariable.
Oncethislinearreconstructionofthenonconservedvariableismade,theconvectivefluxiscalculated
ateachedgeofthecomputationalmeshwith:
C h u n ij u n ij C h ij un,ij C h ij
u n ,ij u n ,i (1 ) u n ,j
C h ij C h Ij si u n ,ij 0
C h ij C h iJ si u n ,ij 0
Where C h I j and C h i J thevaluesofthevariableatedgeLijtakenfromthelinearreconstruction
ofcellsCiyCjrespectively.
5.4. DiscretizationoftheExnerConservationofSedimentEquation
TheExnerequationofconservationofsedimentcanbeas:
Page 45 of 56
z b q sb,x q sb,y
(1 p) DE
t x y
The Exner equation is discretized in a similar way as the suspended load sediment conservation
equation.TheintegraloftheExnerequationineachtwodimensionalcellcanbewrittenas:
n 1
zb,i zb,i
n
q q sb,y
(1 p) A i sb,x + dA D - E i A i
t Ai
x y
Applyingthedivergencetheoremtothesecondtermoftheequation,weget:
n 1
zb,i
n
q n L D - E A
zb,i
(1 p) Ai *
sb ij i i
t jK i
ij
Thevalueofthebedloadineachedgeofthemesh q *sb iscalculatedinanupwindformas:

ij
q
*
sb ij q sb,i si q*sb n 0
ij

q
*
sb ij q sb,j si q n 0
*
sb ij
Theerosion/depositionsourcetermisdiscretizedinacentredformineachmeshcell.
5.4.1. Nonerodiblebedlayer
Thereisthepossibilitytoincludearocklayer,meaninganonerodiblebedlayer.Theimplementationin
thecodeofthisnonerodiblebedlayerhasbeendoneusinga2stagecalculation(predictorcorrector)
ofthenumericalbedloadfluxesthroughtheelementscontours.
Predictor:searchesanewbedlayerelevationduemainlytotheoutflowinanelement.
max 0, q
t
z ip = z in + *
sb n Lij
(1 p) A i jK i
ij
Corrector: if at the predictor the element has eroded below the rock layer, the bedload outflow is
correctedeneachoftheedges,andafterwardsitisupdatedtothenewbedlayerelevation.
p
Si (zi < z roca ) y ( q*sbn > 0) :
ij
Page 46 of 56
z in -z roca
q = q
*,c
sb ij
*
sb ij
z in -z ip
q *,c
sb ji = - q*sb
ij
n 1
zb,i
n
q n L D - E A
z b,i
(1 p) Ai *,c
sb ij i i
t jK i
ij
5.5. ManagingWetDryFronts
Modeling flooding zones, as is the movement of wave fronts in estuaries and coastal zones is
fundamental in environmental hydraulics. In Iber wetdry fronts, both stationary and nonstationary,
thatmayappearinthedomainaremodelledwithafixedfinitevolumemesh,allowingeachelementto
havewaterornotdependingontheflowconditions.Amongthevolumeswhichdonothavewaterand
thosewhichdohave,thereisadrywetfrontwhichisnecessarytotreatappropriatelyfromanumerical
point of view to avoid instabilities and nonphysical oscillations in the results. For this wetdry front,
eitherinafloodfrontorawavefront,thereisawetdrytolerancewd,soifthewaterdepthislower
thanwditisconsideredthatthecellisdryandisnotincludedinthecalculations.Thewetdrytolerance
canbechangedbytheusertoavalueclosetozero,eventhoughinirregularbathymetryproblems,asis
very usual in river an coastal hydraulics, it is advisable to use values in the order of 1cm or 1mm to
increasethestabilityoftheresultswithoutspoilingtheirprecision.Inanycase,thewaterheightisnever
forced to be equal to zero, with the purpose of avoiding mass losses in the calculation domain. The
numericalschemeusedtosolvethewetdryfrontisstableandnondiffusive.
ThemanagementthewetdryfrontsinIberisstable,conservativeandnondiffusive,meaningthatthe
fronts are adequately solved without instabilities of the numerical type, even when these occur in
strongbedslopes.Everyfiniteelementhasanassociatedbedlevel.SchematicallyshowninFigure9.
Page 47 of 56
Figure9.Schematicrepresentationofthewetdrybedlevelmanagement.
Betweentwovolumeswithdifferentbedleveloneofthefollowingsituationsmayoccur,aspresented
inFigure10:

Figure10.Severalsituationsofthewaterlevelsbetweentwoadjacentcells.
In the first figure both volumes have water, so no front is produced and thus no special treatment is
necessary.Intheothertwocasesthereisawetdryfront.Thedifferenceisthatinthesecondcase,the
freewaterlevelinthewetcellishigherthanthebedlevelofthedrycell,whileinthethirdcaseitis
lower.Onlyinthethirdcaseitisnecessarytouseaspecialtreatment,whichconsistsofredefiningthe
bedslopetoimposeareflexiveconditioninthefront.Inthiscasethebedslopeisredefinedas:
h i - h j si h j z b,j - z b,i
z b,ij =
z b,j - z b,i si h j > z b,j - z b,i
Thereflexiveconditionisimposedas:
qn,ij = q x,ij n x,ij + q y,ij n y,ij 0
The use of the mentioned conditions provides can reproduce the hydrostatic solution for any
bathymetry,withoutstreatchingthefrontandwithoutgeneratinganyspuriousoscillationsinthefree
surface. This type of wetdry front treatment has been successfully used in the modeling of both
stationary and nonstationary processes, being particularly useful in the simulation of flood zones in
riversandcoastalzones,andincalculatingtheevolutionofwavefronts.
AmoredetaileddescriptionoftheimplementationinIbercanbefoundinthereferencesprovidedin
section6.
Page 48 of 56
IBER
Shallowwatersturbulentflow
REFERENCES
Page 49 of 56
6. REFERENCES
Apsley,D.,&Stansby,P.(2008).BedLoadSedimentTransportonlargeslopes:ModelFormulationand
ImplementationwithinaRANSsolver.JournalofHydraulicEngineering,ASCE,vol134(10).
Ariathurai,R.,&Arulanandan,K.(1978).Erosionrateofcohesivesoils.JournaloftheHydraulicsDivision.
ASCE,104(HY2),279283.
Bermudez, A., Dervieux, A., Desideri, J., & VazquezCendon, M. (1998). Upwind schemes for the two
dimensional shallow water equations with variable depth using unstructured meshes. Computational
MethodsAppliedforMechanicalEngineering,Vol155.
Blad, E., & GmezValentn, M. (2006). Modelacin del flujo en lmina libre sobre cauces naturales.
Anlisisintegradoenunaydosdimensiones.Monografa.Barcelona:CIMNE.
Blad, E., GmezValentn, M., SnchezJuny, M., & Dolz, J. (2008). Preserving steady state in Finite
VolumeComputationsofRiverFlow.JournalofHydraulicEngineering,ASCE,Vol134(9).
Cea,L.(2005).Unstructuredfinitevolumemodelforunsteadyturbulentshallowwaterflowwithwetdry
flows:numericalsolveransexperimentalvalidation.Phdthesis.ACorua:UniversidaddeACorua.
Cea, L., French, J., & VzquezCendn, M. (2006). Numerical modelling of tidal flows in complex
estuaries including turbulence: An unstructured finite volume solver and experimental validation.
InternationalJournalforNumericalMethodsinEngineering,vol67(13).
Cea,L.,Puertas,J.,&VzquezCendn,M.(2007).Depthaveragedmodellingofturbulentshallowwater
flowwithwetdryfronts.ArchivesofComputationalMethodsinEngineering,Stateoftheartreviews,
vol14(3).
Chiew, Y., & Parker, G. (1994). Incipient Sediment Motion on NonHorizontal Slopes. Journal of
HydraulicsResearch,649660.
Chow,V.,Maidment,D.,&Mays,L.(1988).AppliedHydrology.McGrawHill.
CIMNE.(2009).GiDThepersonalpreandpostprocessor.Barcelona:www.gidhome.com.
Corestein,G.,Blad,E.,GmezValentn,M.,Dolz,J.,Oate,E.,&Piazzese,J.(2004).NewGiDinterface
forRamfloodDSSprojecthydraulicsimulationcode.Proceedingofthe2ndConferenceonadvancesand
applicationsofGiD.Barcelona:CIMNE.
Dey. (2003). Threshold of sediment motion on combined transverse and longitudinal sloping beds.
JournalofHydraulicResearch,#4405415.
Page 50 of 56
Garca,M.,&Parker,G.(1991).EntrainmentofBedSedimentintosuspension.ASCEJournalofHydraulic
Engineering,#4414435.
MeyerPeter,E.,&Muller,R.(1948).Formulasforbedloadtransport.Proceedingofthe2ndConference
(pp.3964).Stockholm:IAHR.
Rastogi, A., & Rodi, W. (1978). Predictions of heat and mass transfer in open channels. Journal of
HydraulicsDivision,104(3)397420.
Smith, J., & McLean, S. (1977). Spatially averaged flow over a wavy surface. Journal of Geophysical
Research,82(12)17351746.
Smith, J.D. (1977). Modelling of sediment transport on continental shelves. The Sea: ideas and
observationsonprogressinthestudyoftheseas.NewYork:JohnWileyandSons.
VanDorn,W.(1953).WindStressonanartificialpond.JournalofMarineResearch,Vol12.
van Rijn, L. (1987). Mathematical modelling of morpholigical processes in the case of suspended
sedimenttransport.Delft:DelftHydraulicsCommunications#382.
van Rijn, L.C. (1984). Sediment Transport, part I: Bedload transport. Journal of Hydraulic Engineering ,
vol110(10).
Wong,M.(2003).DoesthebedloadequationofMeyerPeterandMullerfititsowndata?30thCongress
(pp.7380).Thessaloniki:InternationalAssociationofHydraulicResearch.
Wong,M.,&Parker,G.(2006).ReanalysisandCorrectionofBedloadrelationofMeyerPeterandMuller
usingtheirowndatabase.JournalofHydraulicEngineering,132(11)11591168.
Wu,W.(2004).Depthaveragedtwodimensionalnumericalmodelingofunsteadyflowandnonuniform
sedimenttransportinOpenChannels.JournalofHydraulicEngineering,130(10)10131024.
7.
Page 51 of 56
NOMENCLATURE
awidthofthelayerinwhichbedloadtransportoccurs
c a , c*a instantconcentrationandequilibriumconcentrationofsuspendedsolidsataheightz=aabove
theriverbedlevel
Cdepthaveragedsuspendedload
C* Depthaveragedsuspendedloadinequilibriumconditions
Cddischargecoefficientatagateoraweir
Cfbedfrictioncoefficient
Cmadjustcoefficientfortheturbulentviscosityinaparabolicmodel
dwallwalldistance
Dijlateraldispersionterms
Dssedimentdiameter
D50meansedimentdiameter
D* Nondimensionalsedimentdiameter
Dsx, Dsy Suspended load dispersion due to the nonhomogeneity of the velocity profiles and of the
verticaldirectionsedimentconcentration
DEbalancebetweenthebedloadandsuspendedload
fpotentialinfiltrationrate
ggravityacceleration
hdepth
irealinfiltrationrate
kturbulentkineticenergy
kssoilssaturatedpermeability
Ksgrainsroughnessheight
kvdwinddragcoefficient
Page 52 of 56
L0Initialdepthofthesoilssaturatedregion
LDepthofthesoilssaturatedregion
M2parameterthatmeasurestheerosionratebecauseofresuspension
Msmasssource/drainterm
Mx,Mymomentumsource/drainterms
nManningcoefficient
nsequivalentManningcoefficientduetothegrain
pporosityofthebedloadsediments
qnunitnormalflowattheinletcontour
qsb,x,qsb,ycomponentsofthebedloaddischarge
Qdischarge/Flow
Rhhydraulicradius
Rnondimensionalsubmergedspecificweight
Sc,tSchmidtnumber
SeEffectiveinitialsoilsaturation
Sijdeformationtensor
Tnondimensionalparameterwhichmeasurestheexcessoffrictionatthebedlayerabovethecritical
valuethatdefinesthestartofmovement
uffrictionvelocityduetobedroughness
u'i u' j TurbulentstressesorReynoldsstresses
Ux,Uydepthaveragedhorizontalvelocities
|U|depthaveragedmeanvelocity
V10Windvelocityat10metersheight
Wssedimentationvelocityofsolidparticles
ZsFreelayerlevel
ZbBedlayerlevel
Page 53 of 56
coefficient thatrelatesthemean concentrationofsuspendedsolidsandtheconcentrationnearthe

riverbed
Bedslope
Moleculardiffusionofsuspendedsolidscoefficient
ijKroneckerdelta
turbulencedissipationrate
changesinthesoilshumidcontentasthesaturationadvances
iinitialsoilshumiditycontent
esoilseffectiveporosity(drainable)
rretentioncapacity(irreduciblehumidityornondrainable)
=0.41vonKarmanconstant
latitudeoftheconsideredpoint
cbedmaterialsinternalfrictionangle
fluidscinematicviscosity
tturbulentviscosity
waterdensity
s Sedimentdensity
b Totalstressduetobedfriction

bs Bedstressduetograin
c Criticalbedstress
*b , *bs Totalgrainstresses(nondimensional)
*c Criticalbedstress(nonadimensional)
sfreesurfacefrictionduetowindinducedroughness
exx,exy,eyytangentialeffectivehorizontalstresses
ijv Viscousstresses
Page 54 of 56
Soiltotalporosity
soilnonsaturatedreginsuction
earthsrotationangularvelocity
Wssedimentvelocity
Page 55 of 56

Hydraulic Reference Manual Iber

Uploaded by

Copyright:

Available Formats

Hydraulic Reference Manual Iber

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Hydraulic Reference Manual Iber

Uploaded by

Copyright:

Available Formats

Twodimensionalmodellingoffreesurfaceshallow

Two-dimensional modelling of free surface shallow water flow

2. HYDRODYNAMIC MODULE ...................................................................... 7

2.1. Introduction ........................................................................................7

2.2. Hydrodynamic Equations .......................................................................7

2.3. Bed Friction .........................................................................................8

2.4. Superficial Friction due to Wind ..............................................................9

2.5. Effective Stresses.................................................................................9

2.6. Hydrodynamic boundary Conditions ...................................................... 11

2.6.1. Closed Boundaries ........................................................................ 11

2.6.2. Open Boundaries .......................................................................... 13

2.7. Internal boundary conditions ............................................................... 14

2.7.1. Gates .......................................................................................... 15

2.7.2. Weirs .......................................................................................... 15

2.7.3. Gate-Weir Combination.................................................................. 16

2.7.4. Local Losses ................................................................................. 16

2.8. Infiltration ........................................................................................ 17

2.8.1. Green-Ampt ................................................................................. 17

2.8.2. Horton ........................................................................................ 18

2.8.3. Linear ......................................................................................... 19

2.9. Initial Abstraction............................................................................... 19

3. TURBULENCE MODULE .......................................................................... 21

3.1. Introduction ...................................................................................... 21

3.2. Scales of Turbulence in Shallow Waters ................................................. 22

3.3. Constant Turbulent Viscosity................................................................ 23

3.4. Parabolic Profile of turbulent Viscosity ................................................... 23

3.5. Mixing Length Model ........................................................................... 24

4. NON-STATIONARY SEDIMENT TRANSPORT MODULE ............................ 28

4.1. Sediment Conservation Equations......................................................... 28

4.2. Bed Transport.................................................................................... 29

4.2.1. Stress Partition ............................................................................. 29

4.2.2. Bedload Discharge ........................................................................ 29

4.2.3. Bed slope correction ...................................................................... 31

4.2.4. Avalanche model .......................................................................... 32

4.2.5. Considerations on a non-erodible layer ............................................ 32

4.3. 2D Turbulent Suspended Load Module................................................... 32

4.3.1. Turbulent Suspended Load Equation ................................................ 32

4.3.2. Entrainment/Deposition (E-D) term Calculation ................................. 33

4.3.3. Sedimentation Velocity .................................................................. 35

5. NUMERICAL SCHEMES .......................................................................... 37

5.1. Calculation Mesh ................................................................................ 37

5.2. Discretization in finite volumes for the 2D-SWE Equations ........................ 38

5.2.1. Discretization of the convective flux terms ........................................ 39

5.2.2. Discretization of the bed slope source term ...................................... 42

5.3.1. Depth-averaged Transport Equation ................................................ 43

5.4. Discretization of the Exner Conservation of Sediment Equation ................. 45

5.4.1. Non-erodible bed layer .................................................................. 46

5.5. Managing Wet-Dry Fronts .................................................................... 47

ije ijv u'i u' j D ij

WhereCdis theweirsdischarge coefficient,Zsisthewaterlevelupstream oftheweirandZwisthe

Boundary Regime Boundaryconditions

Total discharge in normal direction to the

Inlet Specific discharge in normal direction to

(ZD ZB ) / (ZU ZB ) DischargeEquations

Transition 0.670.80 Q = Cd B h 6g (ZU ZD )

SubmergedFlow 0.801.00 Q = Cd B h 2g (ZU ZD )

(ZD ZW ) / (ZU ZW ) DischargeEquation

FreeFlow <0.67 Q = Cd B (ZU ZW )1.5

SubmergedFlow >0.67 Q = 2.598 Cd B (ZD ZW ) (ZU ZW )0.5

qsb = 8 bs - c = 8 bs - *c

qsb = 3.97 bs - *c

bs,eff = bs + Do sin b

eff,crit = c,0 cos