A regularization strategy for modeling mixed-sediment river morphodynamics

(1)

Delft University of Technology

A regularization strategy for modeling mixed-sediment river morphodynamics

Chavarrías, Víctor; Stecca, Guglielmo; Siviglia, Annunziato; Blom, Astrid

DOI

10.1016/j.advwatres.2019.04.001

Publication date

2019

Document Version

Final published version

Published in

Advances in Water Resources

Citation (APA)

Chavarrías, V., Stecca, G., Siviglia, A., & Blom, A. (2019). A regularization strategy for modeling

mixed-sediment river morphodynamics. Advances in Water Resources, 127, 291-309.

https://doi.org/10.1016/j.advwatres.2019.04.001

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

(2)

ContentslistsavailableatScienceDirect

Advances

in

Water

Resources

journalhomepage:www.elsevier.com/locate/advwatres

A

regularization

strategy

for

modeling

mixed-sediment

river

morphodynamics

Víctor

Chavarrías

a,∗

_,

_Guglielmo

_Stecca

b

_,

_Annunziato

_Siviglia

c

_,

_Astrid

_Blom

a a Faculty of Civil Engineering and Geosciences, Delft University of Technology, Delft, The Netherlands

b National Institute of Water and Atmospheric Research (NIWA), Christchurch, New Zealand

c Laboratory of Hydraulics, Hydrology and Glaciology, Swiss Federal Institute of Technology (ETH), Zurich, Switzerland

a b s t r a c t

Anotabledrawbackinmixed-sizesedimentmorphodynamicmodelingisthefactthatthemostcommonlyusedmathematicalmodelinthisfield(i.e.,theactivelayer modelHirano,1971)canbeill-posedundercertaincircumstances.Undertheseconditionsthemodellosesitspredictivecapabilities,asnegligibleperturbationsin theinitialorboundaryconditionsproducesignificantdifferencesinthesolution.Inthispaperweproposeapreconditioningmethodthatregularizesthemodelto recoverwell-posednessbyalteringthetimescaleofthesedimentmixingprocesses.Wecompareresultsoftheregularizedmodeltodatafromfournewlaboratory experimentsconductedunderconditionsinwhichtheactivelayermodelisill-posed.Theregularizedactivelayermodelcapturesthechangeofbedelevationand surfacetextureaveragedoverthepassageofseveralbedforms.Neithertheactivelayermodelnortheregularizedoneaccountforsmallscalechangesduetoindividual bedforms.

1. Introduction

Thepresenceofmixed-sizesedimentisakeyfeatureofrivers. Sed-iment sortingpatternsdevelopin thestreamwisedirection (e.g.,the characteristicdownstreamfiningprofileSternberg,1875),inthe trans-versedirection(e.g.,bendsortingAllen,1970),andinthevertical di-rection(e.g.,bedarmoringParkerandKlingeman,1982anddune sort-ingBlometal.,2003).Modelingapplicationsinwhichthemixed-size characterofrivermorphodynamicsisnotnegligiblemandatetheuseof asuitablecontinuitymodelaccountingformassconservationofeach oftheconsideredsedimentsizefractions.Hirano(1971)wasthefirst todevelopa mixed-sedimentcontinuitymodel. Heassumed thatthe river bed can be vertically dividedinto an active top part(the ac-tivelayer),which interacts withtheflow,andaninactivesubstrate. Inthemodel,sedimenttransport andfrictiondependonthetexture oftheactivelayer,whereas thesedimentin thesubstrateonlyplays aroleifnetaggradationcreatesnewsubstratesedimentornet degra-dationleadstotheentrainmentofsubstratesedimentintotheactive layer.

Althoughithasbeenfruitfullyusedtorepresentphysicalphenomena relatedtomixed-sedimentfornearlyhalfacentury(seeChavarríasetal., 2018), the activelayer model suﬀers from a drawback. Under cer-tainconditionsitbecomesill-posed(Chavarríasetal.,2018;Ribberink, 1987; Steccaet al., 2014).Amodel isill-posed ifaunique solution doesnotexist,orifthesolutiondoesnotdependcontinuouslyonthe initialandboundaryconditions(Hadamard,1923).Ifa modelis ill-posed, inﬁnitesimal variations in theinitial or boundary conditions

∗_{Corresponding}_author.

E-mailaddress:v.chavarriasborras@tudelft.nl(V.Chavarrías).

yieldasigniﬁcantdeviationofthesolutionwithinaninﬁnitesimaltime (Hadamard,1923).Whensolvingthemathematicalmodelby numeri-calapproximations,perturbationsintheinitialandboundaryconditions simplyarisebytruncationerrors.Thismakesanill-posedmodel unsuit-ableinpractice.

Theproblemofill-posednessarisesfromaninaccurate representa-tionof thephysical processes(JosephandSaut,1990).Forinstance, atwo-fluid modelforincompressibleandinviscidflow intwolayers withavelocitydiscontinuityisill-posed(vonHelmholtz,1868;Kelvin, 1871).Itisregularized(i.e.,becomeswell-posed)ifviscouseffectsare takenintoaccount(JosephandSaut,1990).Fromthisperspective,the preferredapproachtoregularizetheactivelayermodelwouldbe the developmentofanewmodelthatincludesthosephysicalmechanisms thatarenotaccountedforbytheactivelayermodel.

There exist alternatives to the active layer that typically aim to improve the physical description of sediment mixing process.

Ribberink(1987)introducedasecondlayertoaccountforthemixing duetodunesexceptionallylargerthantheaverageduneheight.Besides producingaverticalsortingproﬁlethatbetterreproducestheresults ofalaboratoryexperiment(Blom,2008),Ribberink’stwo-layermodel makestheoccurrenceofill-posednesslesslikely,althoughitdoesnot completelyavoidit(Sieben,1994).Luuetal.(2006,2004)proposeda modelinwhichtheactivelayerisreplacedbythesedimenttransport layerrepresentingthesedimentintransportratherthanthesediment atthebedsurface.Thethicknessofthesedimenttransportlayeris esti-matedwithaclosurerelationsuchastheonedevelopedbyEgashira and Ashida (1992). Although conceptually diﬀerent, the model by

https://doi.org/10.1016/j.advwatres.2019.04.001

Received12June2018;Receivedinrevisedform20February2019;Accepted2April2019 Availableonline5April2019

(3)

Luuetal.(2006,2004)ismathematicallyequivalenttotheactivelayer model,whichimpliesthatitcanalsobeill-posed.

BlomandParker(2004)andBlometal.(2006,2008)developeda modelinwhichbothbedelevationandbedgrainsizedistributionare treatedusingaverticallycontinuousformulation(Parkeretal.,2000). Thisimpliesthatthereisnodistinctionbetweentheactiveand inac-tivepartofthebed.ThemodelbyBlomandcoauthorssatisfactorily describestheverticalstratigraphyduetodunesatlaboratoryscale,but itrequiresatimesteptoosmalltobeapplicableatlargescale. More-over,itswell-posednesshasnotbeenstudied.Simplifyingthe contin-uousframeworkproposedbyParkeretal.(2000),thevertically con-tinuousmodelbyViparellietal.(2017)overcomestheneedforasmall timestep.Althoughapplicableatlargespatialandtemporalscales,their modeldoesnotsolvetheproblemofill-posedness(Chavarríasetal., 2018).

Giventhefacts that:(a)There isnotyeta practicallyfeasible al-ternativetotheactivelayermodel,(b)theactivelayermodelremains well-posedoveralargerangeofapplications(Chavarríasetal.,2018), and(c)itisacomputationallycheapmodelimplementedinseveral soft-warepackages,hereourobjectiveistodevelopastrategytoavoid ill-posednesswhilemaintainingtheconceptualframeworkoftheactive layermodel.Tothisend,wedeveloparegularizationstrategythat re-coverswell-posednessoftheactivelayermodelandweconduct4 labo-ratoryexperimentsunderconditionsinwhichtheactivelayermodelis ill-posedtoobtaindatatowhichwecomparetheresultsofour regular-izedmodel.

Thepaperisorganized asfollows.InSection2we review strate-giesforregularizingill-posedmodels.InSection3wepresentthe reg-ularizationstrategy.Section4presentsthelaboratoryexperimentsand

Section5focusesonthenumericalrunstoreproducethe experimen-talresults.InSection6wediscussthelimitationsoftheregularization strategy,aswellasotherpossiblemodelingstrategies.

2. Overviewofregularizationtechniques

Inthis section we review techniquesused toregularize ill-posed problems.Propagationproblemsaremost completelymathematically representedby a setof partialdiﬀerential equationsconstituting an initialvalueproblem. Intheseproblemsaninitialstatechangeswith timesubjecttoconditionsattheboundariesofthedomain.The matrix-vectorformulationprovidesacompactexpressionofthesetof equa-tions(e.g.,CourantandHilbert,1989;LynandGoodwin,1987;Toro, 2001):

𝜕𝐐

𝜕𝑡 +𝐀𝜕𝐐𝜕𝑥 =𝐒, (1)

whereQisthevectorofdependentvariables,Aisthesystemmatrix,and

Sisthevectorofsourceterms.Thevelocityatwhichsmallwaves prop-agatethroughoutthedomain(i.e.,theeigenvaluesofmatrixA)must berealfortheproblemtobewell-posed(e.g.,Hadamard,1923;Ivrii andPetkov,1974;Kabanikhin,2008;Lax,1957;1958;1980;Mizohata, 1961).Whentheeigenvaluesarereal,theproblemishyperbolic.Ifthe eigenvalueshaveanimaginarycomponent(theproblembeingelliptic orofmixed-type),aninitialvalueproblemisill-posed.

