Ill-posedness in modeling mixed sediment river morphodynamics

Delft University of Technology

Ill-posedness in modeling mixed sediment river morphodynamics

Chavarrias Borras, V.; Stecca, Guglielmo; Blom, Astrid

DOI

10.1016/j.advwatres.2018.02.011

Publication date

2018

Document Version

Final published version

Published in

Advances in Water Resources

Citation (APA)

Chavarrias Borras, V., Stecca, G., & Blom, A. (2018). Ill-posedness in modeling mixed sediment river

morphodynamics. Advances in Water Resources, 114, 219-235.

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

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.

ContentslistsavailableatScienceDirect

Advances

in

Water

Resources

journalhomepage:www.elsevier.com/locate/advwatres

Ill-posedness

in

modeling

mixed

sediment

river

morphodynamics

Víctor

Chavarrías

a,∗

_,

_Guglielmo

_Stecca

b,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, 10 Kyle St, Riccarton, Christchurch 8011, New Zealand.

c Department of Civil, Environmental and Mechanical Engineering, University of Trento, Via Mesiano 77, 38123 Trento, Italy/

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 15 March 2017 Revised 12 February 2018 Accepted 13 February 2018 Available online 16 February 2018

Keywords:

River morphodynamics Ill-posedness Mixed sediment

a

b

s

t

r

a

c

t

Inthispaperweanalyzethe Hiranoactive layermodel used inmixedsedimentriver morphodynam-icsconcerningitsill-posedness.Ill-posedness causesthe solutiontobeunstable toshort-wave pertur-bations.Thisimplies thatthe solutionpresents spuriousoscillations, theamplitudeofwhichdepends onthedomaindiscretization.Ill-posednessnotonlyproducesphysicallyunrealisticresultsbutmayalso causefailureofnumericalsimulations. Byconsideringatwo-fractionsediment mixtureweobtain an-alyticalexpressionsforthemathematicalcharacterizationofthemodel.Usingtheseweshow thatthe ill-poseddomainislargerthanwhatwasfoundinpreviousanalyses,notonlycomprisingcasesofbed degradationintoasubstratefinerthantheactivelayerbutalsoinaggradationalcases.Furthermore,by analyzingathree-fraction modelweobserveill-posednessunder conditionsofbed degradationintoa coarsesubstrate.We observethat oscillationsinthe numericalsolution ofill-posedsimulations grow untilthemodelbecomeswell-posed,asthespuriousmixingoftheactivelayersedimentandsubstrate sedimentactsasaregularizationmechanism.Finallyweconductaneigenstructureanalysisofa simpli-fiedverticallycontinuous modelformixedsedimentforwhichweshowthatill-posednessoccursina widerrangeofconditionsthantheactivelayermodel.

