Comparison between experimental and numerical stratigraphy emplaced by a prograding delta

A one-dimensional model that is able to store the stratigraphy emplaced by a prograding delta is validated against experimental results. The laboratory experiment describes the migration of a Gilbert delta on a sloping basement into standing water, i.e., a condition in which the stratigraphy emplaced by the delta front is entirely stored in the deposit. The migration of the delta front and the deposition on the delta top are modeled with total and grain-size-based mass conservation models. The vertical sorting on the delta front is modeled with a lee-face-sorting model as a function of the grain size distribution of the sediment deposited at the brinkpoint, i.e., at the downstream end of the delta top. Notwithstanding the errors associated with the grain-size-specific bedload transport formulation, the comparison between numerical and experimental results shows that the model is able to reasonably describe the progradation of the delta front, the frictional resistances on the delta top, and the overall grain size distribution of the delta top and delta front deposits. Further validation of the model in the case of variable base level is currently in progress to allow for future studies, at field and laboratory scale, on how the delta stratigraphy is affected by different changes of relative base level.


Introduction
A fluvio-deltaic deposit can be considered composed of two parts: the delta top, where sediment is transported, eroded, and deposited by fluvial-type processes, and the delta front, where sedimentary processes are characterized by avalanches, deposition of sediment from suspension, and particle entrainment and deposition by submarine currents (Swenson et al., 2000).Depending on the dominant submarine processes, delta front deposits can either be characterized by upward-fining, upward-coarsening, or more complex vertical depositional patterns (e.g., Rohais et al., 2008, for the Gilbert deltas in the Gulf of Corinth, Greece).This paper focuses on relatively coarse grained Gilbert deltas, with steep fronts and a thick delta front deposit compared to the delta top deposit (Edmonds et al., 2011).Angles of the delta front of about 20-35 • are observed in the Serra Ciciniello section of the Potenza Basin, Italy (Longhitano, 2008).
In this study we present the validation against experimental observations of a one-dimensional delta migration model that is able to predict and record the spatial, i.e., vertical and streamwise, variation of the grain size distribution within the deltaic deposit, under the assumption that grain flows are the dominant delta front sedimentation process.In particular, the model describes (i) grain flow deposition on the delta front and the emplacement of upward-fining units, as well as (ii) fluvial deposition on the delta top.It can thus be of aid in the interpretation of ancient and modern deltaic deposits in which grain flows are the dominant delta front deposition process.
The work of Kleinhans (2005), Blom et al. (2003Blom et al. ( , 2013)), and Blom and Kleinhans (2006) on grain flow deposits emplaced on the lee face of dunes, bars, and on Gilbert delta fronts demonstrates that the vertical sorting pattern of the grain flow deposits does not significantly change from Published by Copernicus Publications on behalf of the European Geosciences Union.
migrating bedforms to prograding Gilbert deltas.It can indeed be modeled with the same mathematical relations as a function of the height of the lee face and of the sediment transport on the topmost part of the lee face, the brinkpoint.Thus, the procedure for the storage of stratigraphy, validated herein for the case of a single downstream-migrating lee face (or delta front), may be extended in the future to multiple migrating lee faces, such as those observed in the cases of stacked deltaic complexes and downstream-migrating dunes or bars.
The upward-fining pattern of grain flow deposits is the result of several processes that can be summarized as follows (Fig. 1).Bedload sediment is deposited in a wedge on the topmost part of the slip face until the static angle of repose of the material is reached.This mechanism is termed grain fall, and it is characterized by preferential deposition of coarser grains in the upstream part of the grain fall deposit.When the static angle of repose is exceeded, the wedge collapses, a grain flow is initiated, and the remobilized sediment avalanches down the lee face.During the grain flow, sediment sorting takes place, and coarse sediment is deposited over the lower part of the lee face, while fine sediment remains trapped in its upper portion.The formation and emplacement of grain-fall-grain-flow deposits is schematically represented in Fig. 1a and b, respectively.In Fig. 1 the delta is assumed to migrate over a horizontal substrate for illustration purposes only.
A numerical model that needs to reproduce the stratigraphy emplaced by a prograding delta is composed of three sub-models (types 1-3) that respectively describe (1) the total sediment mass conservation in the system (e.g., Wright and Parker, 2005a, b), (2) the mass conservation of sediment in each grain size range (e.g., Hirano, 1971), and (3) the sorting process on the lee face (e.g., Blom and Parker, 2004;Blom et al., 2006).Each sub-model has its purpose.In particular, total sediment mass conservation models, type 1, predict the rates of channel bed aggradation and delta front migration.Type 2 models account for the mobility of sediment particles of different sizes on the delta top.Finally, lee-facesorting models, type 3, synthetically describe the grain-fallgrain-flow mechanism that occurs on the lee face.Bookkeeping procedures to store stratigraphy are implemented to record the spatial variation of the characteristics of the deposited material.
The comparison between experimental measurements and one-dimensional numerical predictions of stratigraphy, defined as the vertical and streamwise variation of grain size distribution within the deposited sediment (Viparelli et al., 2010a), is performed for the case of an experimental Gilbert delta prograding into standing water (Ferrer-Boix et al., 2013).The experiment is characterized by the following conditions: (i) the system is always net depositional; (ii) the stratigraphy emplaced by the delta front is entirely stored in the deeper portion of the deposit; (iii) the stratigraphy emplaced on the delta top is stored in the upper portion of the deposit; (iv) the mode of sediment transport on the delta top is lower plane bed bedload regime, i.e., migrating dunes and bars that would partially rework the upper portion of the delta top are not present; and (v) grain flows are the dominant deposition process on the delta front.
The mass conservation sub-model (type 1) is validated by comparing measured and predicted longitudinal profiles of delta elevation and delta front migration rates.Due to the different depositional processes on the delta front and on the delta top (Fig. 1), these experimental results allow for the validation of both the grain-size-specific mass conservation submodel (type 2) and the lee-face-sorting sub-model (type 3).In particular, the type 2 sub-model is validated with the comparison between measured and predicted grain size distributions of the topmost layer of the delta top deposit at the final state of the experiment, while the lee-face-sorting sub-model (type 3) is validated by comparing the grain size distribution of the front deposits, i.e., the deposit below the delta top.
The paper is organized as follows: the relevant characteristics of the laboratory experiment and the numerical model are presented in Sects. 2 and 3, respectively.The comparison between measured and numerical results is discussed in Sect. 4. The results of the study and the plans for future work are summarized in the last section of the manuscript.A brief discussion of the model sensitivity to the grain-size-specific bedload model is presented in the Appendix.

