Measuring Subaqueous Progradation of the Wax Lake Delta with a Model of Flow Direction Divergence

Remotely sensed flow patterns can reveal the location of the subaqueous distal tip of a distributary channel on a prograding river delta. Morphodynamic feedbacks produce distributary channel tips that become shallower over their final reaches before becoming deeper over the unchannelized foreset. The flow direction field over this morphology tends to diverge 10 and then converge providing a diagnostic signature that can be captured in flow or remote sensing data. Twenty-one measurements from the Wax Lake Delta (WLD) in coastal Louisiana, and 317 measurements from numerically simulated deltas show that the transition from divergence to convergence occurs in a distribution that is centered just downstream of the channel tip, on average 132 m in the case of the WLD. With these data we validate the Flow Direction to Channel tips (FD2C) inverse model for remotely estimating subaqueous channel tip location. We apply this model to 33 remotely sensed images of 15 the WLD between its initiation in 1974 and 2016. We find that the distributaries grew unevenly, 6 of the primary channels grew at rates of 60-80 m/yr while one grew at 116 m/yr. We also estimate the growth rate of the total area enclosed by the subaqueous delta platform to be 1.83 km/yr with no obvious rate changes over time.


Introduction
River deltas host productive ecosystems and hundreds of millions of people worldwide.Over the past century, river deltas have changed rapidly, putting these large human populations at risk (Barras et al., 2008;Erban et al., 2014;Wilson et al., 2017;Wu et al., 2017).Monitoring morphologic change on river deltas is key to their sustainable management (Peyronnin et al., 2017).Existing remote sensing techniques provide synoptic monitoring of deltas, but are generally limited to monitoring subaerial or very shallow regions (Couvillion et al., 2011;Li and Damen, 2010;Rahman et al., 2011;Rangoonwala et al., 2016).However, most deltas are far larger than their subaerial portions.For instance, the subaerial area of the Wax Lake Delta in 2015 was 50 km 2 (Olliver and Edmonds, 2017), while the subaqueous area was an additional 82 km 2 (Shaw et al., 2016a).The difference arises because the Wax Lake Delta has an extensive delta front that lies below low tide.While the subaerial and shallow land areas are where marshes are established (Johnson et al., 1985), the subaqueous delta forms the platform upon which subaerial islands grow (Cahoon et al., 2011;Shaw et al., 2018).Hence, the subaqueous platform extent is important as a leading indicator of future marsh growth, necessary data for navigation, and the key area metric for estimating delta volume and volume change (Geleynse et al., 2015).Unfortunately, only a small fraction of global river deltas has been directly surveyed in a manner that resolves their subaqueous portion.This is partly due to the vast area of deltas and partly because year-round turbidity fundamentally limits bathymetric lidar or multispectral remote sensing techniques (Gao, 2009).Many shallow regions along the US coast far from navigation corridors have not been officially surveyed since the 1930s (e.g., NOAA, 2017).Here we make progress in subaqueous delta monitoring by recognizing key connections between delta front bathymetry and the flow field that organizes over it.We then exploit this connection by remotely sensing the flow direction using streak lines on the water surface visible on some deltas (Fig. 1).We use the Wax Lake Delta as a field site due to available bathymetric maps and because it frequently exhibits streak lines that resolve delta front flow directions.In Sect.2, we review the coupling of emergent delta front bathymetry and flow patterns.In Sect.3, we present the flow direction to channel (FD2C) model of estimating the location of channel tips using the remotely sensed flow direction field.In Sect.4, the model is validated on the Wax Lake Delta and with four numerical models of deltas.The model is applied to 33 images of the Wax Lake Delta spanning its development from 1974 to 2016 in Sect. 5 in order to estimate progradation rates of individual channels and the growth rate of total delta area.The strengths and limitations of the model and results from its application are discussed in Sect.6.