Thetwo-fluidshallowflowmodel(i.e.,amodeloftheflowoftwo layers of superimposedfluids at different velocities) is knowntobe ill-posedwhenthedifferenceinflowvelocitybetween theupperand lowerlayers exceedsa certainthreshold(Ardron,1980; Armi,1986; Lawrence,1990;Long,1956;Pelantietal.,2008).Ingeneralterms ill-posednessarisesinmultiphasemodels(e.g.,bubblesinafluid)(Harlow andAmsden,1975;KumbaroandNdjinga,2011;Murray,1965; Stew-art,1979;StewartandWendroff,1984).Multiphasemodelsare regu-larizedbyaccountingfortheforcesattheinterfacebetweenthetwo fluids(AbgrallandKarni,2009;Drewetal.,1979;Liskaetal.,1995; Lyczkowskiet al., 1978;Ramshaw andTrapp, 1978;Stewart, 1979;

Stuhmiller,1977;TiseljandPetelin,1997;Travisetal.,1976).Although thephysicsofmultiphaseproblemsisbetterrepresentedwhen includ-ingtheeﬀectsoftheinterfaceforces,thisapproachdoesnotcompletely eliminatethepossibilityoftheproblembeingill-posed.

Fernández Nieto (2003), Castro Díaz et al. (2011), and

Sarno et al. (2017) introduce an additional term in the momen-tumequations toaccount forfrictionbetweentheﬂuidlayers.Their regularization strategyyieldsa well-posedmodelandhasaphysical origin.However,theadditionalphysicaltermdependsonthetimestep ofthenumericalsolution,whichimpliesthatitcannotbeconsidereda fullyphysically-basedsolution.

Thenumericalsolutionofamathematicallyill-posedmodelcanbe well-posed(ChenandPeng,2006;Chenetal.,2007;SavaryandZech, 2007;Spinewineetal.,2011)ifthenumericalsolutionneglects informa-tioninthephysicalequations(Grecoetal.,2008).Wordeddiﬀerently, insuchacasethephysicalequationsareill-posed,butthenumerical equationsthatweactuallysolvearewell-posed.Inparticular,when us-ingtheHLLsolver(acommonapproximateRiemannsolverproposed byHartenetal.,1983,seeToro,2009),oneonlyusesthefastestand slowesteigenvaluesofthesystem,whichimpliesthatthedynamicsdue totheintermediateceleritiesarenotresolved.Thishidestheproblem ofill-posednessratherthansolvesit.

Indetermining thesteady(equilibrium)stateofa fluiddynamics problem,acommonlyadoptedstrategytoachievefastconvergenceis tomodifytheceleritiesatwhichinformationpropagates(i.e.,the sys-temeigenvalues)(Chorin,1967;GrabowskiandBerger,1976;Plows, 1968;SohandBerger,1984).Forinstance,inaerodynamics,thespeed of soundmaydiffersignificantlyfrom theair velocity,which causes aslowconvergencetosteadystate(ChoiandMerkle,1993;Fengand Merkle,1990;Godfreyetal.,1993;vanLeeretal.,1991). Precondition-ingmethods(Turkel,1987;1993;1999)aimatbringingthe eigenval-ues ofthesystemclosertoeachothersuchthatalargertimestepis allowed.

Analogously,the“bedcelerity” (i.e.,thespeedofthewaverelatedto changesinbedelevation(LynandAltinakar,2002;MorrisandWilliams, 1996;Steccaetal.,2014;DeVries,1965))isgenerallyslowcompared totheceleritiesassociatedwithperturbationsoftheﬂow.Thisfacthas encouragedtheuseofa“morphodynamicaccelerationfactor” in mor-phodynamicmodelingtoreducethecomputationaltime(Latteux,1995; Lesseretal.,2004;Ranasingheetal.,2011;Roelvink,2006). Mathemat-ically,aswewillshowlater,theuseofamorphodynamicacceleration factorisequivalenttotheapplication ofaparticularpreconditioning method.

Byalteringthecelerityatwhichinformationpropagates,the tran-sientstateofthepreconditionedproblemisalteredwithrespecttothe originalproblem,butbothproblemsconvergetothesamesteadystate solutioniftheboundaryconditionsaresteady.Adrawbackof precondi-tioningisthefactthat,whentheproblemissubjecttounsteady bound-aryconditions,preconditioningmethodsmodifythesteadystate,asthey indirectlymodifythetimingoftheboundaryconditions(Turkel,1999). Forthisreason,theboundaryconditionsofapreconditionedmodelneed tobeadjustedifthesevarywithtime.

Thefactthatapreconditioningmethodaltersthetransientstatewas usedbyZanottietal.(2007)toregularizethetwo-fluidmodel.They modifiedthesystemofequationsbyintroducingtwoparameters.One parametermodifiesthecontinuityequation,whichaffectstheimaginary partoftheeigenvalues.Dependingontherelationsbetweenvelocities anddensitiesofthetwofluids,aspecificvalueofthisparametermakes theimaginarypartequaltozero.Apartfrommodifyingtheimaginary component,theparameteralsomodifiestherealpartoftheeigenvalues. Theyintroduceasecondparameterthataffectsallequationstorecover theoriginalrealpartoftheeigenvalues.Theycomparethesolutionof theregularizedmodeltoanalyticalsolutionsandtheyshow thatthe regularizedtwo-fluidmodelisstable.Inthenextsectionwewillfollow asimilarapproachtoderivearegularizationstrategyfortheactivelayer model.

(4)

3. Regularizationstrategyfortheactivelayermodel

Inthissectionweproposeastrategyforrecoveringthewell-posed character of the system of equations for modeling mixed-sediment rivermorphodynamics.Themodiﬁedsetofequationsispresentedin

Section3.1.InSection3.2-wederivetheparametersusedtorecover thewell-posedcharacterofthemodelconsideringasimplifiedcasewith twosedimentsizefractionsandsteadyflow,whichallowsustoobtain analyticalexpressions.Wethenextendthevaliditytounsteadyflow con-ditions(Section3.3)andtoconditionswithmorethan2sedimentsize fractions(Section3.4).InSection3.5wediscusstheimplementationof thestrategy.

3.1. Modiﬁedsystemofequations

Weconsiderone-dimensionalhydrostaticflowoverabedcomposed ofNnon-cohesive sizefractions.Theflow isdescribed bythe Saint-Venant(1871)equations.We assumea Chézy-typefrictionin which thenondimensionalfrictioncoefficientisindependentoftheflowand bed parameters. The sediment transport rateis considered toadapt instantaneouslytochanges in thebed shear stress(Belland Suther-land,1983).Themassconservationofthebedsedimentisdescribed bythe Exner (1920)equation,andthe 𝑁−1 activelayerequations (Hirano,1971)accountfortheconservationofthemassofeachgrain sizefractionwithinadiscretetoplayerofthebedsurface(i.e.,the ac-tivelayer).Giventheflow,friction,andsedimenttransportassumptions, themodelcannotrepresentsmall-scaleprocesses(i.e.,processesatthe scaleof bedelevationfluctuationsdue thestochasticnatureof sedi-menttransport,ripples,dunes,orbedloadsheets).Inotherwords,the variablesrepresentparametersaveragedoveraperiodlargerthanthe characteristictimeofsmall-scalebedelevationfluctuations(Armanini anddiSilvio,1988;Blometal.,2008;Parkeretal.,2000;Ribberink, 1987;WongandParker,2006).WerefertoAppendixAforthemodel equationsandthematrix-vectorformulationofthesystem.

AnalogoustoZanottietal.(2007)(Section2),thesystemof equa-tionsinEquation(1)ismodiﬁedmultiplyingthetimederivativeterm byadiagonalmatrixMtoregularizetheproblem:

𝐌𝜕𝐐_𝜕𝑡 +𝐀𝜕𝐐

𝜕𝑥 =𝐒. (2)

MatrixMmodiﬁesthetransientstateonly.Thepreconditioning tech-niquedoesnotaﬀectthesolutionofthesteadystate(i.e.,𝜕∕𝜕𝑡=0).

Themorphodynamicmodelunderunisizeconditionswasanalyzed byCordieretal.(2011).TheyfoundthattheSaint-Venant-Exnermodel isalwayswell-posedassumingaChézy-typefriction.Thisconﬁrmsthat theill-posedcharacterofthemixed-sizesedimentmodelresultsfromthe inappropriaterepresentationofthemixingprocessesbytheactivelayer model(Chavarríasetal.,2018).Forthisreason,weproposea regulariza-tionstrategythatrecoversthewell-posedcharactermodifyingthe celeri-tiesatwhichmixedsedimentprocessesoccur.Thisisdonebymeansofa setofparameters𝛼k[– ]for1≤𝑘≤𝑁−1thatmultiplythetime

deriva-tiveofeachactivelayerequation.SimilarlytoZanottietal.(2007),we consideraparameter𝛽 [– ]thatcanbeusedtorescalethecelerities af-terbeingmodiﬁedby𝛼k.Westipulatethatthisparameter𝛽 aﬀectsonly

thesedimentprocesses(includingtheExner(1920)equation)butnot theﬂow.

Themodiﬁedsystemofequationsmustbemassconservativewith respecttothesediment.Thisimpliesthat𝛼kcannotbegrainsize

de-pendent(i.e.,𝛼𝑘=𝛼 ∀k)andthatthepreconditioningtechniqueisonly

applicablewhentheactivelayerthicknessisconstant(AppendixB). 3.2. Derivationoftheregularizationcoeﬃcients

Inthissectionwederivethevaluesofthecoefficients𝛼 and𝛽 that enableregularizationoftheactivelayermodel.Weconsidera simpli-fiedcasewithtwosedimentsizefractionsundersteadyflowconditions,

asthisallowsustoobtainanalyticalexpressionsoftheregularization parameters.