The laboratory experiment
The laboratory experiment was performed in the 12 m long and 0.60 m wide tilting flume at the Hydrosystems Laboratory, University of Illinois, Urbana-Champaign (Ferrer-Boix et al., 2013).The parent material was a mixture of sand and pea gravel with a geometric mean diameter, D g , of 3.43 mm, and geometric standard deviation of 1.75, shown in Fig. 2, where the blue line represents the cumulative grain size distribution.The yellow diamonds denote the fractions of sediment finer than the bound diameters, D bi , used in the numerical runs presented in Sect. 4. The red line connects yellow squares denoting the fractions of sediment contained in each characteristic grain size range, i.e., the grain size range bounded by two consecutive bound diameters, D bi and D bi+1 .The sediment in each characteristic grain size range is modeled as uniform and with characteristic grain size D i equal to the geometric mean of the bound diameters (e.g., Parker, 2004).
The laboratory flume operated in sediment feed mode; that is, water and sediment were fed at a constant rate at the upstream end of the flume.In particular, the flow rate of 47 L s −1 was controlled with an electromagnetic flow meter, and the feed rate of parent material was set at 800 g min −1 with a screw-type feeder.
The flume was tilted with a bottom slope of 2 %.An initial 10 cm thick layer of parent material was placed on the bottom of the flume for the entire length of the experimental section.The downstream water elevation was set at 26 cm above the initial deposit by means of a transverse wall located 9 m downstream of the flume entrance.The experiment started when the flow and the sediment feeder were simultaneously turned on.A downstream-migrating Gilbert delta formed and migrated downstream.The experiment terminated when the Gilbert delta reached the transverse wall, i.e., after 10.42 h.
Longitudinal profiles of delta elevation were periodically recorded (i.e., after 0.15, 0.67, 1.05, 2.30, 3.97, 5.98 ,and 8.50 h from the beginning of the experiment) with four ultrasonic transducer probes (Wong et al., 2007).At the end of the experimental run, a last longitudinal profile was measured and core samples were collected in six locations of the deposit (i.e., 3.5, 4.5, 5.5, 6.5, 7.5, and 8.5 m downstream of the flume entrance) with the metallic box described by Blom et al. (2003).Each core sample was then sliced into 2 cm thick layers.Each layer was oven dried and its grain size distribution was measured to characterize the emplaced x is a streamwise coordinate, x b and x t are the coordinates of the brinkpoint and of the delta toe.H is the water depth on the delta top.q is the input water discharge per unit channel width, q bTfeed is the total (i.e., summed over all the grain sizes) volumetric sediment feed rate per unit channel width, p ifeed characterizes the grain size distribution of the fed material, and ξ b is the elevation of the downstream standing water above the datum.η represents the elevation of the delta top above the datum, η b is the elevation of the brinkpoint, and L a denotes the active layer thickness.c b denotes the migration rate of the brinkpoint, and c t is the migration rate of the delta toe.
stratigraphy.Further details on the experiment are reported in Ferrer-Boix et al. (2013).