Bathymetry and flow patterns on river-dominated deltas
There is virtually no limit to the paths that a parcel of water can trace across a domain with arbitrary bathymetry and boundary conditions.This seemingly unlimited degree of freedom limits the skill of inverse models predicting bathymetry that place no constraints on the possible morphology (Alpers et al., 2004;Romeiser and Alpers, 1997).
However, direct study of river deltas reveals that emergent patterns can be found within their bathymetry and flow patterns.If the bathymetry and flow take on predictable patterns, this greatly reduces the degrees of freedom in a system, thus improving predictability.The postulation of emergent flow patterns has a long tradition in coastal geomorphology (e.g., Edmonds and Slingerland, 2007;Wright, 1977).Extensive work has been done to predict initial flow and sedimentation patterns associated with turbulent jets entering basins with simple initial bed morphology (Fagherazzi et al., 2015).Our work seeks to extend this approach to systems with complex emergent topography and multiple interacting channels.If a certain bed morphology produces distinct flow patterns visible on remotely sensed imagery, then that pattern can be used to predict the underlying morphology.We make this case for the flow patterns on the delta front of a prograding delta.
For well-developed, prograding river deltas, the bed morphology of a channel terminus can be idealized as an adverse bed slope along the thalweg (shallowing with distance downstream) and a basinward slope along the levee (deepening with distance).Together, these produce a gradual loss of channel confinement (Fig. 2).Channels lose definition at the channel tip where the thalweg elevation equals the levee elevation.This transition occurs gradually over more than seven channel widths for the Wax Lake Delta (Shaw and Mohrig, 2014).Beyond the location where channel definition is lost, the unchannelized delta front grows gradually deeper with distance away from the channel.Although dimensions vary, this general morphology has been observed on the Wax Lake Delta (Shaw et al., 2016b;Shaw and Mohrig, 2014), Brant's Pass crevasse on the bird-foot delta of the Mississippi River (Esposito et al., 2013), the Mobile and Apalachicola River deltas (Edmonds et al., 2011b), and the St. Clair River Delta (Fig. 1b; NOAA, 2017).Additionally, numerical models often produce this morphology (Caldwell and Edmonds, 2014;Geleynse et al., 2010;Liang et al., 2016).These deltas can be qualitatively classified as river dominated (Galloway, 1975), both by their large fluvial sources and relatively small waves and tides.Waves and tides can alter this morphology (Leonardi et al., 2013;Nardin and Fagherazzi, 2012).Hence, we limit ourselves to river-dominated conditions here.
Recently, flow patterns have been measured across channel tips on the Wax Lake Delta.Various techniques have been used to show that on Gadwall Pass on the Wax Lake Delta (Figs. 1a, 3), roughly 50 % of water discharge leaves channels laterally (Hiatt and Passalacqua, 2015;Shaw et al., 2016b).This is due to hydrological connectivity between the distributary channels and interdistributary bays across subaqueous levees and due to reduction of channel crosssectional area over the final reach of the distributary channel (Coffey and Shaw, 2017;Hiatt and Passalacqua, 2017).One way to track the flow field in this transitional zone is through streak lines on the water surface.In many coastal settings, slicks of naturally occurring oil and biogenic debris accumulate on the air-water interface (Alpers and Espedal, 2004;Espedal et al., 1996;Garabetian et al., 1993).Despite thicknesses on the order of nanometers, streaks produced by this material are readily observed from boats (Espedal et al., 1996), in near-infrared aerial and satellite imagery, and from synthetic-aperture radar backscatter (Hühnerfuss et al., 1994).Shaw et al. (2016b) showed that the tangent of such streak lines is similar to direct measurements of flow direction, even when streak lines were mapped and measurements were made months apart.Similar streak line patterns have been observed on other delta fronts as well (Fig. 1) and should indicate flow direction where three-dimensional flow patterns and unsteady changes to flow are minimal.Our cursory analysis suggests that streak lines form mostly on deltas with established marshes that are building into freshwater basins or basins for which river discharge is enough to make the proximal receiving basin fresh.The latter condition characterizes the Wax Lake Delta and Atchafalaya Bay (Holm and Sasser, 2001;Li et al., 2011).
While streak lines record the depth-averaged flow direction field, they provide no information about the flow velocity magnitude (speed).However, if d is the unit vector field aligned with the local flow direction (dimensionless), h is the flow depth field (L), |U | is the velocity magnitude field (L/T), and temporal variations (dh/dt) are minimal, then conservation of fluid mass can be manipulated to produce a set of equations that relate spatial velocity change ( Ǎ), vertical constriction ( B), and lateral divergence ( Ď; Shaw et al., 2016b).Analyzing flow patterns on the delta front downstream of Gadwall Pass on the Wax Lake Delta, Shaw et al. (2016b) found that adverse bed slopes ( B > 0 m −1 ) were generally associated with flow direction divergence ( Ď < 0 m −1 ).In contrast, downstream of the channel tip on the basinward sloping delta foreset ( B < 0 m −1 ) the flow direction field converged ( Ď > 0 m −1 ).The transition from negative to positive Ď occurred 400 m (two channel widths) downstream of the channel tips in that study.
A converging flow direction field for delta front flows is counterintuitive: turbulent jets emanating from channel mouths generally expand ( Ď < 0 m −1 ) with distance downstream (Kundu et al., 2011) due to lateral shear with still water or bed friction.Qualitatively, jets expand because deceleration promotes an increase in the cross-sectional area of the jet core.When the flow depth is constant or decreasing, this increase in area is accomplished by widening and flow direction divergence.However, when the depth increases rapidly compared to the increasing cross-sectional area (and the jet does not detach from the bed), then the jet width must contract and flow directions must converge.A scaling of shallow water jets by Özsoy and Ünlüata (1982) showed that jets can converge or contract in width when the basinward bed slope (∇h•d > 0) exceeds the dimensionless Darcy-Weisbach friction factor divided by 8 (equivalent to the commonly used friction factor C f ).Recent numerical modeling by Jiménez-Robles et al. ( 2016) also shows that jets can exhibit flow direction convergence when basinward slopes exceed ∼ 1 %.The maximum foreset slopes on the Wax Lake Delta are about 2 × 10 −3 , which is slightly too gradual to produce flow direction convergence from either of these studies.However, we note that there is a physical basis for flow direction convergence on delta fronts that supports the convergence we observe in streak lines.

