A reduced-complexity model for river delta formation – Part 1: Modeling deltas with channel dynamics

We develop a reduced-complexity model (RCM) delta formation model, in contrast to reductionist models based on high-resolution computational ﬂuid dynamics. The basic framework of this model (referred in this paper as “DeltaRCM”) consists of stochastic parcel-based cellular routing schemes for water and sediment and a set of phenomeno-5 logical rules for sediment deposition and erosion. The outputs of the model include ﬂow ﬁeld, water surface topography and bed topography that evolves in time. Results show that DeltaRCM is able to: (1) resolve a wide range of channel dynamics, including elongation, bifurcation, avulsion and migration; (2) in response to the changes in input parameters, produce di ﬀ erent types of deltas such as alluvial fan deltas at experimental 10 scale. We also identify three key areas of particular model sensitivity, even at the RCM level: (1) avulsion dynamics is sensitive to dynamic free-surface topography; (2) channel network structure is sensitive to instability at channel mouths which creates bars; and (3) out-of-channel sedimentation is sensitive to water surface slope along channel margins. We also demonstrate a simple stratigraphy tracking component which can 15 display the structure of the deposit in terms of distribution of coarse and ﬁne materials along with the age of the deposit. DeltaRCM is a useful tool for understanding the dynamics of river deltas within a relatively simple cellular representation of water and sediment transport.


Introduction
Home to hundreds of millions of people, major coastal cities and infrastructure, immensely productive wetlands, and some of the most compelling and diverse landscapes on Earth -yet low-lying and vulnerable to storms and rising sea levelsdeltas are emerging as among the most critical environments in a changing world (Syvitski et al., 2009).They are also immensely complex.The science of deltas comprises, in roughly equal parts, geomorphology, ecology, hydrology, organic and Introduction

Conclusions References
Tables Figures

Back Close
Full microbial geochemistry, and human dynamics.The physical dynamics alone would present a formidable challenge if they were restricted to just turbulent flow interacting with sand; but most natural deltas involve major additional complications such as fine-grained cohesive sediment (mud) and strong, two-way interactions with biota.A fundamental debate is developing across the sciences as to the best way to model and understand such complexity (e.g., Murray, 2003;Overeem et al., 2005;Paola and Leeder, 2011;Paola et al., 2011;Hajek and Wolinsky, 2012).Should we try to capture everything, creating models that simulate the processes in as much detail as current knowledge and computing power allow; or should we simplify, even at the risk of losing connections with reality?Modeling of deltas in recent years has produced excellent examples of both approaches, which we review below.Our aim here is to present a model that resides in the middle ground between detailed simulation and abstract simplification.We use a method based on weighted random walks, where the random walks are constrained by rules based on a hybrid of simplified governing equations for fluid motion and phenomenological representation of sediment transport processes.With suitable rules, DeltaRCM is able to produce delta morphologies that compare well with those produced by more complex models such as Delft3D and with the morphology of deltas in the field.We believe that the availability of abundant computing power strengthens rather than weakens the case for so-called reduced-complexity models such as the one we propose here.Understanding -as opposed to simulating -complex natural phenomena requires a spectrum of approaches and a clear understanding of the advantages and disadvantages of each.
The paper begins with a review (Sect.2) of current approaches to modeling deltas, emphasizing previous reduced-complexity models.The detailed implementation of our model is presented in Sect.3, and results from it in Sects.4 and 5.In Sect.6 we discuss the meaning of the model results to date.Introduction

Conclusions References
Tables Figures

Back Close
Full As with any morphodynamic model, the most direct delta formation model would solve the governing equations for water flow and sediment particles based on first principles, i.e., the conservation of mass and momentum or energy, in detail, given all the necessary initial and boundary conditions.However, this is still not practical, not only because of limits of computational power, but also because of the potential error accumulation in such complex "full physics" models (Hajek and Wolinsky, 2012).Existing models for delta formation cover a wide range of scales and complexity (Fagherazzi and Overeem, 2007;Paola et al., 2011).
On the simple side, models based on spatially averaged delta surface topography can predict average delta dynamics, such as laterally averaged surface profile, position of the shoreline, and position of the alluvial-bedrock transition (Parker et al., 2008;Kim et al., 2009;Lorenzo-Trueba et al., 2013).These models do not attempt to provide detailed structure, such as topography and channel networks.On the more complex side, to date the most inclusive physics-based delta formation model is Delft3D, which solves a depth-integrated version of the Reynolds-averaged Navier-Stokes equations (shallow water equations) with a turbulence closure term for horizontal Reynolds stresses, and coupled with empirical sediment transport formulas based on bed shear stress (Lesser et al., 2004;Edmonds and Slingerland, 2007).Delft3D can resolve deltaic processes from smaller, engineering scales such as river mouth bar formation and bifurcation (Edmonds and Slingerland, 2007) to larger, geological scales such as the whole delta morphodynamics controlled by sediment cohesion (Edmonds and Slingerland, 2009), wave, tides and antecedent stratigraphy (Geleynse et al., 2010).Delft3D is widely considered the best high-resolution delta model available to the research community, and its utility is greatly enhanced by the release of an open-source version in 2012.In the middle ground of the model complexity spectrum are the so-called reduced-complexity models (RCMs).These models feature descriptive constructions and intuitive simplifications over the hierarchy of natural processes, in contrast to highly Introduction