The numerical model
The numerical model to predict the stratigraphy emplaced by a downstream-migrating Gilbert delta is built by means of coupling (i) a delta progradation model, type 1 (Wright and Parker, 2005a, b); (ii) an active layer model for mass conservation of nonuniform sediment on the delta top, type 2 (Hirano, 1971;Parker, 1991a, b); (iii) a lee-face-sorting model for the delta front, type 3 (Blom et al., 2013); and (iv) a procedure for the storage of stratigraphy (Viparelli et al., 2010a).The role of each model is schematically represented in Fig. 3 with the definition of the relevant model parameters and boundary conditions.
As discussed by Blom (2008) different modeling approaches can be used to couple grain-size-specific mass conservation models (type 2) and lee-face-sorting models (type 3).The active layer approximation (Hirano, 1971, as modified by Parker, 1991a) is used herein because (i) it can be implemented with reasonably large spatial and temporal steps, allowing for future laboratory-and field-scale applications (Blom, 2008), and (ii) it reasonably reproduces the stratigraphy emplaced under lower plane bed bedload transport conditions on the delta top (Viparelli et al., 2010a, b).
To simplify the schematic representation of the system, in Fig. 3 the delta progrades on a horizontal basement.The model, however, is designed to handle an arbitrarily sloping basement with constant slope S b .
The model parameters represented in Fig. 3 are the streamwise coordinate x and the streamwise locations of the brinkpoint and the delta toe, x b and x t , respectively.As further discussed in the remainder of this section, the model boundary conditions are expressed in terms of volumetric feed rate of water, q, and sediment, q bT , and elevation of the downstream standing water above the datum, ξ b .The grain size distribution of the sediment feed is characterized in terms of fraction of sediment in the generic grain size range i, p ifeed .The water depth on the delta top is denoted by H , the elevation of the deltaic deposit above the datum with η, and the thickness of the active layer with L a .The elevation of the brinkpoint above the datum is η b , and its migration rate is denoted by c b .Similarly, the migration rate of the delta toe is denoted by c t .
The numerical model is a one-dimensional (laterally averaged) model of delta growth based on the standard shallow water equations of open channel flow and on the equation of sediment conservation (e.g., Parker, 2004).Before outlining the governing equations and the numerical scheme for the storage of grain size stratigraphy, the simplifying assumptions are listed below.Some of these assumptions are introduced to apply the model at laboratory scale and can be relatively easily relaxed for field scale applications: 1. the volume bedload transport rate is orders of magnitude smaller than the flow discharge, so that the quasi-steady approximation (De Vries, 1965) holds and the bed elevation profile can be considered as unchanging in the hydraulic calculations; 2. the channel cross section is rectangular, with constant width B and vertical smooth sidewalls; 3. the flow is Froude-subcritical and the shallow water equations are reduced to the equation for a backwater curve, so that the equations can be integrated upstream starting from the brinkpoint, i.e., the downstream end of the delta top; 4. the laboratory flume is long enough that entrance effects can be reasonably neglected; 5. grain flows are the dominant depositional process on the delta front.

Calculation of the flow
As in Viparelli et al. (2010a), the shallow water momentum equation is modified to account for the different shear stresses acting on the smooth flume sidewalls and on the rough delta top.This correction is necessary to properly model bedload transport in a laboratory flume (Vanoni, 1975).The water mass and momentum conservation equations for open channel flow take the form respectively, where H denotes the water depth, U is the cross-sectionally averaged flow velocity, g is the acceleration of gravity, S denotes the slope of the delta top, ρ is the water density, τ e is the effective shear stress associated with the resistances of the smooth sidewalls and of the rough bed, and t and x are a temporal and a streamwise coordinate, respectively.
The simplified version of the Vanoni-Brooks decomposition (e.g., Francalanci et al., 2008) of the shear stress into a sidewall and a bed component to estimate the effective shear stress, τ e , implemented in the previous versions of the model is substituted with the complete Vanoni and Brooks (Vanoni, 1975) decomposition, as modified by Chiew and Parker (1994).These formulations are based on the assumption that the cross section can be decomposed into two noninteracting regions, the bed region and the sidewall region, where the mean flow velocity and the energy gradient are respectively equal to the mean flow velocity, U , and the energy gradient, S f , of the entire cross section.The continuity equation thus takes the form where A is the cross-sectional area, and A w and A b respectively denote the area of the wall and of the bed region.
The main difference between the Chiew-Parker and the simplified Vanoni-Brooks decomposition (e.g., Francalanci et al., 2008) is related to the partition of the cross section between the smooth sidewall region and the rough bed region.The underlying assumption of the simplified Vanoni and Brooks formulation is that the boundary between the sidewall region and the bed region is a 45 • straight line.In the complete formulation (Vanoni, 1975;Chiew and Parker, 1994) the areas of the sidewall and of the bed regions are computed from the flow characteristics, with a better estimate of the shear stress on the rough boundary (Chiew and Parker, 1994).
The Chiew-Parker decomposition is based on the following form of the momentum balance for the cross section: where τ b and τ w respectively denote the shear stresses on the rough bed and on the smooth sidewalls, and P , P b , and P w are the wetted perimeters of the entire cross section (B + 2 H ), the bed (B), and the sidewall region (2 H ).
In the case of turbulent open channel flow, the shear stresses can be expressed as the product of the water density, the mean velocity squared, and a nondimensional friction coefficient C f , τ = ρC f U 2 .Since the mean flow velocity is assumed to be the same in the bed region, in the sidewall region, and in the entire cross section, Eq. ( 4) can be rewritten as where C fe is an effective nondimensional friction coefficient associated with the resistances on the sidewalls and on the bed, and C fb and C fw denote the nondimensional friction coefficients for the bed and the sidewall region, respectively.Under the assumption that the Darcy-Weisbach relation can be applied to the entire cross section, to the bed, and to the sidewall region, the energy gradient is given as where r, r b , and r w denote the hydraulic radii (i.e., the ratios between the cross-sectional areas and the wetter perimeters) for the entire cross section, for the bed, and for the sidewall region, respectively.Recalling that the Reynolds number of the cross section is defined as Re = rU/ν, with ν denoting the kinematic viscosity of the fluid, Eq. ( 6) can be rewritten as where Re b and Re w are the Reynolds numbers of the bed and the sidewall region, respectively.The unknowns in Eqs.
(3), (5), and ( 7) are the friction coefficients, C fe , C fb , and C fw ; the area of the bed region, A b ; and the area of the wall region, A w .Two closure relations are thus needed to solve the problem.
The first closure relation expresses C fb as a function of the hydraulic radius of the bed region, r b , and of the roughness height of the delta top, k s , as Figure 3 in Viparelli et al. (2010a) shows that this bed resistance model is appropriate to describe flow resistances in the bed region with Ferrer-Boix et al. ( 2013) sediment mixture and flow conditions if (i) the roughness height is assumed to be equal to 1.5 D s90 and (ii) the active layer thickness is assumed to be equal to D s90 .Here D s90 denotes the diameter of the sediment stored in the active layer such that 90 % of the active layer sediment is finer.
The second closure is the relation for hydraulically smooth walls given by Vanoni (1975) to compute the Darcy-Weisbach sidewall friction coefficient f w = 8 C fw as a function of the Reynolds number of the wall region Equations ( 1) and ( 2) are reduced to the classical backwater form using (i) the quasi-steady approximation (De Vries, 1965) to drop the time dependence as well as (ii) the definition of effective shear stress as a product of water density, mean flow velocity square, and effective friction coefficient, τ e = ρ C fe U 2 .The backwater equation thus takes the form where Fr denotes the Froude number defined as U/(gH ) 0.5 .In the numerical run described below, Eq. ( 10) is integrated in the upstream direction with the downstream boundary condition ξ = ξ b = 0.26 m, with ξ denoting the water surface elevation above the datum and the subscript b indicating the downstream end of the delta top, i.e., the brinkpoint (see Fig. 3).