The FD2C model
If a channel tip's location controls the flow direction field, we seek an inverse method of estimating the channel tip location from the flow direction field that can be used with remotely sensed imagery.Previous analysis (see Sect. 2) showed that the transition from adverse bed slopes to basinward bed slopes at channel tips is coupled with the transition from flow divergence to flow convergence as tracked by streak lines.We name this model of coupled bathymetry and flow the C2FD model (Channel to Flow Divergence).If this correlation is persistent and predictable, then the location of the channel tips can be related to a critical point in the Ď field x Ď .Analysis of the Wax Lake Delta and numerical delta simulations show that x Ď is where Ď = 0 m −1 and Ď is changing from negative to positive in the downstream direction (Fig. 2).We name this inverse model FD2C for Flow Divergence to Channel Tips.

Application
The Ď field can be calculated using streak lines (Figs. 3,4).First, we trace the curvilinear shape of all streak lines manually in ArcGIS.Streak lines are also mapped down the center of primary distributary channels if there is a streak line or not.This is done because flow direction is generally found to follow the trends of large channels.Assuming that the local flow direction is everywhere tangent to the streak line, we sample each streak line at 25 m increments along the line, noting the local direction of the line.This produces a dataset of points P (x , y , d), where x and y are the easting and northing spatial coordinates (UTM Zone 15N) and d is the unit vector tangent to the mapped streak line.Flow direction d is recorded as a unit vector with components in the x and y directions, for which we use easting and northing coordinates: x +d 2 y = 1).The flow direction field is then constructed by interpolating d x and d y independently from P .We use the biharmonic spline interpolation technique of Sandwell (1987) because of the smooth interpolation results.The resulting fields were again normalized by their magnitude to ensure the field remained unit vectors.Finally, the flow convergence field Ď is calculated on the grid as ∂d y ∂y (Fig. 3b).For numerical models (Fig. 4), the Ď field was calculated directly from the modeled depth-averaged velocity field and thus required no interpretation of streak lines or interpolation, so the Ď field is exact in that case.

Estimating channel tip location
We test the FD2C model by comparing the location of channel tips to the critical divergence point x Ď on the Wax Lake Delta, as well as on a set of numerically modeled deltas.For each primary distributary channel, we draw a transect down the center of the subaerial reach of the distributary channel extending into the basin (Fig. 4) and track bathymetry (η(x)) and flow direction divergence ( Ď(x)) as a function of distance x along it (Fig. 5).The channel tip η is defined as the global maximum elevation along the transect with the location defined as x η.The critical divergence point (x Ď ) is defined as the first downstream location where both Ď = 0 and d Ď(x)/dx > 0 m −2 .Note that this location is interpolated along Ď(x) and therefore may appear to have sub-grid resolution.The difference l = x Ď −x η is defined as the distance downstream of the channel tip where the critical divergence point occurs (Fig. 5).Note that if l < 0 m, then the critical divergence point occurs upstream of the channel tip.
The method is designed to estimate one channel tip location that is along the subaqueously defined distributary channel axis.The benefit is that this channel axis is easily defined in imagery, but it means that the method cannot account for bends or branches in the subaqueous reach.

