Subtidal Flow Reversal Associated With Sediment Accretion in a Delta Channel

In branching delta channel networks, tides intrude into parallel river outlets causing complex tidal behavior. The tidal motion can impact the division of river discharge over distributary channels. Under some circumstances, the discharge averaged over a tidal cycle may even reverse, drawing ocean water into the delta, such as observed in the Yangtze Delta, the Fly Delta and the Colorado Delta. Here, we study the flow reversal associated with sediment accumulation of a distributary channel in terms of tidal propagation and subtidal discharge dynamics, by focusing on the Yangtze Delta. The Yangtze Delta channel configuration represents a deep and wide main channel and a smaller side channel that has rapidly accreted over the past decades. A new mechanism is presented, which results in seawater transport across the shallow side channel, posing a risk to freshwater availability. The shallowest section in the side channel nearly acts as a tidal divide. In this section, the tide has the character of a standing wave, which implies that Stokes transport converges from opposite sides of the channel. This leads to significant variation of the total water storage in the side channel and subtidal water levels in the side channel being elevated above that in the main channel. The bulge of water that is stored during spring tide leaves the side channel toward neap tide, explaining reversal of the tide‐averaged discharge and seawater intrusion when the river discharge is low.


Introduction
Flow division processes at river bifurcations govern the allocation of water and sediment discharge over distributaries and have a significant influence on the evolution of the downstream channel morphologies (Bolla Pittaluga et al., 2003;Burge, 2006;Kleinhans et al., 2008;Zhou et al., 2017). This issue is equally important but more complicated for estuarine bifurcations subject to tides (Hoitink & Jay, 2016;Wolinsky, 2010). Tides propagating toward estuarine bifurcations from the ocean (see Figure 1 for an example) can reduce or enhance inequality in the flow allocation that would happen in absence of the tidal motion (Buschman et al., 2010;Sassi et al., 2011). Tidal impacts on the flow division are dependent on geometrical properties of the distributary channels, including depth, length, roughness, and convergence of the outlets (Zhang et al., 2010). Asymmetry in these geometrical properties of bifurcations leads to tideinduced residual circulation in estuarine bifurcations (Buschman et al., 2010;Guo & Valle-Levinson, 2007;Liu et al., 2009;Zhou et al., 2018). Deepening of a single channel in a delta to improve navigation conditions (Vellinga et al., 2014) or to extract sand  generally enhances the asymmetry in the discharge distribution, favoring the deeper channel such that side channels receive a progressively smaller part of the discharge. In the partially abandoned side channels, the landward sediment transport induced by the flood-dominant tidal motion dominates over the seaward flushing of sediment during high river discharge, resulting in sediment accumulation . Here, we describe a new mechanism explaining how this development may result in seawater being transported through a shallow side channel, into the delta interior.
In delta channel networks, the tidal motion may influence the river discharge distribution by generating Stokes transport, which can be defined as Stokes velocity integrated over the width of a channel. Stokes velocity, in turn, refers to the difference between the Lagrangian velocity and the Eulerian mean velocity, causing a net drift of a fluid element in a propagating wave. Recent research has revealed the importance of Stokes transport in distributary channels (Buschman et al., 2010;Zhang et al., 2017). For a progressive ©2019. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2019WR025945
Key Points: • In a width-convergent delta channel subject to accretion, the shallowest part of the channel may act as a tidal divide • Stokes transport converges at this tidal divide from opposite sides, which causes subtidal variation of water storage in the channel • Reversal of the subtidal surface level gradient causes the residual flow to intrude far upstream in the delta