Calculation of sediment transport and deposition on the delta top
Bedload sediment transport on the delta top is modeled with the version of the Ashida-Michiue bedload relation of Viparelli et al. (2010b).This grain-size-specific bedload relation is derived for mobile bed equilibrium conditions obtained in the same laboratory flume and with the same sediment mixture of Ferrer-Boix et al. (2013).
During the Viparelli et al. (2010b) experiment the flume was operated in sediment-recirculating mode, meaning that the sediment collected in the sediment trap was recirculated to the upstream end of the flume.Thus, during the condition of nonequilibrium, the sediment input rate and its grain size distribution were not constant in time (Viparelli et al., 2010a, b).
In addition, in a sediment-recirculating flume, the total volume of sediment in the system does not change in time, and so the grain size stratigraphy of the bed deposit and the equilibrium conditions are dependent on the initial experimental conditions (Parker and Wilcock, 1993).Ferrer-Boix et al. (2013) operated the laboratory flume in sediment feed mode, i.e., with constant grain-size-specific sediment input rate, and with a volume of sediment in the flume that increased in time.
In a sediment feed flume the conditions of mobile bed equilibrium are independent of the initial condition of the experiment and are dictated by the upstream input of water and sediment only (Parker and Wilcock, 1993).It is thus reasonable to expect that disequilibrium bedload transport conditions in a sediment feed flume, such as those of the Ferrer-Boix et al. ( 2013) experiment, are somewhat different from those observed in a sediment-recirculating flume (Viparelli et al., 2010a).
As shown in Fig. 2, the parent material is divided in M (M = 9 for the numerical run presented herein) grain size ranges with characteristic diameters D i .The volumetric bedload transport rate per unit channel width, q bT , is equal to the Grain-size-specific nondimensional volumetric bedload transport rates per unit width, q * bi , (Einstein parameters) are defined as (Parker, 2008) where F i represents the fraction of active layer sediment in the generic (ith) grain size range, and R denotes the submerged specific gravity of the parent material, i.e., (ρ s − ρ)/ρ, with ρ s denoting the density of the sediment.R = 1.58 for the experiment discussed herein.
The grain-size-specific Einstein parameters are computed as where β is an adjustment coefficient equal to 0.27 for the considered experimental conditions, τ * bi is the grainsize-specific nondimensional shear stress on the bed region (Shields number), and τ * ci is its reference value for significant bedload transport of sediment in the generic (ith) grain size range.
The grain-size-specific Shields number is defined as (Parker, 2008) its reference value for significant bedload transport is estimated with the hiding/exposure function derived by Viparelli et al. (2010b) that is valid for the Ferrer-Boix et al. ( 2013) experimental conditions where D sg is the geometric mean diameter of the active layer sediment, and τ * scg represents the reference Shields number for significant motion in the case of uniform sediment.τ * scg is equal to 0.043.
The Exner equation of conservation of total (i.e., summed over all the grain sizes) sediment mass conservation takes the form (Parker, 2004) where η denotes the elevation of the delta top above the datum (Fig. 3) and λ p is the bulk bed porosity, equal to 0.35 in the numerical run discussed below (Viparelli et al., 2010a).Equation ( 16) is solved to compute the aggradation rate of the delta top and thus update the longitudinal profile of the Gilbert delta top at each time step.The grain-size-specific equation of conservation of sediment in the generic (ith) grain size range takes the form (e.g., Hirano, 1971;Parker, 2004) where L a denotes the thickness of the active layer, F i is the fraction of sediment in the ith grain size range in the active layer, and f Ii represents the fraction of sediment in the generic grain size range exchanged between the active layer and the deposit during channel bed aggradation or degradation.
In the case of delta top erosion, the grain size distribution of the sediment exchanged between the emplaced deposit and the active layer, f Ii , is equal to the grain size distribution of the deposit.However in the case of an aggrading delta top, the grain size distribution of the sediment transferred to the deposit is assumed to be a weighted average between the grain size distribution of the active layer and of the bedload (Hoey and Ferguson, 1994), where p i denotes the fraction of sediment in the generic grain size range in the bedload, i.e., p i = q bi /q bT .Toro-Escobar et al. (1996) discuss the reason why the parameter α should be greater than 0 and smaller than 1.When α = 0 the grain size distribution of the sediment transferred to the substrate during channel bed aggradation is equal to the grain size distribution of the bedload, and the downstream fining observed in gravel-bed rivers cannot be modeled.However, if α = 1 the surface material is directly transferred to the substrate during channel bed aggradation, and the formation of the coarse pavement observed in gravel-bed rivers that regulates the mobility of particles differing in size (Parker et al., 1982;Parker and Klingeman, 1982) cannot be modeled.
Equation ( 17) is solved to compute the time rate of change of F i and thus update the grain size distribution of the active layer at each time step.