Validation
Summary statistics of l (Table 1, Fig. 6) provide a means of testing the FD2C model, which indicates that l is generally small.Data were drawn from the Wax Lake Delta by comparing delta front bathymetry collected in July 2010 (two channel tips), August 2011 (six tips), February 2015 (six tips), and July 2016 (seven tips) to imagery from 14 October 2010, 1 October 2011, 19 April 2015, and 5 April 2016, respectively.Over these 21 measurements, l had a mean of 145 m and a median of 132 m (Fig. 6a).The sample had an interquartile range of 701 m (Table 1, Fig. 6).This supports the claim that Ďcr is generally near η, but also shows that the variance is large.
The distribution of l on the Wax Lake Delta was compared to the unsteady hydrologic conditions present when the aerial image was collected (Fig. 7; data are in the Sup-  plement).The measurements were made over river discharge values spanning low flow to flood discharge upstream of the delta at Calumet, LA (USGS no.07381590).Measurements also spanned very low tide to relatively high tide and characteristic rates of rising and falling tide measured at Amerada Pass in Atchafalaya Bay (NOAA no.11354).In each case, the variation in l was far larger than any correlation with these parameters, and r 2 values were each less than 0.02.Even if the linear fits were statistically significant, they would explain no more than 200 m of variation in l over the common values of discharge and tidal conditions.These analyses suggest that unsteady flow is not an important control on the distribution of l.
In order to achieve some validation independent of the Wax Lake Delta, the FD2C model was also evaluated on four numerical river deltas originally presented by Caldwell and Edmonds (2014).These deltas were modeled using Delft3D on a 25 × 25 m 2 grid.Model runs A1a1, A1e1, D1a1, and D1e1 were used.These runs had an upstream discharge of 1000 m 3 s −1 and no tidal or wave forcing.They differed in incoming median grain diameter between 0.01 and 0.1 mm, the sorting of the sediment distribution, and the fraction of the sediment that was cohesive.Full descriptions of the runs are found in Caldwell and Edmonds (2014).Measurements began at time step 500 to allow a significant deposit to develop and then every five time steps thereafter.At each time step, up to five of the largest distributary channels in terms of flow velocity were measured.Fewer measurements were made if fewer than five channels were present.These analyses yielded a total of 374 samples (Table 1).
For the four modeled deltas, the median l ranged from 12 to 199 m (Fig. 6a).While the ranges of l were up to 2110 m in the case of D1e1, the interquartile ranges were between 156 and 254 m.Some transects drawn on numerical deltas did not yield x Ď because the criteria for x Ď were not met.In these failed cases, the Ď transect was always positive or trended from positive to negative.This meant that l could not be measured and the FD2C method could not be applied.Such cases accounted for 17 % of the transects on delta A1a1, 21 % of delta D1a1 (Table 2), and 8 % of the total transects measured on numerical deltas.
Measurements from the modeled deltas also show that the distribution of l is relatively stable over time (Fig. 6b).A linear regression was fit to l versus time and the slope of the data was not significant by a t test (p > 0.10 for each numerical delta).The slope that was found (1.6 ± 2.6 m yr −1 for D1e1) would introduce a small error to l relative to the uncertainty of l (order 100 m) even if it slowly grew over many decades.This near stationarity suggests that l can be assumed constant in time, even as a delta progrades.We also investigated whether l is a function of upstream channel width and flow depth at the channel tip (Supplement).However, none of these parameters showed predictive power over l.
Taken together, these analyses validate the FD2C method for prograding deltas with several distributary channels.Measurements from the Wax Lake Delta and numerical models all show that the central tendency is for l to be about 100 m with little dependence on unsteady hydrologic conditions, and model data show that the distribution remains stationary over time.We apply this result to the Wax Lake Delta to measure the growth of its subaqueous channel tips.

Methods
The FD2C method was applied to estimate the locations of channel tips over time on the Wax Lake Delta using 33 images between 30 January 1974 and 5 April 2016.Images are near-infrared imagery from Landsat 2, 5, and 8, SPOT, and an overhead photomosaic.See the Supplement for imagery metadata and the Acknowledgements and "Data availability" sections for image availability.For each image, the FD2C method (Sect.3) was applied by mapping a transect starting at the edge of subaerial exposure (delta shoreline) and extending along the seven primary distributary channel axes of the WLD (Fig. 3b) to find x Ď .The first two images (30 January 1974 and 9 February 1979) showed minimal subaerial delta exposure, so transects were mapped over abrupt changes in Ď and grouped to East Pass in the eastern portion of the delta, Gadwall Pass in the central portion of the delta, and Campground Pass in the western portion of the delta.The channel tip location was then estimated as x Ď − l, or an average 145 m upstream of x Ď according to measurements from the Wax Lake Delta itself (Fig. 6, Table 1).Channel tip growth was then tracked as the Euclidian distance between the delta apexes (UTM Zone 15N: 651673 E, 3267186 N).Estimated channel tips were connected to one another and the pre-delta shoreline to measure the area within the delta's subaqueous platform.Each channel tip occurs at a crest in bathymetric elevation and progrades via erosion of the deposit in front of it.By connecting these tips, we enclose an area that has received significant deposition, but not yet enough to become subaerially emergent, even at low tide.
The enclosed area also contains channels and all subaerially emergent regions, but excludes some seaward deposition associated with the delta forest.We name the region the total delta area.Monte Carlo techniques were used to include uncertainty in l in calculating this area.First, x Ď was found for each of the seven primary distributary channels using the technique above.The location of x η was determined by randomly sampling (with replacement) one of the 21 measured values of l that were measured on the Wax Lake Delta (Fig. 6a) and then estimating the location of the channel tip x η = x Ď − l .The seven channel tips were connected by straight lines and then connected to a pre-delta shoreline of Atchafalaya Bay mapped from 1974 imagery (Fig. 3).The pre-delta shoreline extends 10 km up the original Wax Lake estuary (Shlemon, 1972); however, delta area is truncated north of 3269274 N (Zone 15, UTM) in order for area results to be more comparable to existing datasets.The truncated area of the original Wax Lake estuary is 17.2 km 2 and can be added to all area estimates if desired.In order to join the straight lines connecting channel tips to the pre-delta shoreline, the East Pass channel tip was connected to the pre-delta shoreline with a ray of 27 • azimuth (Fig. 3b).The Campground Pass channel tip was connected to the mainland with a ray of 0 • azimuth.These azimuths were chosen to accurately reflect the marginal deposition on the Wax Lake Delta over the imagery used in the study.Total area is not sensitive to these choices.The area of the resulting polygon was calculated 10 4 times with different random sampling to account for the distribution of l (Monte Carlo sampling).The 16th, 50th, and 84th percentiles of area were then recorded for a given image.This process was repeated for each image to track delta area over time.