Supporting Information:
• Data Set S1 • Data Set S10 • Data Set S11 • Data Set S12 • Data Set S13 • Data Set S14 • Data Set S15 • Data Set S16 • Data Set S17 • Data Set S18 • Supporting Information S1 tidal wave in deep water, the phase lag between the water surface level and flow velocity is 0, leading to a maximum Stokes transport. For tides in distributary channels a phase lag occurs, caused, for example, by the effects of bottom friction and the river discharge (Buschman et al., 2009;Godin & Martinez, 1994;Godin, 1999). The Stokes transport reduces to a minimum when the phase difference approaches 90°, which corresponds to a standing tidal wave. In a multiple channel system, the Stokes transport in a particular channel is not necessarily equal to the return flux integrated over the channel cross section. The mass transport generated by (nearly) in-phase variation of the water level and cross-section averaged flow velocity in a tidal channel can be returned to sea through a different channel. This can have significant implications for the division of river discharge over bifurcated channels (Alebregtse & de Swart, 2016).
The division of return currents that compensate for Stokes transport over parallel distributary channels depends on tidal mean water level gradients that occur locally at a channel bifurcation, which depend on the aforementioned geometrical channel properties identified by Buschman et al. (2009). Sassi et al. (2011) focused on the tidal impacts on the flow division over tidal river channels in an estuary, showing that the effect of the tide is generally to equalize the flow division that would occur in the absence of tides. The differences in water level setup induced by river-tide interaction between neighboring channels play a crucial role in this. The enhanced subtidal water level gradient may steer the flow to shallower and longer channels in the Mahakam Delta (Sassi et al., 2011). In a subsequent study, Sassi et al. (2012) quantified the associated tidal signature on delta morphology by adopting the hydraulic geometry concept, which refers to a set of empirically derived power law relations between a river channel's width, mean depth, and mean flow velocity and the discharge conveyed by the channel. The hydraulic geometry of distributary channels revealed that the effect of the tide on flow division at bifurcations in a tidal river estuary may increase with the along channel distance from the delta apex, exerting a progressively larger influence on channel sections located closer to the coast.
Here we focus on the extreme situations in which the tide may act to occasionally reverse the direction of the discharge averaged over a tidal cycle in a distributary, which was observed in the Yangtze Delta in China.
Previous studies have shown that this reversal of the subtidal discharge is the main cause of seawater intrusion in the Yangtze Delta Xue et al., 2009) and can be adequately represented in a barotropic model (Zhang et al., 2017). From these full-complexity modeling studies, the effects of channel depth and width convergence in generating the flow reversal remain unclear. In this manuscript, based on results from a reduced complexity model, we will show how the propagation of tides varies under different configurations of channel depth and width convergence. Details of the channel configuration control the interplay of Stokes transport and the Eulerian mean current, governing the occurrence of subtidal flow reversal.
In a single channel estuary, it is often assumed that the Stokes transport balances the Eulerian mean "return" current, such that the subtidal discharge is equal to 0 (e.g., Hoitink & Jay, 2016). This is not necessarily the case, because the tide-averaged total water storage in a channel can differ for successive tidal periods in a spring-neap cycle. We show that in a two-channel system, subtidal water storage variation can become very important. We demonstrate how water accumulation in a side channel may create a bulge of water resulting in a reversed longitudinal water level gradient, driving a landward tide-averaged Eulerian flow that dominates the total discharge.
Section 2 explains the motivation of this study in more detail. Section 3 presents the reduced complexity model inspired by the bifurcation of the Yangtze Delta. Section 4 analyzes the tidal behavior in the idealized model and describes the mechanisms governing the tide-averaged flow. Section 5 discusses the key insights gained from this study, and section 6 finalizes this paper with conclusions.