Calculation of sediment transport and deposition on the delta front
The bedload transport rate that reaches the brinkpoint is deposited as grain fall deposit on the upper part of the delta front.Thus, the overall grain size distribution of the grain fall deposit is equal to the grain size distribution of the bedload at the brinkpoint at the specific time.When the static angle of repose of the sediment is exceeded, a grain flow is initiated, and sediment is distributed over the delta front.In particular, coarse sediment is deposited more abundantly in the lowermost part of the front and finer sediment is trapped more abundantly in the upper portion of the lee face.Vertical sorting of sediment on the lee face of the delta front is modeled with the lee face model of Blom et al. (2013).In particular, it is described in terms of a sorting function, ω i , defined as where f si represents the volume fraction content of sediment in the ith grain size range on the slip face at elevation z above the datum, and p i,b represents the volume fraction content of sediment in the generic grain size range in the bedload at the brinkpoint.In the Blom et al. (2013) the sorting function is assumed to linearly vary with the nondimensional elevation z * = (z − η ba )/ , where η ba is the average elevation of the slip face and is the height of the slip face deposit relative to a vertical coordinate, δ i is the lee sorting parameter, defined as with σ qbb denoting the standard deviation on the sedimentological ϕ scale of the bedload at the brinkpoint, and τ * bbsg representing the Shields parameter at the brinkpoint evaluated with the geometric mean diameter of the active layer, D sg .ϕ reli is the adjusted relative arithmetic grain size defined as where ϕ i is the characteristic grain size D i on the ϕ scale, ϕ i = − log 2 D i , and ϕ mtop is the adjusted arithmetic mean grain size of the lee face deposit

Grids for the storage of the stratigraphy
The delta growth problem is characterized by a moving boundary at the downstream end of the delta top, the brinkpoint.Thus, Eqs. ( 10), ( 16), and ( 17) could be integrated in a moving-boundary coordinate system, in which the streamwise coordinate, x, is made nondimensional with the coordinate of the brinkpoint, x b (Swenson et al., 2000).In this moving-boundary system the distance between the computational nodes does not change in time but, due to the movement of the brinkpoint, it varies in the dimensioned coordinate system x (e.g., Wright and Parker, 2005a).
In an active layer model, the moving-boundary transformation requires cumbersome interpolations of the size distributions associated with each computational node.The grain size distributions associated with each node change for (i) fluxes of sediment in the streamwise direction due to the changing dimensioned spatial distance between the computational nodes x, as well as for (ii) vertical fluxes of sediment due to aggradation and degradation of the bed deposit.The streamwise fluxes of sediment in the active layer and in the bed deposit are estimated by interpolating the grain size distributions associated with the computational nodes, with a consequent loss of stratigraphic information.
Since the ultimate scope of the numerical model is to store and access the stratigraphy emplaced by the prograding delta, the governing equations are not solved in the movingboundary coordinate system.A grid with a fixed distance between the computational nodes, x, is used to model sediment transport and deposition upstream of the brinkpoint (Eke et al., 2011;Viparelli et al., 2011a).The distance between the brinkpoint and the last grid node is denoted as x brink .As the brinkpoint moves downstream, x brink increases.When x brink > x, a new grid node is added to the fixed grid, as shown in Fig. 4. As in the previous figures, the delta for Fig. 4 progrades on a horizontal basement; however the formulation holds for an arbitrary sloping basement with slope S b .
The migration rate of the brinkpoint, c b , is computed under the assumptions that (i) all the sediment is trapped on the delta front and (ii) the lee face has a constant slope S l (e.g., Wright and Parker, 2005a), where q bbT denotes the total bedload transport rate at the brinkpoint, η b is the elevation of the brinkpoint above the datum, and x b and x t respectively denote the streamwise coordinates of the brinkpoint and of the delta toe.Equation ( 24) is derived by integrating the Exner equation, Eq. ( 16), on the delta front.As sediment is deposited on the delta front, the delta toe migrates downstream with velocity c t .The migration rate of the delta toe is computed by imposing the continuity of bed elevation, i.e., the elevation of the lowermost point of the delta front must be equal to the elevation of the basement.The continuity condition for the movement of the delta toe takes the form (e.g., Wright and Parker, 2005a) where S b denotes the basement slope.Due to the assumption of the constant slope of the delta front, Eqs. ( 24) and ( 25) are solved to update the streamwise coordinates of the brinkpoint and of the delta toe and thus the longitudinal profile of the delta front.
The bookkeeping procedure of stratigraphy in the delta deposit, i.e., upstream of the brinkpoint point, is that of Viparelli et al. (2010a).The bed is divided in two partsthe relatively thin and well-mixed (i.e., no vertical variation of the grain size distribution) active layer and the substrate, whose grain size distribution can vary in the vertical direction.The grain size distribution of the substrate is stored in the grid represented in Fig. 4 at time t.The substrate deposit is divided into horizontal well-mixed layers.The lowermost grid node, node 1, is located on the datum; the uppermost grid node, node N, is at the active-layer-substrate interface, i.e., at elevation η − L a above the datum.The grain size distribution associated with the grid node j is representative of the layer bounded by the grid nodes j and j − 1.The vertical distance between the consecutive grid nodes from node 1 to node N − 1 is L s , equal to 2 cm in the numerical run presented herein.
The vertical distance between node N − 1 and node N is z < L s .As the delta top aggrades, sediment is stored in the topmost part of the substrate and z increases.When z becomes greater than L s , a new grid node is added to the grid (see Fig. 4 at time t + t).The distance between node N and node N − 1 is equal to L s and the new node N + 1 is added to the grid at the active-layer-substrate interface.The grain size distribution of the material stored in each layer is a weighted average over the thicknesses of the topmost layer of the grid and of the sediment deposited at each time step.The grain size distribution of the sediment transferred to the substrate during channel bed aggradation is f Ii , given by Eq. ( 18).
The stratigraphy of the delta front is stored in a moving grid (Fig. 4).The streamwise distance between the N f grid nodes is equal to x front = (x t − x b )/(N f − 1), which may vary in time due to the different migration rates of the brinkpoint and of the delta toe.Horizontal fluxes of sediment from the front to the fixed grid one time step later and between the nodes of the moving gird are estimated by interpolating the grain size distributions of the sediment stored at the same elevation above the datum.Vertical fluxes of sediment are computed to transfer the newly deposited sediment to the existing front substrate with the same averaging procedure implemented for the delta deposit, as shown in Fig. 4 (Viparelli et al., 2011a).