Results
Channel tip progradation rates are shown in Fig. 8. Between 1974 and 2016, each of the seven primary distributary channels extended at least 2 km.In clockwise order, East, Pintail, Greg, Main, Gadwall, Mallard, and Campground Pass had average progradation rates of 74 ± 9, 75 ± 13, 89 ± 13, 73 ± 10, 116 ± 10, 66 ± 15, and 60 ± 20 m yr −1 , respectively (±1 standard error).All primary distributary channels except Gadwall Pass grew at rates between 60 ± 20 (Campground) and 89 ± 13 m yr −1 , which are nearly indistinguishable given the uncertainty.In contrast, Gadwall Pass grew at a significantly faster rate of 116 ± 10 m yr −1 .Looking beyond simple linear regression, we used the "segmented" package in R (Muggeo, 2003) to search for break points or dates with different progradation rates before and after.However, no statistically significant break points were found.
The delta area estimated using the FD2C model is shown in Fig. 9.The delta area shows an apparently linear increase in area over time from 38.6 km 2 in 1974 to 113.4 km 2 in March 2016.The growth rate over this period is 1.72 ± 0.13 km 2 yr −1 and maintains this trend remarkably well over decadal timescales.The data have a root mean square error of about 7.86 km 2 associated with the Monte Carlo sampling of area.This uncertainty is generally smaller than the residuals for which data points departed from the linear trend, which averaged 13.69 km 2 over the dataset.