© 2018TheAuthors.PublishedbyElsevierLtd. ThisisanopenaccessarticleundertheCCBYlicense.(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Themixedcharacterofthesedimentisapropertynecessaryto explainphysicalphenomenasuchasdownstreamfining(Sternberg, 1875; Blom et al., 2016), the gravel sand transition zone (Yatsu, 1955;Blometal.,2017),theformationofbedloadsheets(Seminara et al., 1996), or bed surface armoring (Parker and Klingeman, 1982).Hirano (1971)wasthefirsttodevelopa massconservation modelformixed-sizesediment.The modelassumesthat the top-mostpartofthe bed,i.e.the activelayer, interactswiththeflow andisinstantaneouslymixed. Belowthe activelayerliesthe sub-stratewhichcanhaveverticalstratification.Inthisschematic rep-resentationofthemorphodynamic processesonlytheactivelayer sedimentisaffectedbyentrainmentanddepositionalprocesses.A vertical flux of sedimentoriginates from changes in elevationof theinterfacebetweentheactivelayerandthesubstrate.

Oneofthecriticalaspectsoftheactivelayermodelisthefact thattheverticalextentoftheactivelayer,oractivelayerthickness, shall be a prioriassigned.However, it cannot be physically mea-sured,asitstemsfromtheaboveschematicrepresentation(Siviglia

∗ _{Corresponding author.}

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

et al., 2017; Church and Haschenburger, 2017). The active layer thicknessisrelatedtothetimescaleoftheprocessunder consider-ation(BennettandNordin,1977;Rahueletal.,1989;Sieben,1997; Wu,2007).Inplanebedconditionsandshorttimescalestheactive layerthicknessisassumedtobeproportionaltothesizeofa char-acteristiccoarsefractioninthebed, forinstance, D 84 or D 90 (e.g., Pettsetal.,1989;Rahueletal.,1989;ParkerandSutherland,1990). Ifbedforms arepredominantandthe timescale under consider-ation involves the mixinginduced by the passage of severalbed forms, the active layer thickness is typically related to a charac-teristicbedformheight(e.g.,DeigaardandFredsøe,1978;Leeand Odgaard,1986;ArmaniniandSilvio,1988).Theactive layer thick-nessmayvaryover spaceandtime,although oftenitisassumed tobeauniformconstant.

Theactivelayer modelingframeworkhasproven tobeable to representawide varietyofphysicalphenomenasuchasbed sur-face armoring (e.g.,Park and Jain,1987) and the morphodynam-ics ofgravel-bed rivers (e.g.,Vogel etal., 1992) and tidal basins (e.g.,Carnielloetal., 2012).Moreover,itisimplementedinalarge amountofsoftwarepackagessuchasTelemac(Villaretetal.,2013), Delft3D(Sloff andMosselman,2012),andBASEMENT(Vetschetal., 2006).

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

bed

t=0

t=t

₁

t=0

t=t

₁ water water water water bed sorting a b unisiz e mixed

Fig. 1. Schematic of the effect of a perturbation in bed elevation in ( a ) a unisize sediment case and ( b ) a mixed sediment case. In the latter case, a perturbation in bed elevation introduces another wave, which is mainly related to the bed surface grain size distribution. Yet, each wave perturbs the flow, bed elevation, and bed surface grain size distribution. The arrows indicate the direction of propagation of the perturbations under subcritical flow conditions. The words “water”, “bed”, and “sorting” refer to a perturbation in water flow, bed level, and surface grain size distribution, respectively.

The mathematical representation of river morphodynamics shouldbewell-posed.Thismeansthat themathematicalproblem musthavea unique solutionwhich dependscontinuously onthe data(Hadamard, 1923).If the solutiondoesnot depend continu-ouslyonthedata,themodelisunfittorepresentthe correspond-ingphysics.

Despiteitswidespreaduse,theactivelayermodelhasone ma-jor mathematical shortcoming: the model can change its math-ematical character under some parameter settings. Therefore the mathematical problem that represents the physics of river mor-phodynamicscanbecome ill-posed.Thisfact wasfirstrecognized by Ribberink (1987). To this end he simplified the active layer modelby consideringan equation forthe meangrain sizeof the active layer sediment rather than one active layer equation for each grain size fraction. He found that under aggradational con-ditionstheproblemisunconditionallywell-posedandthesystem maybecome ill-posed underdegradational conditionsifthe sub-strateisfinerthantheactivelayer(i.e.degradationinanarmored river).Ribberink(1987)includeda thirdlayer betweentheactive layerandthesubstratetomodeltheeffectsofdunesexceptionally largerthantheaverageduneheight.Althoughthismodelincludes morephysicalmechanismsandimprovesthe predictionofmixed sedimentprocessesin dune-dominatedcases,it maystillbecome ill-posed(Sieben,1994).

To understandthe conditionsin whichthe active layer model becomesill-posedwe focusonhowinformationpropagatesalong ariver.We firstconsider acertain reachcharacterized by normal flow and immobile sediment. A perturbation of the flow propa-gatesalongtheriverintheformoftwowavestravelingatspeeds equalto u ±



gh where u [m/s]denotesthemeanflowvelocity, h

[m]theflowdepth,and g [m/s2_{]istheaccelerationduetogravity.} Ifsedimentismobile,yetuniform,aperturbationinbedelevation (e.g.,asedimenthump)willpropagatewithaspeedthatistermed the“bedcelerity” (de Vries,1965;LynandAltinakar,2002;Stecca etal., 2014). Asthe bedelevationaffectsthe flow,thebed eleva-tionperturbationalsoinducesaperturbationoftheflow.Thus, un-derunisize sedimentconditions,aperturbation ofthebed eleva-tionleadstothreewaves(Fig.1a).Althougheachofthewaves per-turbbothbedelevationandflow,twoofthewavesperturbmainly theflowwithoutmuchchangeinbedleveliftheFroudenumber

(F r = u/



gh )issufficientlysmall(deVries,1973;Needham,1990; Zanré andNeedham,1994).

The considerationof mixedsediment(of two sizefractions to simplifytheexample)introducesanotherceleritywhichistermed the“sortingcelerity” (Suzuki,1976;Ribberink,1987;Stecca etal., 2014). Thus, under mixed sediment conditions (with two grain sizes),aperturbationofbedelevationcausesfourwaves.Although eachwaveperturbstheflow, bedelevation,andsurfacegrainsize distribution, two of these perturb mainly the flow, one mainly the bedlevel, andone mainlythe surface grain size distribution (Ribberink,1987;Steccaetal.,2014)(Fig.1b).Sieben(1994) identi-fiedaregionofparameterswhere,forasedimentmixture consist-ingoftwograinsizeclassesunderbeddegradationintoasubstrate finerthanthe activelayer,the modelisunconditionallyill-posed. Thisoccurs whenthe “sorting celerity” equals the“bed celerity”. ThiswasconfirmedbySteccaetal.(2014),whoobserved,through numericalcomputationofthesystemeigenvalues,suchmodel be-havioralsoincaseofmorethantwosedimentfractions.

Furthermore, Stecca et al. (2014) analytically confirmed the outcomes of Ribberink’sanalysis using a more realistic unsteady modelfortwosedimentsize fractions.Theyconsideredgrain size selectivity ofthebedload buthiding ina limitedmanner. Hiding accountsforthefactthat grainsizefractions finerthana charac-teristic meangrain size of the mixturehide behind largergrains andsotheyexperiencealargercriticalbedshearstresscompared totheunisizecase(Einstein,1950;Komar, 1987a;1987b).The op-posite happensforcoarse sedimentfractions,whichexperience a larger exposure tothe flow than ina unisize case. In their anal-ysis Stecca et al. (2014) showed that the model can become ill-posedunderdegradationalconditionsifandonlyifthesubstrateis finerthanareferencegrainsizedistributionwhichisrelatedtothe grainsizedistributionofthebedload,insteadoftheactivelayer(as inRibberink’s(1987)analysis).

Toovercometheproblemofsettingtheactive layerthickness, Parkeretal.(2000)developedastochasticframeworkwithoutthe needforadistinctionbetweentheactiveandinactivepartsofthe bed. BlomandParker (2004),Blometal.(2006),andBlomet al. (2008)developedamodelthataccountsfordunesortingandthe variabilityofbedelevationbasedonthestochasticframework de-velopedby Parker etal.(2000).The modelassociates a probabil-ityofgrainsizeselectiveentrainmenttoallelevationswithin the bed, and hence allows for sediment at any elevation to be en-trained and contribute to the bedload discharge. Viparelli et al. (2017)developeda simplifiedverticallycontinuousmodel assum-ingslowchanges inbedelevationandasteadyprobability distri-butionofentrainment,deposition,andbedelevation,whichmake their modelsuitableforlargespaceandtimedomains. Sofarthe well-posednessofthecontinuousmodelhasneverbeenassessed.

Ourmain objectiveis to analyzethe problemofill-posedness oftheactivelayermodelusedformixedsediment morphodynam-ics. The present paperprovides four key improvements with re-specttopresentlyavailableknowledge:(i)weobtainanalytical ex-pressionsto characterizea simplifiedmodel (i.e.,tofindwhether itisill-posedorwell-posed)withtwosedimentfractionsonly,(ii) we study the effect ofmodel parameter choice on ill-posedness, (iii)wefindnew(previouslyneglected)ill-poseddomains,and(iv) we studytheconsequences ofill-posednessin numerical simula-tions. Our second objectiveis to mathematically characterize the verticallycontinuousmodeldevelopedbyViparellietal.(2017).In thenextsectionwepresentthegeneralsetofequationsfor mod-eling mixedsedimentriver morphodynamicsusing (a) the active layermodeland(b)theverticallycontinuousmodeldevelopedby Viparellietal. (2017).The models aresimplified andanalyzed in Section 3. Weanalyze the effectofmodel parameters on the ill-posednessoftheactive layermodelinSection4.InSection 5we

studythe consequences ofill-posednessusing numericalruns. In Section6werelaxandstudythesimplificationsofouranalysis. 2. Modelequations

In this section we presentthe equations used to modelriver morphodynamics.These equationsrepresentone-dimensional hy-drostaticflowoveramobilebedcomposedofanarbitrarynumber

N ofnon-cohesivesedimentfractionscharacterizedbyagrainsize

d k[m],wherethesubscript k identifieseachfractioninincreasing

size(i.e., d 1 < d 2 <...<d N).

Inthefollowingsectionwedescribetheflowequations.As pre-viousresearchhasclarified(Ribberink,1987;Steccaetal.,2014),a key parameter indetermining well-posedness of the active layer modelistheactivelayerthickness.Inthispaperwebothconsider a model with constant, andwith unsteady (time-varying) active layerthickness.Whilethewell-posednessofthemodelwith con-stant active layer thickness has been analyzed in previous work (Ribberink,1987; Stecca et al., 2014), toour knowledge no anal-ysis of the well-posedness of the model with unsteady active layer thicknessisavailable, althoughuse ofsuch amodelis doc-umented inthe literature (e.g.Karimet al., 1983). Theequations oftheadapted activelayermodelarepresentedinSection 2.2.In Section2.3wepresenttheverticallycontinuousmodelderivedby Viparelli et al. (2017). The closure relations for both models are treatedinSection2.4.InSection2.5wepresentacompactmatrix formulationofthemodelequations.

2.1. Flow equations

Theflow isdescribedby the1D ShallowWaterEquations(i.e., theSaint–Venantequations,Saint-Venant(1871))considering con-stantwidth.Assumingsteadyflowconditionsthewaterdischarge is uniform and conservation of momentum reduces to the so-called backwater equation (Appendix A). When assuming steady flow overamovablebedcomposedofsedimentofdifferentsizes we implicitlyassume thattheflowadaptsinstantaneously to per-turbationsinbedelevationandgrainsizedistribution.Worded dif-ferently, we assume that flow perturbations propagate infinitely fast relative to perturbations in bed elevation and surface grain size distribution. This assumption is referred to in literature as the“quasi-steadyflowassumption” (deVries,1965;Zhangand Ka-hawita,1987;1990;CaoandCarling,2002a).Thequasi-steadyflow assumptionisacceptableprovidedthat theFroudenumberis suf-ficiently small,1− Fr 2_=O

₍

₁

₎

_{(de Vries,1973; Sieben,1999; Lyn} and Altinakar, 2002). Note that in this context the term “quasi-steady” has a meaningdifferentfrom, forinstance, its useinthe modelingofflood waveswhere“quasi-steady” refers tonegligible inertiainthemomentumbalance.

2.2. Adapted active layer model equations

The conservation of the total amount ofsediment in the bed is formulated by the Exner equation (Exner, 1920). The active layer equation describes massconservation foreach size fraction (Hirano,1971).AppendixBpresentsthedetailsoftheactivelayer model.

Toanalyzethemodelwithunsteadyactivelayerthickness, we firstneedtosetaclosurerelationexpressingthethicknesschange in time. We consider an empirical empirical power relation be-tween duneheight H [m]andflowdepth h [m](Yalin,1964;Gill, 1971):

H=aLhbL, (1)

where a L[m1−bL]andb L[-]areconstants.Allen(1968a,b)proposed

valuesof a L=0.1∼ 0.2 m1−bL and b L=0.9∼ 1.2(with h in[m]).

Assuming that the active layer thickness L a [m] is equal to the

meanduneheight (Blom, 2008),we relate the activelayer thick-nesstotheflowdepthasfollows:

La=aLhbL. (2)

Toobtainanequationfortheactivelayerthicknessvariationwe differentiatetheconstitutivelaw,Eq.(2),withrespecttotimeand thensubstitutethecontinuityEq.(A.1)init:

∂

La

∂

t =−aLbLh

bL−1

∂

q

∂

x. (3)

SubstitutionofEq.(3)intotheactivelayerEq.(B.2)yieldsthe fol-lowingadaptedactivelayerequation:

∂

Mak

∂

t − f I k

∂

qb

∂

x +f I kaLbLhbL−1

∂

q

∂

x+

∂

qbk

∂

x =0, (4) where t [s] denotes the time coordinate, x [m] the streamwise coordinate, q = uh [m2_{/s] thewater discharge per unit width, q}

b

[m2_{/s]isthesedimenttransportrateperunitwidthmultipliedby} 1/

(

1− p)where p [-]isthebedporosity(i.e.,thesediment trans-portrate q b accountsforpores), q bk[m2/s]isthesediment

trans-portrateper size fraction, M _ak [m]is thevolume ofsedimentof size fraction k in the active layer per unit of bed area, and f I k

[-]is the volume fractioncontent of size fraction k atthe inter-facebetweentheactivelayerandthesubstrate.InFig.2weshow aschematicrepresentationofthemainvariablesoftheactivelayer model.

2.3. Simplified vertically continuous model equations

Theconservedquantityintheverticallycontinuousmodel (sim-ilarto M ak intheactivelayermodel)istheproductofthe

cumu-lativeprobabilityofbedelevation(P e [-])andthevolumefraction

contentofaspecificgrainsize class k (f _k[-])(Parkeretal.,2000; Pelosietal.,2014).Theverticalcoordinateis z [m].Tosimplifythe problem,theprobability distributiondependsonasecond vertical coordinate y ₌ z ₋

η

whichiscenteredatthemeanbedelevation. Assumingslowchangesinmeanbedelevationandaconstant(in timeandspace)probability distributionofbedelevation,Viparelli etal.(2017)obtainan equationforthechangeintimeofthe vol-umefractioncontent, f k:

Pe

∂

fk

∂

t =−pe

∂

qbk

∂

x −

∂

qb

∂

x

∂

fkPe

∂

y , (5)

where p e [m−1]istheprobabilitydensityfunctionofbedelevation

(Fig.2).

As in the active layer model,information is only advected in streamwise direction, i.e., the conservationequation does not in-cludedivergence terms inthe y direction andthe only indepen-dentvariable inspace is x .In contrastto the active layermodel, thereisnoinactivesubstrateandsedimentatallelevationsplays arole.Thisisillustratedbythedependenceoftheprobability func-tiononthe y coordinateandthegradientinthe y directioninEq. (5).Thus,althoughthesystemofequationsisone-dimensional,the mathematicalcharacterofthemodelisa propertydependingnot onlyonthestreamwisecoordinate x butalsoonthevertical coor-dinate y .

2.4. Closure relations

We apply the Chézy lawforthe friction slope. Thus, the fric-tionslopeisproportionaltothesquare ofthemeanflowvelocity dividedbytheflowdepth, S f=C fu 2/

(

gh

)

,where C f[-]isa

nondi-mensionalfrictioncoefficient.Forsimplicity,weassumeaconstant nondimensionalfrictioncoefficientthatisindependentoftheflow andbedparameters.

η

₀

η

h

L

_a

p

_e

F

_ak

f

_kI

f

_sk

(z)

p

_e

f

_k

(z)

Active Layer Model Vertically Continuous Model

q

q

_bk

Fig. 2. Representation of the main variables of the active layer model ( Hirano, 1971 ) and the vertically continuous model proposed by Viparelli et al. (2017) .

Table 1

Values of parameters in Eq. (6) according to several authors. Author A [-] B [-] θc [-]

Meyer-Peter and Müller (1948) 8 1.5 0.047

Engelund and Hansen (1967) 0.05/ C f 2.5 0

Fernandez-Luque and Van Beek (1976) 5.7 1.5 0.037 – 0.0455

Wong and Parker (2006) (1) 4.93 1.6 0.047

Wong and Parker (2006) (2) 3.97 1.5 0.0495

We apply a generalized form of the Meyer-Peter and Müller (1948)transportrelation,in whichthesedimenttransport rateis apowerfunctionoftheexcessbedshearstress:

q∗_bk=A[max

(

θ

k−

ξ

k

θ

c,0

)

]B, (6)

where q ∗_bk [−] is a nondimensional sediment transport rate (Appendix C), A [-] and B [-]are nondimensional parameters,

θ

k

[-]is thenondimensional bedshearstress ofsize fraction k , also knownasShields (1936) parameter,

θ

c [-]isthe nondimensional

criticalbedshearstress,and

ξ

k[-]isthehidingcoefficient.Table1

summarizesappropriatevaluesof A, B ,and

θ

caccordingtoseveral

authors.

A common hiding function is the one due to Egiazaroff (1965)(Appendix C). Asimpler relationwasdevelopedby Parker etal.(1982):

ξ

k=



_Dm dk



b , (7)

where D m [m] is a characteristic mean grain size and b is a

nondimensional parameter. A value of b >1(Dhamotharan etal., 1980; Misri etal., 1984; Kuhnle, 1993) impliesthat hiding is so strongthatthecoarserfraction(s)inthemixtureis(are)more mo-bilethanthe finerone(s),i.e.,reverse mobility(SolariandParker, 2000).

A final closure relation is required (only for the active layer model) forthe volume fractioncontent of sedimentof size frac-tion k attheinterface betweentheactive layerandthesubstrate,

f I

k.Whentheinterfacelowersthetextureattheinterfaceisequal

tothat atthe topmost partof the substrate. When the interface elevationincreasesvariousrelations canbe appliedfor f I

k.Hirano

(1971)proposedthatduringaggradationthegrainsizedistribution attheinterfaceisequaltotheoneoftheactivelayer.Accordingto Parker(1991)alsothebedloadsedimentplays aroleinthe aggra-dationalfluxtothesubstrate.HoeyandFerguson(1994)combined bothconceptsinoneparameter

α

[-]spanningtherange[0,1]that describesthecontributionoftheactivelayerrelativetotheoneof

thebedload: fI k=

⎧

⎪

⎨

⎪

⎩

fsk

(

z=

η

− La

)

if

∂

(

η

− La

)

∂

t <0

α

Fak+

(

1−

α

)

pk if

∂

(

η

− L a

)

∂

t >0 , (8)

where p k= q bk/q b [-]isthefractionofsedimenttransport rateof

sizefraction k .

2.5. Matrix formulation

Inthissection weintroduce amatrixformulationtoasses the well-posednessofthesystemofequations.

Asystemofpartialdifferentialequations(PDEs)canbe mathe-maticallyclassifiedasbeingofahyperbolic,elliptic,ormixedtype (e.g.,CourantandHilbert,1989).Tothisendwewritetheproblem inmatrix-vectorform(e.g.,Toro,2001):

∂

Q

∂

t +A

∂

Q

∂

x =S. (9)

Thisequationistheone-dimensionalquasi-linearnon-conservative formoftheadvectionequation.Qisthevectorofdependent vari-ables,Aisthesystemmatrix,andSisthevectorofsourceterms.

Asystemishyperbolicatapoint(x, t ) ifalltheeigenvaluesof matrixAarereal.Physicalpropagationproblemsaremodeledwith hyperbolicsystemsofequations.Ifalleigenvaluesarecomplex,the systemistermedelliptic.Ellipticsystemsmodelequilibrium phys-icalproblems.IfmatrixAhasbothrealandcomplexeigenvaluesit isamixed-typesystem.

A space-time dependent problem, in which we prescribe boundaryconditionsasafunctionoftimeandaninitialcondition (as is the casein modeling river morphodynamics), governed by an elliptic set of equations is ill-posed (Hadamard, 1923; Joseph andSaut,1990;Kabanikhin,2008). Thisisconfirmedby a pertur-bationanalysisthat showsthat,ifall eigenvaluesofmatrixA are real,perturbations ofa referencestate are bounded (AppendixA in thesupplementary material). However, ifthere is atleast one complexeigenvalue(or,precisely,atleasttwo,becauseofthe com-plex conjugate), perturbations grow exponentially. The exponen-tialgrowth dependsonthe product ofthe imaginarypart ofthe eigenvaluesandthewave numberoftheperturbation,which im-pliesthatthesolutionofanill-posedproblemisunstabletoshort perturbations.Attemptstonumericallyintegrateanill-posed prob-lemthereforeproduceresultsthatcontinuetochangeasthegrid isrefined (Woodhouse etal., 2012;Barker etal.,2015), asin nu-merical solutions perturbations always existdueto atleast trun-cation errors. In numerical simulations the wave number of the shortestpossibleperturbation isinverselyrelated tothe horizon-taldiscretization(



x ).

Byonly usingthe eigenvalues ofA tocharacterize the system of equationswe are neglecting the effectoffriction (AppendixA in the supplementary material). Yet,this suffices hereas friction becomesrelevantforsmallwavenumbersonly(AppendixAinthe supplementary material)andthe most criticalwave numbers, as regardstooscillationgrowth,arethelargeones.

Asasinglecomplexpairofeigenvaluesmakestheproblem ill-posed,wedonotmakeadistinctionbetweenthenumberof com-plex eigenvalues.Weterma problemwithatleastapairof com-plexeigenvaluesaselliptic.

We recast in matrix-vector form the Saint–Venant equations, (A.1)and(A.2),theExnerequation,(B.1),theactivelayerthickness equation,(3),andtheadaptedactivelayerequation,(4).Thevector ofdependentvariablesis:

Qal=

h,q,

η

,La,[

M



ak

] N−1



T , (10)

thesystemmatrixis:

(11) andthevectorofsourcetermsis:

Sal=

0,−ghSf,0,0, []



N−1



T . (12)

Thebrackets([])highlightthosetermsthatarevectorsor ma-trices.

We also recastin matrix-vector formthe Saint–Venant equa-tions,(A.1)and(A.2),theExnerequation,(B.1),andthe conserva-tionequation oftheverticallycontinuousmodel,Eq.(5).The vec-torofdependentvariablesis:

Qvc=

h,q,

η

,



[fk] N−1



T , (13)

thesystemmatrixis:

(14) andthevectorofsourcetermsis:

Svc=

0,−ghSf,0,



[] N−1



T . (15)

3. Characterizationofthemathematicalmodels

In this section we analyze the mathematical character of the modelsdescribedinSection2.Eigenvaluescomputednumerically can be obtained for an unlimited number of fractions. Here we studya simplecaseassuming steadyflow andtwo size fractions toobtainanalyticalexpressionsoftheeigenvalues.

As in our casethe temporal change of the active layer thick-ness depends on the spatialgradient of the water dischargeper unitwidth,Eq.(3),thesteadyflowassumptionimpliesaconstant active layer thickness. Yet, in a numerical simulation where the steadyflowassumptionisusedbuttheupstreamdischargevaries withtime (i.e.,alternating steadyflow), theactive layerthickness mayvarywithtime.However,insuchacasetheperturbationsdue toachangeinactivelayerthicknesspropagateinfinitelyfast rela-tive to the perturbations in bedelevation and surfacegrain size distribution.

Theimplicationsofmorethantwosedimentsizefractionsand anactivelayerthicknessasafunctionoftheflowdepthare stud-iedinSection6.

3.1. Steady active layer model consisting of two size fractions

Substitution of the backwater equation, (A.3), in the Exner equation, (B.1), and the active layer equation, (B.2), allows us to obtaina reducedmodel wherethe vectorof dependentvariables (QalS2)is:

QalS2=[

η

,Ma1]T, (16) thevectorofsourceterms(SalS2)reads:

SalS2=−Sf

u

ψ

1− Fr2[1,

γ

1]

T_{, (17)}

andthesystemmatrix(A_alS2)is:

AalS2=u

⎡

⎢

⎣

ψ

1− Fr2

χ

1

ψ

1− Fr2

γ

1

χ

1

μ

1,1

⎤

⎥

⎦

. (18) Wedefine

ψ

[-]as:

ψ

=

∂

qb

∂

q, (19)

whichis a parameter relatedto the intensityof total bedload in theflowandrangesbetween0(nullsedimentdischarge,i.e.,fixed bed)andO

(

10−2

)

(highsedimentdischarge),(e.g.,deVries,1965; LynandAltinakar,2002;Steccaetal.,2014).

The parameter

γ

1 [-] is a measure of the fraction content of sedimentintransportrelative tothefractioncontentofsediment attheinterfacebetweentheactivelayerandthesubstrate(Stecca etal.,2014):

γ

1=c1− f1I, (20)

where c 1∈[0,1] [-] isa parameter expressing the increase in the sedimenttransportintensityofthefinefractionrelativetothe to-talsedimenttransportintensity(Steccaetal.,2014):

c1= 1

ψ

∂

qb1

∂

q . (21)

Wenowintroducetheparameter

χ

1 [-]whichisa nondimen-sional measure of the derivative of the total sediment transport rate with respect to the volume of fine sediment in the active layer:

χ

1= 1 u

∂

qb

∂

Ma1. (22) Theparameter

μ1,1

[-]isdefinedas:

μ

1,1=d1,1− f1I, (23)

where d 1,1[-]isanondimensionalmeasureofthederivativeofthe sedimenttransportrateofthefinefractionwithrespecttothe vol-umeoffinesedimentintheactivelayer:

d1,1= 1 u

χ

1

∂

qb1

∂

Ma1. (24) We obtain the eigenvalues of the system matrix finding the rootsofitssecond degreecharacteristicpolynomial.The eigenval-uesarenondimensionalizeddividingbytheflowvelocity:

λ

alS2i= 1 2



λ

b+

λ

s1±



Δ

alS2



fori=1,2, (25) wherethediscriminantis:

Δ

alS2=

(

λ

b−

λ

s1

)

2+4

λ

b

λ

s1

_μ

γ

1

1,1.

(26) The eigenvalues of the system carry coupled information on both the bed elevation and the surface grain size distribution whichshowsthataperturbationinbedelevationcausesa pertur-bationinsurfacegrainsizedistributionandviceversa(Section1). Yet, we identify two nondimensional celerities that approximate thechangesinbedelevation(

λ

b) andinsurfacegrainsize

distri-bution(

λ

s1)independently.

Thebedcelerity,whichisindependentoftheactivelayer thick-ness,wasfirstderivedbydeVries(1965)forunisizesediment:

λ

b=₁

ψ

− Fr2. (27)

Wedefinethenondimensionalsortingcelerityas:

λ

s1=

χ

1

μ

1,1. (28)

ThissortingceleritydiffersfromtheoneofRibberink(1987)as he considered a perturbation in the mean grain size while here the sorting celerity relates to a perturbation inthe volume frac-tioncontentofeachgrain sizefractionindividually.Theproposed expressionforthesortingcelerityinEq.(28)isageneralizationof theexpressionproposedbySteccaetal.(2014),aswehaverelaxed Stecca’sassumptionoflimitedhiding.

Themathematical characterofthemodeldependsonthesign ofthediscriminant

ΔalS2

,Eq.(26).If

ΔalS2

> 0thetwoeigenvalues arerealandthesystemishyperbolic. If

ΔalS2

<0 theeigenvalues arecomplexandthesystemiselliptic.Alargedifferencebetween the bed celerity and sorting celerity reduces the likelihood that themodel becomeselliptic.Hyperbolicity isguaranteed if

γ

1>0. If

γ

1<0andthebedandsortingceleritiesareequal,ellipticityis guaranteed(Sieben,1994;Steccaetal.,2014).

Assumingthatreversemobilitydoesnotoccur(Section2.4), c 1 islargerthanthevolumefractioncontentoffinesedimentinthe activelayer (F a1)duetothegrain sizeselectivityofthe sediment transportrelation(Steccaetal., 2014). Ifwe alsoassume thatthe sedimenttransferred to the substratein aggradational conditions has the same grain size distribution as the active layer (Hirano, 1971), then the parameter

γ

1 is always positive in aggradational conditions. Only a substrate finer than the active layer yields a negative value of the parameter

γ

1. Thus, a two-fraction active layermodelcan onlybeill-posed ifthebeddegradesintoa sub-strate that is finer than the active layer (a result also found by Steccaetal.(2014)consideringunsteadyflow).

InSections4.1and4.2weassesstherelaxationofthe assump-tions that reverse mobilitydoes not occur andthat the aggrada-tionalfluxtothesubstratehasthesamegrainsizedistributionas theactivelayer.

3.2. Steady vertically continuous model consisting of two size fractions

We apply the same procedure used to analyze the active layer model to the vertically continuous model (Section 2.3). In this manner we obtain the discriminant of the eigenvalues (AppendixD):

Δ

vcS2=

(

λ

b−

λ

sc1

)

2+4

λ

b

λ

sc1

g1 m1,1,

(29) where

λ

sc1, g 1,and m 1,1 aretheequivalentsto

λ

s1,

γ

1,and

μ1,1

of theactivelayermodel(AppendixD).

Similar tothe activelayer model(Section 3.1), the continuous modelishyperbolicandwell-posedif

ΔvcS2

> 0andviceversa. Al-thoughtheexpressionofthediscriminantofthevertically continu-ousmodel,Eq.(29),issimilartotheoneoftheactivelayermodel, Eq.(26),there isan essential difference betweenthe two. Inthe active layer model the discriminant is a function of the stream-wise position,

ΔalS2

(x ), yetinthecontinuousmodelthe discrimi-nantisalsoafunctionoftheverticalcoordinate,

ΔvcS2

(x, y ).Thus, ellipticityorhyperbolicityisapropertynotonlyofthestreamwise coordinatebutalso oftheelevationin thebed(Section2.3). Hy-perbolicityisguaranteed if g 1>0but,contrarytotheactivelayer model, this parameter can be negativeboth under aggradational anddegradationalconditions.

Dueto grainsize selectivetransportwe canassure that,if re-verse mobility conditions do not prevail, the concentration c 1 is largerthan thevolumefractioncontent representativeofthebed surface F b1,Eq.(C.4).However, F b1isaweightedaverageofall sed-imentandforthisreasonthereisnoguaranteethatforallbed ele-vationstheaveragevolumefractioncontent(F b1)islargerthanthe

localvolume fractioncontentinthe bedsediment(f 1).Moreover, asthereisnodistinctionbetweenaggradationalanddegradational cases,the domain inwhich themodel islikely tobe ill-posed is larger than fortheactive layer model.Thepresence offine sedi-ment atthelocations havinglarger probability ofentrainmentin combination witha “smooth” verticalvariation (small derivative) ofthevolumefractioncontentoffinesedimentreducesthe likeli-hoodofthemodelbecomingelliptic.

4. Activelayermodelparameterstudy

In thissection we assessthe effects ofvarious model param-eterson themathematicalcharacter oftheactivelayer model.To thisendwe studytheanalyticalexpressions oftheeigenvaluesof thesteadymodelconsideringtwosedimentsizefractionsobtained inSection3.1.

4.1. Hiding

Giventhefactthatill-posednessariseswhenconsidering differ-entgrainsizesinthemixtureandthatalargerdifferencebetween grainsizesincreasestheill-poseddomain,intuitivelyhidingshould reducethelikelihoodofill-posedness.Itseffect,however,is oppo-siteaswewillshowhere.

Toexplain thiscounter-intuitiveresultweanalyzethe termin the characteristic polynomial intrinsically related to hiding. This term isthe derivative of the sedimenttransport rateof fine and coarse sediment with respect to the volume of fine sedimentin the activelayer (

∂

q bk/

∂

M a1). Itcan be consideredasthe summa-tionoftwoterms:

∂

qbk

∂

Ma1 = 1 La

⎛

⎜

⎜

⎝

∂

∂

FFaka1 Qbk



presence +Fak

∂

Qbk

∂

Fa1



hiding

⎞

⎟

⎟

⎠

fork=1,2. (30)

We name the first and second terms on the right-hand side the “presenceterm” andthe“hidingterm”,respectively.The“presence term” explainsthat an increasein the volume fraction content of thefine sedimentintheactive layerimpliesboth(1)an increase ofthesedimenttransport rateofthe finefractionasitspresence atthebedsurfaceislarger,and(2)a consequentdecreaseofthe sedimenttransportrateofthecoarsefractionbecauseitspresence atthebedsurfacedecreases.The “hidingterm” indicates thefact that a variation of the volume fraction content of the fine sedi-mentchangesthesedimenttransport rateofbothfineandcoarse fractions dueto achangeinthe meangrainsize ofthesediment mixture.The “presence term” ispositive forthefine fractionand negativeforthecoarsefraction.The“hidingterm” isalways posi-tive.

In a situation where hiding is negligible, an increase of the characteristic size of the coarse fraction or decrease of the fine fraction, which is associated with a larger likelihood that the model is elliptic,causes an increase of the “presence term” and thus of

∂

q _bk/

∂

M a1. With respect to such a situation, hiding de-creasesthe“presenceterm” ofbothfineandcoarsesediment (re-ducingthelikelihoodofellipticity)butintroducesthepositive con-tributionofthe“hidingterm”.Overall,the“hidingterm” may dom-inate,whichincreasesthevalueof

∂

q bk/

∂

M a1andthusofthe like-lihoodofellipticity.

Interestingly, in degradational conditions into a fine substrate (asituationprone tobeelliptic)thehidingtermdominates.Thus, hidingincreasesthelikelihoodthatthemodelisellipticin degra-dationalconditionsintoasubstratefinerthantheactivelayer.

InFig.3a weshow theeffectofhiding onthediscriminant of the steady active layer model considering two size fractions, Eq.

Table 2

Values of the reference case. For these values the active layer thickness can be seen as representative of plane bed conditions ( L a ≈ 2.5 D 90 ) as well as bedform

dominated conditions ( L a ≈ a L h bL ).

Fa1 [-] f1I [-] d1 [m] d2 [m] Cf [-] u [m/s] h [m] La [m]

0 0.6 0.002 0.004 0.015 0.68 0.20 0.01

(26).WeconsiderthereferencecasedescribedinTable2.The sed-iment transport rate is computed using the relation derived by Meyer-Peterand Müller(1948).Toobtain differentvalues of hid-ingwe varyparameter b inthe powerlawhiding functioninEq. (7) between0 and 1 (purple line in Fig. 3a). The yellowline in Fig.3aisobtainedvarying thecharacteristicgrainsizeofthefine fractionbetween0.001m and0.004musingthe Egiazaroff hiding

relation,Eq.(C.5).Thediscriminantdecreasesforincreasinghiding independentfromthehidingfunction.

Besidesthesecasesunderdegradationalconditions,wemay en-counter problems evenunder aggradation. In fact, ifhiding is so strongthat reverse mobilityisinduced, then one ofthe assump-tions ofthe analysisby Stecca etal.(2014),maynot be fulfilled. In detail,it may happen that the referencecontent c 1 related to thefinesedimentinthebedload,Eq.(21),isnot greaterthanthe contentoffinesintheactivelayer F a1,whichwastheirassumption undergrainsizeselectivetransport.Whenreversemobilityinstead determinesconditionssuchthat c 1<F a1,thenthediscriminant,Eq. (26), may be negative, and the model maybecome elliptic even underaggradationalconditions.

4.2. Aggradational flux to the substrate

The sediment transferred to the substrate under aggra-dational conditions using the model by Hoey and Ferguson (1994)(Section 2.4) isalways finer than the sedimentinthe ac-tivelayer.Thisisbecausethebedloadisfinerthanthebedsurface duetograinsizeselectiveprocesses(providedthatreverse mobil-itydoesnotdominate).Thus,applicationofthemodelbyHoeyand Ferguson (1994) implies that under aggradational conditions the interface betweentheactive layerandthe substrateisfiner than theactivelayer.Thismeansthatthecondition

γ

1>0(Section3.1) maynotbefulfilledunderaggradationalconditions,whichimplies that the modelmay become ill-posed. Therefore,a larger contri-bution of the bedload to the aggradational flux to the substrate (smallervalueoftheparameter

α

inEq.(8))impliesalarger like-lihoodofthemodelbecomingelliptic(AppendixBinthe supple-mentarymaterial).

However, in a hypotheticalaggrading casein which thegrain sizedistribution transferred tothe substrateisfullycomposed of bedloadsediment(

α

=0), therelative contentofthefinefraction intheverticalsedimentflux,

γ

1(Eq.(20)),thatcontrolsthesizeof theill-poseddomain(Section3.1),isstillnotassmallasitcanbe found under degradationalcases (Appendix B inthe supplemen-tarymaterial).Thus,ill-posedcasesareexpectedtooccurprimarily underdegradationalconditionsintoafinesubstrate.

4.3. Prefactor in a sediment transport relation and morphodynamic factor

The discriminant (

ΔalS2

) of the steady active layer model for two size fractions, Eq. (26), can be written as

ΔalS2

=A 2

ΔalS2

_, where 

alS2

is the discriminant for a unit prefactor (i.e., A ₌1). The prefactor A increases ordecreases the discriminant butdoes notchangeitssign,andsoitdoesnotchangethecharacterofthe mathematical system. This is confirmed by Fig. 3b,which shows theeffectofvaryingtheprefactor A inthereferencecasedescribed inSection4.1.

1 2 3 hiding factor of the fine fraction

1 [-] -5 0 5 discriminant alS2 [-] 10-4 a hyperbolic elliptic Egiazaroff Power Law 4 6 8 10 prefactor A [-] b hyperbolic elliptic MPM EH 1 2 3 power B [-] c hyperbolic elliptic reference 0 0.02 0.04 0.06 critical Shields stress

c [-]

d hyperbolic

elliptic

Fig. 3. Discriminant alS2 , Eq. (26) , as a function of: ( a ) the hiding for the fine fraction ξ1 , ( b ) the prefactor of the sediment transport formula A , ( c ) the power B , and ( d )

the critical Shields stress θc . The blue line is obtained varying the parameters from the reference state using the Meyer-Peter and Müller (1948) (MPM) sediment transport

relation. The red line is obtained using the Engelund and Hansen (1967) (EH) relation. The yellow line is obtained by varying the characteristic grain size of the fine fraction using the hiding relation by Egiazaroff (1965) , Eq. (C.5) , in combination with the MPM relation, the purple line by varying the coefficient b of the power law function by

Parker et al. (1982) , Eq. (7) , in combination with the MPM relation. The dots represent the reference situation described in Table 2 . Note that there is a different reference situation depending on the sediment transport relation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Sincemorphodynamic timescales areusuallyseveralordersof magnitudelarger than the time scales of the flow (Section 2.1), computations usually cover a significant number of years. The computationaltimeissometimesreducedusingamorphodynamic factorthatmultipliesthedivergenceofthesedimenttransportrate (Latteux,1995;Roelvink,2006;Ranasingheetal.,2011).Thisfactor canalsobe considered asa multiplicationofthe sediment trans-portrateandthereforehasthesameeffectastheprefactorA .Thus, theuse of a morphodynamic factor doesnot changethe mathe-maticalcharacterofthemodel.

This result is obtainedassuming quasi-steady flow. While the prefactorinthe sedimenttransport relationrarelyvariesby more than an order of magnitude, simulations may be run with mor-phodynamicaccelerationfactors O

(

102

₎

_{.Intheselattercases,the} quasi-steadyflowassumptionmaynotbeacceptable,whichlimits theextensionofouranalysis.

4.4. Exponent and critical Shields stress in a sediment transport relation

Thediscriminant,Eq.(26),tendsto0−withincreasingvaluesof

B iftheeffectiveShieldsstressforallsedimentfractionsissmaller than1,orto∞iftheeffectiveShieldsstressislargerthan1forat leastonefraction.Thus,itisdifficulttogeneralizetheeffectofthe exponent.Its variation froma reference situationcan both make thesystemhyperbolic ifthe referencesituationis ellipticorvice versa.InFig. 3cwe show thediscriminant asa functionof B for thesamereferencecasesasinSection4.3.Thehyperbolicsituation whenusingMeyer-PeterandMüller(1948)becomesellipticifthe valueof theexponent B increasestowards thevalue in Engelund andHansen(1967).

TheeffectofthecriticalShieldsstress,

θ

cinEq.(6),onthe

dis-criminantissimilar totheeffect oftheexponent,asitsvariation canboth makea previouslyhyperbolic caseellipticorvice versa. Fig. 3d shows how a decrease of the critical shear stress when usingthe sedimenttransport relationby Meyer-Peter and Müller

(1948)increasesthe discriminantreducing the likelihoodof ellip-ticbehavior.

4.5. Active layer thickness

Thediscriminantoftheeigenvalues,Eq.(26),canbewrittenas aseconddegreepolynomialoftheinverseoftheactivelayer thick-ness, i.e.,

ΔalS2

₌a 1

(

1/La

)

2+a 2

(

1/La

)

+a 3 where a 1>0, a 2, and

a 3 arecoefficients independentofthe active layerthickness.This impliesthat:(1)themodeliswell-posedforasufficientlythin ac-tive layer, (2)the modelis well-posedfor asufficiently thick ac-tive layer,and(3)there existsoneill-posed domainonly (regard-ing the active layer thickness).These resultsof the two-fractions modelconfirmpreviousresultsbasedonthesimplifiedactivelayer model(Ribberink,1987;Sieben,1994).

The inverse ofthe roots ofthe second degree polynomial are thelimit valuesoftheactive layerthicknessthat ensurethat the modeliswell-posed: L±_a =

λ

s1 

λ

b

⎛

⎝

1− 2

γ

1 

μ

1,1±

!

2

γ

1 

μ

1,1− 1

"

2 − 1

⎞

⎠

−1 , (31)

wherewehaveusedthenotation()forthevariableswithunit ac-tivelayerthickness(i.e., L a=1m).Giventhefactsthat theactive

layerthicknessisoneofthemostempiricalparametersofthe sys-temofequations(Section1)andthatrivermorphodynamicmodels oftenrequirecalibration(e.g.,CaoandCarling,2002b),Eq.(31)can beappliedtoselectacertainvaluefortheactivelayerthicknessto avoidasituationthatispronetobeill-posed.

5. Consequencesofill-posedness

Inthissectionweanalyzetheconsequencesofill-posedness us-ing numerical simulations. Ouraim is to provide modellerswith thetoolstodetectoccurrenceofill-posednessintheir resultsand understandhowthe observedunrealisticmodelbehavior changes

Table 3

Overview of the simulations. Only the parameters that are different between simulations are shown.

Simulation La [m] hiding Mf [-] Fs1 z Math. character

1 0.010 No 1 0.6 0.10 Hyperbolic 2 0 . 025 No 1 0.6 0.10 Elliptic 3 0.010 Yes 1 0.6 0.10 Elliptic 4 0.010 Yes 2 0.6 0.10 Elliptic 5 0.010 No 1 0.6 0 . 01 Hyperbolic 6 0.010 No 1 1 . 0 0.10 Elliptic 7 0.010 No 1 1 . 0 0 . 01 Elliptic

withparameterandmodelchoices.Firstwe makenumericalruns toqualitativelyobservetheconsequencesofthenon-lineareffects that are neglected in the perturbation analysis (Section 5.1). In Section5.2weconductasensitivityanalysistogeneralizethe con-sequencesobservedintheprevioussection.

5.1. Numerical examples

ThelinearperturbationanalysisshowninSection2.5indicates that perturbations grow unboundedly ifthe model is elliptic. In this section we run four numericalsimulations at flumescale to analyze the effects of the neglected non-linear terms. The simu-lations areone-dimensional fortheflow andarecomputedusing theDelft3Dsoftwarepackage(Lesseretal.,2004),whichsolvesthe unsteadyShallowWaterEquationsincombinationwiththeactive layer model. Forsimplicity the active layer thickness is assumed constant.Underaggradationalconditionsthesedimenttransferred tothesubstrateiscomposedofonlytheactivelayersediment(i.e.,

α

=1 in Eq. (8)). Substrate stratigraphy is stored using a book-keepingsystem(Viparellietal.,2010;Steccaetal.,2016). All sim-ulations start fromequilibriumconditions undercoarse sediment feeding. A lowering of the base level is imposed, which causes degradation into a fine substrate. We consider a well-posed ref-erencecase(Simulation1)that(initially)hasthesameparameters asthereferencecaseoftheparameterstudyofSection4(Table2). Then,theactivelayerthicknessischanged(Simulation2)and hid-ing isconsidered(Simulation 3).Simulation 4isequalto Simula-tion 3 except for its morphodynamic factor. Table 3 summarizes thedifferencesbetweenthefoursimulations.Theboundary condi-tionsthatareinequilibriumwiththeinitialcondition(Blometal., 2016)aswellasotherparametersaredescribedinTable4.

In Fig. 4a we plot the discriminant of the eigenvalues of the quasi-steady active layer model,Eq. (26), at theinitial time asa function oftheactive layerthickness.Notethat theconditionsof Simulation 1yielda well-posed model(

ΔalS2

> 0) whilethe con-ditions of Simulation 2 yield an ill-posed model (

ΔalS2

<0). In Fig. 4d-e we show the evolution of the bed elevationfor Simu-lations1 and2,respectively.Theill-posed Simulation2showsan oscillatory behaviorthat isnot presentinthewell-posed Simula-tion 1. We havenot imposed anyinitial perturbation, which im-plies that numerical noise is sufficient to trigger the oscillatory behavior. Simulation 3is thesame asSimulation 1,yetthe sedi-menttransportratenowaccountsforhidingusingtheParkeretal. (1982)function, Eq.(7),withexponent b equalto 0.8.Simulation 3isill-posed(Fig.4b)andthesolutionshowsoscillations(Fig.4f) justasintheill-posedSimulation2.

Table 5

Values of the physical and numerical parameters varied in the sensitivity analysis.

Parameter Values Physical d1 [m] 0.0 0 05, 0.0 01, 0.0 02 d2 [m] 0.0 02, 0.0 03, 0.0 04 La [m] 0.010, 0.015, 0.020, 0.050 fs1 [-] 0.6, 0.8, 1.0 Numerical x [m] 0.1, 0.2 z [m] 0.01, 0.02, 0.05, 0.10

Simulation4isthesameasSimulation3exceptforits morpho-dynamicfactorequaltotwo,whichdecreasesthevalueofthe dis-criminant(Fig.4c),causingoscillationstodevelopfaster(Fig.4g).

Forall casesoscillations do not occur atthe upstream endof thedomain.Thisisbecauseoscillationsrequiretime(andsospace) togrow.Inallcasesoscillationsgrowuntilamaximumamplitude isreachedandthenpropagatedownstream.Thismaximum ampli-tudeis such that theconditionsare atthe brink ofill-posedness andwell-posedness.Worded differently,downstreamfromthe lo-cation where the amplitude is maximum the model is ill-posed anditiswell-posedupstreamfromit.

Theoscillationsareassociatedwithdegradationandsubsequent aggradation.The depositedsedimenthasthesamegrainsize dis-tributionastheactivelayer(whichiscoarserthantheinitial sub-strate), so the overall effect of an oscillation is a coarsening of the topmostpart ofthe substrate. This coarsening actsasa reg-ularizationmechanism, which not onlyrestoreshyperbolicity but alsodampensoscillationsthatarrivefromupstreambylimitingthe sourceoffinematerial.

Itislikelythat,becauseoftheregularizationmechanism, com-putations do not crash. As a result the extent and likelihood of ellipticity mayin practicebe underestimated. Yet,the resultsare physicallyunrealisticandimplementinganautomatedcheckofthe eigenvalues would be good practice for software developers and users.

5.2. Sensitivity analysis

Theprevious section hasshownthat, dueto thenon-linearity ofthesystem,anill-posedsimulationgeneratesnon-physical oscil-lationsthatpropagatedownstreamandgrowuntilacertain maxi-mumamplitudeatwhichthemathematicalproblemisatthebrink ofill-posednessandwell-posedness.Inthissection weruna sen-sitivityanalysistogeneralizethoseresults.

Tothis end, we vary 4physical parameters and 2parameters relatedtothe domaindiscretization, usingSimulation 1asa ref-erencecase.Table5summarizestheparametervaluesusedinthe sensitivityanalysis.

As we are interested in studying the behavior of simulations underill-posedconditionsweexcludefromtheanalysisthose sim-ulations in which the combination of parameters yield a well-posedmodel.Aswehaveobservedthattheoscillationsneedspace to grow until a maximum value (Section 5.1), we exclude those simulationsinwhichthedomainisnotlongenoughtodevelopan oscillationthattravelswithaconstantamplitude.Asetof173out of256simulationsfulfillsthesetworequirements.

Table 4

Domain definition, boundary conditions, and numerical parameters. The symbols not defined in the text are: reach length ( L ), channel width ( B ), simulation time ( T ), lowering rate of the downstream water level (low. rate), horizontal discretization length ( x ), and time step ( t ).

L [m] B [m] qb1 [m 2 /s] q b2 [m 2 /s] q [m 2 /s] T [h] low. rate [m/h] x [m] t [s]

0.02 0.04 0.06 0.08 0.1 active layer thickness L_a [m] -1 0 1 discriminant alS2 [-] 10-3 Sim. 1 Sim. 2 a hyperbolic elliptic 1 1.5 2

hiding factor of the fine fraction ₁ [-]

Sim. 1 Sim. 3 b hyperbolic elliptic 1 2 3 4 morphodynamic factor M_f [-] Sim. 3 Sim. 4 c hyperbolic elliptic 0 0.1 0.2 0.3 bed elevation [m] Sim. 1 d 0 0.5 1 1.5 time [h] Sim. 2 e 0 20 40 60 80 streamwise position [m] 0 0.1 0.2 0.3 bed elevation [m] Sim. 3 f 0 20 40 60 80 streamwise position [m] Sim. 4 g

Fig. 4. Discriminant ΔalS2 , Eq. (26) , as a function of ( a ) the active layer thickness L a , ( b ) the hiding function of the fine fraction ξ1 , and ( c ) the morphological factor M f .

Numerical solutions of ( d ) Simulation 1, ( e ) Simulation 2, ( f ) Simulation 3, and ( g ) Simulation 4. See Table 3 for the parameters definition.

In Fig.5awe plotthemaximumflowdepth(h max)

nondimen-sionalizedwith thenormal flow depth (h n), as a function of the

discriminant,Eq.(26).Forthesakeofclarityweplottheresultsof thesimulationswitha horizontaldiscretizationlength (



x ) equal to0.1mandathicknessofthesubstratelayersequalto0.01mand 0.10m(seeAppendixCofthesupplementarymaterialforthe re-sultsofall simulations).The parametersusedtoevaluatethe dis-criminantare those at the start ofthe simulation (normal flow). Theverticalblacklinesjoinsimulationswiththesamephysical pa-rameters(i.e.,theyonlydifferregardingnumericalparameters)and thecolor ofeach dot is related tothe thickness ofthe substrate layer.Thelinearanalysishasshownthatthegrowthratedepends onthediscriminant (Section2.5),howeverthereisonlymild cor-relationbetweenthediscriminantandthemaximumamplitudeof theoscillations.

Athinlydiscretizedsubstrateisassociatedwithalarger ampli-tudeoftheoscillations(Fig.5a).Thiseffectcanbeseenonly em-piricallysinceitisnotaparameterofthesystemofequationsnor doesitappearinthelinearstabilityanalysis.Forallsimulationswe computetheflowdepththatyields avalueofthediscriminantat theinitialconditionequalto0(i.e.,atthebrinkbetween elliptic-ityandhyperbolicity).Thisisdonenumericallyfindingtherootof

ΔalS2

(h ), Eq.(26),considering thewaterdischargeandvolume of

sedimentintheactivelayerandatthesubstrateoftheinitial con-dition.We termthisflowdepth thehyperbolic flowdepth (h _hyp), whichisindependentofthenumericalparametersofthe simula-tionanddependsonphysicalparametersonly.InFig.5bwe

com-parethe measuredmaximumflowdepthandthehyperbolic flow depth.Thegreylinerepresentsthesituationinwhich h max= h hyp.

Weseethatthehyperbolicflowdepthcanbeusedasarough es-timate ofthe maximum flow depththat willoccur in an elliptic simulation.Oneimportantsourceofscatteristhefactthatthe hy-perbolicflowdepthdependsontheinitialconditiononly,whereas themaximumflowdepthalsodependsontheevolutionofthe so-lutionastheoscillationsinteractwitheachother.

InFig.5cweplotthenondimensionalmaximumflowdepthas a function ofthe streamwise location where the maximum flow depthoccurs(nondimensionalizedwiththetotallengthofthe do-main).Forthesakeofclarityweplottheresultsofthesimulations witha thicknessofthe bookkeepinglayers (



z ) equal to0.01 m (see AppendixC ofthe supplementary materialforthe resultsof all simulations).Inthe thinlydiscretized simulations we findthe maximumamplitudemoreupstreamcomparedtothecoarsely dis-cretizedsimulations.Thelocation wherethemaximumamplitude of the oscillations is found is related to its growth rate since a fastergrowingoscillationdevelopsitsmaximumamplitudeinless distancethanaslowerone.Thisresultconfirmsthefindingsof lin-earstabilityanalysis(Section2.5).

Also the maximum flow depth and the domain discretization aremildlycorrelated.Similartotheverticaldiscretization,smaller cellsyieldalargermaximumflowdepthbutthiseffectcanbeseen onlyempirically.

The discretizationof thesubstrate affectsalso theduration of the ellipticbehavior. Fig.6 showsthe longitudinalprofile offour

-8 -7 -6 -5 -4 -3 -2 -1 0 alS2 [-] 10-4 1.2 1.4 1.6 1.8 h max /h n [-] a z=0.01 m z=0.10 m 1.25 1.3 1.35 1.4 1.45 1.5 1.55 1.6 1.65 h_hyp/h_n [-] 1.2 1.4 1.6 1.8 h max /h n [-] b 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 x_max/L [-] 1.2 1.4 1.6 1.8 h max /h n [-] c x=0.1 m x=0.2 m

Fig. 5. Maximum flow depth h max (nondimensionalized by the normal flow depth h n ) that develops as a consequence of ellipticity. In ( a ) the maximum flow depth is plotted

against the discriminant ΔalS2 for a horizontal discretization length ( x ) equal to 0.1 m and for a thickness of the bookkeeping layers ( z ) equal to 0.01 m (red dots) and

0.10 m (blue dots). The vertical black lines connect two simulations in which all parameters but z are the same. In ( b ) the maximum flow depth is plotted against the

hyperbolic flow depth ( h hyp ) nondimensionalized with the normal flow depth. Each black dot is the result of a simulation. The black lines connect simulations with the same

physical parameters and the grey line is the perfect agreement. In ( c ) the maximum flow depth is plotted against the distance from upstream at which the maximum flow depth occurs ( x max ) nondimensionalized with the length of the domain ( L ) for a thickness of the bookkeeping layers ( z ) equal to 0.01 m and for a horizontal discretization

length equal to 0.1 m (orange dots) and 0.2 m (green dots). The black lines connect two simulations in which all parameters but x are the same. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

simulations from the sensitivity analysis at the end of the run (Simulations1,5,6,and7,seeTable3).Thesubstrateofthe well-posed simulations at the end of the runs is unaltered, whereas it hascoarsened in theill-posed casesdueto the oscillatory be-havior, whichacts asaregularization mechanism(Section 5.1).A thinlydiscretizedsubstrateenhancestheregularizationmechanism asone fulllayer iscreatedwiththegrainsize distributionofthe (coarse)activelayerduringtheaggradingphaseoftheoscillation. Ifthebedisdiscretizedintothicklayersthematerialtransferredto thesubstratewill beaveraged withthesedimentalreadypresent in thetop substratelayer. The resultinggrain size distributionof thesubstratemaynotbesufficientlycoarse topreventthemodel frombeingill-posed.

6. Implicationsofconsideringmorethantwosizefractionsor anunsteadyactivelayerthickness

Toobtain an analytical expression ofthe eigenvalueswe have restricted ouranalysisto mixtures ofsediment composedof two sizefractionsandsteadyflow. Inthissectionwe explorethe con-sequencesofrelaxingtheseassumptions.

6.1. Ill-posed domain of a three-size-fractions case

Amodelforthreesedimentsizefractionsistoocomplexto ob-tainanalyticalexpressionsoftheeigenvalues.Wethereforefirst at-tempt to provideinsight addressing a specific casebased on the resultsofthecasefortwosizefractions.

The concept of a finer or coarser active layer relative to the substrateisunequivocallyapplicable inthecaseoftwo size frac-tions.However,thisconceptisnotasstraightforwardforthreesize fractions,asitrequiresthedefinitionofa meangrain size.Asan extensionofthe resultsforthetwosize fractionscasewherethe modelcanonlybeill-posedindegradationalconditionsintoa sub-strate finerthan theactive layer (assumingcertain conditions on the closure relations, Section 3.1), we consider a situation with three size fractions which, regardless of the method to compute themean grain size, isgoverned by degradation into a substrate coarserthantheactivelayer.Thishappens,forinstance,ifthe vol-umefractioncontentsintheactivelayerofthefine,medium,and coarse size fractions are 0.5, 0.5, and 0,respectively; and at the interfaceare0.5,0,and0.5.

Weconsider a sedimentmixturewiththe abovevolume frac-tion contents and characteristic grain sizes of the fine, medium,

0.05 0.1 0.15 0.2 elevation [m] Simulation 1 (well-posed) a Simulation 5 (well-posed) b 2 2.5 3 3.5 4

mean grain size [m]

10-3 35 40 45 streamwise position [m] 0.05 0.1 0.15 0.2 elevation [m] Simulation 6 (ill-posed) c 35 40 45 streamwise position [m] Simulation 7 (ill-posed) d

Fig. 6. Grain size stratification at the end of: ( a ) Simulation 1, ( b ) Simulation 5, ( c ) Simulation 6, and ( d ) Simulation 7. The substrate of the ill-posed simulations coarsens unrealistically.

Fig. 7. Results of an ill-posed simulation with 3 grain size fractions under degradational conditions into a substrate coarser than the active layer: ( a ) bed elevation at selected times, ( b ) surface mean grain size with time, and ( c ) grain size stratification at the end of the simulation.

and coarse fraction equal to 0.001m, 0.003m, 0.005m. All the otherparametersareequaltothereferencecase(Table2).This sit-uationisellipticastwooftheeigenvaluesofthesystemmatrixare complex.Thus,inathreesizefractionscasethemeangrainsizeof thesedimentintheactive layerrelativetothatattheinterface is notavaliddiscriminantofthemathematicalcharacterofthe sys-temofequations.

Fig.7showstheresultsofanumericalsimulationbasedonthe aboveparameters.Thesolutionpresentsoscillationsasinthe pre-viousill-posedcases.However,theamplitudeoftheseoscillations isnowsignificantlysmaller(compareFig.7atoFig.4e).Arelatively large oscillation appears after approximately 3 h which entrains coarsesedimentfromthesubstrate(Fig.7b).Duringthe aggrada-tionalphaseoftheoscillationsfinesedimentfromtheactivelayer

istransferred tothesubstrate. Thus, attheendofthesimulation thetop partofthesubstrateisfinerthan initially(Fig.7c).To il-lustratetheimplicationsofthisresultwestudytheeffectsof dis-cretizing thesame sedimentmixture intotwo or threesediment fractions.Todiscretizethesampleintothreesedimentfractionswe usecharacteristicgrainsizesequalto0.001m,0.003mand0.005m andtodiscretizeitintotwosedimentfractionsweuse0.002mand 0.004 m.The volumefraction contentin themedium size ofthe threefractionmixtureisequallysplitbetweenthefineandcoarse binsofthetwofractionmixture.Wevaryallvolumefraction con-tentsbetween0and1 toobtaindifferentsedimentmixturesand theflow depth (keepingthewater dischargeperunit width con-stant)between0.15mand1.5mtoobtaindifferentflowconditions. Allotherparametersareequaltothereferencecase.

0 0.2 0.4 0.6 0.8 1 Fr [-] -4 -3 -2 -1 0 1 2 3 4 D ma -D mI [m] 10-3 2 fractions a 0 0.2 0.4 0.6 0.8 1 Fr [-] 3 fractions b hyperbolic elliptic

Fig. 8. Ill-posed domain for degradational cases in which the sediment mixture is discretized into ( a ) two and ( b ) three size fractions as a function of the Froude number ( Fr ) and the difference between the mean grain size of the sediment in the active layer ( D ma ) and at the interface between the active layer and the substrate ( D mI ). When

the mixture is discretized into three size fractions the model may be ill-posed under degradational conditions into a substrate coarser than the active layer.

Fig. 8a-b shows the ellipticdomain when the mixture is dis-cretizedintotwoandthreesizefractions,respectively.Onthe ver-tical axiswe plot the difference betweenthe meangrain sizeof thesedimentintheactivelayer(D ma[m])andattheinterface

be-tweentheactivelayerandthesubstrate(D mI [m]).Notethatsome

situationsthatarewell-posedwhenthemixtureisdiscretizedinto two fractions are ill-posed when it isdiscretized into three frac-tions. Wecannotprove thattheeigenvalues ofthe systemmatrix forthreesizefractionsarealwaysrealunderaggradational condi-tions due to the complexityofthe expressions. Nevertheless, we have not obtaineda singlecomplexvalue in anyofthe aggrada-tionaltestswehaveconducted.

6.2. Effect of an unsteady active layer thickness in the ill-posed domain

In this section we analyze the implications of considering a variableactivelayerthicknesswithrespecttotheill-posednessof thesystemofequations.Theflowneedstobeconsideredunsteady tostudythevariabilityoftheactivelayerthicknessifitisrelated todunegrowth(Section3).Wesimplifythesystemassumingtwo sediment size fractions and negligible hiding (Stecca etal., 2014 andAppendixDofthesupplementarymaterial).

We obtain thecharacteristicpolynomial of thesystemmatrix, Eq. (11). We prove that in aggradational and degradational con-ditions into a substrate coarser than the active layer the charac-teristic polynomial has 5 real roots (Appendix D of the supple-mentarymaterial).Thereforethemodeliswell-posedregardlessof the unsteady activelayer thickness.Regarding degradational con-ditions into asubstrate finerthan the active layerwe prove that if

λ

b>

λ

s1F a1/ f 1I, considering a variable active layerthickness in-creases the likelihood of the model becoming elliptic. Note that, assuming a similar order of magnitude of the bed and sorting celerities, inconditionsprone to be elliptic(i.e., degradationinto a substratesignificantlyfiner thantheactive layer) avariable ac-tivelayerthicknessincreasesthedomaininwhichtheactivelayer model is elliptic. We numerically test severalsets of parameters andwe find no casewherethe modelis hyperbolic ifthe active

layeris unsteady butellipticif itis constant (AppendixDof the supplementary material). This suggests that, although we do not providea formalproof, an unsteadyactive layerthicknessalways increasesthelikelihoodofthemodelbeingill-posed.

7. Conclusions

Wehaveassessedthewell-posednessoftheequationsusedto model mixed sediment river morphodynamics. In particular we have studied the system formed by the flow equations ( Saint-Venant,1871)together withtheactive layermodel(Hirano,1971) andasimplifiedverticallycontinuousmodel(Viparellietal.,2017). Ourfindingsarethefollowing:

• Considering two size fractions and the quasi-steady flow as-sumption we obtainan analytical expression forthe discrimi-nant thatdetermines whethertheactive layerandcontinuous modelsareill-posed.

• Assuming(i)twosizefractions,(ii)steadyflow,(iii)noreverse mobility,and(iv)thatunderaggradationalconditionsthe depo-sitionalflux ofsedimentto thesubstrateisentirelycomposed ofactivelayersediment,theactivelayermodelcanbeill-posed under degradationalconditionsinto asubstrate finerthan the activelayeronly.

• The use of a hiding factor increases the likelihood of ill-posedness. Strong hiding that causes reverse mobility may causeill-posednessalsoinaggradationalconditions.

• Aggradationalcasesmaybeill-posedifthedepositionalfluxof sedimenttothesubstrateincludesbedloadsediment(Hoeyand Ferguson,1994).

• Theactivelayermodelmaybeill-posedindegradational condi-tionsintoasubstratecoarserthantheactivelayerifmorethan twosizefractionsareconsidered.

• Considering a variable active layer thickness associated with dune growth increases the likelihood that the active layer modelisill-posed.

• Thesimplifiedverticallycontinuousmodelcanbeill-posed un-der both aggradational and degradational conditions. A small