Motivation: The bifurcation of the Yangtze Delta
The Yangtze Delta is divided into the North Branch (NB) and the South Branch (SB) by Chongming Island. The North Branch is a funnel-shaped channel with its upstream part perpendicular to the South Branch. The length is about 80 km, and the width is narrowing from 10 km at the downstream outlet to 2 km at its upstream entrance. (Lu et al., 2015;Xue et al., 2009). The South Branch is about 120 km long between Baimao and Hengsha, and about 12 and 30 km wide at these two transects, respectively (Wang et al., 2013). The South Branch is divided into the South Channel and North Channel and the South Channel is further bifurcated into the South Passage and the North Passage. In comparison to the North Branch, the South Branch and its downstream distributary channels are wide and deep, with smaller spatial variation in bathymetry. During the past one and a half century, the South Branch has become narrower and deeper. The mean water volume in the channel increased by 15% in this period (Wang et al., 2013). As a result of tideinduced sedimentation and land reclamation over the past half century, the upstream area of the North Branch has gradually become narrower and shallower, with current depths ranging between 1 and 8 m (Figure 1).
The planar configuration of the North Branch leads to deformation of the tidal wave as it propagates inland (Chen et al., 2003). The convergence effect of the North Branch increases the tidal range significantly, such that it exceeds the tidal range in the South Branch. Observations revealed that the tidal prism entering at the mouth of the North Branch accounts for 25% of the total tidal prism in the Yangtze Delta.
Previous studies have reported that the year-averaged percentage of the river discharge issued to the North Branch continued to decline from~25% to~1% over the past century (Chen et al., 2003;Li et al., 2010;Zhang et al., 2011), which implies that the North Branch becomes abandoned. Since there is little influence of the Three Gorges Dam on the annual river discharge change in the Yangtze River (Guo et al., 2018;Wang et al., 2016), this dramatic decrease can be attributed to the channel evolution of the delta channels, and the upstream part of the North Channel in particular. Due to the process of sediment accumulation, the river bed in the upstream entrance of the North Branch was elevated, reducing the allocation of river discharge to this side branch. The enhanced flood-dominant tidal currents carry large amounts of sediment from the offshore into the North Branch, explaining a sustained decrease in year-averaged water volume in the branch, from 11.7 × 10 8 to 5.4 × 10 8 m 3 in the past 50 years. Those values were calculated with respect to mean water level. On the contrary, the area of the North Branch's mouth has a trend of expansion (Dai et al., 2016). During the dry season, the tide-averaged flow reverses in the North Branch (Wu et al., 2006;Wu & Zhu, 2007), which causes seawater to be transported far upstream into this branch. Under extreme circumstances, seawater can even reach the upstream area of the South Branch, beyond the bifurcation.
Reversed (i.e., landward) residual flow occurs in many other estuarine channel networks where a river splits into a river-dominated branch and a smaller tide-dominated side branch. Herein, residual flow refers to flow velocity averaged over a tidal cycle. The Fly (Harris et al., 2004) and Colorado (Carriquiry & Sánchez, 1999) deltas are prime examples. The planar configuration of these systems is funnel shaped, leading to a progressive asymmetry of the tidal wave along its propagation inland. The Fly River bifurcates into three primary downstream channels (the Southern, Northern, and Far Northern Entrances) in the delta area.
Observations show that the Southern Entrance, the main route, transports 60-80% of the river discharge to the Gulf of Papua. The Far North Entrance, which is effectively abandoned at present (Dalrymple et al., 2003), delivers less than 10% of the discharge (Wolanski et al., 1997). A reversal in the direction of the residual current is observed here. The Colorado Delta contains two main distributary branches (the Baja Californian and the Sonoran channel) that branch from a common point. Similarly, results derived from observations show that two opposing net flows occur along these two branches (Carriquiry & Sánchez, 1999), complicating the residual circulation. The direction of sediment transport along the Baja Californian channel is seaward and landward along the Sonoran channel. The interaction of geometry and hydrodynamics in the shallow side channel influences the flow division and the associated seawater intrusion. The purpose of this paper is to present an idealized model that captures the key mechanism governing the tide-averaged flow reversal in an accreting distributary of a multichannel delta system, taking the Yangtze Delta as a starting point.