The FD2C model
The FD2C conceptual model assumes that water leaving a self-formed distributary channel will have a flow direction field that first diverges ( Ď < 0) and then converges ( Ď > 0) with the transition between the two fields occurring near the channel tip where the bathymetric elevation peaks and begins to slope basinward.This model was supported by measure-  ments from Wax Lake Delta and four Delft3D model runs (Figs. 5, 6).Analysis of l using modeled deltas also confirms that temporal trends in l are insignificant.We use this to assume that the distribution of l is stationary and the modern distribution can be applied to the delta in the past.This is important because field measurements on the Wax Lake Delta from the past 5 years do not provide the time span to confirm this in a field setting.Each of the l distributions have significant standard deviations (Table 1).This suggests that although the FD2C conceptual model is accurate to first order, other processes also affect l.Analysis of unsteady conditions showed that they had little predictive power over l measured on Wax Lake Delta (Fig. 7).Wind shear could also play a role, but local wind measurements were only available for two of the four images in which l was measured, and in each case the wind was light (< 4 m s −1 ) and from the east.Hence, we have insufficient data to test the effect of wind setup at this time.However, the wide spread of l even for a single image makes wind control unlikely.
There is also uncertainty associated with the use and interpolation of the streak lines in calculating the Ď field in field cases.Analysis of l on the Delft3D deltas showed standard deviations that were 35 %-81 % of the standard deviation on Wax Lake, suggesting that even when unsteadiness and interpolation errors are neglected, the variation of l remains significant.Given that unsteadiness and interpolation errors explain only portions of the l distribution, we hypothesize that the remainder stems from channel properties.If two channels are near one another, their outflows and the unchannelized flow between them would be constricted compared to a channel that is far from its neighbors.This may be particularly important in places like the eastern portion of the Wax Lake Delta, where Main, Greg, Pintail, and East passes enter Atchafalaya Bay over about 8 km (Fig. 3a).Channel dimensions and input discharges also affect hydrodynamics, including the Ď field.Finally, aspects of channel tips that have been less studied, such as the bed slope of the channels or number of branches in the subaqueous delta front, could have important effects.We expect that further study of flow patterns and bed morphology on complex deltas will shed more light on what sets l.
The variation of l prevents confident interpretations of changes in channel tip location on seasonal or annual timescales.For example, extension and back-stepping of channel tip location on the order of several hundred meters were directly measured on Gadwall Pass between July 2010 and February 2012 (Shaw and Mohrig, 2014).The 512 m standard deviation of l prevent these changes from being estimated with confidence.Even with this limitation, cer-Earth Surf.Dynam., 6, 1155Dynam., 6, -1168Dynam., 6, , 2018 www.earth-surf-dynam.net/6/1155/2018/tainty in measuring change increases with time.The ∼ 60-116 m yr −1 growth rates observed on distributary channels grow larger than the standard deviation after 4-9 years.For platform area growth analysis, the standard deviation produced by Monte Carlo sampling of l measurements is about 8.2 km 2 , and progradation rates are calculated to be 1.72 km 2 yr −1 .Hence, for estimating changes in Wax Lake Delta total area, we expect the method to be able to perform on timescales greater than 4-5 years.Clearly, more research on the effects of channel and delta morphological characteristics and unsteady flows on l are warranted and could increase the sensitivity of the FD2C method for monitoring change.However, the method already detects clear changes in channel tip and delta area at the decadal scale or better for the Wax Lake Delta.
Where can the FD2C model be applied?The model was validated on deltas that were prograding with several active channels in which lateral channel migration was minimal, and there are many tens of deltas globally that have these characteristics.Streak lines must also be present if the Ď field is to be estimated from remote sensing, and streak lines are somewhat rarer, but our cursory search has revealed at least 10 deltas globally (e.g., Fig. 1) with detectable streak lines under certain conditions.In particular, streak lines are common around the bird-foot delta of the Mississippi River's main stem where many coastal restoration projects are planned or are currently operational.One example is the West Bay Diversion, for which the progradation of a delta is an explicit goal (Allison et al., 2017;Andrus and Bentley, 2007;Kolker et al., 2012).The method can also be used with the decades of remote sensing imagery that already exist.Our analysis of the Wax Lake Delta found that streak lines were sometimes visible in synthetic-aperture radar backscatter, Landsat 1, and CORONA imagery, making monitoring from the 1960s or earlier potentially feasible.The FD2C model could also be used as a hypothesis for the hydrodynamic study of other similar deltaic systems, even if they do not support streak lines.For example, flow direction convergence occurred on nearly all delta foresets analyzed here, so it is possible that this pattern occurs but remains unmeasured on many of the world's river deltas.Converging flow patterns could inform the study of delta front sedimentology (e.g., Enge et al., 2010) or the hydrodynamics of a plunging river plume (Lamb et al., 2010).Whether applied to monitor growth or understand delta front hydrodynamics, we present the FD2C model as an advance in understanding the coupled morphology and flow field at distributary channel tips.
5.2 Subaqueous growth of the Wax Lake Delta

Channel tip progradation
The FD2C model allows the progradation rates of individual subaqueous channel tips on the Wax Lake Delta to be measured and provides new insight into decadal growth patterns of the Wax Lake Delta from its initiation to present.The hypothesis of radially symmetric growth often applied to the Wax Lake Delta (Kim et al., 2009;Paola et al., 2011) is largely supported: six of the seven channels have prograded at rates between 60 and 90 m yr −1 .However, the consistently larger progradation rate of Gadwall Pass (116 ± 10 m yr −1 ) also suggests that the delta is becoming more asymmetric over time.Future evolution may correct for the dominance of Gadwall Pass, possibly by a soft avulsion (sensu Edmonds et al., 2011a) whereby Gadwall Pass's progradation decelerates and another channel accelerates.However, the consistently dominant growth rates since 1983 (Fig. 8) and the fact that Gadwall Pass is presently the widest channel (Figs. 1a,3a) suggest that this dominance will continue and delta asymmetry will continue to grow.