Inthiscase,thedependentvariablesofthesystemarethebed eleva-tion𝜂 [m]andthevolumeofﬁnesedimentintheactivelayerperunit ofbedarea,M_a1[– ](Chavarríasetal.,2018andAppendixA):

𝐐s2= [

𝜂,𝑀a1 ]_⊺

. (3)

Thesystemmatrixis:

𝐀s2=𝑢 ⎡ ⎢ ⎢ ⎣ 𝜆b 𝜆s1 𝜇1,1 𝜆b𝛾1 𝜆s1 ⎤ ⎥ ⎥ ⎦, (4)

wheretheparameters𝜆b [−]and𝜆s1[−]arethenondimensional ap-proximated bed andsorting celerities,which (approximately) repre-senttheceleritiesatwhichinﬁnitesimalperturbationsinbedleveland grain size distributionof thebed surfacepropagatethrough the do-main(Chavarríasetal.,2018;Steccaetal.,2014;DeVries,1965and

AppendixA.5),andu[m/s]isthemeanﬂowvelocity.Theparameters 𝛾1 [−]and𝜇1,1[−]relate thechangesin thesedimenttransport rate tothepropertiesofthebed(AppendixA.5).Subscriptsindicatesthat themodelissteadyandsubscript2highlightsthatitaccountsfortwo sedimentsizefractionsonly.

Thepreconditioningmatrixis:

𝐌s2=𝛽 [ 1 0 0 𝛼 ] . (5)

Notethat𝛽 doesnotaffectthemathematicalcharacterofthesystem, asitmodifiesallequationsequally.Wordeddifferently,theparameter 𝛽 changesthemagnitudeoftheeigenvaluesbutnotthetype(realor complex). Wecomputetheeigenvalues(𝜆k for𝑘=1,2)of the

modi-ﬁedsystemofequationsastherootsofthecharacteristicpolynomial det(𝐌s2𝜆 −𝐀s2)=0: 𝜆𝑘=₂𝑢_𝛽 ( 𝜆b+ 𝜆s1 𝛼 ± √ Δ 𝛼 ) for 𝑘=1,2, (6)

wherethediscriminantΔisaseconddegreepolynomialon𝛼 equalto: Δ =𝜆2_b𝛼2+2𝜆b𝜆s1 (₂_𝛾 1 𝜇1,1 −1 ) 𝛼 +𝜆2 s1. (7)

Weconsiderasituationwhichisill-posediftheregularization strat-egyisnotapplied.Thisimpliesthatwhen𝛼 =1(theregularization strat-egyisnotapplied),Δ<0(theeigenvaluesarecomplex).Weaimto mod-ifythesystemofequationsaslittleaspossibleinregularizingit.Worded differently,weaimatchanging𝛼 aslittleaspossiblefrom1.The mini-mummodificationisobtainedwhenthediscriminantisequalto0(i.e., theeigenvaluesareinthelimitforhavinganimaginarypartdifferent than0).Thethresholdvalues𝛼cthatmodifythesystemofequationsas littleaspossiblearefoundbyequating(7)tozero:

𝛼c= 𝜆s1 𝜆b ( 1−2 𝛾1 𝜇1,1 ±2 √ 𝛾1 𝜇1,1 ( _𝛾 1 𝜇1,1 −1 )) . (8)

Therearetwopossiblevaluesof𝛼cthatyieldrealeigenvalues.The discriminant(Eq.(7))asafunctionof𝛼 isaconcaveparabolaas𝜆2

b>0. Moreover,when𝛼 =0,Δ =𝜆2_s1>0andwhen𝛼 =1,Δ<0.Thisshows thatonecriticalvalueofparameter𝛼 isbetween0and1,andthesecond valueislargerthan1(i.e.,0<𝛼c1<1<𝛼c2).

Wecomputethevalueofparameter𝛽 that,assuming𝛼 =𝛼c, recov-erstherealpartthattheeigenvalueswouldhaveiftheyhadnotbeen modiﬁedusingtheparameter𝛼:

𝛽 = 𝜆b+𝜆s1∕𝛼c 𝜆b+𝜆s1

. (9)

Inthiscase,irrespectiveof thevalueof 𝛼,theeigenvaluesof the regularizedsystemareequalto:

𝜆𝑘=𝑢₂(𝜆b+𝜆s1 )

(5)

Ifwedonotuse𝛽 torecovertheoriginalrealpartoftheeigenvalues (i.e.,if𝛽 =1),theeigenvaluesoftheregularizedsystemareequalto 𝜆𝑘= 𝑢₂ ( 𝜆b+ 𝜆s1 𝛼 ) for 𝑘=1,2. (11)

Parameter𝛼 canbeselectedtobelargerorsmallerthan1andif wechoosetouse𝛽 (i.e.,if𝛽 ≠ 1)theeigenvaluesareindependentof𝛼. Summarizing,weﬁndthreepossibleregularizationstrategies:

1. 𝛼 ≠ 1and𝛽 ≠ 1 2. 𝛼 <1and𝛽 =1 3. 𝛼 >1and𝛽 =1

Ingeneralterms,theapproximatedsortingceleritiesarepositive, andundersubcriticalﬂowconditions(i.e,Fr<1)theapproximatedbed celerityisalsopositive.However,duetohidinginthesediment trans-portrelation,underconditionsinwhichill-posednesslikelyoccurs,𝜆s1 maybenegativeregardless oftheFroudenumber(Chavarríasetal., 2018).Inthiscase,Strategies1and2donotguaranteethatthe eigen-values𝜆k>0.Weconsiderthatitisphysicallyunrealisticthat

morpho-dynamicinformationtravelsintheupstreamdirectionundersubcritical ﬂowconditions.Anegativeeigenvaluewouldimplythattheboundary conditionformorphologyneedstobeimposedatthedownstreamend toyieldawell-posedmodel,andthisiscontradictorytothefactthat themorphodynamicstateundersubcriticalﬂowconditionsdependson theloadcomingfromupstream(Blometal.,2017a;2016).Ontheother hand,Strategy3guaranteesthat𝜆k>0(AppendixAofthe

supplemen-tarymaterial).Thus,weconsiderthattheonlypossibleregularization strategyistheoneinwhich𝛼 >1and𝛽 =1.

Weneedtoguaranteethattheceleritiesofthesystemofequations modifiedbytheregularizationstrategyarenotphysicallyunrealistic.In particular,underasufficientlysmallFroudenumber,themodifiedbed andsortingceleritiesmustbesignificantlysmallerthanthecelerities oftheflow.Theregularizationtechniquedoesnotmodifythe approx-imatedcelerityassociatedwithbedelevationchanges(i.e.,𝛽 =1)and decreasesthecelerityassociatedwithmixingprocesses(i.e.,𝛼 >1,we willdiscussthispointinSection6.1).Forthisreason,theregularization techniquedoesnotcausetheceleritiestobephysicallyunrealistic.

Theregularizationstrategyisnotlimitedtoaparticularrangeof pa-rametersettings.Yet,whenusingthevalueof𝛼 derivedinthissection, theFroudenumbercannotbeinthetranscriticalregion,asinthiscase thequasi-steadyapproximationisnotvalid(CaoandCarling,2002b; Caoetal.,2002;ColombiniandStocchino,2005;Lyn,1987;Lynand Altinakar,2002;Sieben, 1999).Inthefollowing sectionweconsider unsteadyﬂow,whichextendstheregularizationtechniquetothe trans-criticalregion.

3.3. Validityunderunsteadyﬂowconditions

Inthissectionweextendthevalidityoftheregularizationparameter 𝛼 foundforsteadyﬂowcases(Section3.2)tounsteadyﬂowconditions.

Fig.1.Maximumimaginarypartofalltheeigenvaluesofthereferencecase (Table1)asafunctionof𝛼.Inthiscase𝛼c=16.1isthesmallestvalueof𝛼 >1

thatyieldsawell-posedmodel(i.e.,alleigenvaluesarereal). Table1

Referencevaluesinthecomparisonofthevalueof𝛼ccomputedanalyticallyand

numerically.

u [m/s] h [m] Cf [ − ] La [m] Fa1 [ − ] 𝑓 1I [ − ] d1 [m] d2 [m]

1 1 0.01 0.20 0 1 0.001 0.005

Whenconsideringunsteadyflowconditions,wecannotobtainan an-alyticalexpressionof𝛼cforregularizingthesystemofequations. Nev-erthelesswecannumericallyfindthesmallestvalueof𝛼 >1forwhich therootsofthecharacteristicpolynomialofdet(𝐌_u𝜆 −𝐀_u)=0arereal values(i.e.,theeigenvaluesarereal),wheresubscriptuindicatesthat themodelisunsteady.MatricesMuandAuarelistedinAppendixA.5. Thisprocedureisnonethelessexpensivecomputationallyincomparison withanalgebraiccalculation.Fig.1showsthemaximumimaginarypart ofalleigenvaluesofareferenceill-posedcase(Table1)considering un-steadyflowforvarying𝛼.Thesedimenttransportrateiscomputedusing afractionalversionoftheEngelundandHansen(1967)sediment trans-portrelation(Blometal.,2017a).Avalue𝛼 >16.1yieldsawell-posed model(i.e.,alleigenvaluesarereal).

Totestthevalidityofthealgebraicvalueof𝛼c obtainedassuming steadyflowweconsiderthesamereferencecase(Table1)andwevary theflowvelocitytoobtainarangeofconditions.InFig.2(a)wepresent thevalueof𝛼cnecessarytoobtainawell-posedmodelcomputed assum-ingsteadyflow(Eq.(8))andnumericallyconsideringunsteadyflow.We concludethatforaFroudenumberbelowapproximately0.6,thereisno significantdifferencebetweenthevaluesforsteadyandunsteadyflow. Thisimpliesthat,forFr<0.6,thevalueof𝛼cobtainedanalytically as-sumingsteadyflowisagoodapproximationoftheactualvalue.

