Morphodynamic Modeling of River-Dominated Deltas: A Review and Future Perspectives

DeltaRCM ( Liang et al., 2016b ) and Delft3D ( Van de Lageweg and Slangen, 2017 ) have been used to show that relative sea level rise (RSLR) can also influence channel dynamics. When a “filling index,” defined as the ratio of accommodation space creation (area times RSLR) to sediment discharge, falls below 0.6, distributary channels tend to reduce in lateral mobility (Liang et al., 2016a,b). Deltas must eventually become inundated for a filling index of less than one. Hence, it is these “drowning” deltas that can be expected to have channels with reduced mobility. Van de Lageweg and Slangen (2017) model the co-evolution of sea level rise and morphodynamics in Delft3D. During sea level rise they allow the network to dynamically adjust to compensate for the accommodation space created. On average a dynamic network—that is, one that is allowed to adjust—lost 20% less land than a static network ( Van de Lageweg and Slangen, 2017 ). They also found that the type of delta morphodynamics, i.e

1 Introduction and scope of the review 1.1 Why simulate river-dominated deltas and their processes?
River deltas are a compelling target for numerical simulation because they contain seemingly organized patterns and shapes at a variety of scales. For instance, most river-dominated deltas, regardless of size, have triangular to semi-circular planform shapes, channel networks, and channel bifurcations. The common presence of these features among most deltas in the world (Caldwell et al., 2019;Nienhuis et al., 2020) suggests there are consistent underlying physical processes controlling delta form and behavior. In this review, we discuss how numerical modeling, and more specifically a type of modeling focused on the morphodynamic feedback, has helped explore some of these key physical processes over the last 15 years.
Morphodynamic modeling of river deltas requires a breadth of approaches because of the range of time and space scales operating in deltaic environments. During the last 15 years, a suite of deltaic models has emerged, ranging from 1D morphodynamic models to idealized models abstracting flow-routing processes, to detailed channel-resolving models that create networks. This suite of models (presented in Section 2) allows for diverse use in many scenarios and has led to new ideas for physical phenomenon that control, for example, the position of the avulsion node and controls on delta progradation (Sections 3.1 and 3.2), the origin of deltaic channel networks (Sections 3.2 and 3.4), and the asymmetry and stability of bifurcations (Section 3.3).
Simulating river deltas not only reveals essential information about their origin, but it also supports multi-disciplinary science. Morphodynamic models have supported the burgeoning multidisciplinary field of studying deltas as social-ecological systems (Brondizio et al., 2016b). Social-ecological systems consist of biophysical (climate, geology, hydrology, ecology) and social factors that have coupled interactions. River deltas are some of the most heavily populated landforms on Earth Szabo et al., 2016), and understanding how these social-ecological systems evolve depends on sedimentation processes required for the delta to grow. These efforts are important because the stress of human occupation coupled with climate change has pushed many deltas to the brink of sustainability (Syvitski et al., 2009;Brondizio et al., 2016a). Morphodynamic modeling has also led to new insights about the stratigraphic architecture of deltaic deposits, and understanding their reservoirs of oil and gas.

