Fully coupled approach to modeling shallow water flow, sediment transport, and bed evolution in

20 Pages·2011·4.43 MB·English

WATERRESOURCESRESEARCH,VOL.47,W03508,doi:10.1029/2010WR009751,2011
Fully coupled approach to modeling shallow water flow, sediment
transport, and bed evolution in rivers
ShuangcaiLi1andChristopherJ.Duffy1
Received12July2010;revised6December2010;accepted21December2010;published2March2011.
[1] Ourabilitytopredictcomplexenvironmentalfluidflowandtransporthingesonaccurate
andefficientsimulationsofmultiplephysicalphenomenonoperatingsimultaneouslyovera
widerangeofspatialandtemporalscales,includingoverbankfloods,coastalstormsurge
events,dryingandwettingbedconditions,andsimultaneousbedformevolution.This
researchimplementsafullycoupledstrategyforsolvingshallowwaterhydrodynamics,
sedimenttransport,andmorphologicalbedevolutioninriversandfloodplains
(PIHM_Hydro)andappliesthemodeltofieldandlaboratoryexperimentsthatcoverawide
rangeofspatialandtemporalscales.Themodelusesastandardupwindfinitevolume
methodandRoe’sapproximateRiemannsolverforunstructuredgrids.Amultidimensional
linearreconstructionandslopelimiterareimplemented,achievingsecond-orderspatial
accuracy.Modelefficiencyandstabilityaretreatedusinganexplicit-implicitmethodfor
temporaldiscretizationwithoperatorsplitting.Laboratory-andfield-scaleexperimentswere
compiledwherecoupledprocessesacrossarangeofscaleswereobservedandwherehigher-
orderspatialandtemporalaccuracymightbeneededforaccurateandefficientsolutions.
Theseexperimentsdemonstratetheabilityofthefullycoupledstrategyincapturing
dynamicsoffield-scalefloodwavesandsmall-scaledrying-wettingprocesses.
Citation: Li, S., and C. J. Duffy (2011), Fully coupled approach to modeling shallow water flow, sediment transport, and bed
evolutionin rivers,WaterResour.Res., 47,W03508, doi:10.1029/2010WR009751.
1. Introduction [3] Theshallowwaterequationsaretypicallyusedtorep-
resentthehydrodynamicsofriverfloods,stormsurges,tidal
[2] Water flow and sediment transport are simultaneous
fluctuations, tsunami waves, and forces acting on offshore
and interactive processes in rivers, floodplains, and coastal
structures [Aizinger and Dawson, 2002]. Methods for solv-
areas. The interaction among these processes is influenced
ingtheshallowwaterequationsincludethemethodofchar-
by both human activities and extreme natural events,
acteristics [e.g., Katopodes and Strelkoff, 1978], finite
resulting in aggradation and degradation in channels and
difference [e.g., Molls and Chaudhry, 1995], finite element
harbors, deterioration of water quality and fisheries among
[e.g.,Hervouet,2000],andfinitevolume[e.g.,Alcrudoand
other environmental effects, and many other forms of eco-
Garcia-Navarro, 1993; Zhao et al., 1994; Anastasiou and
logical disturbance. Examples include dam removal, dam
Chan, 1997; Sleigh et al., 1998; Toro, 2001; Bradford
break, and extreme storm events that induce rapidly vary-
and Sanders, 2002; Valiani et al., 2002; Yoon and Kang,
ingflowandsedimentflushing.Thedisturbanceiscomplex
2004; Begnudelli and Sanders, 2006]. Although each
because of the uneven and changing bottom topography,
methodhasitsownstrengthsandlimitations,itisgenerally
irregular boundaries, rapid and strong erosion with abrupt
true that unstructured grids have advantages for represent-
bed and flow variations, and complicated and uncertain
ing natural channels. An algorithm for ‘‘optimal’’ unstruc-
flow-sediment transport mechanisms. Under these condi-
tured grids was proposed by Shewchuk [1997], which is
tions, one-dimensional, uncoupled strategies are generally
able to provide an ‘‘optimal’’ representation of the domain
not sufficient, and two-dimensional approaches capable of
withtheleastnumberofelementswhilestillconformingto
handling complicated geometry, rapidly varying flow, and
alimitedsetofphysicalandgeometricconstraintsparticular
fully coupled physics are necessary. This research builds
tothephysicaldomain.
onrecentadvancesfornumericalsolutionstofullycoupled
multiphysics problems in engineering and computational [4] With respect to the numerical method, the finite vol-
umemethodallowsforlocaland global mass conservation,
fluid dynamics to flow, sediment, and bed morphology
can be applied to structured or unstructured grids, and
interaction in rivers and tests the model over a range of
requires less memory for explicit calculations as compared
scalesinlaboratoryandfieldexperiments.
to finite difference or finite element methods [Loukili and
Soulaimani, 2007]. Several investigators have solved the
shallow water equations on unstructured grids using finite
1Department of Civil and Environmental Engineering, Pennsylvania
volume methods [Zhao et al., 1994; Anastasiou and Chan,
StateUniversity,UniversityPark,Pennsylvania,USA.
1997; Sleigh et al., 1998; Yoon and Kang, 2004], although
sedimenttransportwasnotconsideredinthosemodels.The
Copyright2011bytheAmericanGeophysicalUnion.
0043-1397/11/2010WR009751 coupled behavior of sediment transport and bed elevation
W03508 1of20
W03508 LIANDDUFFY:SHALLOWWATERFLOWANDSEDIMENTTRANSPORTINRIVERS W03508
changes was experimentally studied by Capart and Young and an incompressible fluid appropriate for vertically well-
[1998], with implications for the dynamics of the flow re- mixed water bodies. The system in conservative form is
gimeaswell.Forcoupledsedimenttransportandbedevolu- writtenas
tion the assumption of nonequilibrium conditions enforces
@ð(cid:2)hÞ @ð(cid:2)uhÞ @ð(cid:2)vhÞ @ð(cid:2)zÞ
a dynamic exchange between sediment deposition and þ þ ¼(cid:2) z þ(cid:2) S ; ð1Þ
@t @x @y @t w p
entrainment,andthiscouplingisexploredinthispaper.
are[5a]vaRielalabtlievetolystfuedwym2-oDdeclosuapnleddfhewyderrodfiyenldamoibcseflrovwatiaonnds @ð(cid:2)uhÞ @h(cid:2)(cid:2)u2hþg2h2(cid:3)i @ð(cid:2)uvhÞ (cid:4) (cid:5)
þ þ ¼(cid:2)(cid:2)gh S þS ; ð2Þ
sediment transport. An example is a large flood event or @t @x @y 0x fx
dam break on an initially dry surface where full water–
sediment–bed form coupling is likely to be important. The h (cid:2) (cid:3)i
dry-to-wet transition followed by a wet-to-dry transition @ð(cid:2)vhÞ @ð(cid:2)uvhÞ @ (cid:2) v2hþg2h2 (cid:4) (cid:5)
þ þ ¼(cid:2)(cid:2)gh S þS : ð3Þ
with an evolving bed surface during postevent relaxation @t @x @y 0y fy
produces interesting multiscale behavior. Recently, several
1-D models were developed to simulate the dam-break- [8] Massconservationequationsareused to describe the
induced sediment transport or high-concentration sediment sediment transport and morphological evolution process.
transportasinhyperconcentratedflowanddebrisflow[Bel- There are two approaches to coupled sediment routing and
los and Hrissanthou, 1998; Fraccarollo et al., 2003; Cao bed evolution, i.e., noncapacity and capacity models (or,
et al., 2004; Ottevanger, 2005; Rosatti and Fraccarollo, customarily, nonequilibrium and equilibrium). The nonca-
2006; WuandWang,2007].HudsonandSweby[2003]and pacity models represent the sediment in a single mode as
Castro Diaz et al. [2008] discussed 1-D bed load transport the total load. Compared to capacity models, noncapacity
models coupled withshallowwater equations byfinite vol- models treat entrainment and deposition as independent
ume methods. A few studies were found in the 2-D case. processes,thedifferencebetweenwhichinfluencesthesedi-
Hudson and Sweby [2005] and Simpson and Castelltort ment discharge and morphological evolution. The nonca-
[2006] extended the 1-D models of Hudson and Sweby pacitymodelsfacilitatethenumericalformulationsincethe
[2003] and Cao et al. [2004] to 2-D on structured grids, empirical entrainment and deposition functions can be
although the models were not tested in laboratory experi- treated as source terms. On the basisof this discussion, the
mentsorrealflowfieldsinthefield.Liuetal.[2008]devel- noncapacity model is adopted here. The conservation of
opeda 2-Dmodelof shallowwater equationsandbedload sedimentsuspendedinthewatercolumnisgivenby
transport. Delft3D is capable of modeling 2-D and 3-D
hydrodynamics andsedimenttransport usinga finitediffer- @ð hÞ @ð uhÞ @ð vhÞ
ence method (http://delftsoftware.wldelft.nl/). Neither of @t þ @x þ @y ¼E(cid:2)DþSs: ð4Þ
these models considered the effects of suspended sediment
onthehydrodynamics. [9] Amassbalanceforlocalvariationinthebedelevation
[6] Inthispaperwepresentastrategyfornumericalsolu- asafunctionofsedimentremovedoraccumulatedisgivenby
tionofthesystemoffullycoupledpartialdifferentialequa-
@z
tions for 2-D shallow water flow, sediment transport, and ð1(cid:2)pÞ ¼D(cid:2)E; ð5Þ
bed evolution. We test the code against published labora- @t
tory, field, and numerical experiments to demonstrate the wheretistime,xandyarehorizontalcoordinates,hisflow
multiscale performance of the model. The model, referred depth,uandvarehorizontalvelocities,zisbedelevation,
to as PIHM_Hydro, is based on a cell-centered upwind fi- isflux-averagedvolumetricsedimentconcentration(L3/L3),
nite volume method using Roe’s approximate Riemann g is gravitational acceleration (L/T2), p is bed sediment po-
solver on an unstructured triangular grid. A multidimen- rosity, (cid:2) is density of saturated bed, (cid:2) is density of the
z
sionallinearreconstructiontechniqueandmultidimensional water-sediment mixture, S and S are the bed slopes in
0x 0y
slopelimiter[JawaharandKamath,2000]areimplemented the x and y directions (L/L), S and S are x and y friction
fx fy
to achieve a second-order spatial accuracy. For model effi- slopes (L/L), S are sources and sinks (e.g., precipitation,
p
ciency and stability, an explicit-implicit method is used in infiltration, etc.), S are sediment sources and sinks, and D
s
temporaldiscretizationwithoperatorsplittingwhereadvec- andEaresedimentdepositionandentrainmentfluxesacross
tion and nonstiff source terms are solved via an explicit theriverbed(L/T),respectively.Thedensitiesofthewater-
scheme with the stiff source terms handled by a fully sedimentmixtureandsaturatedbedaregivenby
implicitscheme.Anumberoftestcasesoverarangeofspa-
tialscalesandhydrologicaleventsareusedtotestthemodel (cid:2)¼(cid:2)wð1(cid:2) Þþ(cid:2)s ; ð6Þ
anddemonstratethepotentialapplication.Thecodeisopen
sourceandavailablefromtheauthors. (cid:2) ¼(cid:2) pþ(cid:2)ð1(cid:2)pÞ: ð7Þ
z w s
2. Methodology [10] We note in equation (4) that diffusive transport is
treatedasnegligible[Bennett,1974].
2.1. MathematicalFormulation
[11] Equation (1) differs from the traditional mass con-
[7] Our system entails two-dimensional shallow water servation equation for shallow water flow since the right-
equations coupled with sediment mass conservation and handtermaccountsforthemorphologicalchange.Wealso
bedtopographyevolution.The2-Dshallowwaterequations note that variable fluid density in (1)–(3) allows fluvial
imply a negligible vertical velocity, hydrostatic pressure, processestocarryconcentratedsedimentflows.
2of20
W03508 LIANDDUFFY:SHALLOWWATERFLOWANDSEDIMENTTRANSPORTINRIVERS W03508
[12] Following the approach by Cao et al. [2004], the (cid:5)¼min½2;ð1(cid:2)pÞ= (cid:4): ð15Þ
systemofequations(1)–(5)canbemanipulatedsothatthe
mixture density appears on the right-hand side (i.e., in the [16] And!iscalculatedusing
sourceterms)togive
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
(cid:2) (cid:6)(cid:3)2 (cid:6)
!¼ 13:95 þ1:09gsd(cid:2)13:95 : ð16Þ
@h @ðuhÞ @ðvhÞ E(cid:2)D ð(cid:2) (cid:2)(cid:2)ÞS d d
þ þ ¼ þS þ w s s ; ð8Þ
@t @x @y 1(cid:2)p p (cid:2)
w 2.2. DomainDecomposition
@ðuhÞþ@ðu2hþgh2=2Þþ@ðuvhÞ¼(cid:2)gh(cid:4)S þS (cid:5) kn[o1w7]nQTuraialintygleunastlrguocrtiuthremd gorfidSshaerwecgheunker[a2te0d07b]y. Tthreiawngellel
@t @x @y 0x fx
allows the user to decompose the domain subject to con-
ð(cid:2) (cid:2)(cid:2) Þgh2@ ð(cid:2) (cid:2)(cid:2)ÞðE(cid:2)DÞu straints,suchasboundaries,observationpoints,and/ornested
(cid:2) s w (cid:2) z ð9Þ
2(cid:2) @x (cid:2)ð1(cid:2)pÞ multiresolution bathymetry [Qu and Duffy, 2007; Kumar
(cid:6) S S (cid:7) et al., 2008]. We have incorporated Triangle into an open
þ p(cid:2) s ð(cid:2) (cid:2)(cid:2) Þu; source GIS (PIHM_GIS [Bhatt et al., 2008]) to facilitate the
(cid:2) (cid:2) s w
w generationofunstructuredmeshesusingGISfeatureobjects.
@ðvhÞ @ðuvhÞ @ðv2hþgh2=2Þ (cid:4) (cid:5) 2.3. NumericalModel
þ þ ¼(cid:2)gh S þS
@t @x @y 0y fy [18] Thesystemofequations(8),(9),(10),(4)and(5)is
ð(cid:2) (cid:2)(cid:2) Þgh2@ ð(cid:2) (cid:2)(cid:2)ÞðE(cid:2)DÞv hyperbolic and nonlinear and subject to discontinuities
(cid:2) s w (cid:2) z ð10Þ (shocks). Extending the 1-D formulation by Fagherazzi
2(cid:2) @y (cid:2)ð1(cid:2)pÞ
andSun[2003],(5)isrewrittenas
(cid:6) S S (cid:7)
þ p(cid:2) s ð(cid:2) (cid:2)(cid:2) Þv:
(cid:2) (cid:2) s w
w @(cid:7) @ð uhÞ @ð vhÞ
þ þ ¼S ; ð17Þ
[13] In equations (9) and (10), there are two additional @t @x @y s
sourceterms.Thesecondandthirdtermsontheright-hand
with
side account for the spatial variations in sediment column
concentrationandthemomentumtransfer,becauseofsedi-
ment exchange between the water and the erodible bottom (cid:7)¼ð1(cid:2)pÞzþ h: ð18Þ
boundary.
[14] Auxiliary equations for the bottom slope are given [19] Thesystemisnowconvenientlyexpressedinvector
by S ¼@z=@x and S ¼@z=@y, and the friction slope is form:
0x 0y
estimatedbytheManningequation:
@U @E @G
n2upffiuffiffi2ffiffiffiþffiffiffiffiffivffiffi2ffiffi @t þ@xþ @y ¼S; ð19Þ
S ¼ ; ð11Þ
fx h4=3
p where Uisa vector ofthe conservativevariables, Eand G
n2v ffiuffiffi2ffiffiffiþffiffiffiffiffivffiffi2ffiffi are the flux vectors in the x and y directions, and S is the
S ¼ ; ð12Þ
fy h4=3 vectorofsourceterms.
wherenisManning’scoefficient. 0h 1 0uh 1 0vh 1
[15] Forthesedimentflux,thereexistsanextensiveliter- BuhC Bu2hþgh2C Buv C
B C B 2 C B C
ature of empirical formulae [e.g., Fagherazzi and Sun, B C B C B C
U¼Bvh C; E¼Buvh C; G¼Bv2hþgh2C;
2003; Capart and Young, 1998; Cao et al., 2004; Wu and B C B C B 2 C
B C B C B C
Wang,2007].Aparsimoniousformwhichcapturestheem- B hC B uh C B vh C
@ A @ A @ A
pirical physics for entrainment and deposition in a mini-
(cid:7) uh vh
mum of parameters is proposed by the authors based on 0 1 ð20Þ
thoseexistingworks: ðE(cid:2)DÞþS þð(cid:2)w(cid:2)(cid:2)sÞSs
B 1(cid:2)p p (cid:2)w C
E¼(cid:3)ð(cid:4)(cid:2)(cid:4)cÞhpffiuffiffi2ffiffiffiþffiffiffiffiffivffiffi2ffiffi; ð13Þ BBSx CC
B C
S¼BS C:
D¼(cid:5)! ; ð14Þ BB y CC
where(cid:3)isacalibratedconstant, (cid:4)¼u2=sgd istheShields B@E(cid:2)DþSs CA
(cid:3) S
parameter,(cid:4) isthecriticalShieldsparameterforinitiation s
c
ofsedimentmovement, (cid:5)isaparameterwhichdependson with
the distribution of the sediment in water column, ! is
the settling velocity of sediment particles in water, S ¼(cid:2)gh(cid:4)S þS (cid:5)(cid:2)ð(cid:2)s(cid:2)(cid:2)wÞgh2@ (cid:2)ð(cid:2)z(cid:2)(cid:2)ÞðE(cid:2)DÞu
u ¼pffigffiffihffiffiqffiSffiffi2ffiffiffiffiþffiffiffiffiSffiffiffi2ffiffiffi is the friction velocity, d is the sedi- x 0x fx 2(cid:2) @x (cid:2)ð1(cid:2)pÞ
(cid:3) fx fy (cid:6) S S (cid:7)
ment diameter, (cid:6) is the kinematic viscosity of water, and þ p(cid:2) s ð(cid:2) (cid:2)(cid:2) Þu;
(cid:2) (cid:2) s w
s¼(cid:2) =(cid:2) (cid:2)1. In this paper, following Cao et al. [2004], w
s w
(cid:5)issetas ð21Þ
3of20
W03508 LIANDDUFFY:SHALLOWWATERFLOWANDSEDIMENTTRANSPORTINRIVERS W03508
where U again refers to average values over the control
S ¼(cid:2)gh(cid:4)S þS (cid:5)(cid:2)ð(cid:2)s(cid:2)(cid:2)wÞgh2@ (cid:2)ð(cid:2)z(cid:2)(cid:2)ÞðE(cid:2)DÞv volume Vi , S ¼1=V R SdV is the numerical approxima-
y 0y fy 2(cid:2) @y (cid:2)ð1(cid:2)pÞ tion of thei soiurce termi ,Vni is the unit outward normal vec-
þ(cid:6) Sp(cid:2)Ss(cid:7)ð(cid:2) (cid:2)(cid:2) Þv: tor to the edge j, (cid:2)j is theijn length of edge j, and Fij is the
(cid:2) (cid:2) s w numerical flux vector through the edge j, which is calcu-
w
ð22Þ lated using an approximate Riemann solver [LeVeque,
2002]. Details of the flux evaluation using an approximate
[20] The source term vector now has five parts: bed Riemann solver are found in Appendix A, including a dis-
slope S , friction slope S, sediment concentration varia- cussion of the strategy for handling discontinuities and
0 f
tions S , sediment exchange S , and the additional source- Roe’sformulationandthenumericalfluxevaluationatele-
c e
sinktermS : mentboundaries[Sleighetal.,1998].
p
2.4. LinearReconstructionandMultidimensional
S ¼ð0 (cid:2)ghS (cid:2)ghS 0 0ÞT; ð23Þ SlopeLimiter
0 0x 0y
[24] In PIHM_Hydro we attempt to preserve high spa-
S ¼ð0 (cid:2)ghS (cid:2)ghS 0 0ÞT; ð24Þ tial accuracy in the flow simulation by adopting a second-
f fx fy
order piecewise linear reconstruction. Many second-order
numerical schemes have been implemented for shallow
Sc ¼(cid:2)0 (cid:2)ð(cid:2)s(cid:2)2(cid:2)(cid:2)wÞgh2@@ x (cid:2)ð(cid:2)s(cid:2)2(cid:2)(cid:2)wÞgh2@@ y 0 0(cid:3)T; ð25Þ water equations on unstructured triangular grids [e.g.,
AnastasiouandChan,1997; Sleighetal.,1998; Hubbard,
1999; Wang and Liu, 2000; Yoon and Kang, 2004] and
S ¼(cid:2)E(cid:2)D (cid:2)ð(cid:2)z(cid:2)(cid:2)ÞðE(cid:2)DÞu (cid:2)ð(cid:2)z(cid:2)(cid:2)ÞðE(cid:2)DÞv E(cid:2)D 0(cid:3)T; other types of grids as well [e.g., Alcrudo and Garcia-
e 1(cid:2)p (cid:2)ð1(cid:2)pÞ (cid:2)ð1(cid:2)pÞ Navarro, 1993; Ambrosi, 1995; Valiani et al., 2002; Cal-
ð26Þ effi et al., 2003]. Extension of structured techniques to
unstructured grids such as the MUSCL approach have
(cid:2) ð(cid:2) (cid:2)(cid:2)ÞS(cid:6) S S (cid:7) achieved only partial success because of the pronounced
Sp¼ Spþ w (cid:2) s s (cid:2)p(cid:2)(cid:2)s ð(cid:2)s(cid:2)(cid:2)wÞu grid sensitivity [Jawahar and Kamath, 2000], with poor
w w
ð27Þ results typically obtained on highly distorted grids. The
(cid:6) S S (cid:7) (cid:3)T
(cid:5) p(cid:2) s ð(cid:2) (cid:2)(cid:2) Þv S S : reconstruction technique proposed by Jawahar and
(cid:2) (cid:2)w s w s s Kamath[2000] was successfullyappliedwith the Harten–
Lax–van Leer (HLL) approximate Riemann solver for
[21] Itisnowconvenienttowritethesystemas shallow water equations [Yoon and Kang, 2004]. Com-
pared to other multidimensional linear reconstruction
@U
þr(cid:5)F¼S; ð28Þ techniques, this one uses a wide computational stencil and
@t
does not strongly depend on vertex values. It has also
where F ¼ (E, G)T. The system is then integrated over an beenshownthathigh-orderschemesmayleadtononphys-
arbitrarycontrolvolumeV (atriangularelementhere): ical oscillatorysolutions near discontinuities [Toro,1999].
i
To avoid oscillations, we limit the solution slope during
Z @U Z Z the linear reconstruction. The multidimensional slope lim-
dV þ r(cid:5)FdV ¼ SdV: ð29Þ iter proposed by Jawahar and Kamath [2000] has the
@t
Vi Vi Vi advantages that (1) the limiter is inherently multidimen-
sional which fits unstructured grids and (2) it is continu-
[22] Application of the Gauss theorem leads to the inte-
ouslydifferentiable.
gralformof(29):
2.5. SourceTerms
Z @UdV þI F(cid:5)nd(cid:2)¼Z SdV; ð30Þ [25] It was found that the source terms in (23)–(27)
@t require careful treatment. The friction slope S and sedi-
Vi (cid:2)i Vi f
ment exchange S were discretized in a pointwise manner
e
where(cid:2)istheboundaryofthecontrolvolumeand n¼(n and evaluated at the element centroid. The sediment con-
x
n )Tistheunitoutwardvectornormaltotheboundary. centration S terms are also discretized at the centroid, and
y c
[23] Acell-centeredfinitevolumeisusedtoapproximate the linear reconstruction procedure readily provides the
(30) by vertically projecting a Delaunay triangle to form gradientofsedimentconcentration ð Þ:
theprismaticcontrolvolume(Figure1).Thestatevariables
of the system are stored at the centroid of the control vol- rðh Þ (cid:2) rh
ume and represented as piecewise constant over the do- r i ¼ ih i : ð32Þ
main. We note that by solving the system for the centroid
ofeachfinitevolumeenablestheimplementationofahigh-
2.6. TimeIntegration
order interpolation scheme [Sleigh et al., 1998], discussed
insection2.4.Finally,(30)isrewrittenas [26] A critical problem in our coupled system is the nu-
merical instabilities related to friction slope and sediment
@@Uti ¼(cid:2)V1 X3 (cid:4)Fij(cid:5)nij(cid:5)(cid:2)jþSi; ð31Þ euxsecdh,anwgheerfeortheshsayllsotewmd(e3p1t)hiss.sAplisteimntio-itmwpoliscyitstemmesthoofdoirs-
i j¼1 dinarydifferentialequations:
4of20
W03508 LIANDDUFFY:SHALLOWWATERFLOWANDSEDIMENTTRANSPORTINRIVERS W03508
Figure 1. (a) Plan view, (b) computational mesh, and (c) close-up view of the near-dam area of the
converging-divergingchannelinBellosetal.’s[1992]experiment.
@U 1 X3 (cid:2) (cid:3) 2.7. BoundaryConditions
i ¼(cid:2) F (cid:5)n (cid:2) þS þS þS ; ð33Þ
@t Vi j¼1 ij ij j 0i ci pi [30] Open andsolidwallboundaryconditionshave been
implemented.Theleft-hand(slip)solidwallboundarycon-
ditionisgivenby
@U
@ti ¼SfiþSei: ð34Þ 0h(cid:3) 1 0hL 1
BBu(cid:3)(cid:5)nCC BB(cid:2)uL(cid:5)nCC
[27] In the first step, advection and source terms for bed BBu(cid:3)(cid:5)t CC¼BBuL(cid:5)t CC; ð37Þ
slope and sediment concentration are solved using an B@ (cid:3) CA B@ L CA
explicitmethod.Then,thevaluesobtainedfromthefirststep
z z
define initial conditions for the remaining equations, which (cid:3) L
use an implicit method (backward differentiation formula). where u denotes (u, v)T and the subscripts L and asterisk
PIHM_Hydro uses the advanced ODE solver CVODE arethevariablesattheleftsideandboundary,respectively.
[HindmarshandSerban,2005]. Thevelocitycomponentscanbecalculatedby
[28] Explicit time integration is performed by the first-
orderEulermethodoratotalvariationdiminishingRunge- (cid:6)u (cid:7) (cid:6)n (cid:2)n (cid:7)(cid:6)u (cid:5)n(cid:7)
(cid:3) ¼ x y L : ð38Þ
Kuttamethod[ShuandOsher,1988]givenby v n n u (cid:5)t
(cid:3) y x L
U1¼Unþ(cid:3)tfðUnÞ; [31] The open boundary conditions are more compli-
3 1 1 cated. A simple one is the free outfall condition, where the
U2¼4Unþ4U1þ4(cid:3)tfðU1Þ; ð35Þ waves pass the boundary without reflection. It can be
described as U ¼ U . For other cases, it is sometimes
1 1 2 (cid:3) L
Unþ1¼3Unþ3U2þ3(cid:3)tfðU2Þ; found that the physical boundary condition(s) is not suffi-
cient,andthetheoryofcharacteristicsisusedtoderivesuf-
wherefistheright-handsideof(33).Themethodhasbeen ficient information at the boundaries. For more detailed
shown to improve stability and preserves high-order accu- discussion about the open boundary conditions, refer to
racy(third)[ShuandOsher,1988]. Zhaoetal.[1994],AnastasiouandChan[1997],andSleigh
[29] It is well known that the explicit scheme has a sta- etal.[1998].
bilityrestrictionontheCourant-Friedrichs-Lewycondition. [32] A very small flow depth commonly occurs near the
Anadaptive(cid:3)tisusedinthemodelusingthefollowing: wetting-drying boundaries, which introduces numerical
instabilities with unreliably high velocities that make the
minðdÞ
(cid:3)t(cid:6) h(cid:4)pffiffiffiffiffiffiffiffiffiiffiffiffiffiffiffi (cid:5)i ; ð36Þ friction slope and sediment source terms (equations (24)
2max u2þv2þc and(26))verystiff,asmentioned.Tosolvethisproblem,in
i
addition to employing the semiexplicit time integration
where i is the cell index and d represents the whole set of scheme,a simplewetting-dryingalgorithmwasalsoimple-
i
distances between the ith centroid and those of its neigh- mented. In this approach, a tolerance depth (h ) was cho-
tol
boringcells. sen in that the flow velocities were set to zero when the
5of20
W03508 LIANDDUFFY:SHALLOWWATERFLOWANDSEDIMENTTRANSPORTINRIVERS W03508
flowdepthwaslessthanh .Ineachexampleh wascho- 3.2. Case2:Two-DimensionalRainfall-Runoff
tol tol
senempirically. NumericalExperimentWithWetting-DryingSurface
[37] A rainfall-runoff numerical experiment by diGiam-
marco et al. [1996] was used to compare simulations of
3. FieldandLaboratoryImplementation
rainfall-driven flow with rapid wetting-drying of the sur-
[33] The first objective of this research was to develop a face. In this example, a tilted V-shaped catchment (Figure
robust fully coupled numerical code suitable for applica-
3) is generated by introducing a single, 90 min duration,
tions over a wide range of spatial and temporal scales and
10.8 mm/h intensity rainfall event. The catchment is com-
riverineprocessesforflow,sediment,andbedevolution.A
posedoftwo1000m(cid:5) 800mplanesconnectedbya1000
second objective was to bring together experimental data
m long (cid:5) 20 m wide channel. Infiltration is zero over the
sets to test and compare PIHM_HYDRO over a useful
entire domain. In this case the channel depth is set to zero
range of spatial and temporal scales and also to compare
which means the planes and channel are continuously con-
the model with other published models. For the laboratory
nected so that a backwater effect could be considered. The
scale, we found three examples that represent small-scale
bed slopes are 0.05 and 0.02, perpendicular to and parallel
hydrodynamics and sediment transport, and two examples
to the channel, respectively. Manning roughness coeffi-
were developed for the field-scale application. Each of
cients are0.015 s/m1/3 fortheplane and0.15 s/m1/3forthe
the four experimental results is available at http://www.
channel. The simulation was run for 180 min with a free
pihm.psu.edu for the purpose of future experimental and
boundary condition at the channel outlet and no-flow
modeltesting.
boundary conditionselsewhere. Initially, the water depth is
zeroovertheentiredomain.
3.1. Two-DimensionalLaboratory-ScaleDamBreak
WithaDry-WetTransitionandConverging-Diverging [38] Thesimulationused740trianglesandwascompared
withintegratedfinitedifference(IFD)[diGiammarcoetal.,
Channel
1996], MIKE-SHE [Abbot et al., 1986], IHM [VanderK-
[34] Bellosetal.[1992]performedsimulationsofinstan- waak, 1999], and MODHMS [Panday and Huyakorn,
taneous dam failure for a range of initial conditions. This
2004], all of which use either the kinematic or diffusive
experimentisusedtotestPIHM_Hydroforsimulatingflow
approximations. In these models, IFD, MIKE-SHE, and
dynamics in an irregular flow domain with an initially dry
MODHMS weakly couple the 2-D flow on the plane with
bed with nonzero bed slope and roughness. The 21.2 m (cid:5)
the1-Dflowinthechannel,whileIHMfullycouplesthe2-
1.4mwideexperimentalflumehadarectangular,converg-
D flows on the plane and in the channel. The hydrograph
ing-diverging cross section (Figure 4). A movable gate
predicted by PIHM_Hydro is compared with the multimo-
(x ¼ 0 m) was used to simulate the instantaneous dam
del average in Figure 4a. The multimodel average is the
break. Eight probes were installed along the centerline of
meanvalueofallmodeloutputs.PIHM_Hydroisalsocom-
thechanneltomeasuretheflowdepths.
pared individually with IFD, MIKE-SHE, IHM, and
[35] The initial water level is 0.15 m upstream of the MODHMS in Figure 4a. The simulation data from these
damandzerodownstreamfromthedamwiththebedslope
models were obtained by digitizing the figures from the
of 0.002. A solid wall boundary condition was applied at
publications,whichintroducedsomeerror.
the upstream end and sidewalls, and a free out-fall condi-
[39] Excellent agreement is illustrated between PIHM_
tion was applied at the downstream boundary. The domain
Hydroandthemultimodelaverage(usingNash-Sutcliffcri-
was discretized into 3886 triangles (Figure 1). The Nash-
teria) as with the other models, especially for the peak dis-
Sutcliffemodelefficiencycoefficient(NSE)[NashandSut-
charge and receding limb. Although there is a slight
cliffe, 1970] was used to quantitatively evaluate the accu-
discrepancyintherisinglimbforthesemodels,thetimesto
racy of the model for observed and predicted quantities.
peak discharge are almost identical. According to Figure
The predicted and measured flow depths at the different
4a, PIHM_Hydro predicts a constant and stable plateau of
positionsareshowninFigure2. thepeakdischarge(4.86m3/s),whichisanalyticallycorrect
[36] Upstream of the dam the predictions match the andshowsthemodeltobestable.Amongthesesimulations,
measurements very well with minimal error. At down-
PIHM_Hydro had the closest agreement with IHM, fol-
streamlocationswherecriticaltosupercriticalflowandthe
lowed by MODHMS, MIKE-SHE, and IFD. respectively,
wetting-dryingprocessisobserved,theflowdepthsarestill
again using Nash-Sutcliff criteria. One possibility for the
well reproduced. However, it was found that downstream
similarity is that only PIHM_Hydro and IHM fully couple
locations (x ¼ þ5.0 m) were quite sensitive to Manning’s
the 2-D overland and channel flow to consider the back-
roughness coefficient, possibly because of depth averaging
water effect. Figure 4b shows the accumulative mass bal-
and lack of a vertical velocity term in the transition from
anceerrorscalculatedusingtheequation
subcritical to supercritical flow [Martin and Gorelick,
2005]. Figure 2 reveals that the arrival times of the wave- ðprecipitation(cid:2)discharge(cid:2)storageÞ
err ¼ : ð39Þ
front are accurately simulated, and for this small-scale test h precipitation
casethemodelwasfoundtobebothaccurateandstablefor
the wetting-drying process as well as for supercritical [40] Themassbalanceerrorisverysmall,withthemaxi-
flows.Themodelisabletosimulatetheearlytimeandlate mumvaluelessthan0.8%at t¼45minandalmostzeroat
timebehaviorofthefloodwavesaswellasthesharppropa- t ¼ 180 min. For this benchmark case, PIHM-Hydro per-
gationofthedryingfront.Sedimenttransportwasnotcon- formanceisstable andmassconservativeforcasesofrain-
sidered in the original experiment, and other model fall-driven overland channel flow and intensive wetting-
simulationswerenotfoundforthisexperiment. dryingprocesses.
6of20
W03508 LIANDDUFFY:SHALLOWWATERFLOWANDSEDIMENTTRANSPORTINRIVERS W03508
Figure 2. Two-dimensionallaboratory-scaledambreakwithdry-wetandconverging-divergingchannel:
measured(circles)andpredicted(solidline)temporalvariationsofflowdepthatdifferentlocationsupstream
ofthedam,atthedam,anddownstreamfromthedam.Notalllocationsshownhaveobserveddata.
7of20
W03508 LIANDDUFFY:SHALLOWWATERFLOWANDSEDIMENTTRANSPORTINRIVERS W03508
Figure 3. Two-dimensional rainfall-runoff dry to wet numerical experiment: tilted V-shaped catch-
ment (not to scale). The red lines represent the no-flow boundaries. The blue line denotes the channel
outlet, where the free outfall boundary is applied [diGiammarco et al., 1996]. The black lines are the
break in slope between side slope planes and the channel. The side slope is 0.05, and the bed slope is
0.02.
3.3. Case3:Two-DimensionalLaboratory-ScaleFlow model reproduced the flow dynamics over movable beds
andSedimentTransportFollowingDamBreak verywell.Althoughthepredictionofbedsurfaceelevations
isnotasgoodasthewatersurface,itcanbeconsideredsat-
[41] In this test case, model flow dynamics with fully
isfactory considering the complexity of the bed evolution,
coupled sediment transport is implemented for two labora-
the small timescale, the small magnitude of the bed eleva-
tory-scale dam break experiments with flow over movable
tion change, and the measurement uncertainty. According
beds. These are described in the literature as the Taipei
to Figure 6, the model predictions of the wavefront loca-
experiment [Capart and Young, 1998] and the Louvain
tions and the erosion magnitude were in fairly good agree-
experiment [Fraccarollo and Capart, 2002]. In the Taipei
ment withthe measurement in the experiment.A hydraulic
experiment,sedimentparticleswereartificialsphericalbeads
jumpformedneartheinitialdamsitesbecauseofrapidbed
of uniform size diameter of 6.1 mm, specific gravity of
erosion. In the Louvain experiment, the location of the hy-
1.048, and settling velocity of 7.6 cm/s. In the Louvain
draulic jump was well predicted by the model. It can be
experiment,thesedimentparticleswerereplacedbycylindri-
observed from both the predictions and measurements that
cal PVC pellets with equivalent spherical diameter of 3.5
the hydraulic jump propagated upstream in the Louvain
mm,specificgravityof1.54,andsettlingvelocityof18cm/s.
experiment.
Horizontal prismatic flumes were used in both experiments.
Thetestreachwas1.2mlong,70cmhigh,and20cmwide [44] In the Taipei experiment, the measured data at 3t0
were used to calibrate themodel. Itis evident that thecali-
fortheTaipeiexperimentand2.5mlong,25cmhigh,and10
brationisnotasgoodasthatintheLouvainexperiment.By
cm wide for the Louvain experiment. The sluice gates were
using the calibrated parameters, Figure 7 shows that the
placedinthe middleoftheflumesinbothcasestorepresent
model predicted the wavefront location and the erosion
the dam break. The initial water depth was 10 cm upstream
magnitudequitewell.However,theagreementbetweenthe
and 0 cm downstream for both experiments. The solid wall
predictionandmeasurementwasnotasgoodforthemagni-
boundaryconditionwasappliedtoalltheboundaries.
tude and the locations of the hydraulic jump (Figure 10).
[42] The 2-D computational meshes are shown in Figure
The hydraulic jump moved upstream in the prediction but
5, with 745 and 764 triangles for the Taipei and Louvain
remainedstationaryintheexperiment.Thecauseneedsfur-
experiments,respectively.Inthismodel,severalparameters
therinvestigation.
need to be determined. The sediment porosity (p) and
Manning’sroughnesscoefficient(n)wereobtainedfromthe [45] The mass balance errors of sediment transport were
calculatedusing
literature[WuandWang,2007].Herepwassetto0.28and
0.3 for the Taipei and Louvain experiments, respectively,
while n was set as 0.025 s/m1/3 for both experiments. For Pð(cid:3)zÞi(cid:5)Aið1(cid:2)pÞþP(cid:3)ðhcÞi(cid:5)Ati¼5t0
the parameters (cid:3) and (cid:4)c, calibration was done using meas- errz¼ i Pð(cid:3)zÞ (cid:5)iA ; ð40Þ
ured data at 5t and 3t (t ¼ 0.101 s) for the Taipei and i i
0 0 0 i
Louvain experiments, respectively, and the calibrated pa-
rameters were then used to predict the flow and sediment where ð(cid:3)zÞ and (cid:3)ðhcÞ are the change of bed elevation
i i
dynamicsfortheothertwotimes.Parameters(cid:3)and(cid:4) were and sediment load per area at each grid i. In the Louvain
c
calibrated to 2.2 and 0.15 for the Taipei example and 5.0 experiment, the err values are 0.05%, 0.007%, and 0.05%
z
and0.05fortheLouvaincase. at5t ,7t ,and10t ,respectively.IntheTaipeiexperiment,
0 0 0
[43] In Figures 6 and 7, the predicted water (H) and bed the errz values are 0.02%, 0.005%, and 0.001% at 3t0, 4t0,
(z)surfaceelevationsalongthecenterlineofthechannelare and5t ,respectively.Massbalanceerrorsareverysmallin
0
compared with the measured values. Nash-Sutcliff criteria eachcase.
wereusedtocomparemodeledandobservedsurfaceeleva- [46] Figure8showsthepredictedsedimentconcentration
tions for the two experiments. In the Louvain experiment, profiles in both experiments at different times. The model
the measured data at 5t were used to calibrate the model predictssharpforefrontsofthesedimentconcentrationpro-
0
forthesameunknownparametersasbefore.Ingeneral,the file, which were caused by extremely high sediment
8of20
W03508 LIANDDUFFY:SHALLOWWATERFLOWANDSEDIMENTTRANSPORTINRIVERS W03508
Figure 4. Evaluation of model performance for the 2-D rainfall-runoff dry to wet numerical experi-
ment. (a) Comparison of simulated hydrograph with MODHMS [Panday and Huyakorn, 2004], IHM
[VanderKwaak, 1999], MIKE-SHE [Abbot et al., 1986], and IFD [diGiammarco et al., 1996] and (b)
massbalanceerrors.
concentrationsinthedambreakwavefronts.Similarpredic- settling velocity of sediment particles in the Taipei case
tions were also made by several 1-D models [Fraccarollo which make erosion easier and deposition slower. It is also
andCapart, 2002; WuandWang, 2007] in theTaipei case noted that the peak sediment concentration decreased with
at4t .IntheWuandWangmodel,however,asmallinitial timeintheLouvaincasewhileitincreasedwithtimeinthe
0
downstreamflowdepthwasspecifiedratherthantheactual Taipei case during their time scales. This difference is
dry bed condition, which may have caused errors in mass partly related to the different time scales of the two cases,
balance and wave arrival time. By comparing Figures 8a where the Taipei experiment is run for a shorter duration
and8b,itisnotedthatthepeakconcentrationsintheTaipei thantheLouvainexperiment.
caseweremuchhigherthanthoseintheLouvaincase.The [47] Figure 9 illustrates the temporal variations of flow
difference is because of the smaller density and lower depths, flow discharge, and sediment discharge predicted
9of20
W03508 LIANDDUFFY:SHALLOWWATERFLOWANDSEDIMENTTRANSPORTINRIVERS W03508
Figure 5. Plan view of the computational mesh for (a) the Louvain experiment and (b) the Taipei
experiment.Thethickverticallinesdenotethedamineachcase.
bythismodelatx¼0.03,0.2,and0.4m.Weobservedthat s/m1/3. A uniform Manning coefficient of 0.033 s/m1/3
the flow depth and flow and sediment discharges approach was advised by other researchers [e.g., Goutal, 1999; Her-
equilibrium values after a period of time. In Figures 8 and vouet, 2000; Valiani et al., 2002]. The sensitivity of the
9 it is also found that volumetric sediment concentration, model to the Manning’s coefficient was first studied on
flow depth, and sediment discharge were very sensitive a 26,000-element triangle mesh (Table 1). It shows
to particle size and density while discharge was not. Like- that PIHM_Hydro does work best at n ¼ 0.033 s/m1/3,
wise, sediment discharge is also affected by the sediment which is in agreement with the previous studies as well
particle properties indirectly. Overall erosion and deposi- as the recommendation for winding natural rivers by
tionalterthebedsurface,whichinturninfluencestheflow Henderson[1966].
depthprofile. [50] Table 2shows how thegrid resolution or terrain re-
solution influences the numerical solutions. Because of the
3.4. Case4:RiverandFloodplainDynamicsDuring largemagnitudeofwatersurfaceelevations, bothNSEand
MalpassetDamBreakEvent root-mean-square error (RMSE) are used here to quantita-
[48] Case4involvesanactualdambreakandfloodwave tivelyassessthepredictionerrors.Itisevidentthatthefiner
propagating over a dry bed with complex topography and spatial resolutions produce better solutions, up to the point
boundaries. The Malpasset dam was located in a narrow when the grid resolution reaches (cid:7)38,208 in this case, and
gorge of the Reyran river valley, about 12 km upstream of accuracydoesnotimprovemuchevenifthegridsgetmuch
FrejusinsouthernFrance.Themaximumreservoircapacity finer. By comparing the two metrics at the gauge stations
was 55,106 m3. In December 1959, the dam failed explo- andthepolicesurveypointsinTable2,wecanseethatthe
sively at night, partly because of exceptionally heavy rain. maximum water surface elevations (H) at the gauge points
The flood wave propagated along the Reyran valley to the are more sensitive to the grid and terrain resolutions than
city of Frejus, and 433 casualties were reported. Details of those at the police survey points. This is because the mag-
thiscatastrophearegivenbySoaresFraz~aoetal.[1999]. nitudes of the water depths (h) in the river are larger than
[49] The topography of the computation domain is the bed elevation (z), while along the banks the bed eleva-
shown in Figure 10 with the location of dams, transform- tionsarelarger. Themaximumwater surfaceelevation isa
ers,gauges,andpolicesurveypoints[SoaresFraz~aoetal., reasonable metric to assess the model accuracy. Tables 1
1999], with the solution domain encroaching on the shal- and 2 illustrate that when the number of grids reaches
lowtidalestuary.Thebottomelevationrangesfrom (cid:2)20m 26,000, the simulation becomes satisfactory for all mea-
below sea level to 100 m above sea level. The domain surementlocationsaccordingtoourNSEcriterionandthat
was discretized into different resolution meshes using the increasing the resolution to 38,208 grids provides some
domainboundaryand the13,541 topography survey points improvement, while finer-resolution grids do not signifi-
as the constraints. The 38,208-triangle mesh and a close- cantlyimproveresults.
up view of the computational mesh near the dam are pre- [51] Toexaminetheconvergenceofthemodelwithsuc-
sented in Figure 11. The mesh is much finer immediately cessive grid refinement, we evaluated Roache’s grid con-
downstreamfromthedamandalongtheriver,whererapid vergence index [Roache, 1994] at each gauge point. From
andabruptflowoccurred.Theinitialwaterlevelintheres- thisanalysisa38,208-gridmeshwithn¼0.033s/m1/3was
ervoir is 100 m above sea level and the channel bottom found to provide the best result and was used to compare
downstream from the dam is set dry in the model. The PIHM_Hydro with other models of the Malpasset dam
solid wall condition is imposed along all boundaries. The break.
Strickler coefficient ranges from 30 to 40 m1/3/s1 [e.g., [52] The predicted evolution of flood inundation is
SoaresFraz~aoetal.,1999; HervouetandPetitjean,1999], shown with the police-surveyed points in Figure 12. The
correspondingtoManning’scoefficientfrom0.033to0.025 cellswithawaterdepthsmallerthan0.01mareconsidered
10of20