Fig.2. Comparisonbetween(a)thesteadyandunsteady val-uesoftheregularizationparameter𝛼c,and(b)theexactand

(6)

3.4. Validityundermultiplesizefractionsconditions

Inamodelwithmorethan2sizefractions,wecannotanalytically obtainthevalueof𝛼cthatregularizestheactivelayermodel.Similar totheunsteadycase,itispossibletonumericallyobtainthesmallest valueof𝛼 >1thatyieldsrealeigenvaluescomputedastherootsofthe characteristicpolynomialdet(𝐌s𝜆 −𝐀s)=0(matricesMsandAsare pre-sentedinAppendixA.5).Again,thisprocessisrelativelyexpensivein computationalterms.Inthissectionweproposeamethodtoobtainan approximatevalueof𝛼cforsuchcasesandcompareittotheexactvalue obtainednumerically.

Assumingsteadyﬂow,asystemthatmodelsNsedimentsize frac-tionshasNequations(AppendixA).WereducethesystemofN equa-tionstoanapproximatesystemof2equationsfollowingtheapproach ofRibberink(1987).WesumtheNactivelayerequationstoobtainone equationthatmodelsthechangesofthemeangrainsizeofthebed sur-facesediment(AppendixA.3).Subsequently,weapplythesame tech-niqueastheonewehaveusedinthecaseof2sizefractionstoobtaina criticalvalueof𝛼 thatguaranteesthattheapproximatemodelis well-posed: 𝛼cm= 𝜆m 𝜆b ⎛ ⎜ ⎜ ⎝ 1−2𝛾m 𝜇m ±2 √ 𝛾m 𝜇m (_𝛾 m 𝜇m −1 )⎞ ⎟ ⎟ ⎠, (12)

wherethesymbolsaretheequivalentofthecasefortwosizefractions intheapproximatemodel(AppendixA.5)

Weconsideracasewith3sedimentsizefractions,wherethefine andcoarsefractionshavecharacteristicsizesequalto𝑑1=0.001mand 𝑑3=0.005m,respectively.Thevolumefractioncontentsofthe3size fractionsintheactivelayerare𝐹a1=0,𝐹a2=0.9,and𝐹a3=0.1.The substrateisfullycomposedoffinesediment.Wevarythemediumgrain size(d2)toobtainarangeofconditions.Theremainingparametersare thesameastheonespresentedinTable1.InFig.2(b)wecomparethe exactvalueof𝛼c(computednumerically)totheapproximatedone (com-putedusingEq.(12)).Theapproximatedvalueof𝛼cfollowsthesame trendastheexactone.However,theapproximatedvalueisbothlarger andsmallerthantheexactonedependingonthesedimentconditions. Thisimpliesthatthecurrentapproximateapproachmaybeinsufficient toregularizetheactivelayermodelinthecaseofmorethan2sediment sizefractions.

Theapproximatesystemofequationscanbeill-posedunder degra-dationalconditionsintoaﬁnesubstrateonly(Chavarríasetal.,2018; Ribberink,1987).However,a3sizefractionscasecanbeill-posed un-derdegradationalconditionsintoacoarsesubstrate(Chavarríasetal., 2018),whichfurtherlimitstheapplicabilityoftheapproximatesolution forthethresholdvalueof𝛼.

3.5. Implementation

Inthissectionwedescribeourapproachfornumericallysolvingthe systemofequationsandapplytheregularizationstrategy.

WehavedevelopedthenumericalresearchcodeElvtomodel mixed-sizesedimentrivermorphodynamics(Blometal.,2017a;2017b)which solvestheequationsforﬂow,bedelevation,andthebedsurfacegrain sizedistributioninadecoupledmanner(i.e.,inseriesandnotasa cou-pledsystemofequations).Thus,ourcodeisnotappropriatefor solv-ingtranscriticalsituations(Lyn,1987;LynandAltinakar,2002;Sieben, 1999)orcaseswithahighsedimentconcentration(CaoandCarling, 2002a;MorrisandWilliams,1996).

Theone-dimensionalspatialdomainis discretizedusing an equis-pacedgrid.Allvariablesarecomputedatthecellcentersandare con-sideredconstantineachtimestep.Hereweassumesteadyflow,whichis representedbythebackwaterequation(Eq.(16)).Thisordinary differ-entialequationisintegratedusingthestandardfourth-orderfinite differ-enceRunge–Kuttamethod(RK4).TheExner(1920)equation(Eq.(17)) andactivelayerequation(Eq.(19))aresolvedinconservativeform

us-ingafirstorderupwindschemeincombinationwithforwardEulerto integrateintime.Wediscretizetheverticaldomaininafinitenumber ofcellshavingacertainthicknesstoaccountforstratigraphicchanges inthesubstrate.Ourschemeisbalancedfortheverticalfluxesbetween theactivelayerandthesubstrate(Steccaetal.,2016).Thismeansthat massconservationisguaranteedindependentofthesubstrate discretiza-tion.ThetimestepvarieswithtimeandiscomputedsuchthattheCFL number(Courantetal.,1928)isconstantandequalto0.9(Toro,2009; Wu,2007).Thedetailsofthenumericalimplementationaredescribed inAppendixBofthesupplementarymaterial.

Whentheregularizationstrategyisapplied,weﬁrstdeterminethe mathematicalcharacter ofthemodel(i.e.,well-posedorill-posed)at eachnodeusingtheapproachproposedbyChavarríasetal.(2018).For thecaseof2sizefractions,thisisdoneevaluatinganalgebraicequation, andformorethan2sizefractionswenumericallycomputethe eigen-valuesofthesystemmatrix.Atcontinuation,foreachnodewecompute thethreshold value𝛼c that guaranteesthatthe modeliswell-posed. Again,this isdoneevaluatinganalgebraicexpression(Eq.(8)) for2 sizefractionsanditisdonenumericallyformorethan2sizefractions (Section3.4).

Theregularizationstrategyyieldsequaleigenvalues(i.e.,inatwo sizefractionscase𝜆1=𝜆2,Eq.(11)).Thisimpliesthattheproblemis hyperbolicbutnotstrictlyhyperbolic(Cordieretal.,2011;Lax,1980; Toro,2009).Inanon-strictlyhyperbolicproblem,thesolutionmaynot be uniqueandresonancemayoccur, whichgivesrisetostrong non-linearinteractions(IsaacsonandTemple,1992;Liu,1987).Inavoidinga non-strictlyhyperbolic-problem,wemodifythevalueof𝛼cusingasmall parameter𝜖 >0[– ]suchthat𝛼∗₌_𝛼

c(1+𝜖),where𝛼∗[– ]isthevalue usedforupdatingthebedsurfacegrainsizedistribution.Forthecases wehavestudiedavalueof𝜖 =0.005issuﬃcienttoavoidtheproblems associatedwithnon-stricthyperbolicity.

Ill-posednesscausesshort-waveinstability(Chavarríasetal.,2018; JosephandSaut,1990andSection5.1.2)meaningthatperturbations willgrowunstableatratesdependingontheinverseoftheirlength. Dif-fusioncounteractstheseeffectsbydampeningperturbations(Grayand Ancey, 2011). Regularization of the problem can be provided by numericaldiffusionifafirst-order(diffusive)methodisused.However, iftheunderlying problemis ill-posed,cellrefinementwillbeableto revealitsill-posedcharacterevenifafirst-ordermethodisusedinits solution,aswedointhispaper.Thisisbecause,withdecreasingcellsize, thenumericaldiffusioncoefficientofafirst-ordermethodwillgenerally decrease,whileatthesametimeshorter(moreunstable)perturbations willbe solved.Therefore,anill-posedproblemwillshowno conver-genceduetoitsinherentinstabilitywhenthemeshisprogressively re-fined,regardlessofthelow-ordermethodinuse.

Weobservesuchabehaviorin Section5.1.2whereweshowthat ourlow-ordernumericalschemesufficestocapturetheconsequences ofill-posednessbyrevealinginstabilityandnon-convergingcharacter in simulationsconductedwithintheill-posedrange. Itislikely that, withahigher-order(non-diffusive)method,thesefeatureswouldhave becomeapparentevenatlowermeshresolutionduetoabsenceof spuri-ousdiffusiondampeningperturbations.However,itmustbeconsidered thatourupwindschemeischaracterizedbysmallnumericaldiffusion coefficient,andthatSteccaetal.(2016)andSivigliaetal.(2017)have shownthatafirst-orderupwindschemewithafinegridresolutionis sufficienttocapturethemainfeaturesofmixed-sizesediment morpho-dynamicsimulationssuchastheonesweconduct.

4. Laboratoryexperiments

Inthis sectionwe describethelaboratoryexperiments conducted underconditionsinwhichtheactivelayermodel(Hirano,1971)is ill-posedinordertoobtainadatasettowhichwecancomparetheresults oftheproposedregularizationstrategy.Wedescribetheexperimental plan, materials,andmeasurementsin Section4.1.InSection4.2we presenttheexperimentalresults.

(7)

4.1. Experimentalplanandmeasurements

Weconducted4laboratoryexperiments(I1,I2,I3,andI4).The ex-periments reproduceddegradational conditionsinto aﬁne substrate, whichareconditionspronetobeill-posed(Chavarríasetal.,2018; Rib-berink,1987;Steccaetal.,2014).Theexperimentswereconductedin a14mlong,0.40mwide,and0.45mhightiltingﬂumeintheWater LaboratoryoftheFacultyofCivilEngineeringandGeosciencesofDelft UniversityofTechnology.Attheupstreamend,aturbulencedissipation devicewasinstalled(item(a)inFig.3).Aninclinedplanewasplaced downstreamfromtheturbulencedissipationdevice(item(b)inFig.3) toallowforanalluvialbed(item(c)inFig.3).Thestructurewas cov-eredwithgluedsedimentsuchthatfrictionwassimilartotheoneofthe alluvialbed.Itselevationcouldbeadjusted.