What will this review cover?
We limit our focus to morphodynamic models that simulate river-dominated deltas. We consider models morphodynamic when they have coupling among changes in bed elevation, fluid flow, and sediment transport. We consider deltas river-dominated when their dynamics and morphology are strongly linked to river processes at the coastline. River-dominated deltas are common (Nienhuis et al., 2020;Caldwell et al., 2019), and they often are among the larger and more densely populated deltas Szabo et al., 2016). We review studies that seek to understand the subaerial expression of river-dominated deltas, which is usually the result of basinward progradation as fluvial sediment accumulates nearshore and of the creation of distributary networks without progradation. River-dominated deltas can also be influenced by waves and tides. We discuss some of those influences, but do not consider cases where wave-driven processes, primarily involving alongshore sediment transport (e.g., Ashton and Giosan, 2011), exert a strong influence on delta morphodynamics. We do not cover steeper-sloped, coarse-grained fans and fan-deltas, and we do not cover estuarine deltas, or deltas generated from tidally reversing flows, such as ebb and flood-tidal deltas (e.g. Our review focuses on deltaic processes that are relevant to timescales ranging from Ο(1) years to Ο(1000) years. This includes processes from mouth bar building to deltaic lobe switching but does not include longer term processes such as transgression and regression of the deltaic shoreline during cycles of sea-level changes; the interested reader is directed to Overeem et al. (2005), for a review on models typically used for longer term, stratigraphic purposes. We also do not look at subsidence caused by either shallow compaction or deep-seated faulting. However, human-driven compaction can occur over short time-scales, e.g., decades (e.g., Higgins et al., 2013).
Even given our restrictions in scope, there is an impressive range of morphodynamic models used to understand river-dominated deltas. In our view, this range is a natural outcome of how the community uses models. As a guiding principle, most researchers seek to U N C O R R E C T E D P R O O F use or create models that are of enough complexity to address the science question of interest. Complexity refers to the level of detail, the dimensionality of the governing equations (1D, 2D, or 3D), or the number of processes in the model. Because interests in deltaic environments range over huge spatial and temporal scales, the processes and complexity required to answer the corresponding questions is impressively broad. Because of this, the community has developed a suite of models with different levels of assumptions and approximations in how they represent deltaic processes. Those approximations are embodied in models ranging from simplified 1D representations of a deltaic prism, to 2D planform models that have an abstraction of flow laws, to 3D models that solve the shallow water equations. These models have been used to both explain the origin of natural phenomenon and simulate the natural environment as accurately as possible. We focus our review on studies that seek to explain the origin of natural phenomenon (e.g., bifurcation stability or avulsion) in river-dominated deltas. delta, there are two coeval channels at the avulsion node because human engineering has stopped the steeper-sloped Atchafalaya River from capturing the Mississippi River. Some other deltas have multiple active channels at the avulsion node, e.g., the Nile, but many do not.

Organization of this review paper
Our review consists of three main components. Section 2 describes the ways researchers simulate deltaic environments. In Section 2 and its subsections there are detailed descriptions of models, their governing equations, and their limiting assumptions. In Section 3, we focus on what has been learned from these different modeling approaches. Each subsection in Section 3 ends with a summative statement of the key remaining knowledge gaps. In most cases, the results of multiple modeling approaches from Section 2 will contribute to a given subsection in Section 3. Finally, in Section 4 we discuss the future needs for morphodynamic modeling of deltaic environments.
2 Approaches commonly used to simulate river deltas

1D morphodynamic profile models
How the main deltaic channel aggrades and degrades over the topset sets the pace of many deltaic processes, particularly avulsion and delta lobe progradation. Because of this, there is a class of 1D morphodynamics models that focus on the aggradation and degradation of the main channel centerline of river deltas. These profile models simulate the main channel and usually extend from the backwater zone through the channel network and into the receiving basin. On the Mississippi River delta, this domain stretches for nearly 500 km along the river centerline from the avulsion node to the shoreline (Fig. 1). These 1D morphodynamic profile models generally come in two styles that depend on whether flow is represented as steady and uniform or steady and non-uniform. Steady-uniform profile models describe the delta as a triangular prism superimposed upon a relatively planar basement profile and the topset slope adjusts in response to gradients in sediment transport (Swenson et (Fig. 2). The second style of 1D morphodynamic profile models capture non-uniform flow and subsequent bed adjustment (De Vries, 1973;Cunge, 1980;Garde, 2000;Parker, 2004;Garcia, 2008). Both model types tend to represent foreset progradation using a moving boundary formulation (discussed below). In the following sections, we first review how 1D morphodynamic models with steady, uniform flow are implemented and solved, and then we review what processes can be included in 1D profile models when non-uniform flow is represented.

1D profile models with steady, uniform flow and a moving boundary framework
The simpler 1D morphodynamic profile models of deltaic dynamics assume steady, uniform flow, and calculate the fluvial topset slope with conservation of mass equations. Other aspects of the delta are modeled geometrically. These models usually include a low gradient depositional fluvial topset, and a subaqueous offshore region or foreset of equilibrium slope ψ. These interconnected environments move over a relatively planar basement slope of magnitude β (Fig. 2). These environments are in turn delimited by three moving boundaries: the alluvial-bedrock transition located at x = r(t), which separates the bedrock (or basement) from the topset, the shoreline located at x = s(t), which separates the topset from the foreset, and the delta toe located at x = p(t), where the subaqueous sediment wedge intersects with the bedrock (Fig. 2). To determine the trajectories of these boundaries, their location at any point in time needs to be determined internally as a part of mass balance in the solution to the overall problem (Swenson et  (1) where η is the riverbed elevation with respect to the current sea level, x is the distance downstream with x = 0 located at the intersection between the initial sea level and the basement (Fig. 2), and θ is a nonlinear exponent. v is defined as the fluvial diffusivity and calculated as a function of water discharge and grain size characteristics (Paola et al., 1992;Paola, 2000;Swenson et al., 2000;Postma et al., 2008;Lorenzo-Trueba et al., 2009). When θ = 0, the relationship between sediment flux and riverbed slope is linear. The rate of bed elevation change from mass balance, typically written assuming dilute suspended load is (2) where λ p is the bed porosity. When this equation is combined with Eq. (1), it forms the classic fluvial diffusivity equation. Neglecting the backwater effect (see Section 2.1.2 for more detail) and assuming flow is steady and uniform, then conservation of flow mass and momentum are where q w is volumetric water discharge per unit width, h is the average water depth, U is average velocity, C f is a friction coefficient, u ⁎ is the shear velocity given as , S is bed slope, g is acceleration due to gravity, ρ is fluid density, and τ b is the bed shear stress. C f is often assumed to be a constant (Paola et al., 1992 (3a-c), and assuming only a fraction δ of the fluvial topset area is occupied by active channel water flow, we obtain the following expression (4) Next, we consider a general volumetric bed-material load (including bed load and suspended bed material (not washload)) equation per unit channel width as follows: (5) where R = (ρ s − ρ)/ρ, ρ s is the sediment density, D is the particle diameter, τ * = is the Shields number (τ c * is the value at the threshold of motion), and A and m are empirical constants. Eq. (5) can be simplified under two cases. First, for braided gravel-bed rivers, the channel width self-adjusts to keep the bed shear stress slightly above critical, i.e., τ * = (1 + ϵ)τ c * (Parker, 1978;Parker et al., 2007 As a second order differential equation in space, solving Eq. (7) generally requires two boundary conditions. However, when the domain boundaries change over time, two additional moving boundary conditions (i.e., r and s) are needed to determine their location. First, researchers often impose a given sediment flux q s0 at the alluvial-bedrock transition: Second, the fluvial topset elevation is set equal to the alluvial-basement transition to the basement elevation: Eq. (9) is only required when the alluvial-basement transition is part of the solution. This can be avoided by assuming the alluvial basement transition is fixed at the origin (i.e., x = 0).
Third, the elevation of the fluvial topset at the shoreline is set equal to the current sea level Z: The fourth condition relies on the analogy between heat and sediment diffusive transport; for more details see Swenson et al. (2000). Given the geometry presented in Fig. 2, and assuming the shoreface toe only migrates seawards (i.e., dw/dt > 0), we can define the offshore water depth to be filled by sediment as ψ/(ψ − β) · (sβ + Z) (Lorenzo-Trueba et al., 2013). But the foreset slope ψ is generally orders of magnitude larger than the basement slope β, so it is often assumed ψ/(ψ − β)~1. Under these conditions, the rates of migration of the shoreline and the delta toe are equal (i.e., dw/dt = ds/dt), and we can relate the sediment flux that reaches the shoreline with the rate of migration of the shoreline as follows: (11) With Equations (7-11), and a given initial geometry: s(0) = r(0) = h(x, 0) = 0, we can fully describe the dynamics of the fluvial topset under a range of sea-level change scenarios as long as the system maintains the wedge geometry depicted in Fig. 2 (i.e., ds/dt ≥ 0). Recent work by Anderson et al. (2019) has relaxed this geometric constraint to account for scenarios beyond just shoreline seaward migration, such as cycles of shoreline retreat and progradation.

1D profile models with non-uniform flow: Simulating the backwater zone
The backwater zone is a long reach in coastal rivers where the influence of the downstream standing body of water causes non-uniform flow. These flow dynamics impact bed elevation changes and are often incorporated into more complex morphodynamic models (compared to Section 2.1.1). In this section we extend the formulation presented in Section 2.1.1 to incorporate non-uniform flow. These types of models allow for additional processes to be represented, including lateral flow spreading in the offshore river plume, variable flood discharges, and deltaic lobe formation.
Non-uniform flow in deltaic channels can extend for hundreds of kilometers upstream from the river mouth, typically across the entire delta topset (Jerolmack and Swenson, 2007) and is referred to as the backwater zone (e.g., Lamb et al., 2012). In models that incorporate non-uniform flow, Eq. (3) is replaced by (12) and can be solved with an equation for conservation of sediment on the river bed (e.g., Eq. 2). For rectangular channels Fr 2 = U 2 /gh and w is flow width (Chow, 1959). The assumption of uniform flow is valid only when the left-hand side of Eq. (12) is small compared to terms on the right-hand side, which is unlikely to be the case in low-gradient deltaic rivers that have small Fr and large due to lateral flow spreading offshore or through deltaic bifurcations. Some models have implicitly incorporated lateral flow spreading in the offshore river plume by holding the water surface elevation constant at the river mouth, under the assumption that rapid spreading would minimize differences in water surface elevation Karadogan et al., 2009;Moran et al., 2017). But not all river plumes spread rapidly, potentially invalidating this assumption (Falcini and Jerolmack, 2010;Mariotti et al., 2013;Canestrelli et al., 2014). Lamb et al. (2012) addressed this issue by specifying a flow-spreading angle (i.e., a flare) in the 1-D model, which allowed a small amount of water-surface elevation change during changing discharge regimes. Their results showed lateral spreading enhances deposition at the river mouth, which in turn, creates an adverse bed slope such that the water depth shallows towards the topset-foreset break (Chatanantavet and Lamb, 2014) (Fig. 3).  incorporated lateral spreading at the river mouth into a moving boundary and shock-capturing framework (Hotchkiss and Parker, 1991;Swenson et al., 2000) that allowed for simulations with progradation and transgression. Despite the addition of lateral flow spreading, these channel-width averaged formulations cannot produce bars or bifurcations at river mouths (see Section 2.3 on planform models).
Most 1D profile models for evolution of riverbeds and deltas have adopted the characteristic flood hypothesis, which assumes that morphodynamics are adequately represented by a single characteristic discharge. This discharge is often assumed to be the bankfull discharge, and this is captured by multiplying Eq. (12) by an intermittency factor that accounts for the effective reoccurrence interval of that flood event (e.g., Paola et al., 1992). But the assumption of a characteristic discharge potentially misses some important behaviors. As first envisioned by Lane (1957), non-uniform flow combined with variable-discharge floods produces patterns of erosion and deposition in the backwater zone that do not occur in constant discharge simulations. For example, at low discharges non-uniform flow models have concave-up water surface elevations, causing flow deceleration and sediment accumulation in the backwater zone Nittrouer et al., 2012) (Fig. 3). In contrast, during high flows, the water surface becomes less concave, which, combined with the horizontal-to-adverse bed slope, causes a downstream reduction in flow cross-sectional area, and an increase in flow velocity and bed erosion in the lower part of the backwater zone Chatanantavet and Lamb, 2014), similar to observations of bed scour on delta topsets (Nittrouer et al., 2011;Shaw and Mohrig, 2014). The net effect of multiple cycles of high and low flow events in morphodynamic models is a river bed with a horizontal to convex-up geometry, as compared to models that assume a characteristic discharge (with or without non-uniform flow) (Chatanantavet et Chadwick et al. (2019Chadwick et al. ( , 2020 accounted for sediment lost to lateral channel migration and overbank flows on deltaic lobes in a 1D model by assuming floodplain sediment is evenly spread across the floodplain. They modeled several distinct lobes using 1D representations for each lobe and concatenated the profiles to construct a 2D delta that grew through growth and abandonment of multiple lobes. Moodie et al. (2019) took a different approach to merge 1D deltaic river model with a radial averaged delta floodplain model to better simulate the long-term evolution of the Yellow River delta, China, in the presence of coastal reworking processes. These hybrid models incorporate the 1D, Cartoon of the modeling framework of the backwater zone showing lateral spreading of the river plume at the shoreline in planview that leads to adverse bedslopes (A) and non-uniform flows in cross section (B). Non-uniform flows drive a change from deposition during small floods to erosion during large floods that results in an upward convexity of the river profile. These patterns can be imprinted on trends of net aggradation from backfilling due to progradation or relative sea level rise, for example. After Lamb MP, Nittrouer JA, Mohrig D and Shaw J (2012) Backwater and river plume controls on scour upstream of river mouths: Implications for fluvio-deltaic morphodynamics. Journal of Geophysical Research, 117: F01002.
process-based physics for river channels, including backwater hydrodynamics, flow spreading, and delta-front progradation, but use geometric assumptions and mass balance to extend the results to 2D landforms. They are important links between 1D and 2D modeling by providing a framework that is computationally tractable and straightforward to generalize.

1D flux partitioning models
As described in Section 2.1, 1D profile models can be useful for understanding how delta topography evolves during delta growth. However, river deltas commonly have channel networks with multiple simultaneously active branches. Understanding how deltas partition fluxes, such as water and sediment, in their networks is important for predicting how the planform shape of a river delta will change. Recognizing this, there is a class of 1D models for river deltas with the goal of understanding how fluxes of water and sediment are partitioned through the distributary network and how those fluxes lead to bifurcation stability or instability (Section 3.4 covers applications of these models).
In deltaic networks, there are two primary channel junctions: confluences, where two channels join to form a larger channel, and bifurcations, where a channel splits into two smaller downstream distributaries (Fig. 4A). From a 1D modeling perspective, bifurcations pose the greater challenge because the partitioning of sediment and water between the two downstream branches is not known a priori. Here, we review 1D models for the flux partitioning of river bifurcations. The first 1D model flux-partitioning model was developed by Wang et al. (1995). In this model, water discharges into each downstream branch are assumed to satisfy the equations of steady and uniform flow (i.e., Eq. 3). Additionally, the water discharge through each branch must satisfy continuity, so that: where Q is the volumetric water discharge, and subscripts a, b, and c refer to the upstream and two downstream branches, respectively. Finally, the water depth in Eq. (3) is found assuming a constant water level across the entrance to the two downstream branches: where η is the channel bed elevation, and h is the channel depth. In order to compute the sediment fluxes into each branch, Wang et al. (1995) assumed that volumetric sediment fluxes follow a similar continuity relationship as the one Eq. (13), and specified the following relationship for the ratio of sediment fluxes between the two downstream branches: where Q s is the volumetric sediment flux, w is the width, and k is an empirical parameter. Finally, sediment flux exported out of the downstream branches was computed from a sediment transport relationship in the form of Eq. (5), with some exponent m, and τ c * set to zero. The balance between sediment supply and export in each downstream branch controls how the depths and flux partitioning evolves through time. Based on linear stability analysis, Wang et al. (1995) found that if , a symmetric bifurcation is unstable, implying that any slight perturbation in the discharge partitioning grows through time, leading to the complete capture of the flow by a single downstream branch. The downside of the above approach is that the physical mechanism controlling k is unclear, and estimates obtained from different settings differ drastically (Kleinhans et al., 2013). Slingerland and Smith (1998) took a different approach by relaxing the assumption of normal flow and assuming sediment is transported entirely in suspension. The model consists of a crevasse channel that branches off from a primary channel. The bottom elevation of the crevasse channel is offset from the bottom elevation of the primary channel. As a result of the bed elevation offset, the flow into the crevasse channel was derived from the upper portions of the water column, where the sediment concentration is relatively dilute. The mean concentration of the incoming flow to the crevasse channel is obtained by integrating the vertical concentration profile from the height of the crevasse channel bottom to the water surface. By assuming constancy of water level at the bifurcation point (Eq. 2) and gradually varied flow in each branch (Eq. 10), the water and sediment discharge into the crevasse channel is obtained.
A third approach that is most suitable for bedload-dominated bifurcations was developed by Bolla Pittaluga et al. (2003) (hereafter BRT). The model consists of an upstream branch a and downstream branches b and c (Fig. 4B). Under the assumption of Eq. (14), steady and uniform flow into each branch is computed via Eqs. (3a-c). The region just upstream of the bifurcation is divided into two laterally interacting cells with bed elevations η b and η c . Each cell has its width set equal to that of the branch immediately downstream, and its length set to αw a , where α is an order-one parameter. The transverse water discharge Q y between the two cells is obtained from continuity, assuming the flow entering the divided reach is laterally uniform. The combination of the transverse water flow Q y and the presence of a transverse bed slope between the two cells sets up a transverse volumetric sediment flux Q sy between the two cells, given by: where r is a parameter that sets the strength of bed load steering due to a transverse bed slope (Ikeda, 1982), τ a * , is the upstream channel Shields stress, and ∂ η/∂y is the transverse slope evaluated from η b and η c . Finally, mass conservation equations for cells b and c is: where Q sa is the total volumetric sediment flux into the two upstream cells, Q si is the sediment flux at the downstream end of cell i. Depending on model parameters, notably width-to-depth ratio and τ a * , a bifurcation with symmetric flux partitioning may be stable or unstable (Bolla Pittaluga et al., 2015). When it is unstable, the bifurcation will tend to evolve towards one of two reciprocal stable asymmetric equilibria.
The BRT model captures the dynamics of flow partitioning and stability of bifurcations with minimal complexity, and it has been extended in a variety of ways.

2D morphodynamic planform models of delta dynamics
Two-dimensional planform models are designed to resolve the shape of the two-dimensional deltaic surface and are primarily used to explore controls on delta growth, morphology, and channel network dynamics (Fig. 5). Some planform models are based on mass balance and laterally average the topset topography (Kim et al., 2009a,b). These models use many of the same governing equations and constraints specified in Section 2.1.1., and will not be discussed here. Instead, we focus on those planform models that resolve the channel network on the topset. These models are advantageous because they resolve the dynamics of channel formation and do not require models to determine flux distribution among those channels. However, that advantage comes at a steep computational cost. There is a rich history of planform deltaic modeling that was summarized by Overeem et al. (2005), and we review models developed since that time. In particular, Delft3D has emerged as the standard tool for simulating 2D deltaic evolution, there have been new developments in rule-based or "synthesis" modeling, and moving boundary models have been adapted for deltaic growth in 2D.

Conservation-based models
Conservation-based models solve differential equations that describe the motion and conservation of mass, momentum, and sediment. The most widely used conservation-based model for deltaic studies is the open-source model, Delft3D (https://oss.deltares.nl/). Delft3D has been applied to simulate fluvial, deltaic, tidal, oceanographic, and coastal environments. Here we cover those aspects of Delft3D most relevant to deltaic modeling. Examples of deltas created in Delft3D are shown in Fig. 5C and D. Delft3D solves the shallow water equations, which are composed of conservation of momentum in two orthogonal directions (x and y) and conservation of fluid mass: where U, V are the depth-averaged velocity vectors in the x, y directions, g is acceleration due to gravity, ƞ is bed elevation, h is flow depth, τ b is basal shear stress vector (τ bx , τ by ), ρ is fluid density, and ν t is the horizontal eddy viscosity. Delft3D solves these equations to compute flow characteristics (such as water depth, turbulence characteristics, flow velocities and directions) dynamically in time over a three-dimensional computational grid. Even though Delft3D can compute a velocity field that varies with depth, it is still essentially a 2D simulation since the vertical pressure is assumed to be hydrostatic (i.e., the shallow water assumption). There is a non-hydrostatic pressure option that is fully 3D, but it is not usually applied to deltaic simulations. Eq. 18(a-c) are solved with a finite-difference solution, and when in three-dimensional mode there are a variety of turbulence closures model (e.g., kappa-epsilon model). In two-dimensional mode there is also a horizontal large eddy model that captures the transport of mass and momentum by large eddies.
In its sediment transport and morphology modules, Delft3D calculates the amount of bed load and suspended load transport of non-cohesive and cohesive sediment, considering the interchange of sediment between the bed and water column. Bedload fluxes can be calculated from a variety of equations, and Eq. (5) serves as one example. Bedload fluxes and directions are modified by the bed slope and the user can choose the specific formulation. Inappropriate values of the bed slope effect can have drastic results on morphodynamic predictions from models like Delft3D (Baar et al., 2019). Suspended sediment fluxes are modeled by the two or three-dimensional advection-diffusion equation. For suspended load, the entrainment or deposition of sediment is determined by the differences among the actual sediment concentration and the equilibrium concentration. The suspended load entrainment makes distinctions between cohesive and noncohesive particles. For cohesive sediment particles, during erosion and deposition the exchange between the water column and the bed are calculated with Partheniades-Krone equations considering user-defined thresholds for erosion and deposition (Partheniades, 1965). For the non-cohesive fraction, sediment entrainment is commonly modeled (but other choices are available) with the van Rijn method (van Rijn, 1993), which depends on the settling velocity of the particle and the shear velocity.
Morphological adjustment in Delft3D models is typically slow. To speed up computation time a morphological scale factor is applied to the bed elevation change from one time step (Ranasinghe et al., 2011). Thus, for a morphological scale factor of 10, this effectively assumes that the bed elevation change is constant and does not affect the flow for 10 time steps. This is a reasonable assumption when time steps and bed elevation changes are small. Ranasinghe et al. (2011) provide guidance in how to set this factor. That said, newer approaches adaptively solve for the maximum acceleration factor, which reduces computational time (Carraro et al., 2018).
In Delft3D users can specify the effects of vegetation on flow and sediment transport. One method to model vegetation is using Baptist's formulation (Baptist, 2005)

Modeling with synthetic flow-routing
In the previous Section (2.3.1), the models of river deltas simulated hydrodynamics by solving an approximation to the Navier-Stokes equations, and represented sediment processes with relatively high resolution in time and space. In this section, we review models using alternative parameterizations-called "rules" in some contexts-for routing flow and sediment across the delta at bifurcations or avulsions. Rule-based models have been called "cellular models" (e.g., Murray and Paola, 1994) because of the similarities with cellular automata.
These rule-based models attempt to explore processes and interactions on scales similar to the phenomena of interest (Murray, 2003). Because such models can involve parameterizations meant to synthesize the effects of dynamics on much smaller scales-which typically involve many more degrees of freedom than the phenomena of interest-they have also been called "reduced complexity" models (e.g., Murray, 2007), or "appropriate complexity" (French et al., 2016;Larsen et al., 2016), or "synthesist" models (Paola, 2000). These models have been contrasted with "process-based" models, which in practice typically means models based on explicitly representing fluid dynamics by solving Eqs. (18a-c), for example. However, since all the models we include in this review treat processes (represented in various ways and resolving various scales), we avoid this terminology. In addition, all scientific models reduce complexity, relative to the myriad processes on many scales in natural systems, so we minimize the use of the term "reduced complexity," preferring to use "synthesis." For modeling river deltas, synthesis models are typically motivated by the hypothesis that the essential aspects of fluid flow, for a particular context, can be represented in simple ways, like basing flow routing on bed slopes, or with flow relations like Eqs. (3a-c). In the context of deltas, these models often combine synthetic parameterizations with 1-D profile models (like those presented in Section 2.1), and channels are typically sub-grid scale, which means flow and sediment transport are not resolved within channels and instead are represented by synthetic parameterizations. The synthesis approach has two distinct advantages: it is computationally efficient, and it has explanatory power (Murray, 2003). Researchers have used synthesis approaches to investigate delta evolution, on the scales of individual lobe evolution to multiple lobes, for both reasons.  (5) to represent normal flow and sediment transport within channels, respectively, while also using a cellular, rule-based approach for establishing flow paths in the channel network. They used a flow-routing scheme introduced by Murray and Paola (1994), based only on local bed slopes, to determine how discharge is distributed where channels branch. Channel branches develop at avulsions, which are triggered when the difference between the elevation in a channel cell and the elevation in a neighboring cell becomes sufficiently large, relative to channel depth. After an avulsion occurs, a new channel follows a path of steepest-descent, and both branches of the channel initially remain open. The evolution of the longitudinal elevation profiles within both branches of the network, combined with the water routing scheme at the branching point, determines whether one of the branches later becomes inactive.
In a more abstract exploration of delta-lobe morphodynamics, Seybold et al. (2007) uses relatively simple parameterizations for flow and sediment routing, and a separate one for erosion and deposition (see example output in Fig. 5A). Their approach to flow routing builds on that of Murray and Paola (1994) by including water-surface elevation as a variable in each cell, so that water distribution from one cell into neighboring cells depends on water surface heights and slopes, as well as on bed topography. Water surface heights in the model evolve, along with the flow paths, according to a phenomenological treatment. Another phenomenological parameterization determines erosion and deposition, as a function of water discharge and surface slope, with an adjustable relative weight. The relatively simple model framework allowed the authors to explore basic controls on lobe evolution, showing for example that: (1) simple topographic smoothing can mimic the effects of wave-driven alongshore transport; and (2) when erosion and deposition depend mostly on discharge (rather than water surface slope), lobes resembling the Mississippi River Delta birds foot result (Seybold et Liang et al., 2015a,b), water parcels move from cell to cell through a random walk weighted by the depths in neighboring cells and by a parcel-averaged local flow direction. Water surface elevations are updated using a parameterization based on a normal flow assumption. Bedload and suspended-sediment transport are treated separately, using standard equations combined with random-walk approaches like that used for water parcels. DeltaRCM captures enough of the essential physics, without explicitly representing conservation of momentum, to reproduce deltaic morphologies that resemble other deltas (Fig. 5B).
Lauzon and Murray (2018) used DeltaRCM to investigate interactions between vegetation dynamics and morphodynamics in delta evolution, adding parameterizations for: (1) vegetation growth and mortality, as a function of morphology and erosion/deposition rates; (2) the effects of vegetation on flow routing; and (3)  Several other delta models treat channel-scale morphodynamics as sub-grid-scale processes, and have components that are synthesis based, including those addressing delta evolution over the timescales of multiple avulsion and lobe switching cycles described in

Moving boundary models of channel network growth
The 2D models discussed in Sections 2.3.1 and 2.3.2 rely on discretization of equations or phenomena that are solved on a grid. The choice of grid-cell size and how the equations are discretized can influence the solution and leave the imprint of the grid in the result. Considering this, distributary channel network growth may also be modeled using moving boundary techniques. As described in Section 2.1.1, environments that can be separated into two distinct domains (i.e., delta and basin) can be modeled through the evolution of the boundary separating them (Swenson et assumption that topographic variability is negligible within the channelized and unchannelized domains, allowing interface behavior to be studied in detail. Recently, this moving boundary approach has been adapted to study river delta networks (Ke et al., 2019). Boundary movement in the moving boundary delta channel network model (MB_DCN) is a function of sediment transport across the channel network, which is a nonlinear function of τ * . An appropriate hydrodynamic model must be chosen for the channelized and unchannelized domains to determine τ * , such as, where Λ = C f *U, where C f is a nondimensional friction coefficient, U and H are characteristic velocity and flow depth scales, h is local flow depth. Eq. (19) is a simplification of shallow water flow, which follows a quadratic friction model (e.g., Eq. 4). However, Eq. (19) performs remarkably well for shallow, slow moving unchannelized flows (Van Oyen et al., 2012Coffey and Shaw, 2017). When the flow equation is linearized in this way, mathematical techniques, such as the boundary element method used in MB_DCN (Brebbia et al., 2012), can be used.
In MB_DCN, water surface slopes in channels are assumed to be negligible relative to those in the neighboring unchannelized region. This simplifying assumption has not been extensively validated, but flow paths away from channels are generally steeper than those along the channel (Hiatt and Passalacqua, 2017). These hydrodynamic constraints prevent MB_DCN from resolving deposition within channels or water surface gradients between channels. While these simplifications appear appropriate for rapidly prograding distributary channel networks, they cannot reproduce channel avulsion within the network (Edmonds et al., 2009).
From an initial condition of a semicircular network boundary, the MB_DCN model produces a branching prograding channel network (Ke et al., 2019). Water surface slopes and associated progradation rates are largest at channel tips and smallest in concave interdistributary bays (Fig. 5E). This causes channel extension as channel tips extend into the basin and further reduces shear stresses within interdistributary bays and along channel margins. The prograding network boundary is unstable, meaning that small variations in progradation rate are amplified, eventually producing a channel network boundary with significant perturbations. Hence, channel tips with slightly reduced convexity have significantly reduced progradation rates, and quickly become non-prograding "islands" or "interdistributary bays" between two new channels (Fig. 5E). The source of this instability has not yet been derived for MB_DCN.

What do morphodynamic models tell us about how deltas grow?
Models have played an essential role in untangling the causes and effects of river delta morphology and growth. With the suite of models reviewed here, scientists have tested preexisting ideas and generated new hypotheses for the origin of deltaic features that range from millennial-scale simulations aimed at understanding deltaic lobe switching to flood-scale simulations aimed at determining fluxes through a bifurcation. Because of this, we organize Sections 3.1-3.5 from large-scale, long time-frame topics to small-scale, short time-frame topics. In each section, we provide an overview of recent developments and conclude with key points summarizing the state of the research.

Delta lobe formation and backwater dynamics
One of the most obvious and important features of river deltas is that they are composed of discrete lobes (Figs. 1 and 6). Delta lobes generally fan out from a central apex, or avulsion node, which gives deltas their characteristic size and shape ( Fig. 6A and B). Because of this, the origin of the avulsion node is an important topic. Alluvial fans and fan deltas have similar shapes as deltas, but the avulsion node or fan apex is typically determined by a change in confinement, such as a canyon-fan transition, or a change in bed slope that can be tectonically controlled (Ganti et al., 2014). The change in slope forces sediment deposition at the fan apex, resulting in channel filling and repeated avulsions about that location (Reitz et al., 2010;Reitz and Jerolmack, 2012). But many deltas are part of low relief landscapes with no obvious change in confinement or tectonic control and still have an avulsion node (though see Hartley et al., (2017) for a different opinion) . In these cases, lobe sizes can range from decimeters at laboratory scale (Ganti et al., 2016a) to hundreds of kilometers at the scale of the Mississippi Delta (Fisk, 1944) (Fig. 6A and B). Lobe switching occurs via avulsion, and this tends to occur at a statistically regular interval for a given delta, which can range from decades to centuries ( Fig. 6A and B).
River avulsions are thought to occur when the riverbed aggrades to some critical elevation relative to a reference elevation rendering it unstable (Slingerland and Smith, 2004). Once setup for an avulsion, a triggering event shifts flow to the steeper, shorter path to the shoreline, which abandons the old channel while developing a new channel Ganti et al., 2016b;Salter et al., 2018). Therefore, river avulsion frequency depends on the channel filling timescale, given as T A~h /v a , where v a is the aggradation rate (Jerolmack and Mohrig, 2007). A testable hypothesis then is that h/v a is minimized at the delta apex, which creates repeated avulsions about a persistent node (Chatanantavet et al., 2012), and consequently that the thickness of aggradation (i.e., Δη = v a T A ) must be maximized there (Ganti et al., 2016b).
Morphodynamic models have been used to explore two main ideas for what controls avulsion location on deltas in the absence of forced confinement: backwater hydrodynamics and backfilling due to progradation. Although backwater hydrodynamics and backfilling we emphasize that they are not mutually exclusive. Indeed, backfilling of a channel network typically occurs for rivers undergoing progradation, just like a zone of hydrodynamic backwater occurs for all subcritical lowland rivers with floods. Backfilling avulsions are driven by aggradation due to river mouth progradation, or base-level rise. As a deltaic lobe progrades, the river channel lengthens, and therefore the channel must also aggrade to maintain a sediment-transport slope (Schumm, 1993;Edmonds et al., 2009). This aggradation process, which tends to start near the river mouth and propagate upstream, was termed "morphodynamic backwater" by Hoyal and Sheets ( 6C).
The possibility that backwater hydrodynamics control the avulsion node location was suggested by Jerolmack and Swenson (2007), who showed that the avulsion length, L A , measured as the distance from the shoreline to the avulsion node, tends to scale with the backwater length, L b = h/S (Fig. 6C). The backwater length is the characteristic length scale over which non-uniform flow occurs, and can be derived through a scaling analysis of Eq. (12) by comparing the spatial acceleration term (left-hand-side) to the pressure gradient term (first term on the right-hand) in the limit of small Fr (Paola and Mohrig, 1996). Using a 1D profile model, Lamb et al. (2012) showed that spatial acceleration during drawdown conditions typically extends some fraction of L b upstream from the river mouth, whereas spatial deceleration during low flow can extend multiples of L b , and both depend on Fr and the depth of water in the backwater zone as compared to the normal flow depth.
Chatanantavet et al. (2012) modeled bed level adjustments in the backwater zone using a 1D profile model and proposed that floods of variable discharge created an avulsion node in deltas. Their model included lateral spreading of flow at the river mouth to drive non-uniform flow under a series of floods scaled after the Mississippi River. They showed that the alternation of low flow events, which tend to deposit sediment in the upper part of the backwater zone, and high flow events, which tend to scour in the lower part of the backwater zone, led to a zone of preferential sediment accumulation midway through the backwater zone. They inferred that the avulsion timescale, T A , would be minimized at this location leading to a preferential avulsion node location that scales with L b . Their model has been used to explain the location of the avulsion node on the Mississippi River (Chatanantavet et al., 2012) and Yellow River (Ganti et al., 2014). Importantly, their model showed that under constant-discharge conditions, the river tended towards a state of near-uniform flow without a preferential zone of accumulation, similar to previous models (Parker et al., 2008). This result occurred because constant discharge flow tends to level the bed towards a graded state by depositing sediment in locations with deep and slow flow, and eroding the bed in areas with shallow and fast flow, until near constant depths and uniform flow occur everywhere. Overall net aggradation, which is a necessary condition for avulsions, occurred in their model due to backfilling from progradation and/or due to relative sea-level rise, but variable flood discharges were needed to focus bed aggradation at a persistent location that scaled with the backwater length. These necessary conditions were later confirmed in physical experiments (Ganti et al., 2016a,b) and in subsequent morphodynamic models (Chadwick et al., , 2020 To explore this contradiction, Chadwick et al. (2019) used a 1-D morphodynamic model to show that model initial conditions and the reference profile from which relative aggradation is measured to assess avulsion likelihood can have a major effect on model results (Fig. 7). The reference profile is not consistently defined; some have proposed it is the elevation of the channel levees, elevation of the floodplain, or neighboring abandoned-lobe topography (Bryant et al., 1995 slope that extends from the upstream boundary of the model to the downstream shoreline of the delta. As progradation occurs, and new land is created at approximately sea level, a kinked reference profile emerges and aggradation of a concave-up river profile is maximized at the kink (Fig. 7A). Avulsions occur near this kink because, for the river to aggrade by h due to backfilling, it must prograde by~L b , thus producing avulsion node locations that scale with L b even without variable discharge floods. Chadwick et al. (2019) sought to test the influence of the kinked-reference profile, and modeled multiple deltaic lobes through a spin-up phase, and then used riverbed profiles of abandoned lobes as a reference profile. They showed, like previous work, that under these conditions variable discharge floods are needed to drive a persistent avulsion node location, except for cases of extreme relative sea-level rise ( Fig. 7B and C). While the assumptions that determine how the reference profile evolves-i.e., depositional processes on the surrounding floodplain-is an important choice, there are still other important differences between the 1D profile models (Chadwick et

Delta lobe formation: Knowledge gaps and future research directions
The most recent models that produce backwater-scaled avulsion nodes now include a moving boundary at the topset-foreset break, non-uniform flow hydrodynamics, parameterized flow spreading at the river mouth, variable flood discharges, and cycles of multiple lobe development and abandonment that allow models to evolve beyond the memory of initial conditions (e.g., Chadwick et al., 2019; Moodie et al., 2019). These models provide a reasonably simple framework that incorporate physical processes that can produce avulsion locations and delta lobe sizes that compare to natural observations . Yet, there are several important areas for future model development and model intercomparison that are needed. Perhaps the most important need is to further explore the consequences of different model assumptions regarding processes that create the reference profile used to assess avulsion likelihood (i.e., the topographic evolution of levees, floodplain and abandoned lobes). The evolution of reference profile has as much control on the model prediction for avulsion location as the does the evolution of the channel bed (Fig. 7). Existing models use end member assumptions about floodplain deposition that lead to contrasting reference profiles. In one approach (Ratliff et al., 2018), floodplain deposition, via crevasse splays, is minimal, although coastal wetlands are assumed to keep up with sea-level rise, leading to reference profile with localized, high curvature-i.e., a kink. In the other approach (e.g. Looking forward, there is surprisingly little work done on lobe-scale dynamics in two-dimensional planform models. One could use synthesis models (see Section 2.3.2); for instance, DeltaRCM (Liang et al., 2015a) seems promising because it includes rules to represent non-uniform flow dynamics and captures backwater and drawdown dynamics (Liang et al., 2015b). Similarly, fully 2D morphodynamics models, like Delft3D, include the physics that result in backwater-scaled deltaic lobes, but delta-lobe resolving simulations have yet to be evaluated in terms of lobe size or avulsion location. Both Delft-3D model runs (e.g., Nijhuis et al., 2015) and DeltaRCM runs (e.g., Liang et al., 2016b) usually have small domains compared to the backwater length, and this may explain why major avulsions occurred near the upstream domain boundary in those simulations rather than at the backwater length.

Controls on delta lobe progradation
During conditions of modest seal-level change and continued sediment supply from the river, deltas typically prograde basinward. The dynamics, speed, and controls on progradation are often investigated with 1D profile models and 2D planform models.
The advantage of simpler 1D profile models is that the shoreline position through time (i.e., the progradation rate) can be solved for analytically under specific scenarios of sediment supply, sea-level variations, and geometric configurations. An important outcome of all these closed-form analytical solutions is that the sedimentary wedge maintains a self-similar shape, and consequently the alluvial-bedrock transition and shoreline trajectories can be described as follows: (20) where λ ab and λ sh are constants calculated through the solution of two algebraic equations (Lorenzo-Trueba et al., 2009;

Lorenzo-Trueba and Voller, 2010).
A key parameter to determine λ ab and λ sh , and therefore the rate of delta progradation, is the ratio of the fluvial to bedrock slopes at the alluvial-bedrock transition, i.e., R ab = q s /βv. This ratio is physically bounded within the range 0 ≤ R ab ≤ 1, and field estimates put it in the range 0.2 ≤ R ab ≤ 0.8 (Lorenzo-Trueba et al., 2009). With an appropriate definition of the diffusivity (see Eq. 6), R ab can be written as a function of the characteristic bed sediment size (gravel or sand), q w and q s . A low R ab value results in a faster rate of shoreline migration than the alluvial-basement transition, which in turn leads to a low subaerial sediment volume respect to the subaqueous portion (Fig. 8). A high R ab value, however, causes the alluvial-bedrock transition to migrate more rapidly than the shoreline, leading to a delta with more subaerial sediment volume than subaqueous (Fig. 8).
Another important parameter controlling the sediment volume partitioning between the topset and the foreset is θ, which is determined by the choice of sediment transport formulation (see Eq. 1). Similarly as θ increases, the curvature of the fluvial topset decreases, which leads to an even higher increase in the subaerial volume, particularly for R ab > 0.4 (Fig. 9). In the limit θ → ∞, the curvature of the fluvial topset profile is zero and we recover the linear geometric solution (Kim and Muto, 2007), which presents the highest rate of alluvial-bedrock migration relative to the rate of shoreline migration. Under the assumption of a linear basement slope as depicted in Fig. 2, the solution predicts a decelerating progradation rate as delta area grows and/or basin depth increases. But, delta progradation can accelerate, for instance, in "back-tilted" basins where the water depth shallows basinward (Hajek et al., 2014; Zhao et al., 2019).
These analytical predictions generally match laterally averaged progradation rates from more complex 2D planform models. Overeem et al. (2005) reviewed a range of synthesis 2D planform models of delta growth and showed that similarly to 1D models, progradation rate monotonically decreases over time as delta size increases. Of course, this is not surprising since both 1D and 2D models are based on mass conservation. But, 2D models allow spatially variable progradation rates (Syvitski and Daughney, 1992), and river avulsion (Hutton and Syvitski, 2008). More recently using Delft3D, Gao et al. (2018; 2019) established that intermittent discharge creates a similar progradation rate as more simplified models, and at low discharge periods the delta area actually shrinks because the accommodation space grows too large for the sediment supply.
Using 2D planform models, workers have also shown that parts of the delta shoreline can prograde faster or slower than the laterally averaged value. This happens when subsidence steers delta growth into areas of high accommodation, creating localized progradation (Liang et al., 2016a). Or this can occur when the delta has a distributary channel network because sediment flux is unevenly distributed at each bifurcation.

Delta lobe progradation: Knowledge gaps and future research directions
Over the past two decades, models of shoreline progradation have improved our quantitative understanding of the factors that affect delta growth, and more generally the evolution of sedimentary basins (Swenson et (Gao et al., 2019), and how parts of the shoreline might deviate from the laterally averaged rate (Caldwell and Edmonds, 2014).
Through these models we have gained a better understanding of a number of key controls on shoreline progradation such as sediment supply, water discharge, basement geometry, and relative sea-level changes. Despite these advances, however, there are significant knowledge gaps that will require additional effort moving forward. For instance, the modeling framework presented here neglects the interplay of physical and biotic processes of sedimentation that take place in the topset region, even though some 2D planform models account for this (Sections 2.3 and 3.3). In particular, variations in organic sediment dynamics in the fluvial topset can have a significant effect on the rate of shoreline progradation and retreat (Paola et al., 2011; Lorenzo-Trueba et al., 2012). The modeling framework presented here also neglects the effect of wave energy on shaping the morphology of the foreset, which can in turn affect the sediment partitioning between the topset and foreset regions (Swenson et al., 2005; Voller et al., 2020).
Moving boundary problems only have analytical solutions in a limited range of scenarios. For instance, closed-formed analytical solutions only can account for sea-level variations proportional to the square root of time, which are not necessarily representative of what happens in natural systems. Looking forward, two new categories of numerical methods have been developed for the dual alluvial-basement transition and shoreline moving boundary problem: fixed grid methods and deforming grid methods. Fixed grid methods employ a discrete domain that remains fixed in time and space, and require an auxiliary variable to track the boundaries. In particular, the enthalpy method tracks the shoreline based on the fraction of the water column filled by sediment (Anderson et al., 2019). In contrast, in Fig. 9 Solution space for λ ab and λ sh against 0 ≤ R ab ≤ 1 for different values of Θ. λ ab and λ sh are constants that determine the rates of migration of the alluvial-bedrock transition and the shoreline respectively. In the limit Θ → ∞, the modeling framework recovers the linear geometric solution presented by Kim and Muto (2007).
deforming grid methods the discrete domain evolves to ensure that there is always a line of nodes located at the moving boundary ).

Deltaic channel network formation
River delta channel networks are an obvious target for morphodynamic planform models because of the variety of channel networks that appear on modern deltas Syvitski and Saito, 2007;Tejedor et al., 2016Tejedor et al., , 2017. Network models are often an extension of river mouth bar models (see Section 3.5), and are generally sensitive to similar hydrologic and sedimentologic controls. Since 2005, there has been a proliferation of models that resolve network growth, in turn documenting the controls on network formation. This section will document the many similarities and the few key differences observed between these studies.
Planform delta models that resolve the channel network almost exclusively begin with a water and sediment source from a channel mouth directed into a water-filled basin. This condition produces rapid progradation of a deposit because of sediment mass balance (see Section 2.3.1). During this phase of rapid progradation, numerical models frequently produce a channel network that branches to varying degrees (Fig. 5). In a counter example, the MB_DCN (moving boundary delta channel network model) only resolves erosion along the network boundary, but also produces a branching network due to rapid initial progradation. The correlation of rapid initial progradation and network formation among diverse models is quite strong. In fact, the classic delta morphology ternary diagram was originally conceived as just two end members: constructional deltas with branching networks associated with significant sediment delivery, and destructional deltas where waves and tides removed significant material (Fisher, 1969;Galloway, 1975).
In most planform delta models, channel branching starts at the delta apex, where several coeval branches are maintained for long time periods (Fig. 5C and D). This morphology can be repeated as distributary channels branch again as they reach parts of the basin with little sediment accumulation (Fig. 5B). However, recent modeling shows that branching can arise from different processes. For example, Delft3D-based models and the DeltaRCM suggests channels branch when they have stable turbulent jets that favor mouth bar deposition (Edmonds and Slingerland, 2007;Falcini and Jerolmack, 2010;Caldwell and Edmonds, 2014;Canestrelli et al., 2014;Liang et al., 2015b). The recent moving boundary model MB_DCN does not resolve deposition nor jet instability and still produces branching channels (Ke et al., 2019). In that case, network extension is driven by the distribution of erosion rates at the channel network boundary (that includes the subaqueous channel margins as well as the tip ; Fig. 5E). The maximum erosion rate is often at the channel tip (location of maximum curvature), allowing a single channel to extend. However, this interface may become unstable with two or more maxima, in which case multiple channels extend. Decreased curvature of the network boundary, increased water flux per unit width, decreased critical shear stress of sediment motion, and non-local effects of the network shape appear to promote this instability. Once channels branch and form a network, planform delta models consistently show that network morphology is sensitive to sediment characteristics. For channels carrying sand and coarser sediment, the delta usually has many channels that cut and fill, and migrate laterally, whereas channels transporting a finer load, create deltas with fewer, more elongated channels ( , 2015b). This control is tied primarily to the dimensionless critical stress of erosion τ cr * , which is largest for cohesive, fine grained sediments. In MB_DCN, the same effect is observed but for a different reason. In that case, high τ cr * prevents branching because it suppresses the creation of multiple maxima of τ * (Ke et al., 2019). Van der Vegt et al. (2016) showed that the proportion of bedload within the incoming sediment load can control network formation and morphology. They observed that morphological patterns created by increasing sediment cohesion also can be generated by decreasing the proportion of bedload sediment.
Network dynamics are also controlled by vegetation. Models have shown that a certain vegetation density maximizes floodplain sand deposition rates. This is because the presence of vegetation along channel levees and in islands increases friction, and the proportion of water and sediment in the channel increases (  . Using DeltaRCM with a vegetation routine, they found that cohesive sediment and vegetation have similar effects on channel and shoreline shape, and that for deltas built from sand-rich inputs, the presence vegetation was comparable to increasing the input cohesive sediment by approximately 25%. In contrast, vegetation exerted a stronger control than cohesive sediment on channel-network stability. In addition, Lauzon and Murray (2018) found that the presence of vegetation led to lower average delta elevations because vegetation forced more deposition at the shoreline-a finding consistent with Nardin et al. (2016) (Fig. 10). In a recent study, Wright et al. (2018) show that patchy vegetation that is below the disconnectivity threshold would still allow sediment onto island. Importantly though, their model does not include sediment. However, they suggest that not accounting for patch dynamics may underestimate the amount of sediment transported onto deltaic islands.  (Liang et al., 2016b) and Delft3D (Van de Lageweg and Slangen, 2017) have been used to show that relative sea level rise (RSLR) can also influence channel dynamics. When a "filling index," defined as the ratio of accommodation space creation (area times RSLR) to sediment discharge, falls below 0.6, distributary channels tend to reduce in lateral mobility (Liang et al., 2016a,b). Deltas must eventually become inundated for a filling index of less than one. Hence, it is these "drowning" deltas that can be expected to have channels with reduced mobility. Van de Lageweg and Slangen (2017) model the co-evolution of sea level rise and morphodynamics in Delft3D. During sea level rise they allow the network to dynamically adjust to compensate for the accommodation space created. On average a dynamic network-that is, one that is allowed to adjust-lost 20% less land than a static network (Van de Lageweg and Slangen, 2017). They also found that the type of delta morphodynamics, i.e., river, wave, or tide-dominated, has no effect on land loss from sea level rise.

Channel network formation: Knowledge gaps and future research directions
Simulations of channel network formation have progressed considerably since the first channel-resolving simulations of planform delta growth (Overeem et al., 2005). The proliferation of models has allowed researchers to investigate how channel network formation depends on sediment characteristics   (Lauzon et al., 2019).
From these models, we know that the controls of network growth under steady conditions (constant sea level, no subsidence, constant discharge inputs) appear to be basin shape, sediment characteristics, and water discharge. Controls on the network evolution under unsteady and non-uniform conditions are far less known, which is a major knowledge gap. From an unsteady flow perspective, it is likely that the trends of channel formation under increasing flow and channel siltation under decreasing flow will be borne out . However, the details of this response likely depend on the scale and connectivity of an existing network's structure. A myriad of conceivable changes to sediment supply (e.g., incoming grain size changes) and human modification (e.g., dredging and artificial levees) could significantly affect an evolving delta, and could be systematically explored with numerical models. Vegetation patchiness is also a non-uniform condition that has not received much attention. Hydrodynamic models show that patchy vegetation creates a stress-divergence feedback that leads to channelization of surfaces (Temmerman et al., 2007;Yamasaki et al., 2019). An important challenge is to develop ecogeomorphic models capable of simulating patch-scale vegetation dynamics and delta evolution. Early in delta development, vegetation may enhance network formation via the stress-divergence feedback, whereas, on older surfaces, vegetation may limit deposition by forcing water and sediment to stay in channels.
Another major knowledge gap and a promising future research direction is determining what controls channel branching (Edmonds and Slingerland, 2007;Ke et al., 2019). For example, it is unknown if branching is controlled by the shape and progradation rate of a channel tip ("network controlled") or by a depocenter that splits flow ("deposition controlled"), or if different sources of control exist for different boundary conditions. Such unknowns could be addressed with systematic numerical experiments and field observations.

Stability of river channel bifurcations
Once a channel network develops, the longevity of that network (at least on timescales shorter than lobe switching) depends on the morphodynamic evolution of its bifurcations. Most modeling studies of bifurcations predict which configurations of flux partitioning among the downstream branches are stable versus unstable. The existence of stable configurations implies that bifurcations with a range of different initial conditions evolve towards a similar state. Bifurcations in nature are typically asymmetrical, wherein one channel is wider and deeper than the other. Asymmetric bifurcations can be remarkably long-lived, with some transmitting water and sediment for nearly 3000 years , suggesting they might be a stable configuration.
A majority of studies on flux partitioning and stability of bifurcations adopt a 1D approach (see Section 2.2 for more detail). Bolla Pittaluga et al. (2003) showed that the number of equilibrium solutions and the stability of those equilibria for a bifurcation with symmetric geometry and forcing conditions depends on the upstream channel width-to-depth ratio w a /h a . For low w a /h a , the only equilibrium solution is a stable symmetric partitioning (Fig. 11A). In contrast, for high w a /h a , the symmetric equilibrium is unstable, meaning that any small perturbation away from symmetry tends to grow. However, this growth does not continue indefinitely; eventually the transverse sediment flux becomes high enough to prevent complete flow capture by a single branch, and the bifurcation settles on one of two stable reciprocal asymmetric equilibria. Bolla Pittaluga et al. (2015) showed that the role of the Shields stress in the upstream channel τ a * is more complex than that of w a /h a . Using sediment transport formulae derived for bedload-dominated rivers, the symmetric equilibrium changes from stable to unstable as τ a * is decreased. However, when using formulae for high-mobility rivers, the relationship is reversed; increasing τ a * causes the symmetric equilibrium to transition from stable to unstable, in agreement with Delft3D modeling by Edmonds and Slingerland (2008). Finally, Salter et al. (2018) showed that downstream deposition tends to stabilize bifurcations, thereby introducing a dependence on downstream channel length and bypass. One challenge in using the Bolla Pittaluga et al. (2003) model to predict bifurcation dynamics is the dependence on the parameter α, which sets the length of the divided reach immediately upstream of the bifurcation (see Fig. 4B for definition of divided reach). Conceptually, setting α is a way of parametrizing how bed elevation information from the downstream branches propagates upstream into the feeder channel, suggesting a connection to the theory of morphodynamic influence developed by Zolezzi and Seminara (2001). According to this theory, the aspect ratio ( ) sets the dominant direction of morphodynamic influence: when φ is below the proposed setting the value α such that the neutrally-stable aspect ratio coincides with φ R , thereby eliminating the need to set α through calibration. Their results show that resonance plays an underappreciated role in determining bifurcation stability. Next, we consider how asymmetric forcing affects bifurcation dynamics. One such forcing is a difference in downstream branch slopes. Bolla Pittaluga et al. (2003) explored changing the downstream slope ratio S b /S c . They found that for high τ a * , there is a single stable equilibrium solution, with a discharge ratio Q b /Q c that increases smoothly as S b /S c increases. In contrast, for low τ a * , they found three equilibrium solutions for S b /S c close to one, but for high S b /S c the only equilibrium solution was one where branch b captured most of the flow. Another type of forcing is an upstream meander bend. Considering this, Kleinhans et al. (2008) argued that spiral flow in a meander bend steers near-bed sediment from the outer bank towards the inner bank. As a result of this steering, the downstream branch along the inner bank receives more sediment relative to its discharge, causing it to silt up, whereas the bifurcate adjacent to the outer bank tends to deepen . A slope advantage for the inner bifurcate can compensate for the meander bend, allowing the bifurcation to become "quasi-balanced." A third type of forcing is a difference in the angles of the downstream branches with respect to the upstream branch. Classic experiments by Bulle (1926) showed that bed load is steered disproportionately towards the high-angle branch.
In a particularly noteworthy advance, Redolfi et al. (2019) developed a unified framework for bifurcations with both asymmetric and symmetric forcing. Under sub-resonant conditions (low width-to-depth ratio) the flux partitioning is completely controlled by the asymmetric forcing conditions (Fig. 11B). There is a single stable solution that becomes increasingly biased as the forcing strength increases. However, under super-resonant conditions the symmetric and asymmetric responses coexist. For modest forcing strength, three equilibria are present, each of which changes smoothly as the forcing is increased. However, if the forcing becomes strong enough, two of these equilibria merge and then disappear, and the only remaining equilibrium is one where the branch favored by the forcing dominates. This behavior introduces the possibility of catastrophic jumps and hysteresis as the forcing conditions are varied slowly back and forth from favoring one branch to favoring the other.
How showed that migrating meander bends upstream of a bifurcation introduce oscillations in the flow partitioning. In general, for an unforced bifurcation with a single symmetric stable equilibrium, adding periodic forcing introduces small oscillations about the equilibrium state. On the other hand, when the unforced bifurcation has two stable asymmetric equilibria, small forcing results in small oscillations about one of two stable equilibria, but large forcing is enough to overcome the potential barrier between the stable equilibria, causing large oscillations between reciprocal asymmetric states.
Most bifurcation studies have assumed a constant bed or water level as the downstream boundary condition, which causes the bifurcation to evolve towards a bypass state where incoming and outgoing sediment fluxes are in balance. However, river deltas are typically net-depositional. To explore the role of deposition on bifurcation dynamics, Salter et al. (2018) used two types of downstream boundary conditions: progradation with a constant sea level, or steady vertical sedimentation with fixed branch length and a specified bypass fraction. In both cases, even without time-varying external forcings, autogenic oscillations in the flux partitioning occurred due to coupling between upstream and downstream feedbacks; the upstream feedback is the tendency for a bifurcation to evolve towards an asymmetric flow partitioning, and the downstream feedback is that the slope of the branch receiving more sediment becomes lower and lower relative to the other branch, eventually causing an avulsion to occur (Fig. 12). When one branch receives more flow, it progrades faster, and eventually far enough that the flow switches to the opposite branch. This cycle repeats, and causes oscillations in the discharge partitioning. These oscillations continue until the branches arrive at the downstream edge of the domain, where the final cell of the domain is an "infinite sink," e.g., a shelf edge that absorbs all incoming sediment. At that point, the system evolves towards a bypass state and attains a "frozen" asymmetric discharge partitioning. Depending on model parameters (e.g., w/h a and downstream branch length), also explored adding asymmetric forcing to the model in the form of differential subsidence and found that small forcing introduced bias towards one side during an avulsion cycle, whereas sufficiently large forcing changed the dynamics qualitatively, causing autogenic oscillations to shut off. suspension. An open problem is to integrate these two models into a common unified framework for modeling bifurcations with both bed load and suspended load.

Bifurcation stability: Knowledge gaps and future research directions
The two-way feedbacks between bifurcation flow partitioning and planform geometry are a promising area of future research. As mentioned in Section 2.2, Miori et al. (2006) provides a starting point for including width-adjustment within 1D bifurcation models, but it relies on a poorly understood timescale for width adjustment. The timescale of width adjustment could significantly impact flux partitioning through a distributary network and may ultimately contribute to the kind of soft avulsions created by Salter et al. (2018) (Fig. 12). These results suggest that fluxes through a channel network dynamically evolve through time, without channel closure, because channels lengthen, and widths adjust. Outside of the single bifurcation case, this has never been shown for a network. Presently the modeling tools exist to create a new class of 1D models that address this problem. This would likely require 1D bifurcation models to

River mouth bar dynamics
The distal tips of distributary channels are sites of morphodynamic change due to accelerating and decelerating sediment-laden flow fields and the corresponding bed level changes. Many modeling studies focus on simulating the dynamics at this key point because this is where delta growth originates. Indeed, at the mouth of most distributary channels there is a triangular shaped deposit that is variably referred to as a river mouth bar, middle ground bar, or distributary mouth bar. Numerical models based on the solution of fluid and sediment transport equations have helped unravel how flow hydrodynamics affect the formation of river mouth bars (Fagherazzi et al., 2015).
The shapes of mouth bars are related to the stability of the turbulent jet (Fig. 13). A turbulent jet forms at river mouths when the river effluent transitions from confined to unconfined and turbulence generated at the jet margins exchanges momentum with the ambient water that results in lateral jet spreading. During turbulence generation, some of the energy can be transferred to larger eddies (rather than the inertial subrange) leading to an unstable, meandering jet Giger et al., 1991;Jirka, 1994 showed that the river jet may be temporally stable or unstable (meandering), depending upon the jet stability parameter and inlet Reynolds number. Whether stable or unstable, the jet determines sediment particle flow paths and bed deposition patterns in front of the river mouth. Using a 2D depth-integrated model, Canestrelli et al. (2014) found that unstable jets tend to focus deposition on the jet margins preferentially creating levees, whereas stable jets preferentially create mouth bars similar to the central bar in Fig. 14.  Fig. 13 Examples of stable (top) and unstable (bottom) jets created in Delft3D. Instability is created by decreases the aspect ratio and friction of the channel. Colorbar shows fluid velocity (m/s). After Canestrelli A, Nardin W, Edmonds D, Fagherazzi S and Slingerland R (2014) Importance of frictional effects and jet instability on the morphodynamics of river mouth bars and levees. Journal of Geophysical Research: Oceans. doi: 10.1002/2013JC009312. Using Delft3D, Edmonds and Slingerland (2007) suggested once a mouth bar forms there are four stages of bar growth: (a) parallel subaqueous levees are deposited along the edges of the jet; (b) subaqueous levees lengthen basinward and the river mouth bar begins to prograde; (c) the subaqueous river banks continue extending basinward and flare laterally around the now stagnant bar; d) the river mouth bar emerges creating a deltaic island which might be colonized by vegetation. River mouth bars develop within the river jet where the negative gradient of longitudinal sediment flux is steepest. These gradients originate from a decreased velocity with a subsequent drop in sediment transport capacity. Canestrelli et al. (2014) showed in numerical experiments that this occurs at a distance within two channel widths from the river mouth for stable jets and at higher distances for unstable jets.
Modeling has also been used to study how waves and tides influence the spreading of turbulent jets and subsequent mouth bar development (Fig. 14, wave-dominated). Waves can suppress mouth bar formation or change the direction of progradation, resulting in largescale differences in bar shape (Nardin and Fagherazzi, 2012;Nardin et al., 2013). This happens most drastically when the waves approach the shoreline at a high angle. They showed that wave-driven alongshore sediment transport can change bar deposition, and instead create a spit that spans the width of the river. Spits form when the fraction of sediment bypassing the river mouth is large, and the bypassing fraction depends on the relative momentum flux between the jet and waves. In contrast, when waves are locally generated and perpendicular to the shoreline they stabilize the jet and favor deposition near the river mouth (Nardin et al., 2013).
The morphologic impact of tidal processes on mouth bars also has been investigated in Delft3d simulations (Leonardi et al., 2013). Tides enhance jet spreading and lead to faster bar aggradation. Leonardi et al. (2013), showed that depending on the relative strength of river inertia with respect to tidal energy, different hydrodynamic processes dominate the sediment deposits with consequent development of distinct morphologies. Tidally dominated systems are characterized by a predominance of coast-normal, elongate tidal bars, while wave-generated, coast-parallel barriers and/or beaches are rare. When the tidal discharge is much higher than the fluvial discharge, elongated deposits form due to tidally induced bidirectional sediment transport, parallel to the river outflow. The fluvial-dominated case is typical of systems having a negligible tidal prism or during floods when the river flow is high enough to prevent flow reversal. Under these conditions the presence of tides increases jet spreading and results in a progressive water wave near the river mouth, with maximum velocity occurring at low tide.
Numerical experiments also have addressed how bed slope, unsteady discharge, and vegetation influence mouth bar growth. Jiménez-Robles et al. (2016) investigated the role of the basin slope and the river discharge angle (relative to the shoreline) on the hydrodynamics of turbulent jet and a mouth bar morphodynamics using Delft3D. They found that when varying the discharge angle and the receiving basin bottom slope, the turbulent jet changes its direction and spreading with consequences for the river mouth bar morphology. Their model created three distinct bar morphologies: central-bar, side-bar and sand-spit. The central bar develops in front of rivers debouching with low discharge angles. Higher values of this angle produced a lateral bar or a sandy spit. The impact of the basin slope was strongest on depositional style, rather than morphology. Steeply-sloped basin created more aggradation, whereas shallower basins results in bar progradation.
As mouth bars aggrade to the water surface, they can be colonized by vegetation. Vegetation seemingly should influence mouth bar growth. In case of submerged aquatic vegetation ( SAV meadows on river mouth bars morphodynamics. The effects of SAV on the river hydrodynamic is mainly derived by the water flow diverting from the bar to channels. Subsequently, sediment deposition occurs at the bar ridge margins promoting upstream bar migration. Seasons with high SAV density enhance the development of sandy bars closer to the river mouth, also, with implications on the bar morphology. Sediment characteristics also impact those river mouth bar formations. Sediment size is the main control because settling velocity or mode of transport influence where particles leaving the river mouth are deposited. Delft3d experiments by Caldwell and Edmonds (2014) showed that incoming sediment mixtures that are coarse in size, promote mouth bar growth because those sediments are preferentially deposited closer to the river mouth. If the incoming mixture is finer-grained, the sediment was carried into the basin by the jet rather than forming mouth bars.

Mouth bar dynamics: Key knowledge gaps and future research directions
Numerical studies reveal that river mouths are dominant sites of sediment deposition (Mikhailov, 1966;Edmonds and Slingerland, 2007;Fagherazzi et al., 2015). One topic that deserves more attention in numerical simulations on mouth bars is the role of erosion (Shaw and Mohrig, 2014; Ke et al., 2019). Afterall, erosion of previously deposited sediment at the channel mouth may be responsible for initiating branching and channel formation, rather than flow splitting around depositional mouth bars. Understanding the relative roles of deposition and erosion in setting mouth bar dynamics and channelization at the delta front is an exciting future direction.
The nature of the flow field at the river mouth remains uncertain. Most theory and modeling predict a turbulent expanding jet (Bates, 1953). But recent work suggests the flow may be diffusive with minimal turbulent mixing (Coffey and Shaw, 2017;Ke et al., 2019). The notion that river mouths have turbulent jets stems from engineering experiments  as well as field data (Wright, 1977;Ernst et al., 1996). The engineering experiments, however, often do not simulate natural deltaic river mouth conditions. For example, they commonly consider channels with unnatural aspect ratios. Also, most simulations aimed at understanding deltaic deposition do not account for buoyancy effects on the turbulent jet. At flood stage, the saltwater wedge is pushed out past the river mouth, but at lower discharges buoyancy may influence deposition and erosion patterns on the mouth bar. A key knowledge gap to address is collecting field and numerical data that simulate river mouth hydrodynamics under realistic initial and boundary conditions that include buoyancy. This would allow one to determine whether turbulent jets or diffusive flow dominates hydrodynamics on the delta front under natural conditions.

Future needs in morphodynamic modeling of river deltas
Over the last few decades, the collective interest in deltaic modeling has increased significantly. This interest has spurred model development to the point where state-of-the-art models can simulate dynamic deltas that create channel networks through avulsion and bifurcation, and respond to forcings like relative sea-level rise. Most impressively, this level of simulation can be accomplished through a range of models-from simple moving boundary formulations, to synthesis models, to more detailed models like Delft3D.
We have purposefully focused on models aimed at understanding how deltas work from a generic perspective; that is, most of the studies we reviewed did not model any particular delta. These types of generic modeling studies advance new, field-testable hypotheses for how deltas work. In each subsection of Section 3, we summarized the key knowledge gaps and future research directions relevant to those topics. Here we focus on two major questions or needs that we think are required to catalyze the next generation of deltaic models.

New models to bridge the time and space scale gaps
Deltaic models are presently capable of simulating long time-scale processes, such as lobe switching and avulsion in 1D, and short time-scale processes, such as network formation and mouth bar dynamics in 2D. However, we presently lack 2D models that simultaneously capture long and short time-scale processes. These kinds of models could fill a much-needed gap that would also allow researchers to address new problems. For instance, a 2D model that simulates channel network formation and lobe switching, would have the advantage of directly representing the reference profile (see Sections 2.1.2 and 3.1), which would allow avulsions to emerge from the morphodynamics rather than having to parameterize them, as it currently done in 1D models. These 2D models would also be a powerful addition to the stratigraphic community because they would have long time-scale processes important for sedimentation, like lobe switching, but also resolve the detailed facies architecture created by channels. For example, we know from experiments and field data that deltaic lobes stack and shingle in space (Blum and Roberts, 2012;Ganti et al., 2016a), but that behavior has never been reproduced in a 2D model that also resolves channels.
As computing power grows, it may be possible that current planform models, like Delft3D, can be adapted for this purpose. But we believe these goals can be more effectively achieved by building on community models that have already been developed. In this sense, modeling tools and the basic model interface created by the Community Surface Dynamics Modeling System (CSDMS) may be useful. CSDMS allow users to link models together, possibly creating models of intermediate complexity. For example, CSDMS could allow researchers to marry the essential physics from a Delft3D-style simulation with a simpler synthesis model for avulsion and lobe emplacement.

An integrated understanding of how fluvial and marine forces affect delta morphology
Since the seminal paper of Galloway (1975) it is often claimed that relative influences of rivers, waves, and tides determine delta morphology. We are unaware of any single model that can reproduce the diversity of river, wave, and tide-dominated features that Galloway tried to explain. This is in part because most modeling has focused on simulated end-member behavior. Though we did not review studies that incorporate waves and tides, those that do exists also rely on simplified end-member models (Ashton and Giosan, 2011;Nienhuis et al., 2015Nienhuis et al., , 2016Nienhuis et al., , 2018aHoitink et al., 2017). We see the lack of mixed process deltaic modeling that accounts for waves, tides, and rivers as an important future direction of research.
There are two main problems that prevent progress on this issue. First, we currently lack quantitative definitions of what river, wave, and tide dominance mean. The Galloway (1975) diagram does not define those axes and this is a problem that still prevents it from being used to predict delta morphology. Admittedly, some of the studies reviewed here include waves and tides (e.g., Geleynse et al., 2011), but until we quantitatively define what river, wave, or tide dominance means it is difficult to place those studies in the proper context. Recently, Nienhuis et al. (2020) made an important advance on this problem and presented a process-based framework to quantify delta morphology. Their field data compilation demonstrates that delta morphology varies as a function of river, wave, and tide-dominance. Second, no single model has reproduced the variety of river-influenced features (e.g., mouth bars, avulsions, bifurcations), wave-influenced features (e.g., spits, beach ridges, and longshore bars), and tide-influenced features (e.g., funnel shaped river mouths, looping channel networks, and tidal bars). We are not sure if this is because existing models are incapable of producing these landforms or if searching through an undefined parameter space has slowed down progress. Likely, it is some combination of both. Either way, an important next step is to use the framework of Nienhuis et al. (2020) to explore the morphologies of mixed process deltas.

Conclusions
Numerical modeling has emerged as a useful tool to understand the evolution of river deltas. Because rivers deltas contain patterns and processes at a variety of scales in time and space, a suite of models have emerged to address these issues. Deltaic models of shoreline progradation and sediment erosion and deposition through the backwater zone can be accomplished with relatively simple one-dimensional profile models. These models have helped characterize how the aggradation or degradation of the profile is influenced by shoreline progradation rate, water and sediment supply, and offshore basin geometry. Moreover, these profile models been used to test compelling hypotheses for the origin of the avulsion node.
More complex two-dimensional planform models are capable of simulating mouth bar deposition, delta lobe growth and channel network formation. These 2D planform models arise from solving the detailed hydrodynamics and sediment transport (i.e., Delft3D), synthesis approaches, or through a moving boundary formulation. These models have shown that channel branching and the larger network arise from a combination of deposition and erosion at the channel mouth. For the river-dominated case considered here, the shape of the network and the delta depend on sediment properties (e.g., size, distribution, type), vegetation, and flow conditions including intermittent flow and flow strength.
The channel network is a key feature of river-dominated deltas and models have been developed specifically to understand how the key network nodes-the bifurcations-remain stable for a given flux partitioning. One and two-dimensional models are typically used to predict what partitioning of fluxes into the downstream branches results in a stable or unstable bifurcation. For simplified conditions of a straight channel bifurcating into two channels, models can predict a range of stable solutions that depend on the upstream flow strength, the sediment transport regime, and the asymmetry of the downstream forcings.
We see 2D models that simultaneously represent a variety of scales, from lobe switching to channel network formation, as a promising area for model development. Developing models that harness the power of detailed Delft3D simulations with the explanatory power of synthesis models should provide new insight. While these types of models in principle can be generated with the Community Surface Dynamics Modeling Systems tools, they have yet to be applied to river delta simulations.
It is worth noting that in almost all modeling approaches and studies we describe here is limited field data available to test model predictions or constrain model input. Because of this, advancing model-field comparisons will certainly lead to new insights about deltaic dynamics and inspire new modeling approaches.