Delta area growth
Previous studies of delta growth have focused on the emergence of subaerial land using Landsat imagery.Allen et al. (2012) investigated the area of subaerial land growth over the entire Wax Lake Delta.They determined a growth rate in Landsat imagery as a function of time, water discharge, and tide level.They found that the subaerial delta grew at a rate of 1.1 km 2 yr −1 between 1983 and 2002 and then reduced growth to near zero afterward.Olliver and Edmonds (2017) focused on the emergence of just the central islands of the WLD, neglecting some marginal areas of the delta included by Allen et al. (2012).Analyses were based on two images per year selected for minimum and maximum biomass, which mitigated the large swings in area shown by Allen et al. (2012).The authors also interpreted a break in growth rate at about 1999, with a growth rate from 1984-1999 of 1.88 ± 0.42 km 2 yr −1 and a growth rate from 1999-2015 of 0.78±0.44km 2 yr −1 .The total delta area as measured by the FD2C method grew at a rate of 1.72 ± 0.13 km 2 yr −1 from 1974 to 2016 without any break points in growth rate.
In all cases, the FD2C method produces delta area estimates that are > 40 km 2 larger than estimates of subaerial land.The FD2C method is designed to track the location of subaqueous channel tips significantly below any water level datum.Furthermore, the FD2C method includes distributary channels and subaerial land as part of the delta area.Therefore, it stands to reason that the area estimates would be far larger than subaerial methods.However, we also note that the total delta area as it is defined here neglects some foreset deposition, so is not necessarily an upper bound on delta area.The focus on the Wax Lake Delta also ignores the ∼ 10 000 km 2 deposit of silt and clay accumulating on the Atchafalaya Shelf that is fed by the Wax Lake and Atchafalaya Delta (Draut et al., 2005;Neill and Allison, 2005).
Despite the larger absolute area, the FD2C method yielded long-term growth rates that were broadly similar to subaerial www.earth-surf-dynam.net/6/1155/2018/Earth Surf.Dynam., 6, 1155-1168, 2018 growth rates from initial emergence until 1999 (Fig. 9).During that time, the subaqueous platform was always ∼ 40 km 2 larger than the subaerial platform, suggesting that the subaqueous platform was being converted to intertidal or subaerial land at a rate similar to the rate of subaqueous platform production through progradation.
From 1999 to present, the subaqueous total delta area creation continued unabated, but the processes converting subaqueous platform to intertidal or subaerial land were reduced by 59 % or more.We suggest that this is a consequence of channels prograding radially and becoming further apart with distance from the delta apex.Field and remote sensing studies of the Wax Lake Delta have confirmed that island aggradation rates decrease and timescales of emergence increase with distance from the edges of the primary channels (Bevington and Twilley, 2018;Olliver and Edmonds, 2017;Wagner et al., 2017).We reason that increased island width limits vertical accretion because suspended sediment from the channels must travel farther through a vegetated island.Despite this reduction in subaerial growth rate, the total delta area growth rate was roughly constant, consistent with a roughly constant sediment supply.The growing disconnect between subaqueous progradation and island aggradation suggests that these processes begin to decouple as deltas become more mature.
This has interesting implications for the Wax Lake Delta and coastal restoration initiatives.Ecology (Carle et al., 2015;Olliver and Edmonds, 2017) and carbon sequestration (Shields et al., 2017) are highly dependent on island elevation, and the ratio of subaerial land to total delta area appears to decrease with time.The coupled subaerial-subaqueous monitoring scheme used to discover this transition is likely transferrable to other large-scale coastal restoration efforts in Louisiana (streak lines are frequently observed across coastal Louisiana), allowing sediment accumulation and marsh formation to be independently tracked.

Conclusions
The morphodynamic evolution of channel mouths can produce flow patterns and bed morphology that are closely coupled.The delta front morphology consists of subaqueous channels that grow shallower in the downstream direction, subaqueous levees that allow water to exit the channel laterally, and a sloping delta foreset.The flow direction field over this morphology diverges in the final reach of the channel and converges on the delta foreset.The location of the transition from divergence to convergence relative to the channel tip ( l) varies by many hundreds of meters, but is on average 0-200 m downstream of the channel tip in both field data from the Wax Lake Delta and numerically modeled deltas.This distribution of l appears to be independent of hydrodynamic unsteadiness.It also appears stationary, allowing it to be applied through time.We present the FD2C method to relate the divergence of flow direction estimated using remote sensing of streak lines on the water surface to channel tip location with quantitative uncertainty.
The FD2C method provides a means of estimating the progradation of channel tips and total delta area of the Wax Lake Delta from its initiation in 1974 through 2016.The method involves uncertainties associated with flow field characterization and channel tip location estimation.However, the method allows key aspects of the Wax Lake Delta's progradation to be characterized for the first time, such as individual subaqueous channel tip progradation rates and the growth of the total delta area.Channel tips grow at rates ranging from 69-116 m yr −1 that appeared constant at the decadal scale.The subaqueous total delta area also grew steadily between 1974 and 2016 at a rate of 1.72 ± 0.13 km 2 yr −1 .The reduction of subaerial growth rates below this rate around 1999 suggests that the Wax Lake Delta has become less efficient at building subaerial land as islands have grown wider, even as subaqueous channel extension continued unabated.This monitoring technique furthers our understanding of the Wax Lake Delta, which is a example of an uncontrolled river diversion being investigated as a possible land-building strategy in Louisiana.The FD2C model can be applied to similar deltas in which direct field measurements are impossible or scarce.
Data availability.Imagery metadata, channel tip location estimates, and time series of area are all available in the Supplement of this paper.Interpolated images are available in an online repository (Shaw and Haynes, 2018).Landsat imagery used in this study was downloaded from Google Earth Engine.Bathymetric maps of Wax Lake Delta are available in Shaw (2013) and Shaw et al. (2016a).