Weconsiderareferencesystemwithcoordinateoriginatthebottom oftheflumeatthedownstreamendofthemetalstructure.Thez-axis isparalleltogravityandpointingup.Thex-axisfollowsthestreamwise directionoftheflume,beingpositiveinthedirectionoftheflow.The y-axisisperpendiculartotheothertwoaxesformingarighthanded orthonormalbasis.

Weusedtwosedimentsizefractions(fineandcoarse)with charac-teristicgrainsizes(computedasthearithmeticmeanin𝜙 scale)equalto 2.1mmand5.5mm.Thestandarddeviationofthetwosizefractionsis 1.1mmand1.2mm,respectively.Thebedsurfacewasinitiallyflat,with aconstantslope,andcomposedofcoarsesedimentonly.Belowa0.03m thicklayerofcoarsesediment,weinstalledapatchoffinesedimentof varyinglengthLp[m](Fig.3andTable2).Weimposedaconstant wa-terdischargeandaconstantsedimentfeedrateofthecoarsefraction only,whichwasinequilibriumwiththeinitialcondition(Table3and AppendixC ofthesupplementarymaterial).Thesedimentwas intro-ducedusingafeederplacedontopoftheflume(item(d)inFig.3).The downstreamwaterlevelwasloweredatarateof0.01m/hduring8h byadjustingasharp-crestedweirat𝑥=12.60m(item(g)inFig.3).The loweringofthewaterlevelledtobeddegradationandentrainmentof thefinesedimentinthepatch.Wehavetestedthatintheseconditions theactivelayermodelisill-posedregardlessoftheactivelayerthickness andsedimenttransportrelation.

Sedimentwascollectedinasandtrap(item(e)inFig.3)atthe down-streamendoftheflume(𝑥=12.10m).Thesedimentwaspumpedfrom thesandtrap(item(f)inFig.3)intoatankpositionedonaweight bal-ancenexttotheflume.Thissystemallowedustocontinuouslymeasure thesedimenttransportrate.Thewaterinflowwasmeasuredusingan

Table2

Length(Lp)andposition(initialxp0andﬁnalxpfcoordinates)

ofthepatchofﬁnesedimentbelowthecoarsebedsurface. Experiment Lp [m] xp0 [m] xpf [m] I1 0.50 4.70 5.20 I2 1.00 4.49 5.49 I3 2.00 4.47 6.47 I4 4.00 4.47 8.47 Table3

Experimentalconditions,whereqdenoteswaterdischargeperunitwidth,s₀ initialbedslope,qb0sedimentfeedrateperunitwidth,hﬂowdepth,umean

ﬂowvelocity,andFristheFroudenumber. q [m 2 /s] _s

0 [ − ] qb0 [m 2 /s] h [m] u [m/s] Fr [ − ] 0.150 3.50 × 10 −3 _{7.86 × 10}−6 _0.187 _0.799 _0.59

ultrasonicflowmeterandthedownstreamwaterlevelusingaposition sensor.Weobtainedprofilesofthewaterandbedelevationusinglaser sensorsthatwerefixedtoacarriage(item(h)inFig.3).Acamerawas mountedonthecarriagetomeasurethegrainsizedistributionofthebed surfaceusingthetechniquedevelopedbyOrrú etal.(2016a,b)(item(i) inFig.3).Tothisend,thecoarsesedimentwaspaintedredandthefine sedimentblue.Ourexperimentalset-upallowedustomeasureeithera profileofbedandwatersurfaceelevationorthebedsurfacegrainsize distributionatacertainlocationwithtime.

Forthemodelingofthelaboratoryexperiments(Section5),itis im-portanttoobtainturbulentflow conditionsofa relativelydeepflow (i.e.,flowcannotbeaffectedbyindividualgrains),wheresedimentis predominantlytransportedasbedload.Theconcentrationofsediment needstobesosmallthatwecanassumeclearwater.Theseconditions weresatisfied(AppendixCofthesupplementarymaterial).

4.2. Results

Allexperimentsweregovernedbythesameconditionsbefore the ﬁnesedimentinthepatchwasentrained.Weobservedthesuperposition ofbedformsoftwodiﬀerentlengthscales(Fig.4).Secondarybedforms approximately0.5mlongand0.01mhighweresuperimposedon pri-marylongerbedformsoftheorderof3mandtwiceashigh.Theprimary

Fig.3. Sketchoftheﬂumeset-up:(a)Turbulencedissipator,(b)metalplatewithgluedsediment,(c)alluvialbed,(d)feeder,(e)sandtrap,(f)sedimentpump,(g) weir,(h)lasersensorsforwaterandbedsurfaceelevation,and(i)cameraformeasuringthebedsurfacegrainsizedistribution.

(8)

Fig.4. Measuredbedelevationbeforeﬁnesedimentofthepatchesisentrainedshowingthesuperpositionofbedformsoftwodiﬀerentlengthscales(ExperimentI4 at1:51h).

bedformsareinterpretedasincipientgraveldunes(Carling,1999; Car-lingetal.,2005).Thecharacteristicsofthesefeaturesremainedsteady asthebeddegraded.Thesteadinessofthefeatures’characteristicsis conﬁrmedinapreparatoryexperimental runwithoutapatchof ﬁne sediment(AppendixCofthesupplementarymaterial).

Afterapproximately2hthebedhaddegradeduptoapointatwhich thetroughofalongbedformwaslowerthanthetoppartofthepatch (Fig.5(a)).Atthatmomentfinesedimentwasexposed,entrained,and transported.Thelargermobilityofthefinesedimentcreateda down-streammovingdegradationalwave(Fig.5(b)).Aserosionproceeded, theshearstresswasreduced(duetotheincreasedflowdepth),which reducedthedegradationrate.Meanwhile,thesubsequentbedform ad-vancedandstartedtofillthedegradationalpitwithcoarsesediment (Fig.5(c)). Overall,thepassageof bedformsinducedentrainmentof finesedimentandsubsequentcoarseningofthetoppartofthesubstrate. Sincethedegradationalwaveincreasedindepthindownstream direc-tion,alsothethicknessofthecoarsetoplayerincreasedindownstream direction(Fig.5(d)).

Thesubstratecoarseningmechanismcreatedanirregularinterface between coarse andfine sedimentcompared tothe initial situation wheretheinterfacewasparalleltothebedsurface.Asaconsequence, theentrainmentoffinesedimentbecameapseudo-randomprocessin spaceandtime.Degradationalwavesformedatthoselocationswhere finesedimentwasclosesttothebedsurface.Yet,mostofthewavesgrew foronlyalimitedlength,as,duetotheirregularinterface,atsomepoint thesedimentpresentatthetroughwascoarseratherthanfine. Some-timestheinterface wassufficientlyparalleltothebedsurfaceanda largedegradationalwaveformed.Thisisseeninthecontentofcoarse sedimentatthebed surfaceof thepatch (Fig.6(a),(c),(e)and(g)) andinthebedelevation(Fig.7).Oneortwosmalldegradationalwaves formedafterthepassageofalargedegradationalwave,characterized bythefactthatthebedsurfaceiscomposedofmainlyfinesediment andthetroughofabedformreacheselevationssignificantlylowerthan average.

Alongerpatchallowedforthedevelopmentoflonger(inspaceand time)anddeepererosionalwaves(Figs.6and7).Yet,thedecreasein degradationrateasthewaveadvancedactedasasaturationmechanism limitingtheheightofthewave.Thus,theprobabilityoflowerbed el-evationatthepatchzonewasnotsignificantlylargerforanincreasing patchlength(Fig.8).Afterthepatch,wherethesubstratewascomposed ofcoarsesedimentonly,waveheightdecreasedandthebedelevation profiletendedtotheoneupstreamofthepatch(Fig.7).Yet,the pres-enceoffinesedimentdownstreamofthepatchslightlyincreasedthe heightofthebedformswithrespecttothebedformsupstreamofthe patch(Fig.8(a)and(c)).Bedformsdownstreamofthepatchwere char-acterizedbyacoarsefrontandfinetail,andwereapproximately2grain sizesofthecoarsesedimenthigh.Thesecharacteristicsmayindicatethe presenceofbedloadsheets(Dietrichetal.,1989;Reckingetal.,2009; Whitingetal.,1988)orbedformsinatransitionalphasetosmalldunes. Thedomaindownstreamfromthepatchwasnotlongenoughto

pre-Fig.5.Sketchofthecyclicentrainmentofsubstratesediment:(a)Bedforms formedoutofcoarsesedimentsonly,(b)ﬁnesedimentfromthepatchis en-trainedinthetroughofabedform,(c)adegradationalwaveformsandtravels downstream,(d)coarsesedimentfromupstreamﬁllsthepitleftbythe degra-dationalwave.

ciselyconcludeonthetypeofbedforms.Thechangesinvolume frac-tioncontentofcoarsesedimentatthebedsurfacewerelesspronounced downstreamofthepatchcomparedtoatthepatch(Fig.6(b),(d),(f) and(h)).Thisisbecauseﬁnesedimententrainedatthepatchdispersed inthedownstreamdirection.

(9)

Fig.6. MeasuredsurfacefractioncontentofcoarsesedimentasafunctionoftimeforvariouslengthsofthepatchLp:Atthecenterofthepatch(a,c,e,g),andatthe

downstreamend(b,d,f,h).Notethatthestreamwiselocationofthecenterofthepatchvariesforeachexperimentwhilethedownstreampositionisthesameforall cases(𝑥=9.15m).

5. Numericalmodeling