Conclusions References
Tables Figures

Back Close
Full detailed but computationally complex models such as Delft3D, while still evolving the topography and channel network without simplifying to the degree of spatially averaged models.The most common form of models in this category is a rule-based cellular routing scheme, such as the braided river model by Murray andPaola (1994, 1997) and some of the early erosional-landscape models (e.g.Willgoose et al., 1991).In terms of channel-resolving delta formation models, an excellent example is Seybold et al. (2007Seybold et al. ( , 2009Seybold et al. ( , 2010)).In their model, the water flux field is calculated on a lattice grid via a set of simplified hydrodynamic equations which are equivalent to a diffusive-wave form of the shallow water equations with constant diffusivity.A few other examples of delta-related channel-resolving RCMs include an avulsive delta building model by Sun et al. (2002) and a channel-floodplain co-evolution delta building model, "AquaTellUS" by Overeem et al. (2005).RCMs are less computationally intensive than CFD-based high-fidelity models yet still produce morphodynamic features at system scales, such as stream braiding and floodplain aggradation.While computational efficiency is often considered the reason for developing RCMs, their most important advantage is the flexible rule-based framework which allows direct translation of phenomenological observations into the model (as opposed to hoping that they will emerge given a sufficiently detailed description of the underlying mechanics).The challenges of building a RCM for delta formation are the following: (i) the low topographic slope of the majority of river deltas (10 does not provide a strong guide for topographic flow routing, which is a key component in many RCMs for geomorphodynamic systems, (ii) the low slope together with relatively deep, slow channel flow creates a low-Froude-number environment such that the flow senses information from both upstream and downstream directions, making it difficult to design localized rules which are essential for RCMs, (iii) the self-organized distributary channel network includes loops that further complicate flow routing, and (iv) many river deltas have suspended load and wash load as a primary sediment input components, which make sediment routing more complex than in a bedload-only system.In addition, the low-Froude-number flow condition implies, as the Froude number Introduction

Conclusions References
Tables Figures

Back Close
Full tends to zero, a "rigid-lid" condition in which the shape of the free surface is nearly flat.This condition potentially offers computational advantages as the flow depth can be estimated from a fixed surface elevation (usually sea-level or a simple function using backwater equations) and bed elevation, but is almost decoupled from the bed topography.
In this work, we present an RCM delta model using the "weighted random walk" method.The basic goal is to develop a model that includes just enough of the dynamics to tackle the main problems listed above.A detailed model description is given in the next section, followed by results and comparisons with experimental and field deltas along with the results of more detailed delta models, and then a discussion of the strengths and weaknesses of our model approach.

Model construction
DeltaRCM has two components: a cellular flow routing scheme as the hydrodynamic component, and a set of sediment transport rules as the morphodynamic component.The model uses a lattice of square cells for its domain, where water and sediment flux are routed in a cell-by-cell fashion.The model evolves in time by updating the depthaveraged flow field, water surface elevation, sediment flux, and bed elevation at each time step.

Model setup
The physical setting of our delta formation model is simplified to a rectangular basin of constant water depth (h B ) with a short inlet channel on one side (Fig. 1).At the inlet we assume a constant water discharge Q w0 (m 3 s −1 ) and sediment discharge Q s0 (m 3 s −1 ).
The boundary with the inlet channel is a wall boundary such that no water or sediment crosses.The other three boundaries are ocean boundaries with the boundary condition of a fixed sea level, H SL .Introduction

Conclusions References
Tables Figures

Back Close
Full For water and sediment routing, we first define a set of global parameters that remain constant for each model run: (1) a reference water depth h 0 , i.e. a representative flow depth for the system, and (2) a reference slope S 0 , which is a representative overall water surface slope of the system.For example, for a lowland river delta, a typical value of h 0 is from a few meters to tens of meters, with S 0 on the order of 10 −4 to 10 −5 , while for an experimental fan delta, a typical value of h 0 is tens of millimeters and S 0 on the order of 10 −2 .The values are not precise but rather represent scale values, and may require trial and error to validate for each specific system.The depth of the inlet channel is set at h 0 and the inlet flow velocity is calculated as U 0 = Q w0 h 0 W , which will be referred to as a reference velocity of the system.W is the inlet channel width, specified for each model run.
The domain is shown in Fig. 2, with cell size δ c , a value that depends on the target scale of the model run, e.g., in the results section we use 50 m for a field scale delta and 2 cm for a laboratory scale fan delta.The total number of cells along the dip direction (from the inlet, into the basin) is N X and the number of cells along the strike direction (perpendicular to the inlet, across the basin) is N Y .Typically N X and N Y are both on the order of a hundred, with N Y being roughly twice as large as N X to allow a semi-circular delta growth.The inlet has a width of N 0 cells.Typically N 0 is around 5. The primary quantities associated with each cell include: (i) water unit discharge vector q w = (q x , q y ), (ii) water surface elevation H, (iii) bed elevation η.These primary quantities are updated at each time step.Other useful quantities such as velocity vector u = u x , u y and water depth h can be easily calculated from the primary quantities by: h = H − η and u = q w h .Two types of parcels that carry a water or sediment attribute are routed through the domain.A time step is defined by the addition of n w water parcels and n s sediment parcels.This is done through a sequence of water parcels carrying an equal fraction of the total input water discharge during a time step followed by sediment parcels carrying an equal fraction of the total input sediment discharge during a time step.Introduction

Conclusions References
Tables Figures

Back Close
Full Within each model run, the size of the time step ∆t is constant.As is often the case in numerical modeling, the choice of ∆t is a balance between computation efficiency and model stability.In each time step, the total amount of sediment added to the domain is measured by ∆V s = Q s0 ∆t.A smaller ∆V s means less change to the topography, and allows the cellular routing scheme to perform better with a more consistent terrain, but obviously it will take more steps to build the delta to a certain size.Here we introduce a reference volume, which is the volume of a channel inlet cell from the bed to water surface.If we assume that channels on the delta self-organize in scale with the reference depth h 0 , this reference volume gives a good measurement of the characteristic topographic change on the growing delta.Currently we set the time step size so that the sediment volume added in each time step satisfies: Therefore time step size is given by: (3)

Model operation
The operations can best be understood by describing the processes in a single time step.There are four distinct phases: (1) the addition and routing of the water; (2) updating of the water surface elevation; (3) routing the sediment parcels and updating the bed elevation through deposition and erosion; and (4) updating of the routing direction, a vector field that determines the direction of flow through each cell in the domain.
Each of these phases is described in turn.Introduction

Conclusions References
Tables Figures

Back Close
Full To prepare, we divide the upstream water discharge (Q w0 ) and the total sediment input volume during a time step (∆V s ) into parcels.Typically, we use n w = 2000 water parcels and each water parcel carries equal amount of discharge: And we use n s = 2000 sediment parcels and each sediment parcel carries equal amount of sediment volume: (5)

Phase 1: water routing
At the start of a time step we assume that we have a delta with known shape and topography, i.e. at each cell we have a value of the water surface elevation H, bed elevation η, and water depth (difference between the water surface elevation and the bed elevation) h.We also have, at each cell, a unit vector F , referred to as the routing direction, which indicates the average downstream direction of flow through that cell.If the current time step is the first step in the model run, the routing directions are all in line with the inlet channel.
For the purpose of routing water, we define a binary cell state: 0 -dry, 1 -wet.This is done by doing a sweep through the domain and marking cells with a water depth larger than a small threshold value h dry as wet cells.This threshold value is typically a fraction (10 %) of the characteristic depth scale of the environment of interests or 0.1 m, whichever is smaller.For example, for a natural delta, h dry is typically 0.1 m, while for an experimental delta in laboratory, h dry is typically a few millimeters which is 10 % of characteristic flow depth.
The process in the first part of the time step requires us to route, in turn, each of the water parcels through the domain.When the parcel is at a given cell, this requires Introduction

Conclusions References
Tables Figures

Back Close
Full making a choice of which of the 8 neighbor cells it will move to.We achieve this by using a so called "weighted random walk" where the movement is dictated by a predefined probability distribution between the 8 neighbor cells.The specification of the probability distribution is as follows: At a given cell, first we calculate the routing weights for the 8 neighbor cells.With the local routing direction F specified, the routing weights are determined by two factors: (i) the angle between the relative direction of the neighbor cell i and the routing direction, which we will estimate using a dot product method that we describe below, and (ii) the resistance to the flow from each neighbor cell i .In this model we calculate the routing weight for neighbor cell i as: where resistance R i is estimated as an inverse function of local water depth h i , For the current version of flow routing, the exponent θ is set to 1, hence leading to the following relationship of the routing weight: The cellular direction vector, d i , is a unit vector pointing to neighbor i from the given cell.Finally, ∆ i is the cellular distance: 1 for cells in main compass directions and √ 2 for corner cells (Fig. 3).
The weights above are calculated only for wet neighbor cells of the given channel cell.All dry neighbor cells take a weight value of 0. At each channel cell we can then calculate routing probabilities p i : Introduction

Conclusions References
Tables Figures

Back Close
Full To obtain a discharge vector at each cell based on the motion of visiting water parcels, our starting point is to construct, for each visiting parcel, an average vector of the input and output vectors (Fig. 4).So the result is, for each channel cell, a set (size N visit ) of vectors, each expressing the average path of a visiting parcel through that cell.A summation of this set of vectors provides, after appropriate normalization, a representative direction for water parcels through the cell.In this way, a vector with this direction and a magnitude Q cell = N visit Q p_water can be regarded as a discharge vector for the cell (Fig. 4).
Then for the purpose of later sediment transport, we need to estimate local flow unit discharge and velocity.To do this we take the cell discharge vector (m 3 s −1 ) and divide it by the cell size δ c to obtain a unit water discharge vector (m 2 s −1 ):

Phase 2: water surface calculation
Water surface elevation is essential in this model not only in that it participates in the calculation of flow depth but even more importantly, the gradient of water surface plays a major role in determining the routing probabilities, w i (Eq.8), of water parcels.
In this reduced-complexity model, our goal is to obtain a sufficiently accurate surface profile without solving the full 2-D hydrodynamic equations.We propose a method that uses a finite difference scheme along the movement path of individual water parcels, analogous to the simplified surface solver developed by Rinaldo et al. (1999).
To start with the simplest formulation, we assume that water surface slope along a channel streamline can be approximated by the reference slope S 0 , and in the ocean condition H = H SL , ideally along any given streamline we can reconstruct the surface profile with a simple finite difference calculation.In the model, however, instead of tracing a flow streamline, we take advantage of the walking path of water parcels, which can be considered as an approximation to the flow streamlines.The difference between the water parcel paths (the "zigzag" version of streamlines) and the real flow streamlines is illustrated in Fig. 5.In the following we explain how to construct water surface profile along a water parcel path with a given reference slope S 0 .First, we need to locate the part of the path that is on the delta surface, as the part in ocean is considered flat.In general, a water-parcel path starts at one of the inlet cells, moves from one cell to an adjacent cell, and ends at one of the downstream ocean boundary cells.We distinguish the cells along the path on the delta surface and the cells in the open ocean by checking two values at each cell such that either a cell is in the ocean if both of the following criteria are met: 1. Local bed elevation η is lower than a threshold value η shore (set to η shore = H SL − 0.9h ref ) 2. Local flow speed |u| is smaller than a threshold value U shore (set to U shore = 0.5U ref ) or a cell is on the delta otherwise.
With a given water-parcel path, the calculation starts from the end of the path and goes backward towards the inlet.For the kth cell in the direction of calculation, The purpose of the term is to take into account the angle between the parcel path and the streamline.
This calculation gives the surface profile along the path of an individual water parcel and is repeated for all water parcel paths.There are two extra situations to be taken care of: 1.If a cell is visited by multiple water parcels, all the values from each visiting path are recorded and an average value is taken among these stored values in the end to obtain a single value for water surface elevation at each cell.
2. If a cell is not visited by any water parcels, its water surface elevation remains the old value (from the previous time step).
This newly calculated surface profile is recorded as H temp .We then apply a diffuser to smooth the calculated surface profile, which is typically spiky due to the 1-D stepwise method of calculation.The diffusion is applied as: We have used a diffusivity ε = 0.1 and we apply the smoothing calculation in Eq. ( 11) for 10 times in each time step.This number is selected by checking samples of the resulted surface profile until no obvious spikes appear.We will discuss more in detail how sensitive the results are along with other features in calculating the free surface.
In the end, the water surface elevation is updated with an under-relaxation scheme for numerical stability: The under-relaxation coefficient is set to 0.1, which allows the surface profile to transit slowly and smoothly from one time step to another, avoiding numerical instability.Introduction

Conclusions References
Tables Figures

Back Close
Full

Phase 3: sediment transport and bed topography update
Now both the flow field, q w , and water surface elevation, H, are updated.These two variables will remain constant until the next time step.To calculate the changes to the topography in a time step, we propose two sets of rules for the transport, deposition and erosion of sediment.The first set describes the routing of the sediment parcels, and the second set describes the rate of deposition and erosion as the exchange of sediment volume between sediment parcels and the bed.The rules are phenomenological and the goal is to build them via our understanding of macroscopic behavior rather than via fine-scale physical interactions among fluid, sediment and the bed.To this end, we distinguish two types of sediment that have different behaviors in the model: coarse sediment, referred to as "sand", is non-cohesive, and transported as bed load; fine sediment, referred to as "mud", is cohesive, and transported as suspended load.
A sediment parcel is either a "sand" parcel or a "mud" parcel.At the beginning of each run, an input parameter f sand gives the portion of sand in the total upstream sediment input.Therefore a total number of f sand n s parcels are designated as sand parcels and a total number of (1 − f sand )n s parcels are designated as sand parcels for each time step.

Routing of the sediment parcels
For routing sediment parcels we use the same weighted random walk method as for the routing of water parcels (Eq.6) with two modifications: 1.The routing direction F in Eq. ( 6) is replaced with the newly calculated water discharge vector q w at the given cell (from Phase 1 above), assuming that sediment parcels move with the water flow; Introduction

Conclusions References
Tables Figures

Back Close
Full 2. Transport resistance for sediment maintains the inverse function of flow depth but has different exponents.The idea is that sediment flux tends to concentrate in the lower portion of the water column and therefore it is more likely to follow topographic low areas.For now we use an exponent θ = 2 for sand parcels (bed load) which is twice the value for water, and θ = 1 for mud parcels (suspended load) which is equal to the value for water.The physical reason for the values chosen is that the distribution of the concentration of coarse material is skewed towards the lower portion of the water column and the distribution of fine material is more evenly distributed throughout the water column.
Thus the routing weights for sediment parcels are: for sand parcels, and ( 13) And routing probabilities are calculated as: The rate of deposition and erosion Sediment parcels are routed sequentially in a weighted random walk fashion according to the probabilities calculated with Eqs. ( 13 (ii) how much volume should be exchanged between the sediment parcel and the bed.
The rules for sand and mud parcels are different.
For convenience of description, we refer to the initial volume of each sediment parcel V p_sed as the "reference sediment parcel volume", and the remaining volume during the walking process of a sediment parcel as the "residual sediment parcel volume", V p_res .
The amount of deposition at each cell by an individual parcel is referred to as V p_dep .The amount of erosion at each cell by an individual parcel is referred to as V p_ero .Detailed rules are: For the deposition from a sand parcel we do the following: -At each cell in the domain, we calculate a "transport capacity" for sand flux, q s_cap , as the maximum flux per unit width, which is a non-linear function of local flow velocity U loc .The scaling between sediment flux and flow velocity takes the form of the Meyer-Peter and Müller (1948) formula where q s0 is calculated by dividing the upstream sand flux input by the inlet channel width: -Similar to the calculation of water discharge, as the sand parcels are routed sequentially, we track the accumulated total sand flux, q s_loc , which increases with each visiting bed load parcel: -Deposition occurs if a sand parcel visits a cell that has an accumulated local sand flux exceeding the transport capacity: For deposition from a mud parcel we do the following: -Deposition occurs if a mud parcel visits a cell that has a local flow velocity U loc smaller than a threshold velocity, U dep .The amount of deposition is proportional to the residual sediment volume of the mud parcel as well as the relative difference between the squares of U loc and U dep , a simplified representation of standard empirical laws for fine-sediment deposition (van Rijn, 1984) : The idea is that the finer grain size, the slower flow it requires to settle out.
For the erosion of both types of sediment parcel, we do the following: -Erosion occurs if local flow velocity magnitude is larger than a threshold value, U ero , that differs for sand and mud parcels (García and Parker, 1991): -For a sand parcel, U ero = 1.05U ref .
For volume exchange between sediment parcel and the bed: -At each step, the volume of the sediment parcel is updated as: -The elevation of the local bed is updated as: After all sediment parcels finish their random walk, to take into account the influence of topographical slope on sediment flux in an approximation of the Bagnold-Ikeda expressions (García, 2008), we apply a topographic diffuser that assumes where α is a scaling coefficient, by default set to 0.1; |∇η| is bed slope.The total change to the bed elevation by the topographic diffuser is obtained by summing up the inbound and outbound diffusive fluxes at each cell over the time period ∆t.

Phase 4: update routing direction
Before moving to the next time step, we need to update the routing direction, a unit vector at each cell indicating the downstream direction for routing water parcels.In this last phase of the time step, at each cell we calculate the updated value of the unit water discharge vector q w , water surface elevation H, bed elevation η, water depth h, etc.
To achieve this, we combine two physical processes dictating the flow direction: (i) at an instant in time flow has a tendency to continue in the same direction as the direction at the previous instant due to inertia, and (ii) in the absence of any other drivers the flow goes down slope which in our case is indicated by the water-surface slope rather than bed slope.
First, we calculate a unit vector from the downstream direction based on the previous time step: Then, we calculate a unit vector from the water surface gradient (from the previous time step): Then a linear combination of the two vectors is taken with a partitioning coefficient γ: The value of γ is set to a small number, typically 0.05 in the runs reported here.

Model results
In this section we present various morphological features produced by DeltaRCM with different domain setup and input parameters.All simulations assume no effects from wave or tidal energy, i.e. the delta is a classic river-dominated delta (Galloway, 1975).We investigate: (1) the effects of input sediment composition and (2) the model's ability to simulate deltas at field and laboratory scales.The former has been studied via field observation (Orton and Reading, 1993) and numerical simulation (Edmonds and Slingerland, 2009).The latter is based on the availability of data from experimental deltas; also, we believe that if a model can handle both field and experimental scales, it could potentially inform the interpretations and connections of both.Furthermore, we demonstrate using DeltaRCM as a tool for hypothesis testing, through study of the effects of receiving basin depth.
As discussed above, two types of sediment are routed through the system: coarse (sand) and fine (mud).The ratio of the numbers of these two types of parcels at the inlet gives the ratio of sand and mud coming into the system.To set the physical scale of the simulation, domain and grid size are adjusted by changing cell size and physical input parameters, such as total input water and sediment discharge, and also global parameters such as the reference energy slope.
The input parameters include: 1.The portion of sand in the upstream sediment input, f sand ; Introduction

Conclusions References
Tables Figures

Back Close
Full 2. Global parameters: the reference flow depth h 0 , and basin depth h B , and the reference slope, S 0 .
Strictly speaking, the choice of the reference slope S 0 is dependent on the sand-mud ratio as well as the scale of the physical setting.In our model runs for field scale we use 3 × 10 −4 for purely sandy deltas, 1 × 10 −4 on purely muddy deltas and a linear combination of the two for mixed deltas; for laboratory scale, we use values on the order of 10 −2 for S 0 .The magnitude of the reference slope is scaled with the ratio of bedload and water fluxes that come from the inlet channel, that S 0 ∼ Q s0_bed /Q w0 .

Effects of input coarse/fine sediment ratio
In this group, the domain is a lattice grid of 120 by 60 square cells.Cell size is taken to be 50 m.The channel inlet is 5 cells wide (250 m), with a reference flow depth h 0 = 5 m.The total water discharge is 1250 m 3 s −1 .The total sediment discharge is 0.1 % by volume, which is 1.25 m 3 s −1 .We use a time step calculated from Eq. ( 3) of 25 000 s (∼ 7 h).Both water and sediment discharge stay constant and we assume they represent channel-forming conditions.We show three model runs in Fig. 6 with the portion of sand in the upstream sediment discharge set to 90, 50, and 10 %.The resultant deltas differ systematically based on the input mud fraction in the following characteristics, which are consistent with those found in the investigation on the effects of sediment cohesion by Edmonds and Slingerland (2009): -On a sandy delta the channels are relatively shallow and mobile, without welldefined levees.Flow is less confined.There are large areas of sheet flow.The shoreline is smooth and the delta grows roughly in a semi-circular shape.Introduction

Conclusions References
Tables Figures

Back Close
Full -On a muddy delta, channels are deeper and stable, with well-defined levees.
Channels tend to elongate.The shoreline is rugose, and deltas build in different directions by switching lobes.

Experimental fan deltas
Laboratory experiments, numerical modeling and field observation are three important approaches of understanding the formation of deltas.We would like to test our model across as wide a scale range as possible so we include experimental deltas at laboratory scales.To do this, we change the domain to a lattice grid of 80 by 150 cells with a cell size of 0.02 m.The inlet channel is still 5 cells wide but has flow depth of 0.02 m and a water discharge of 0.6 L s −1 .Basin water depth is 0.02 m.The reference slope is set at 0.01.The time step is estimated at 33.3 s.Sediment input is considered to be coarse-grained only.These conditions are representative of laboratory experiments such as those reported by Reitz and Jerolmack (2012).
In Fig. 7 we show a time series of the resultant deltas.Note the alluvial fan delta characteristics, in which a few active channels quickly switch to build a semi-circular shape with a relatively smooth shoreline.To evaluate the details of the channel switching process, we compare a representative switching event from our model results to the avulsion cycle described by Reitz and Jerolmack from the experimental fan delta (Reitz and Jerolmack, 2012).DeltaRCM captures all the steps in the avulsion cycle (Figs.8 and 9).

Effects of basin depth
It has been suggested that the accommodation -the space that a delta can grow into -plays an important role in the architecture and behavior of a growing delta (e.g., Paola et al., 2000;Heller et al., 2001).However, for the case of river deltas with very low-Froude-number flow, it is still unclear how the depth of the basin affects the overall morphology of the delta.Storms et al. (2007)  delta formation from a river effluent discharging constant flow and sediment loads into shallow and deep receiving basins under homopycnal conditions, and shows that the shallow-basin delta is dominated by mouth bar bifurcations and a shoaling channel network, and exhibits significant stratigraphic complexity and sub-aerial development, while the deep basin delta is dominated by unstable bifurcations, levee breaches and avulsions (Storms et al., 2007).The authors suggest that the shallow basin case resembles Wax Lake Delta.In our model runs 6 and 7, we test scenarios with the same inlet channel conditions and discharge, but different basin depths.In run 6, the receiving basin depth is half of the reference depth (defined by the inlet channel which is supposed to be at equilibrium state in terms of sediment transport), while in run 7, the receiving basin depth is twice of the reference depth.In Fig. 9 we show that our results give similar behaviors to the ones modeled by Storms et al. (2007).For the shallow basin the morphological development is very close to the description of Storms et al., while the deep basin delta has similar outcomes but the middle ground bar and avulsion over the levee are not as clear in the RCM results.
The difference between a shallow and deep receiving basin, according to our model results, is that: channels will still try to maintain the same unit power of transporting sediment by maintaining a certain cross-sectional geometry with levees on the side and erosion or deposition on the bottom; in general a distributary channel network shoals up and channels are stable at shallower depth going seaward.With a shallow basin the amount of work is reduced.Also the narrow space promotes the splitting of flow which enhances the growth of a distributary network; a deep basin increases the time scale of establishing a stable channel, therefore introduces stronger competition among channels by allowing larger differences to develop.Introduction

Conclusions References
Tables Figures

Back Close
Full Finally, we note two interesting emergent features from our model that has also been observed in the field at Wax Lake Delta by John Shaw and colleagues (Shaw et al., 2013;Shaw, 2013) (Fig. 11).First, the channels in the shallow basin delta are erosional initially, which carve into the basin bottom.This is consistent with the observations at the Wax Lake Delta (Shaw et al., 2013).Second, the channel network on this delta develops "tributary" sub-networks on islands (highlighted in Fig. 11), which collect flow both from tie channels directly connected to the main channel network and from sheet flow topping the levees into the islands.Whether this sub-network is erosional or depositional, Shaw (2013) points out that at least the channels comprising it are likely not favorable for deposition.In our model results, we notice the following feature that might explain the situation: 1. the sub-network mainly collects fine sediment from the main channel network, which requires much slower flow to settle out; 2. as the tributary sub-network joins into bigger trunk channels, the ability of the flow to carry sediment increases; 3. at the downstream end of the network, where the trunk channel collecting water coming out of the island meets the open water, the sorting of the sediment deposited is very similar to a normal channel that has a coarser bar-like structure at the mouth.

Recording of stratigraphy
A delta writes its own autobiography by preserving deposited sediment underground.
These sedimentary records open a door to understand the past and to use delta deposits to reconstruct their range of natural behavior.Therefore, the ability to record stratigraphy in a delta formation model enables us to directly investigate the connection between surface and sub-surface processes.In this model, we have two methods that Introduction

Conclusions References
Tables Figures

Back Close
Full track the stratigraphy of model-produced deltas: the first method tracks the distribution of coarse and fine sediment by recording the percentage of sand in each deposition event; the second method tracks the age of the deposit by labeling each deposition event with the time that its sediment enters the domain from the inlet channel.Here we present one example from each.(1) We take a sample run of a field scale delta and 30 % sediment input (Run 4).In Fig. 11, we show a stratigraphic slice in the dip direction along the center line of the inlet channel.In Fig. 12, we show the time series of the stratigraphic slice in the strike direction about 20 cells (1 km in this case) away from the inlet channel.In both figures white represents pure sand and dark blue is pure mud, with mixed deposits represented by linear combinations of the two end members.
Generally speaking, coarse sediment (sand) can be found in channel belts and mouth bars, while fine sediment (mud) can be found in distal regions such as the bottom set of the delta, or on the floodplain or in abandoned channels.
(2) In Fig. 13 we show a sample model run for laboratory conditions (Run 8).Note the evolution of the area pointed to by the yellow arrow.The series of images shows the deposition sequence from an individual avulsion event.

Discussion
One of the themes running through this paper is that even in the framework of a reduced-complexity delta model there are a number of important details that must be modeled fairly accurately to achieve even qualitatively correct model results.These include a reasonably accurate representation of the water surface and the inclusion of suspended sediment deposition and entrainment.This is quite striking compared to the success of even fairly radical reduced-complexity approaches in modeling other morphodynamic environments such as erosional landscapes (e.g., Willgoose, 1991), braided rivers (e.g., Murray and Paola, 1994) and eolian bedforms (e.g., Werner, 1995).
So why is it that deltas seem to require more attention to detail?Can we learn anything Figures

Back Close
Full from this experience that might help us better understand what systems are most and least amenable to reduced-complexity approaches?Since deltas and drainage basins share dendritic channel patterns -one is a distributary network while the other is a tributary network -we first look at the differences between these two systems.In modeling the evolution of drainage tributary networks, even highly simplified relations for water flux and sediment transport give quite reasonable drainage networks and elevation changes in long-term evolution of catchments (e.g., Willgoose et al., 1991).The equation describing the evolution of land elevation in Willgoose et al. (1991) includes two transport processes: fluvial transport process and diffusive transport process.The former is dependent on the discharge and the slope in the steepest downhill direction, and the latter, diffusive transport process, is dependent on slope alone.These simple formulations cannot be easily applied to modeling deltas because deltas are low-gradient environments where the transport direction and capacity are to some extent decoupled from bed elevation and slope.To be more specific, (1) bed slope in low-gradient environments is often uncorrelated with flow direction and strength, for example bed slope points opposite to the direction of flow where channels shoal up towards the shoreline; (2) the water surface, which dominates local flow routing, is largely independent of bed topography; (3) the typical low Froude number flow in low-gradient deltaic environments creates strong backwater effects that imply strong non-locality in flow and sediment flux control (Lamb et al., 2012;Nittrouer et al., 2011) -meaning that downstream conditions control upstream flow dynamics (Hoyal and Sheets, 2009); and (4) river mouth and shore processes such as waves and tides also control the overall morphology of deltas, providing additional process complexity.
According to Werner (1995), for a nonlinear and dissipative system, considerable simplification can be applied if the system exhibits the following two properties: (1) it has a finite number of steady states as "attractors", and (2) it has macroscopic emergent behaviors that are self-organized and consistent with but decoupled from microscopic physics.If we compare drainage networks with deltas, the former exhibits a strong generic pattern and scale-invariant properties such as network patterns that can be Introduction

Conclusions References
Tables Figures

Back Close
Full described by Horton's laws (Horton, 1945).In contrast the networks on deltas have many varieties, responding to a wide range of processes, that no universal geometry can describe them as a whole group.Regarding model complexity, the lack of universality in the system pattern indicates the requirement for a more detailed, system-specific approach in modeling them.
So, is the low gradient the main cause of the modeling difficulty, making deltas more "unforgiving" than erosional landscapes in terms of the accuracy of hydrodynamic calculation?For cellular models that use explicit flow routing schemes, the complexity level rises as factors other than topographic slope alone determine water and sediment routing.It also increases with non-locality in the broad sense of the sensitivity of dynamics at one point to conditions far away in the system.Other contributing factors such as water surface gradient and flow inertia weigh in as the overall topographic gradient decreases.For example, dune fields may have very low to zero average topographic slope, but they have locally high steepness meaning that, as in erosional landscapes, the sediment dynamics are dominated by bed topography.In deltas, however, the controlling factor is the relatively subtle water surface topography, therefore simple descriptions relating sediment deposition and erosion to e.g., local elevation and slope give realistic dune field dynamics but do not work in deltas.
Can we be more systematic about evaluating the amount of detail needed to model a geomorphodynamic system?This is an important fundamental question in morphodynamic modeling, and we do not pretend to resolve it here.But our experience with this model suggests the following guidelines as a starting point: -For gravity-driven systems, the overall gradient of the landform is one important index, in the sense that in high-gradient systems the gradient alone is enough to route the flow.Introduction

Conclusions References
Tables Figures

Back Close
Full steepest-path method (Passalacqua et al., 2010) are sufficient to determine the flow path, without the need for simulation of the flow details.
-Froude number (Fr ): as Fr tends to unity, the backwater length tends to zero (Cui and Parker, 1997), so the simplification of local normal-flow assumption provides a satisfactory means of accounting for momentum balance in the flow.
-For systematic behaviors comparable to or beyond the backwater length scale, in-channel scale hydrodynamic details can be resolved at much lower complexity, such as avulsion models that use single-cell-wide threads to represent channel belts (Jerolmack and Paola, 2007).
-Whether the system to be modeled exhibits a strong generic pattern or scaleinvariant (e.g., fractal) properties.The lack of a universal "pattern" in a dynamic system is an indicator of sensitivity to local detail.

Conclusion
In this paper we have introduced a new reduced-complexity model (RCM) for river delta formation.Key techniques include: (1) water and sediment flux is represented as parcels and routed through the domain in a Lagrangian point of view; (2) the movements of parcels are based on a probability field calculated from rules abstracting the governing physics; (3) deposition and erosion are achieved by exchanging the volume of passing sediment parcels and bed sediment columns, and the condition for this exchange depends on a set of rules that distinguish bed load and suspended load; (4) bed sediment columns record the composition of coarse and fine material in layers; (5) a topographic diffusion process that takes into account cross-slope sediment transport and bank erosion.By varying input conditions such as the ratio of coarse and fine sediment, reference slope, and dimensions of the domain, the simulated deltas give a range of different behaviors that compare well to higher fidelity model results and observations of field and experiment deltas.We find that the relatively simple cellular representation of water and sediment transport is able to replicate delta morphology at the scale of channel dynamics, including the emergent channel network with channel extension, bifurcation and avulsion.Here we summarize the basic components sufficient for a RCM to produce major static and dynamic features of river deltas: a depth-averaged flow field that guides sediment transport; a non-trivial water surface profile that accounts for backwater effects at least in the main channels; representation of both bedload and suspended load; topographic steering of sediment transport.
Even at the RCM level of modeling, the following items still require a physically consistent treatment: the instability at channel mouths that creates bars and subsequent bifurcation; the variation in water surface profile associated with lobe extension that causes channel avulsion; water surface slope along channel sides which creates flooding onto the floodplain.
We see the potential of this type of modeling in a similar way as of physical laboratoryscale experiments, which effectiveness does not necessarily come from the classic scaling or a form of detailed solution of the governing equations (Paola et al., 2009).
The strength of RCMs is to serve as (1) exploratory models that allow for direct representation of phenomenological observation; (2) a tool to identify larger scale processes that are not sensitive to the details of smaller processes; and (3) a framework for hybrid modeling in which higher-resolution model results can be integrated where precise description of smaller scale processes is needed even for larger scale dynamics.Introduction

Conclusions References
Tables Figures

Back Close
Full  Full  (Shaw et al., 2013) observed both in the field (Shaw, 2013) and in our numerical model results.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the water surface slope is always zero.With the downstream water surface boundary ESURFD Discussion Paper | Discussion Paper | Discussion Paper | where ∆ is the cellular distance between the kth and (k − 1)th cell, δ c is cell size, and d| k is the parcel step vector from cell k to cell k − 1. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | )-(15).The change to the bed topography is obtained by the exchange of sediment volume between the moving parcel and the local bed at each cell along the path -during deposition a sediment parcel loses part of its volume and this volume is added to the bed, and vice versa for erosion.We use simple phenomenological rules to decide (i) where deposition or erosion happens and Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

)-
The local flow velocity and flow depth are updated in accordance with each event of deposition or erosion: h = H − η and u = q w h .The reason for updating local flow depth and velocity immediately after each event of deposition and erosion is to avoid excess change to the bed.Similarly we add an extra control on the rate of change to the bed by limiting the amount of deposition and erosion by a sediment parcel so that the change to local depth is less than 25 %, so that the change to local flow velocity is less than 33 %.For example, if local flow depth is 4 m, then the maximum deposition or erosion by a single sediment parcel is limited to 1 m change to the bed.
Discussion Paper | Discussion Paper | Discussion Paper | the diffusive flux is proportional to local sand (bed-load) flux and topographical slope: Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | uses Delft3D to model initial Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Figure 10 . 9 Figure 10 .
Figure 10.Flow features on the island of a delta formed in a shallow basin.(a) Model result 3 from Run #6, where basin depth (2.5 meters) is only half of the inlet channel depth (5 meters); 4 (b) Wax Lake Delta, where basin depth (<5 meters) is much less than the inlet channel depth 5 (>20 meters); (c) schematic drawing showing the "tributary" flow feature on the island (Shaw 6 et al., in preparation) observed both in the field (John Shaw, PhD dissertation) and in our 7 numerical model results.8 9

Table 1 .
List of delta model runs and parameter values.