Figure 1 .
Figure 1.Images of river deltas exhibiting streak lines.(a) Wax Lake Delta (Landsat image LT50230402011002CHM01 band 4).(b) The Saint Clair Delta in Michigan, USA (digital orthophoto quads, Saint Clair Flats NE, NW, SW, SE; IR band).(c) Portion of the Volga Delta in Russia (Landsat LC08_L1TP_168028_20170501_01 band 4).In each image, streak lines are mapped in dashed lines.The line is translated in space slightly (see gray arrows) so as not to cover the streak line in the image.

Figure 2 .
Figure 2. Schematic diagram of delta front morphology and streak line behavior.The color map shows topography with dark colors representing deep areas and light colors representing shallow or subaerial areas.Streak lines are shown as black solid lines.The FD2C method takes advantage of flow direction divergence ( Ď < 0) through the shoaling reach of the channel and lateral flow direction convergence ( Ď > 0) on the basinward sloping delta front.The channel tip occurs roughly where Ď transitions from negative to positive.

Figure 3 .
Figure 3. Method for converting imagery into channel tips and delta area.(a) Landsat image displaying streak lines (14 October 2010; Supplement).The seven primary distributary channels are labeled.(b) Streak lines (thin black lines) are mapped manually on the delta front, and lines are also placed down the center of subaerially emergent distributary channels.The Ď field is interpolated from these streak lines (color map).Thick black lines are transects extending from the seven primary distributary channels.The estimated location x η along each transect is connected via the purple line and rays connect channel tips to the pre-delta shoreline to close the area.(c) The interpreted total delta area is shown.A bathymetric map from June 2010 referenced to mean lower low water (MLLW) shows how the interpreted channel tips compare to direct measurements.

Figure 4 .
Figure 4. Method for comparing Ď to bathymetry using a Delft3D numerical simulation (run A1e1).(a) The velocity field and −0.3 m MSL (green) elevation contour are displayed and transects (pink and red lines) are drawn extending from the largest distributary network channels.(b) The bathymetric profile is collected along each transect.(c) Ď is calculated, and transects of Ď are collected along the transects.Transects A-A' and B-B' are shown in Fig. 5.

Figure 5 .
Figure 5.Comparison of bathymetry (η; panels a and b) and divergence of flow direction ( Ď; panels c and d) for transects A-A' (panels a, c) and B-B' (panels b, d).l is the width of the gray box, or the location where Ď changes from positive to negative (x Ď ; red circle) minus the bathymetric maximum of the channel tip (x η; black circle).The distribution of l is shown in Fig. 6.

Figure 6 .
Figure 6.(a) Histogram of l for numerical models A1a1 (dark blue), A1e1 (light blue), D1a1 (seafoam), and D1e1 (yellow) compared with measurements from the Wax Lake Delta (pink).All histograms are binned at 250 m intervals.Descriptive statistics of these populations are shown in Table 1.(b) The location of l as a function of time (model years) for delta run D1e1.A linear fit (black line) with 50 % confidence interval (gray lines) is shown.The trend of this fit is not statistically significant and is small compared to variation within l.

Figure 7 .
Figure 7.Comparison of l measurements on the WLD to unsteady hydrodynamic conditions.In each case, the black line shows the linear trend, which is negligible.(a) Upstream water discharge, with r 2 = 0.01.(b) Tidal elevation relative to mean lower low water with r 2 = 0.001.(c) Rate of change of tide elevation averaged over 30 min with r 2 = 0.01.

Figure 8 .
Figure 8. Growth of individual channels over time.Each series is plotted at the same scale (horizontal lines 1 km, vertical lines 10 years), but shifted vertically for clarity.The uncertainty associated with l measured at Wax Lake Delta (Fig. 6a) is shown with a standard deviation (σ ) and 90 % confidence interval (CI) at the bottom.The primary distributary channels are shown from west to east: Campground Pass (maroon squares), Mallard Pass (turquoise right-pointing triangles), Gadwall Pass (green left-pointing triangles), Main Pass (purple down-pointing triangles), Greg Pass (yellow up-pointing triangles), Pintail Pass (red squares), and East Pass (blue circles).Dashed lines show linear regressions of each dataset.

Figure 9 .
Figure9.Area of the Wax Lake Delta as a function of time.Purple circles show the growth of the total delta area of the Wax Lake Delta using the FD2C method.The gray region shows the 1σ deviation (16th to 84th percentile) of area found from Monte Carlo sampling of l (Sect.4.1).The dashed line shows the linear fit, with a growth rate of 1.72 ± 0.13 km 2 yr −1 .Brown squares and greed triangles show the subaerial area as documented by Olliver and Edmonds (2017) andAllen et al. (2012) over time.

Table 1 .
Statistics describing the distribution of l for the Wax Lake Delta (WLD) and four delta simulations.SD: standard deviation, Iqr: interquartile range, Min and Max: minimum and maximum, n: number of measurements, n misses: the number of measurements for which l could not be computed.