Inthissectionweapplytheregularizationstrategyinmodelingthe laboratoryexperimentsconductedunderconditionsinwhichtheactive layermodelisill-posed(Section5.1).InSection5.2wecomparethe re-sultsoftheregularizedactivelayermodeltotheresultsofthetwo-layer modeldevelopedbyRibberink(1987)byapplyingthemtoathought ex-perimentunderconditionsinwhichtheactivelayermodelisill-posed. 5.1. Modelingofourlaboratoryexperiments

In Section 5.1.1 we calibrate the numerical model. In Section 5.1.2weconductaconvergencetesttoshowtheconsequencesof ill-posednessandthe beneﬁts of theregularization strategy. InSection 5.1.3weapplythenumericalmodeltothelaboratoryexperiments de-scribedintheprevioussection.InSection5.1.4wetestthe regulariza-tionstrategyassumingthreesedimentsizefractions.

5.1.1. Calibration

Modelingthelaboratoryexperimentsrequiresvaluesfortheactive layerthicknessandthefrictioncoefficient,andthechoiceofasediment transportrelation.Tothisendweusetheresultsofasetof prepara-toryexperiments(AppendixCofthesupplementarymaterial).Tochose asedimenttransportrelation,weruntwoexperimentsconducted un-derequilibriumconditions,whilefeedingthefineandthecoarse sed-iment size fractions. Thesediment transport relation byAshida and Michiue(1971) reproduces ourresults reasonably well (Appendix D ofthesupplementarymaterial).Toobtaintheskinfrictioncoefficient (Cfb)forcomputingthesedimenttransportrate(AppendixA.4)we cor-rectthetotalmeasuredfrictioncoefficientC_fforsidewallfrictionwith

themethoddevelopedbyJohnson(1942)(seeGuo,2015).Weobtain thevalues𝐶f=0.0104and𝐶f b=0.0084. Bedformdragwasnegligible duringtheinitialphaseasbedformswerelow.Whenfinesedimentwas entrained,bedformsgrewandbedformdragmayhaveplayedarole.Itis notreasonabletomodelthisadditionalfrictionusingstandardrelations (e.g.EngelundandHansen,1967;HaqueandMahmood,1983;Wright andParker,2004),astheserelationsprovideabedform-averaged fric-tioncoefficient,whileinourcaselargebedformswereisolatedinspace andtime.Wedecidetouseaconstantfrictioncoefficientandwethink thatthemostsensibleapproachistoneglectbedformdrag.

Areasonablevaluefortheactivelayerthicknessis0.01m,which correspondstothedistancebelowthemeanbedelevationwitha proba-bilityofentrainmentbelowapproximately5%(Blom,2008;Ribberink, 1987). This valueis alsoin accordance with1–3 times D90 as pro-posedby,forinstance,Hirano(1971);HoeyandFerguson(1994),and

Seminaraetal.(1996).

Inonepreparatoryexperimentunderequilibriumconditions,wefed coarsesedimentonlyand,fromsomepoint,westartedfeedingtracer sediment(i.e.,sedimentofadiﬀerentcolor).Modelingthepropagation ofthefrontoftracersediment,weconﬁrmthat0.01misareasonable valuefortheactivelayerthickness(AppendixDofthesupplementary material).

5.1.2. Convergencetest

First we aim to show the consequencesof ill-posedness.To this end, we simulate conditions similar tothe ones of the experiments using theactivelayermodel. Intheexperiments,degradationintoa coarse substrate(i.e., under well-posedconditions) occurred for ap-proximately2h, asthepatchof ﬁnesedimentwas placed3cm be-low theinitialbedsurface.Inorder toobtainill-posedconditionsat

(10)

Fig.7. DetrendedbedelevationasafunctionoftimeinExperiment(a)I1,(b)I2,(c)I3,and(d)I4.Thedashedblacklinesindicatetheboundariesofthepatch. Thebedelevationisdetrendedsubtractingthebedslopeofeachprofileindividually,obtainedfittingafirstdegreepolynomial.

(11)

Fig.8. Probabilitydensityofdetrendedbedelevation:(a) Up-streamofthepatch,(b)atthepatch,(c)anddownstreamof thepatch.

thestartofthesimulations,thepatchofﬁnesedimentisplacedright belowtheactivelayer.Inthis manner,300ssimulationssuﬃce for ourpurpose.Moreover,thepatchextendsoveradistanceof8m(from 𝑥=1mto𝑥=9m)tomaximizethedomainoverwhichthemodelis ill-posed.

Weconduct13simulationsusingcellsizesrangingfrom0.1mdown to2.44× 10−5m.Theresultsdonotconvergeandcontinuetochangeas thegridisreﬁned(Fig.9(a)).Wecomputetheerrorasafunctionofthe cellsizetoquantifythe(lackof)convergence.Asthereisnoanalytical solutiontowhich wecancomparetheresults ofthenumericalruns,

Fig.9. Bedelevationat𝑡=300spredictedusingthe(a)activelayermodel(Hirano,1971)and(b)regularizedactivelayermodel.Eachofthe13linespresentsthe resultscomputedusingadiﬀerentcellsize(rangingfrom0.1mdownto2.44× 10−5_m,_where_darker_colors_represent_smaller_cell_sizes)._Panels₍_c₎_and₍_d₎_present

theerroratacertaintimeusingaparticularcellsize(seeEquation(13))whenusingtheactivelayermodelandtheregularizedactivelayermodel,respectively.In panels(b)and(d)onlyonelineisvisible,asitoverlapsallotherlines.

(12)

wecomputetheerrorbetweentheresultsoftwosuccessivesimulations sand𝑠+1(LoveandRider,2013; Roy,2005).Tothisend,ﬁrstwe interpolatethebedelevationresultsofallsimulationsusingthesmallest cellsize.Theinterpolation,ratherthanlinear,takesintoconsideration thateachvalueisconstantinsideacell.Second,wecomputetheerror asthenorm1ofthediﬀerencebetweenbedelevationsoftwosuccessive simulationsatacertaintimet:

error𝑡_𝑠= 1 𝐿𝑁x 𝑁x ∑ 𝑟=1 || |𝜂𝑟𝑡𝑠−𝜂𝑡𝑟𝑠+1|||, (13)

whereNxdenotesthenumberofcellsofthesimulationwiththesmallest cellsize,andL[m]thedomainlength.Fig.9(c)showstheerrorasa func-tionofthecellsizeforseveraltimes.Ifthecellsizeislarge(forinstance, largerthan0.01m),forshortsimulationtimes(forinstance,shorten than10s),theresultsseemtoconverge.Yet,usingthesamecellsize, theresultsdonotconvergeifoneconsidersalongersimulationtime. Similarly,consideringasimulationtimeequalto10s,theresultsdonot convergewhenthecellsizeissmallerthan0.002m.Thisbehavioris characteristicofill-posedsimulations.Thegrowthrateofperturbations increaseswithdecreasing cellsize.Forthis reason,theconsequences ofill-posednessariseearlierforsmallercellsizes.Givenacertaincell size,ifthesimulationisshortenough,perturbationsdonothavetimeto growandthesolutionseemstoconverge.Forafixedtime,simulations seemtoconvergeaftertheerrorgrows(forinstance,for𝑡=120s, sim-ulationsseemtoconvergeforcellsizesbetween0.001mand0.01m). Thisisdue tothefactthat,atthegiventime,perturbations have al-readygrownsignificantlyandhavecoarsenedthebedmaterialcausing thesimulationtobewell-posed.Afurtherdecreaseofthecellsizeoran analysisatadifferenttimeshowsthattheactivelayermodeldoesnot converge.

Werepeatthesamesimulationsapplyingtheregularization strat-egy.The initialvalue of theparameterthat recovers thewell-posed characterofthesystemis𝛼c=11.6.Inthiscasethesolutiondoesnot showoscillations(Fig.9(b)).Moreover,thesolutionconvergesfora de-creasingcellsizeindependentlyfromthetimeatwhichconvergenceis tested(Fig.9(d)).Thissupportsthefactthattheregularizedmodelis well-posed,contrarytotheactivelayermodel.Therateatwhichthe solutionconvergesconﬁrmsthatthenumericalschemeisﬁrst-order ac-curate(Section3.5).

5.1.3. Twosedimentsizefractions

Wereproducealllaboratoryexperimentsusingacellsizeequalto 0.05m.Theregularizedmodelshowsspatialortemporaloscillations innoneofthecases(Fig.10).Forallcasesthebedelevationdecreases smoothlyinthestreamwisedirection(Fig.10(b),(f),(j)and(n)).This contrastswiththemeasuredtemporalchangeofbedelevation,which presentsbedformsandtheformationofdegradationalwavesatthe up-streamendof thepatch (Fig.10(a),(e),(i)and(m)).Themeasured increaseinwaveheightatthepatch(Fig.10(a),(e),(i)and(m))and

Section4.2)isnotcaptured.Theeﬀectofthepatchisobservedinthe modelresultsinthefactthatdegradationoccursfasterforalongpatch (Fig.10)thanforashortone(Fig.10(b)).

Thecontinuousandsmoothpredictedentrainmentofsubstrate sed-imentyieldsanalmoststeadyvolumefractioncontentofsedimentin theactivelayerbothatthepatch(Fig.10(c),(g),(k)and(o))andatthe downstreamend(Fig.10(d),(h),(l)and(p)).Themeasureddatashows, ontheotherhand,avariablevolumefractioncontentatthebedsurface. Themodelcorrectlycapturesthemeanvalueandnicelyreproducesthat alongerpatchcausesanincrease intheamount offine sedimentat thebedsurface.Thefactthatthemodeldoesnotcapturebedformsis notsurprising,asthemechanismsnecessaryforbedformformationare notpresentinthemodel.Forinstance,thefactthattheflowmodelis basedonthehydrostaticpressureassumptionpreventsmodelling pro-cessessuchasflowseparation.Thepossibilityofcapturingthe forma-tionofthedegradationalwavesatthepatchisalsodiscarded,asfrom theanalysisofwell-posednessweseethattheregularizedmodeldoes