Results and discussion
The numerical simulation of the Ferrer-Boix et al. ( 2013) experiment is performed with a fixed distance between the computational nodes on the delta top, x, of 0.1 m; a temporal interval, t, of 10 s; and 40 moving grid nodes on the delta front.
It is important to note here that the parameters of the total and grain-size-based sediment conservation models, i.e., α, k s , β, and L a , are those of the Viparelli et al. (2010b) experiments.The parameters of the lee-face-sorting model, i.e., exponents and coefficient of Eq. ( 21), are those of the Blom et al. (2013) analysis of laboratory and field data on vertical sorting on the lee face of dunes and Gilbert deltas.Thus, it was not necessary to tune or calibrate any model parameter to validate the numerical model against experimental results.
The validation of the delta growth sub-model (type 1) is performed by comparing (i) the brinkpoint location, x b , in time and (ii) the longitudinal profile of delta elevation at the end of the experiment.Measured and numerical temporal variations of brinkpoint location are represented in Fig. 5, where the vertical bars denote a ±5 % error.The comparison between longitudinal profiles of bed elevation is reported in Fig. 6, where the elevation of the initial deposit is represented in grey, the initial condition for the numerical run is in blue, and the final longitudinal profile is in red.Error bars in Fig. 6 show that the numerical delta profile at the end of the numerical run approximates the experimental data within a ±1 cm interval.Figures 5 and 6 show that the total (i.e., summed over all the grain sizes) sediment mass conservation model is able to reasonably capture the temporal evolution of the longitudinal profile and thus the total bedload transport rates on the delta top.In addition, Fig. 6 shows that the bed resistance model reasonably reproduces the experimental conditions, since the slopes of the numerical and the experimental delta top are reasonably similar.
The validation of the grain-size-specific mass conservation model for the delta top (type 2) is presented in Fig. 7, where the error bars denote a ±5 % interval around the measured points.Due to the lack of experimental data on the grain size distribution of the active layer, the comparison is done in terms of measured grain size distributions of the topmost 2 cm of the experimental delta top (diamonds in Fig. 7) and the average grain size distribution of the topmost portion of the numerical deposit, i.e., the active layer and the two uppermost layers of the grid for the storage of grain size stratigraphy (red line in Fig. 7).The numerical results are averaged over a volume thicker than the experimental samples to have a relative robust estimate of the grain size distribution in a well-mixed layer characterizing the grain size distribution of the active layer and of the topmost portion of the substrate.
The comparison in Fig. 7 shows an overall reasonable agreement between measurements and numerical predictions, which is rarely obtained with a 1-D active layer model.In the sampling sections at 3.5, 6.5 ,and 7.5 m from the entrance of the flume the model underestimates the fraction of sediment in the 1.53 mm size range.This is balanced in the 3.5 m section by a slight overestimation of the sediment in the 5.02 and 7.74 mm ranges, and by a more severe overestimation of the sediment in the 2.83 mm range in the sections 6.5 and 7.5 m downstream of the flume entrance.We suspect that the differences between the numerical results and the experimental data are related to the grain-size-specific sediment transport model, i.e., Eqs. ( 13) and ( 15).
As mentioned above, the Viparelli et al. (2010b) model is based on sediment-recirculating flume experiments.In this experimental setting, mobile bed equilibrium is reached through a rotation of the longitudinal profile around the center of the flume.In other words, only the topmost part of the deposit is reworked and the grain size distribution of the transported sediment is constrained by the grain size distribution of the mobilized sediment (Viparelli et al., 2010a).The flume in the Ferrer-Boix et al. (2013) experiment is operated in sediment feed mode, i.e. with a constant input rate of parent material.It is thus reasonable to expect that the sediment mobility in the Ferrer-Boix et al. (2013) experiment is slightly different than in the Viparelli et al. (2010b) experiments.Unfortunately no experimental data are available to further validate the bedload transport model.
The numerical stratigraphy of the bed deposit is represented in Fig. 8 (with x = 5 cm and L s = 1 cm for illustration purposes only), where the dots represent the grid nodes for the storage of stratigraphy and the color scale represents the geometric mean diameter of the substrate layer.The blue oval indicates the portion of the delta deposit whose stratigraphy is affected by the model initial condition, i.e., a wellmixed deposit of parent material with the longitudinal profile represented in Fig. 6.The stratigraphy of this initial profile does not significantly change in time during the numerical runs because the delta front migrates downstream and only a thin layer of sediment is deposited on top of the initial deposit.The black line in Fig. 8 represents the elevation of the  initial layer of parent material placed on the bottom of the flume.The two lines of red dots in the upper part of the delta top denote the active layer thickness.
The color scheme of Fig. 8 shows that the model is able to reproduce the upward-fining profile emplaced by the downstream-migrating lee face.A closer look at the figure reveals that the delta top deposit has a geometric mean diameter similar to the parent material and finer than the active layer, as observed in gravel bed rivers (e.g., Viparelli et al., 2011b).Further, as observed by Ferrer-Boix et al. (2013) during the experimental work, the sediment stored in the lowermost part of the front deposit appears to become coarser in the downstream direction.
Figure 9 shows the temporal variation of the geometric mean diameter of the bedload at the brinkpoint, D gbb .After an initial adjustment due to the development of the mobile armor on the initial deposit, D gbb and thus the grain size distribution of the bedload at the brinkpoint remain reasonably constant in time.Thus, the observed downstream coarsening cannot be the result of an increasingly coarser bedload transport rate at the brinkpoint in time.Instead, we interpret the apparent downstream coarsening as the result of the increasing delta front elevation.As the Gilbert delta progrades on the steep basement, the front becomes higher, and there is more space to sort the bedload material that reaches the brinkpoint, as modeled with Eq. ( 20).Numerical experiments are currently in progress to investigate whether a similar downstream coarsening can be driven by relative baselevel rise.The blue ovals in Figs. 8 and 9 identify the area in which the bedload transport rate at the brinkpoint is affected by the initial model condition, i.e., the development of a coarse active layer on the unarmored initial delta of parent material.
The lee-face-sorting model (type 3) is validated by comparing experimental and numerical grain size distributions of the delta front deposit in the cross sections at 4.5, 5.5, 6.5, 7.5 ,and 8.5 m from the entrance of the flume.Data collected in the measuring section at 3.5 m are not used in the comparison because, as shown in Figs. 8 and 9, the numerical results are affected by the initial condition.
The comparison between experimental and numerical data is represented in Fig. 10 by vertical profiles of sediment fractions in the characteristics grain size ranges 1.53, 2.83, 5.02, and 7.74 mm.The diamonds represent the experimental data, the red lines are the model results, and the horizontal error bars denote a ±5 % error.The vertical elevation of the diamonds, z, corresponds to the elevation of the center of each sampling layer above the datum, and the vertical error bars identify the thickness of the sampling layer, i.e., ±1 cm.
The comparison in Fig. 10 shows that -notwithstanding the uncertainties related to the grain-size-specific bedload transport relation, and thus the grain size distribution of the bedload passing the brinkpoint -the model is able to reasonably capture the overall grain size distribution of the delta front deposit.Significant differences between the fractions of fine sediment, i.e., 1.52 and 2.83 mm, stored in the deposit in the 6.5 and 7.5 m positions confirm what we have previously observed for the grain size distribution of the delta top deposit, i.e., that the bedload transport model may not always be able to properly reproduce the transport of the finer components of the sediment mixture.
The bedload transport model, i.e., the predicted grain size distribution of the sediment at the brinkpoint, is certainly one of the major sources of error in the prediction of the grain size distribution of the delta front deposit.An additional source of error can be hidden in the lee-face-sorting model.As the delta front migrates on a 2 % sloped basement, the delta front height increases (see Fig. 6).The Blom et al. (2013) lee-facesorting model is a linear model, Eq. ( 20), in the nondimensional elevation z * ,; thus nonlinear effects due to an increasing front height are not explicitly accounted for.The study of vertical sorting on an increasingly high lee face goes well beyond the scope of this paper, but we suspect that it may partially explain the differences between the numerical and the experimental stratigraphy of the considered Gilbert delta.