Parameter Setting
To study the effect of channel geometrical properties on the water motion in the delta channel system, we use alternative channel planform configurations. Whereas many of the processes we investigate are essentially 1-D, we use a 2-D depth-averaged model, as it may better resolve the mass and momentum exchange at the channel junction. Two idealized hydrodynamic models of the Yangtze Delta were established in a 2-D setup based on Delft3D, solving the unsteady shallow water equations (Lesser et al., 2004). A 2DH barotropic approach was adopted because the landward net flow is the key process to cause the seawater intrusion (Xue et al., 2009). The two models represent a river that splits in two branches that debouch in the sea and mainly differ in the geometry of the distributary channels. The idealized models, labeled 1 and 2, adopt a similar shape for the distributaries but differ in their convergence at the outlet ( Figure 2).
The 2-D schematized model is constructed by a 600-km-long river channel with the width of the South Branch varying from 8 km at the tidal bifurcation to 2 km at the upstream boundary. In the default setting, the tidal bifurcation is located 100 km from the open sea boundary (Figure 2), corresponding to the distance from the first tidal bifurcation in the Yangtze Delta to the river mouth. The upper part in this model is chosen longer than the tidal wavelength to create a continuous landward tidal decay and to avoid tidal wave reflections at the landward end (Guo et al., 2014(Guo et al., , 2015.
The width of the channel between the upstream boundary and the tidal bifurcation increases exponentially, as in where L W is the e-folding length, defined as the length over with the width decreases by a factor e. We chose a grid that is numerically efficient yet still provides a sufficient resolution. The width of the south branch and the main channel from the river boundary to the bifurcation ( Figure 2) was initially equally divided over 32 grid cells across the channel. For the North Branch, there are eight grid cells along the cross sections in the initial setting. The grid independence of numerical results was tested with halved and doubled grid cell sizes (Hasan et al., 2012). To evaluate the sensitivity of model results to grid cell size, the deviations of water level and current speed were calculated as We apply this equation to simulation results along the central lines in the North Branch and the South Branch (red lines in Figure 2), where α i represents the water level or current speed derived from the simulation with changed grid size (halved or doubled), α 0 denotes the water level or current speed derived from the case with the initial grid size setting, and n represents the length of the input data series, which corresponds to the duration of a spring-neap cycle divided by the time step in the calculations. The unit of deviation (σ) depends on the unit of the sample (α i in equation (2)). In our case, the deviations of surface elevation and current velocity are 6.8 × 10 −7 m and 9.2 × 10 −6 m s −1 , respectively. Values of σ were all less than 1 × 10 −5 , which indicates that the subtidal water levels and velocities were insensitive to the choice of the grid resolution. A time step of 60 s was used in this model.
In Table 1 the default parameter setting is given for the idealized models. To eliminate the effect of varying bed roughness, a time and space invariant Chézy coefficient of 55 m 1/2 s −1 was used and the horizontal eddy viscosity was set to 100 m 2 s −1 . The longitudinal bed profiles are the same in the two models. The bed level increases linearly form the bifurcation to the upstream boundary, where the still water depth is 2 m (Figure 3a). Downstream of the bifurcation, we adopted a simplified bed level profile for the South Branch, obtained by smoothing bathymetric measurements (Figure 3b). For the North Branch, we initially adopt a constant depth of 8 m and carry out additional simulations in which the profile gradually develops toward a simplified profile, obtained by smoothing measurements.  The principal lunar (M 2 ) and the principal solar (S 2 ) tidal constituents are applied to generate the spring-neap cycle in a semidiurnal tidal regime (Pugh, 1996). At the seaward boundary, an M 2 tidal harmonic with a constant amplitude of 1.2 m is imposed, as well as an S 2 tidal harmonic with an amplitude of 0.6 m. These amplitudes were inferred from the tidal harmonic analysis following Lu et al. (2015), using the measured surface elevation from the tide level gauges (Hengsha and Lianxinggang in Figure 1 b) at the outlets of the Yangtze Delta. At the upstream boundary, a steady, uniform discharge of 6,000 m 3 s −1 is issued to the main channel, representing the low flow situation during the dry season at Datong station (i.e., the tidal limit of the Yangtze Estuary). Upstream of the bifurcation, the slope in the river is 1.2 × 10 −5 , which guarantees a sustained water level setup throughout the model governed by the discharge imposed at the upstream boundary (Kong et al., 2016). Figure 4 shows how the M 2 tide propagates into the branches in Model 2 and in the realistic model (Lu et al., 2015;Zhang et al., 2017). The two models feature similar tidal behavior. The amplitudes of the M 2 tidal constituent reduce with increasing distance from the outlets, toward the bifurcation, and its phases change rapidly. The phase in the North Branch peaks at the shallow section, where the tidal motion is partly dissipated. Both shoaling toward the shallow section and convergence of the channel width increase the M 2 amplitude in the lower reaches of these distributary channels, albeit that the effects are subtle. The difference between the amplitudes and phases between the two branches may induce an imbalance of the return flow (Buschman et al., 2010), which can directly change the flow division at the bifurcation. The idealized model captures this behavior.

Model Verification and Basic Tidal Behavior
The phase lag between variation of the water surface elevation and current velocity of the M 2 tide in Model 2 is shown in Figure 5. In the South Branch, the phase lag varies from 30°to 45°. The phase lag is relatively large in the North Branch, reaching nearly 90°when the tide propagates into the area near the shallowest section, which characterizes a standing wave. The increase and subsequent decrease of the phases indicate a tidal wave convergence location. The correspondence of the convergence location with the shallowest section in the channel marks the strong control of the channel bathymetry on tidal propagation.

Discharge Asymmetry Index
The subtidal discharge division at a bifurcation is here quantified as (Buschman et al., 2010;Sassi et al., 2011Sassi et al., , 2012: where brackets 〈〉 indicate a diurnal tidal average and suffixes SB and NB stand for the South Branch and the North Branch, respectively. This metric is referred to as the discharge asymmetry index (DAI), which represents the tide-averaged flow division and is used to analyze processes at subtidal time scales. For an equal flow division, the index is 0. When more flow is directed to the South Branch, the index is positive and vice versa. Usually, the value of the DAI varies from −1 to 1, which indicates that the direction of the tide-averaged flow in both distributary channels is seaward. When the DAI attains values out of this scope, the direction of the tide-averaged flow in one of the two channels is landward.
In order to detect the effect of channel width convergence on flow division, Model 1 and Model 2 are run by applying the same depth profiles in the distributary channels. The depth is constant (8 m) along the North Branch, and for the South Branch the longitudinal bed profile smoothed based on measurements is used ( Figure 3b). The simulation results reveal that values of the DAI in Model 1 ranges between 0.569 and 0.667 over a spring-neap tidal cycle. For Model 2, the DAI values range between 0.657 and 0.720. The higher DAI can thus be attributed to width convergence, which increases the M 2 velocity amplitudes. Because of the stronger convergence in the North Branch, the M 2 velocity amplitude increase in this channel is more significant. As a result, the tidal prism in Model 2 is 125% larger than in Model 1 in the North Branch and 63% larger in the South Branch. The DAI cannot exceed 1 by adjusting convergence only, demonstrating that more factors should be varied to reproduce the subtidal flow reversal in the North Branch.
The sensitivity of the DAI to depth variation in the North Branch was investigated based on a series of simulations with the depth profiles as in Figure 3b. The depth of the South Branch (d 1 ) was kept constant according to the simplified profile, while the minimum depth in the North Branch (d 2 ) was varied from 0.8 to 8.0 m (d 2 /d 1 = 0.1-1.0). By elevating the depth profile in the North Branch, the systematic change in discharge division during growth of the bar (the shallowest section) was investigated. Figure 6 illustrates the changes in DAI as a function of the depth ratio that varies over a complete spring-neap period, which corresponds to the period of the lunisolar synodic fortnightly tidal constituent, MSf (Pugh, 1996). Results calculated by averaging over an MSf period are in between the results of neap tide and spring tide (shaded areas in Figure 5). Both for spring tide and neap tide, the values of Ψ obviously increase when the depth ratio d 2 /d 1 declines. There is no DAI greater than 1 in the Model 1, which means that although the depth effect is significant, it still cannot explain reversed residual flow in the North Branch.
Only the combined effects of width convergence and channel shallowing in Model 2 results in subtidal flow reversal in the North Branch near the junction.
For comparison, the DAI at the first bifurcation derived from the realistic model during the dry season (Zhang et al., 2017) is represented by the green bar. The DAI difference between the realistic model and Model 2  with a smallest depth ratio can primarily be explained by the fact that in the reduced complexity model only the dominant tidal constituents M 2 and S 2 were included. Since tidal constituents within a tidal species act similarly, the neglected semidiurnal constituents would amplify the effect of the dominant M 2 and S 2 tidal constituents. Therefore, in the full complexity model, the flow reversal occurs at smaller amplitudes of the main tidal constituents than in the idealized model.
Additional simulations were performed to establish the minimum channel expansion needed to achieve subtidal flow reversal for a bed elevation profile with a minimum depth of 0.8 m, as shown in Figure 3. In these simulations, we varied the width at the outlet and adjusted the e-folding length for the width convergence such that at the shoal, 65 km from the outlet, the channel width had converged to 2 km. Inland from the shoal, the width remained constant. We find that subtidal flow reversal requires a minimum width at the outlet of only 2.2 km, with an e-folding length of 682 km.
Further, we investigated how sensitive the length of the transect over which subtidal flow reversal occurs responds to width convergence of the North Branch, adjusting the e-folding length of width convergence as detailed in Figure 7. When the width of the outlet in the North Branch increases from 10 to 50 km, the length of the flow reversal region increases only slightly, from 13.50 to 14.25 km. When the funnel shape of the channel widens, the length of the flow reversal region becomes progressively less sensitive to the width convergence rate and does not increase any further beyond an outlet width of 30 km. Finally, we explored how small the minimum depth in the North Branch would have to be to obtain subtidal flow reversal in case of a constant channel width, as in Model 1. The additional simulations showed that reducing the minimum depth to 0.55 m is sufficient. Hence, expansion of the side channel enhances tendency for subtidal flow reversal, but it is not a necessary condition per se.

Discharge Decomposition
The cross-section averaged flow velocity was decomposed according to U = < U>+U′, where the prime denotes the zero mean variation during a diurnal tidal cycle and the angular brackets 〈〉 indicate averaging over a diurnal tidal cycle. The depth can be written as d = h+<ζ>+ζ′, where h represents the still water depth, which indicates the channel depth when the surface elevation is 0. Then, <ζ> denotes the tide-averaged water surface level deviation from the still water level and ζ′ indicates the tidal elevation relative to the tide average elevation. Assuming a time-invariant channel width, subtidal discharge can be rewritten as where Q S denotes Stokes transport and Q N is the Eulerian mean discharge. The distributions of Q N and Q S can be highly nonuniform in convergent channels, being larger closer to the sea, and smaller going landward (Buschman et al., 2010). Figure 8 shows the longitudinal variation of Q t , Q S , and Q N along the North and South Branches over a spring-neap cycle. The shaded areas in this figure represent the extreme values of these terms, which typically occur during spring tide and neap tide. The solid lines indicate the spring-neap tide-averaged terms, which are derived by computing the average value of the results in each diurnal tidal cycle over a complete MSf period. In the seaward part of the North Branch, Q S is largely balanced by Q N and the absolute magnitudes of those variables weakly decrease and then slightly increase with distance from the mouth. The decreases in the magnitudes of Q N and Q S can be attributed to the fact that the river bed was gradually elevated toward the shallow section. Close to the bifurcation, the Stokes transport Q S becomes seaward and changes relatively little over a spring-neap cycle. The Eulerian mean discharge becomes predominantly landward in this region and is more variable in time. Thereby, much of the variation in Q t , being the sum of the subtidal return discharge and the Stokes transport, can be attributed to variation in the Eulerian mean discharge. For the South Branch, the Stokes transport creates an influx that rapidly decreases from the channel mouth inland. The return discharge shows the opposite behavior, such that the total discharge remains constant along the South Channel.

Subtidal Momentum Balance
To gain further insight in the causes of residual flow, the subtidal momentum balance is investigated, following Buschman et al. (2009): The four terms in this equation indicate the influence of temporal and spatial accelerations, pressure gradients and friction in a diurnal tidal cycle, successively. To quantify the terms in equation ((5)), the crosssectional averaged values of U and Q are calculated from the 2-D model results based on 100 cross sections. Figure 9 shows the variation of the four terms over a complete MSf period along the North Branch and the South Branch. The shaded areas indicate the maximum and the minimum magnitudes of the four terms over aspring-neap cycle, whereas the solid line indicates the average. As expected, the advection and friction terms balance the pressure gradient term, whereas the contribution of the temporal variation term in the subtidal momentum balance is minor.  As a result of limited depth variation, the three dominant terms change little along the South Branch, except for a 50-km-long section close to the coast. From the mouth to the shallow section in the North Branch, the terms reduce in magnitude and then change sign repeatedly. The magnitude of the subtidal friction term is strongly dependent on the tidal range and river flow velocity (Figures 10a and 10b), which is in turn altered by width and depth along the North Branch.
Focusing on the section landward of the crest of the bar in the North Branch, the spatial variation of the terms in the subtidal momentum balance is subtle and complex. In this section, the pressure gradient term is dominant and largely balanced by subtidal friction. The advection term, which accounts for spatial accelerations caused by both depth and width variation, is smaller than the pressure term in this region. It is remarkable that despite the idealized, smoothed bed and width profiles, the contributions of the dominant terms switch sign repeatedly. This demonstrates that an Eulerian mean current represents only a very small longitudinal section of the channel. This limits the degree in which a subtidal momentum balance that is set up based on field measurements from one particular cross section informs about the driving mechanism. We attribute this complexity to the tidal regime change across the shallowest part of the channel. Figure 10c shows the mean surface elevation profiles along the North and South Branches averaged over a spring-neap cycle. The subtidal surface water level in the upper reach of North Branch is higher than that in the corresponding reach of the South Branch. The bulge of water in the North Branch thus reverses the subtidal water level gradient in the upper reach of the North Branch, reducing the allocation of river discharge to the North Branch. The tide-averaged surface elevation in the South Branch increases dramatically near the river outlet and then keeps a gentle slope going landward. Although the profiles in the North Branch also increase in the landward direction, the shape is more complex due to the depth and width variation in this channel. There convex mean water surface profiles in the upper reach of the North Branch explain variation in the pressure term shown in Figure 9. In the section landward of the crest of the bar, the predominant adverse subtidal water surface gradient leads to the reversal of tide-averaged flow (Figure 8). Close to the bifurcation, a weak seaward subtidal water level gradient is locally present, consistent with what is illustrated in Figure 10. These results show it is the section-averaged adverse pressure gradient that drives the flow reversal.

Fortnightly Variation of Discharge and Surface Elevation in the North Branch
The subtidal water level is highest at the location where the North Branch is shallowest (Figure 10, black dotted line), due to the water level setup resulting from Stokes transport convergence from opposite sides in the channel (Figure 9). Figure 11a shows the temporal variation of the subtidal water elevation at this location in the North Branch over a springneap cycle. Moreover, Figures 11b and 11c show the subtidal discharge at the upstream and downstream ends of the North Branch, as well as the tide-averaged water volume in the entire branch. The subtidal discharge at the upstream side features only a weak spring-neap oscillation compared to the downstream end, where the river discharge is small compared to the Stokes transport. Due to the long distance as well as the depth and width variation of the North Branch, a remarkable magnitude and phase difference can be found between the upstream and downstream magnitudes of the subtidal discharge. This implies the MSf discharge cannot be regarded constant over a delta branch, neither in phase nor in magnitude. The subtidal water volume that accounts for the total volume of water in the entire North Branch, averaged over a tidal cycle, lags 3 days behind the subtidal discharge at the downstream section. The subtidal water volume and the maximum subtidal water elevation cooscillate (Figure 11), in phase with the tidal amplitude. Figure 12 shows the subtidal surface level variation as a function of the subtidal discharge, revealing hysteresis over a spring-neap tidal cycle. During spring tide, the setup created by Stokes transport reverses the subtidal discharge, creating the highest discharges to the South Branch right after the peak in water level setup that corresponds to spring tide.

Discussion
Deepening of navigation channels often results in sediment accumulation inside channels, because sediment imported during low flows is not anymore all flushed out during high flows (cf. Hoitink et al., 2017). The idealized model results presented herein demonstrate that the shallowest section in a partially abandoned distributary channel has a strong control on the flow division at the bifurcation and that the subtidal flow can reverse direction to become landward (Zhang et al., 2017). The bed elevation acts like a tidal divide, splitting the side channel into upstream and downstream parts. Large subtidal storage variation in the side channel can be characterized as a bulge of water filling the side channel during spring tide, reducing in size during neap tide when the channel empties. Often, fortnightly variation of water levels is attributed to interaction between river discharge and the tidal motion through bottom friction (Hoitink & Jay, 2016), captured in an MSf harmonic (e.g., Gallo & Vinzon, 2005;Jay et al., 2015). Our results show that convergence of Stokes transport can be a second source of fortnightly water level variation, which may occur during low river discharge. The obtained results demonstrate that the extremely shallow section in the North Branch of the Yangtze Delta plays a key role in the fortnightly variation of storage in this branch. Sediment accretion in delta channels may not be a necessary condition for the convergence of Stokes transport. Any circumstance in which an along-channel tidal phase development switches from increasing to decreasing, or vice versa, may potentially cause reversed subtidal flow.
Model results obtained by Xue et al. (2009) reveal that the saltwater flows into the South Branch of the Yangtze Delta, forming an isolated mass of saltwater ( Figure 13). This phenomenon, which has large implications for freshwater availability, is referred to as seawater spilling from the North Branch. Although the driving forces that control the saltwater intrusion in the Yangtze Delta are complex, the mechanism can be partly explained by applying this barotropic model. The fact that seawater intrudes over the elevated channel bed is counterintuitive, because typically, channel deepening causes sea intrusion rather than the opposite. The idealized model helps to explain the seawater intrusion in the Yangtze Delta during the dry season, which is essentially an effect of the tidal rectification (Xue et al., 2009). The local channel bed elevation in the North Branch hampers tidal motion exchange between the two branches of the Yangtze Delta. Seawater accumulates over the shallow section during spring tide because of Stokes transport convergence, which is partly issued to the South Branch during neap tide. Figure 6 shows that for the most realistic depth configuration a substantial increase of the minimum depth in the North Branch is required to avoid seawater spilling, which corresponds to values of Ψ below unity. To avoid reversal of the subtidal discharge throughout the entire spring-neap cycle, the depth ratio d 2 /d 1 would have to increase above 0.4. The case of the Yangtze Delta can be viewed as an extreme case, revealing processes that may occur based on the rapid accretion continues in a delta channel without counter measures. In a long-term stage, the river bed at a certain section in such a channel was gradually elevated, which leads to the significant dissipation of tidal motion. The resulting intake of seawater into the interior of the delta poses a contemporary threat to deltas subject to anthropogenic pressure.
Width convergence is a hallmark for distributary channels in tideinfluenced deltas (e.g., Hoitink et al., 2017). In the example of the Yangtze Delta, the side channel exhibits the most pronounced funnel shape, reflecting tidal dominance in the North Branch and a fluvialdominated geometry of the South Branch. In branching delta channels, the rates of width convergence are rarely identical for parallel branches. Asymmetry in distributary channels can be found in the Pearl Delta (Zhang et al., 2013), the Mekong Delta (Gugliotta et al., 2017), and various deltas in Indonesia (Buschman et al., 2013;Kästner et al., 2017;Sassi et al., 2013). The deeper distributaries conveying the largest share of the river discharge generally show to be less width convergent, but a heterogeneous geological setting may change this general pattern. By focusing on a channel junction in the Berau Delta, Buschman et al. (2013) analyzed the mass exchange between a relatively deep, width convergent channel and a

Water Resources Research
shallower channel of nearly constant width, which conveys a larger share of the discharge. In the shallower channel, reversed salinity gradients were found, caused by seawater transport from the deeper channel. Buschman and coworkers established substantial differences between Stokes transport in each of three cross sections connected by the junction, both during spring tide and during neap tide. Our results inspired on the Yangtze Delta show that Stokes transport and the subtidal residual discharge can vary strongly along branches in between two nodes of a tidal channel network, caused by channel storage variation over a spring-neap cycle. As a consequence, understanding the subtidal exchange processes of mass or seawater at a channel junction requires analyzing the water surface dynamics at a regional scale.

Conclusions
Results from a reduced complexity model demonstrate the tidal regime change associated with sediment accretion and width convergence in a delta channel, revealing a new mechanism which has implications for seawater intrusion. Depth reduction of a secondary channel results in the development of a tidal divide, which occurs at the shallowest section in this channel. The tide propagates into the secondary channel from both sides and converts into a standing wave at the shallowest section that is subject to sediment accretion. The associated Stokes transport converges at the tidal divide from opposite sides of the channel, forming a bulge of water that reverses the water surface gradient during spring tide and low river discharges. As a consequence, the division of river discharge over the two main branches in the model varies over the springneap cycle. From spring tide to neap tide, the Stokes transport relaxes, and the adverse tide-averaged water level gradient landward from the tidal divide reverses the tide-average flow, explaining spilling of seawater into the main branch. A substantial reduction of the difference in depths between the two branches may be needed to avoid this seawater intrusion induced by Stokes transport convergence. The fortnightly variation of water levels simulated with the reduced complexity model reveal that MSf tidal variation may not per se be the result of interaction of river discharge with the tidal motion. The reduced complexity model developed here captures the large-scale delta channel behavior, neglecting the details of smaller-scale processes. In particular, the model helps to anticipate the consequences of sediment accumulation for tidal dynamics in delta channel networks.
Notation d Depth from bottom to water surface level (m).
h Vertical distance from bottom to mean sea surface level (m). ζ Vertical water surface level fluctuation from mean sea level due to tides (m). s Coordinate along channel, positive seaward (m). L Length along the s coordinate from river inflow (m).

Water Resources Research
L W e-folding length for channel width (m). W Channel width (m). Q Total discharge through cross sectional area (m 3 s −1 ). Q S Stokes induced mass flux due to covariation of flow velocity and water level variation (m 3 s −1 ). Q N Net return discharge, compensating for Q S (m 3 s −1 ). U Cross-sectional averaged velocity (m s −1 ). g Gravitational acceleration (m s −2 ). C Chézy coefficient (m 1/2 s −1 ). Ψ Discharge asymmetry index (). G slope in the model from 0 to 500 km. R tidal range (m). 2016B10414), Jiangsu Key Laboratory of Coast Ocean Resources Development and Environment Security (Hohai University) (Project JSCE201509), and by Jiangsu provincial Six Talent Peaks project (XXRJ-008). A. J. F. Hoitink was funded by the Netherlands Organisation for Scientific Research (NWO), within Vici project "Deltas out of shape: regime changes of sediment dynamics in tide-influenced deltas" (Grant NWO-TTW 17062). All data shown in figures in this manuscript are available as supporting information.