not showany instabilitymechanismthatcouldinducewavegrowth. Forthisreason,themodelresultsrepresentvaluesaveragedover the passageofseveralbedformsanddegradationalwaves.Wechoosenot toﬁlterthemeasuredbedelevationdata,asgiventhecharacteristicsof thebedforms,itwouldintroducealargeamountofspuriousinformation (e.g.,thedegradationalwavewouldstartatthewronglocation)andwe wouldloseasigniﬁcantamountofdataatthebeginningandendofthe domain.

Overalltheregularizedmodelyieldsareasonableapproximationof the meantemporalchange of themeasureddata. Thedegradational trendiscapturedandthesurfacegrainsizedistributionapproximates theaveragemeasuredvalues.Thesubstrateisnotunrealisticallyaltered astherearenooscillationsinthesolution.

5.1.4. Threesedimentsizefractions

Totesttheregularizationstrategyformultiplegrainsizes,wemodel ExperimentI4(Table2)using3diﬀerent grainsizesbyapplyingthe exactsolutiontoobtaintheregularizationparameter.Theﬁnesize frac-tionremainsthesameandthepreviouscoarse sizefractionis repre-sentedinthiscasebytwocharacteristicgrainsizesequalto4.895mm and5.895mm.Foraninitialvolumefractioncontentatthebed sur-faceofthemediumsizesedimentequalto0.375,theinitialbedslope isthesameaswhenusingtwocharacteristicsizesandthesumofthe sedimenttransportrateofthemediumandcoarsefractionswhenusing three sizesisequaltothesedimenttransportrateofthecoarse frac-tionwhenusingtwosizes.Inthismannerthesimulationaccountingfor threesedimentfractionsiscomparabletotheoneaccountingfortwosize fractions.

InFig.11wecomparethebedelevationandmeangrainsizeofthe bedsurfacesedimentpredictedbytheregularizedmodelusing2and3 sedimentsizefractions.Theevolutionofthebedelevationshowsonlya weakdependenceonthenumberofsizefractionsusedtodiscretizethe sedimentmixture.Themodelwith3sizefractionspresentsamild coars-ening(0.2%increaseinmeangrainsize)withtimebeforesedimentfrom thepatchisentrained(after2h).Thiscoarseningisnotvisiblewhen us-ing2sizefractions,becauseinthiscase,duringtheinitialstate,thebed surfacesedimentconsistsofonesinglegrainsizeonly.Weconcludethat theregularizationtechniqueisapplicableforageneralcasewithmore than2sizefractions.

5.2. ComparisonbetweenRibberink’s(1987)two-layermodelandthe regularizedmodel

Toourknowledgethereisnootherlaboratorydatasetapartfromthe onepresentedinSection4towhichwecanapplytheregularizedactive layermodeltotestitsperformance.Thisis becauseeitherthe condi-tionsthatotherresearchershavestudiedyieldawell-posedactivelayer model(e.g.AshidaandMichiue,1971)ortheactivelayermodelis ill-posedbuttheactivelayerthicknessvarieswithtimeduetodunegrowth (Blometal.,2003).Thelattercaseisasituationthattheregularization strategycannotdealwith(Section3.1).However,Ribberink(1987) ap-plieshistwo-layermodeltoathoughtexperimentunderconditionsin whichtheactivelayermodelisill-posed.Inthissectionweapplythe regularizedactivelayermodeltohisthoughtexperimentandcompare ittothetwo-layermodel.

Ribberink(1987)conductedalaboratoryexperimentwith mixed-sizesediment,whichwasdominatedbyaggradationafteraperiodof degradation(ExperimentE8-E9).Theinitialbedwascharacterizedbya uniformslope,composedofabimodalmixture(acoarseandﬁne frac-tion),andwellmixedbothinthestreamwiseandverticaldirection.The sedimentsupply wasinitially inequilibrium. Atemporalincrease of theproportionofthecoarsefractioninthesedimentsupplyperturbed theequilibriumconditionandinducedthedownstreampropagationof acoarseningwave.Thedownstreammigrationofthecoarseningfront causedaprecedingandtemporarybeddegradationasaresultof the

(13)

Fig.10. Comparisonbetweenmeasureddataandregularizedmodelresults:ExperimentI1(a-d),ExperimentI2(e-h),ExperimentI3(i-l),andExperimentI4(m-p). Thefirstandsecondcolumnsshowthemeasuredandpredictedbedelevationwithtime,respectively.Theverticaldashedlinesindicatethepositionofthepatchof finesediment.Thethirdandfourthcolumnspresentthesurfacefractioncontentofcoarsesedimentatthecenterofthepatchoffinesedimentandatthedownstream endoftheflume,respectively.

diﬀerenceinsedimentmobilitybetweenthecoarsesedimentforming thewedgeandtheﬁnesedimentdownstreamofthefrontofthewedge. Eventually,thebedaggradedandwascharacterizedbyalargerslope thantheinitialone,soastoallowforthetransportofthecoarserfed sedimentunderequilibriumconditions.

Duringtheshortdegradationalpartoftheexperiment,thebed sur-facewascoarserthanthesubstrate(i.e.,conditionsinwhichtheactive layermodelispronetobeill-posedChavarríasetal.,2018;Ribberink, 1987;Steccaetal.,2014).However,whilereproducingtheexperiment

numerically, Ribberink(1987)foundthattheactivelayermodelwas well-posed.Subsequently,Ribberink(1987)appliedhistwo-layermodel toathoughtexperimentthatwasequaltoE8-E9exceptforthefactthat thesubstratesedimentwasﬁnerthanintheﬂumeexperimentsuchthat theactivelayermodelisill-posed.Anumericalsimulationofthethought experimentusingtheactivelayermodelshowedoscillationsthat even-tuallymadethecodecrash(Ribberink,1987).Thethoughtexperiment was reproduced well bya numericalcode implementing Ribberink’s

(14)

Fig.11. Bedelevation(a)andmeangrainsizeatthebedsurface(b)asa func-tionoftimepredictedinExperimentI4usingtheregularizedactivelayermodel using2and3sedimentsizes.

Herewerunanumericalsimulationofthethoughtexperiment us-ingourregularizedactivelayermodelandcompareittotheresults of Ribberink’s (1987) two-layer model reported in Fig. 7.9 of

Ribberink(1987).SimulationdetailscanbefoundinAppendixC.

Fig.12presentsthetimeseriesofbedelevationandmeangrainsize ofthebedsurfacesedimentatalocation20mdownstreamfromthe inlet.Duringthefirst20htheeffectsofthecoarseningofthefed sedi-mentarenotfelt20mdownstreamfromtheinlet.Whiletheregularized activelayermodelpredictsaconstantbedelevationandgrainsize dis-tributionofthebedsurfacesedimentduringthisperiodof time,the two-layermodelpredictsafiningofthebedsurface(Fig.12(b)).This isduetothefactthattheinitialgrainsizedistributionoftheexchange layerisnotinequilibriumwiththeoneattheactivelayerandcausesa verticalfluxofsediment.However,thebedelevationremainsconstant aspredictedbybothmodels(Fig.12(a)).

Theaggradationalphaseisprecededbyadegradationalwave,which ismuchmorepronouncedintheregularizedactivelayermodelthan inthetwo-layermodel.Thisisbecauseintheregularizedactivelayer modeldegradationcausesentrainmentofthefinesubstratesediment, whereas in thetwo-layer model theexchange layeracts asa buffer thatslows down theprocess. The coarseningof thebed surface be-tweenapproximately 25 hand40 haspredicted byboth modelsis verysimilar.Whileafter40htheregularizedactivelayermodel pre-dictsa constantgrain sizedistribution of thebed surface sediment, thetwo-layer modelpredictsan asymptoticadaptationtoward equi-libriumconditions.Thiseffect isagaincausedbytheexchangelayer thatcoarsensslowly comparedtotheactivelayeron topof it,asit accounts forthe effectsof occasionally large bedforms.The equilib-riumstatediffersbetweenthetwomodels.Webelievethatthisisdue tothe factthat we do not knowexactly what values wereused by

Ribberink(1987)fortheconstantsinthesedimenttransportrelation (AppendixC).

Theregularizedactivelayermodelcapturesthedynamicspredicted bythetwo-layermodelofRibberink(1987).Theadvantageofthe two-layermodelisthatitaccountsforasourceofverticalmixingthatthe regularizedactivelayermodeldoesnottakeintoconsideration(i.e.,the mixingduetooccasionallylargebedforms).Ontheotherhand,the two-layermodelmaybecomeill-posed(Sieben,1994)whiletheregularized activelayermodelisalwayswell-posed.

Fig.12. Bedelevation(a)andmeangrainsizeofthebedsurfacesediment (b)withtimepredictedforthethoughtexperimentbasedonExperiment E8-E9conductedbyRibberink(1987)usingRibberink’stwo-layermodelandthe regularizedactivelayermodel.Theresultsofthetwo-layermodelareextracted fromFig.7.9ofRibberink(1987).

6. Discussion

Inthissectionwediscussthephysicalinterpretationofthe regular-izationstrategy(Section6.1),aswellaspossibleextensionsandfurther development(Section6.2).

6.1. Physicalinterpretationoftheregularizationstrategy