Conclusion and future work
The comparison between numerical and experimental delta stratigraphy was conducted for the case of a Gilbert delta prograding on a sloping basement into standing water.These experimental conditions are appropriate for the validation of a model of delta morphodynamics because the stratigraphy emplaced by the migrating delta front (i.e., the lee face) is entirely stored within the deposit.In other words, a train of migrating bedforms, such as bars or dunes, does not rework the lee face deposit.
The comparison was done in three steps.First, the flow and the total (i.e., summed over all the grain sizes) sediment conservation models were validated against profiles of channel bed elevation and migration rates of the brinkpoint.This comparison shows the following results: 1. Numerical predictions of the streamwise coordinate of the brinkpoint were in reasonable agreement with the experimental measurements (Fig. 5).Since the migration rate of the brinkpoint was computed with an integral shock condition of the equation of total sediment conservation, the model was able to reasonably predict the total bedload transport rates at the brinkpoint.
The results of the grain-size-specific sediment conservation model on the delta top were validated against the grain size distributions of the topmost part of the delta top deposit.The results of Fig. 7 show that the model was able to reasonably capture the overall grain size distribution of the delta top, but it tended to underestimate the fractions of fine material deposited in the topmost part of the experimental Gilbert delta.This was probably due to a failure of the bedload transport model, based on sediment-recirculating flume experiments and applied to simulate a sediment feed flume experiment.Laboratory measurements on the grain size distribution of the active layer were unfortunately not available to further validate the bedload transport model.Finally, the numerical stratigraphy emplaced by the delta front was compared with the laboratory data in Fig. 10.The model reasonably reproduced the upward fining observed in the laboratory, but due to the errors in the predictions of the grain size distributions at the brinkpoint, the differences between the numerical predictions and the experimental measurements were sometimes larger than 5 %.
The results presented herein represent the first step in the validation of the numerical model.Further validation against laboratory experiments is currently in progress to study the stratigraphy of a Gilbert delta under different scenarios of base level change.In the near future we plan to modify the code to not only store but also to access the stratigraphy emplaced by prograding deltas and model (1) Gilbert delta progradation at the field scale and (2) the formation of stacked Gilbert delta complexes at the laboratory and field scale. in the emplacement of finer delta front deposits in the upstream part of the system compared to the case of a single prograding Gilbert delta.More than a sensitivity analysis on the bedload relation, this exercise becomes the comparison between the stratigraphy emplaced by a single Gilbert delta and by stacked deltas, which is an interesting problem that we hope to study in the relatively near future.