Theill-posedsolutionpredictedbytheactivelayermodelis charac-terizedbyoscillationsthattemporarilyfinethebedsurfaceandcoarsen thesubstrate.Thisbehaviorisalsoobservedinourlaboratory experi-ments(Figs.6and7).Onemaybetemptedtoconcludethattheactive layermodel,althoughbeingmathematicallyill-posed,provides reason-ableresults.Thisargumentiswrongfortworeasons.Thefirstreason isthatthenumericalsolutiondoesnotconvergeforadecreasingmesh size.Thesolutionkeepschangingandoscillationsbecomelargerwhen thecellsizeisreduced(Chavarríasetal.,2018;JosephandSaut,1990). Suchasolutioncannotberepresentativeofphysicalphenomena. Sec-ond,thephysicalprocessesresponsibleforthesmallscalevariabilityin bedelevation(i.e.,ripples,bedloadsheets)arenotaccountedforbythe activelayermodel(Section3.1).Anyresemblanceofthemodelresults withbedelevationfluctuationsduetosmallscalebedformsistherefore coincidence.

The frequently used morphodynamic factor(Φ_𝜂) (Latteux, 1995; Ranasingheetal.,2011;Roelvink,2006)isaparticularcaseofa pre-conditioningmatrixwithparameters𝛽 =1∕Φ_𝜂and𝛼𝑘=1∀k.The pro-posedregularizationstrategycanbeconsideredastheuseofa morpho-dynamicfactornotonlyforthechangesinbedelevation(𝜂)butalso forthechangesingrainsizedistributionofthebedsurface(Mak).The

“sortingmorphodynamicfactor” (Φ_s_k)isthendeﬁnedasΦs𝑘=1∕(𝛼𝑘𝛽).

Wehaveseenthattheonlyapplicableregularizationstrategyisthatin which𝛼𝑘=𝛼 >1∀kand𝛽 =1,whichisequivalenttosayingthatthe

regularizationstrategyisbasedona“sortingmorphodynamicfactor” 0<Φ_s<1.Thisimpliesthatthemixingorsortingprocessesassociated withchangesingrainsizedistributionofthebedsurfacesedimentare sloweddownwithrespecttotheceleritypredictedbytheactivelayer model.

(15)

Theeffectofapplyingtheregularizationstrategyisaslowdownof thesedimentmixingprocessesinthemodelcomputations.Thiseffectis similartotheeffectofa(temporary)increaseoftheactivelayer thick-ness.Fromaphysical perspectivethisslowdownofmixingprocesses maybeassociatedwitha(temporary)increaseoftherangeofelevations coveredbythebedlevelfluctuations(Blometal.,2008).Theslowdown ofmixingprocessesresultingfromapplyingtheregularizationstrategy impliesthattheregularizedactivelayermodelcanbeappliedtoawider rangeof physicalproblems (i.e.,alsothose characterizedbya fairly smalltimescaleofmixing)thantheactivelayermodel.

6.2. Alternativestotheregularizationstrategy

Our regularization strategy is applied locally and temporally. Wordeddiﬀerently,onlywhenandwherethemodelisill-posed,we up-datethegrainsizedistributionofthebedsurfacesedimentusingthe pa-rameter𝛼c.Moreover,𝛼cdependsonlyontheconditionsatthelocation underconsideration(notethatthepreconditioningmatrixisdiagonal). Thisisthesimpleststrategybutonecoulddecidetoavoid discontinu-itiesinthevalueof𝛼cthroughoutthedomainbycouplingneighboring nodes.

Carraroetal.(2018)proposeatechniquetodecreasethe computa-tionalcostofmorphodynamicsimulations.Asinourcase,theirstrategy canbeseenasapreconditioningtechnique.Theyconsiderunisize sedi-mentconditionsandmodifynotonlytheExner(1920)equationbutalso thecontinuityequation.Herewemodifytheactivelayerequationbut nottheﬂowequationsortheExner(1920)equation.Acombinationof bothstrategiescouldyieldatechniquethatbothdecreasesthecostof numericalsimulationsandguaranteesthatthemodeliswell-posed.

Wehavefocusedonrestoringthehyperboliccharacterofthesystem ofequationsandtothisendwebasedourstudyonthelinearsolution (i.e.,shortwaves).Thisfocussufficeshere,asshortwavesaremost sen-sitivetoill-posedness(JosephandSaut,1990).However,the regular-izationstrategymodifiesthecelerityandgrowthratenotonlyofshort wavesbutalsooflongones.Forthisreason,wesuggesttofurtherstudy howlongwavesareaffectedandwhethertheresultsofthe regulariza-tionstrategyarephysicallyrealisticbasedonasimilaranalysistothat ofLanzonietal.(2006).

We have assumeda constant activelayerthickness toavoid the addedcomplexityduetoacumbersomeclosurerelationlinkingthe pre-conditioningparameterstothechangeintimeoftheactivelayer thick-ness.Itmaybepossibletoextendourregularizationstrategyto situa-tionsinwhichtheactivelayerthicknesschangeswithtime(e.g.,dueto dunegrowth)byprovidingsuchaclosurerelation.Ontheotherhand,it isreasonablethattheregularizationstrategyrequiresaconstantactive layerthicknessgiventhefactthatmathematicallythestrategyhasthe sameeﬀectasanincreaseintheactivelayerthickness(i.e.,adecrease inthecelerityofthemixingprocesses).

Wehaveconcludedthattheregularizationstrategyneedsto slow-downthemixingprocesses(i.e.,𝛼c>1)toguaranteethatthe eigenval-uesarealwayspositiveregardlessofthevalueofthesortingcelerity 𝜆s1.However,ifthesortingcelerityisguaranteedtobepositive(e.g., becausehidingisnegligible),theaccelerationofthemixingprocesses alsoyieldspositiveeigenvaluesandawell-posedmodel.Theremaybe casesinwhichthelatterstrategyyieldsmorerealisticresults.Moreover, wehavechosentoguaranteethattheregularizedeigenvaluesare pos-itivereasoningthatmorphodynamicinformationtravelsinthe down-streamdirectionundersubcriticalconditions(Lanzonietal.,2006;Lyn andAltinakar,2002;Steccaetal.,2014;Suzuki,1976).Thisstatement ispartiallycontradictorytorecentstudiesthatconsidersediment trans-portasastochasticprocess(AnceyandHeyman,2014;Furbishetal., 2012).Thestochasticnatureofsedimenttransportyieldsan advection-diﬀusionequationthatmodelstheamountofmovingparticlesperunit ofbedarea.Thediﬀusivecharacterimpliesthatinformationalsotravels intheupstreamdirection.Forthisreason,aregularizationstrategyin

whichinformationtravelsintheupstreamdirectionmaybephysically realisticundercertaincircumstances.

Foracasewithmorethantwosedimentsizefractions(Section3.4), theapproximatevalueof theparameter𝛼c isnot (completely) satis-factoryaswell-posednessisnotguaranteed.Wehaveobservedinour teststhatill-posednessoccurswhen(atleast)twoeigenvaluesofthebed andsortingeigenvaluesaresimilarwithrespecttotheotherbedand sortingeigenvalues.Foracaseconsideringtwosedimentsizefractions this isreferred inliteratureasthe“crossing ofeigenvalues” (Sieben, 1997;Steccaetal.,2014).Wordeddiﬀerently,thediﬀerencebetween twoeigenvaluesmustbelargeenoughforthemodeltobewell-posed.A regularizationstrategybasedonguaranteeingaminimumdistance be-tweeneigenvaluescouldyieldaninexpensivesolutionforthecasewith morethantwosedimentsizefractions.

7. Conclusions

Wehavedevelopedapreconditioningmethodforregularizingthe activelayermodel(Hirano,1971)used inmodelingmixed-sediment rivermorphodynamics.Ourmethodrecoversthewell-posedcharacter ofthesystemofequationsbymeansofoneparameterthatmodiﬁesthe celerityofthemixingprocesses.Physicallythismeansthatthemixing processesaresloweddownorthetimescaleofthemixingprocessesis increased.

Weconduct4laboratoryexperimentsunderconditionsinwhichthe activelayermodelisill-posedandwecomparetheobservationstothe predictionsoftheregularizedactivelayermodel.Theregularizedactive layermodelcapturesthemeanbehaviorobservedintheexperiments associatedwithchangesaveragedoverthepassageofseveralbedforms.

Acknowledgments

This research is partof the research programme RiverCare, sup-portedbytheDutchAppliedandEngineeringSciences(AES)domainof theNetherlandsOrganizationforScientificResearch(NWO),andwhich ispartlyfundedbytheMinistryofEconomicAffairsundergrant num-ber P12-14(PerspectiveProgramme).The involvementof Guglielmo SteccawassupportedbyNIWAundertheSustainableWaterAllocation Programme.Thefruitfuldiscussionsandcommentsonthemanuscript ofLiselotArkesteijnaregratefullyacknowledged.WethanktheEditor Prof.Dr.G.C.Sander,AssociateEditorDr.C.Manes,andthree anony-mousreviewersfortheircomments,whichhavesignificantlyimproved themanuscript.

AppendixA. Modelequations

In this section we present the system of equations for model-ingmixed-sedimentrivermorphodynamics.InSectionA.1wepresent the ﬂow equations. In Section A.2 we present the active layer model (Hirano, 1971).A simpliﬁed activelayer model is presented in Section A.3. In Section A.4 we show the closure relations. In

SectionA.5we presentthesystemofequationsin matrix-vector for-mulation.

A.1. Flowequations

Weconsideraone-dimensionalmixtureofwaterandsediment flow-ingoveramobilebed.Asetofpartialdifferentialequationsthat ac-countsfortheinteractionsbetweensedimentandwaterisfoundby ap-plyingmassandmomentumconservationprinciplesforthemixtureof sedimentandwater(e.g.,Garegnanietal.,2011;Grecoetal.,2012). ThecompletesystemofequationsreducestotheSaint–Venant–Exner model (i.e.,clearwater approximation)underlow sediment concen-trations(𝑐=𝑞b∕𝑞<0.006,whereqb[m2/s]andq[m2/s]arethe sedi-menttransportrateandflowdischargeperunitwidth,respectivelye.g.,