Figure 1 .
Figure 1.Schematic representation of the Gilbert delta stratigraphy.(a) Grain flow deposit, (b) grain fall deposit.

Figure 2 .
Figure 2. Grain size distribution of the parent material.D denotes the grain diameter in millimeters.The blue line is the cumulative distribution, and the yellow diamonds denote the bound diameters used in the numerical calculations.The yellow squares indicate the fractions of parent material contained in each characteristic grain size range, i.e., between two bound diameters.The vertical black line is the geometric mean diameter, D g , of the parent material.

Figure 3 .
Figure 3.The numerical sub-models and the relevant model parameters.x is a streamwise coordinate, x b and x t are the coordinates of the brinkpoint and of the delta toe.H is the water depth on the delta top.q is the input water discharge per unit channel width, q bTfeed is the total (i.e., summed over all the grain sizes) volumetric sediment feed rate per unit channel width, p ifeed characterizes the grain size distribution of the fed material, and ξ b is the elevation of the downstream standing water above the datum.η represents the elevation of the delta top above the datum, η b is the elevation of the brinkpoint, and L a denotes the active layer thickness.c b denotes the migration rate of the brinkpoint, and c t is the migration rate of the delta toe.

Figure 4 .
Figure 4. Model grids.GSD stands for grain size distribution.

Figure 6 .
Figure 6.Comparison between measured (diamonds and triangles) and numerical (lines) longitudinal profiles.The profile of the bottom deposit (grey) is a model boundary condition.The delta profile at t = 0.15 h (blue) is the model initial condition.The delta profile at t = 10.42 h (red) is a model result.Error bars of the measured profile at t = 10.42 h denote ±0.01 m.

Figure 7 .
Figure 7.Comparison between predicted and measured grain size distributions of the topmost part of the delta top.The diamonds represent the measurements and the red line is the numerical result.Error bars denote a ±5 % error.

Figure 8 .
Figure 8. Numerical stratigraphy of the deposit.The dots represent the grid nodes, and the color scale is associated with the geometric mean diameter in millimeters of the substrate layers.The black line represents the top of the initial layer of parent material.The blue oval indicates the stratigraphy affected by the initial conditions.

Figure 9 .
Figure 9. Temporal variation of the geometric mean diameter of the bedload at the brinkpoint.The blue oval indicates the initial numerical adjustment of the model, mostly related to the development of a coarse mobile active layer.

Figure 10 .
Figure 10.Comparison between measured and numerical grain size distribution of the front deposit.The diamonds are the experimental data, and the lines are the numerical predictions.Horizontal error bars denote a ±5 % error.The vertical elevation above the datum, z, of the diamonds corresponds to the elevation of the center of the sample.The vertical error bars at ±1 cm denote the thickness of the sampled layer.

Figure A2 .
Figure A2.Numerical stratigraphy of the deposit for the Ashida and Michiue simulation.The dots represent the grid nodes, and the color scale associated with the geometric mean diameter in millimeters of the substrate layers is the same as in Fig. 8.The black line represents the top of the initial layer of parent material.The orange ovals indicate the stratigraphic record emplaced when the stacked small deltas reach the brinkpoint of the initial deltaic deposit.The geometric mean diameter of the parent material is 3.43 mm.The results are obtained using the Ashida and Michiue grain-size-based bedload relation, which is not appropriate for modeling the Ferrer-Boix et al. (2013) experiment.