Abstract

[1] A dynamic model for the morphological evolution of channels and unvegetated tidal flats is proposed. Channels and tidal flats are schematized as two reservoirs that exchange sediment through the tidal dispersion mechanism, which stems from the presence of a tidal exchange flow and spatial gradients in sediment concentration. The reference concentration in each reservoir is determined by the shear stress associated to tidal currents and surface wind waves, which are a function of the geometry of the system. A simplified procedure to compute flow partition between channels and tidal flats is developed and compared to the numerical solution of the shallow water equations, showing good agreement. In the absence of wind waves, tidal flats reach a stable dynamic vertical equilibrium close to mean high water level, resembling a creek-marsh morphology. For intermediate wind conditions, an additional stable dynamic vertical equilibrium, characterized by a channel flanked by tidal flats close to mean low water, arises. Such equilibrium stems from a sediment exchange dynamic balance between current-dominated channels and wave-dominated tidal flats, and it likely represents the morphological configuration of most tidal flats. Waves associated with intense winds suppress the channelization process. The model suggests that tidal flat elevation is primarily controlled by waves and can be decoupled from channels. Channel depth is also indirectly controlled by waves, through the influence of tidal flat elevation on channel hydrodynamics. Finally, the model predicts that variations in environmental parameters, such as sea level, storminess, and sediment availability, can induce catastrophic morphological shifts.

1 Introduction

[2] Despite the extensive knowledge of hydrodynamic and sedimentary processes in coastal tidal environments [e.g., Parker, 1991; Whitehouse et al., 2000], predicting long-term morphological changes remains a challenging task. This limitation stems from the difficulty of identifying and quantifying the relevant morphological processes and their interactions. Additional complications arise when different landforms are produced by different physical processes. For example, tidal channels are generally incised by currents while tidal flats are the product of waves. As a result, most of the morphological models of intertidal landscapes focus either on channels or tidal flats, neglecting the existing feedbacks.

[3] Tidal channels are generally studied with a one-dimensional schematization [Friedrichs, 1995; Lanzoni and Seminara, 2002; Seminara et al., 2009; Toffolon and Lanzoni, 2010]. In these frameworks, tidal flats are considered as water storage areas, without participating in the momentum dynamics. This storage effect promotes ebb dominance and hence a net seaward transport of sediments [Speer and Aubrey, 1985; Friedrichs and Aubrey, 1988]. The channel morphological evolution is based on channel hydrodynamics, computed assuming that the tidal flat is a fixed element that does not evolve in time.

[4] On the other hand, tidal flats are often schematized as a transect parallel to the main flow direction [Le Hir et al., 2000; Roberts et al, 2000; Pritchard et al, 2002; Pritchard and Hogg, 2003; Mariotti and Fagherazzi, 2010], or as a flat surface characterized by a single elevation [Fagherazzi et al., 2006; Fagherazzi et al. 2007; Marani et al., 2007; Marani et al., 2010]. In both one-point and transect models, the tidal flat is coupled with the rest of the system through a sediment boundary condition, and the draining effect of channels and the consequent reduction of flow on the tidal flat are neglected. This boundary condition represents the dynamics of the rest of the system, such as sediment inputs from channels adjacent to the tidal flat. Because this boundary condition is set a priori, the coupling between tidal flat and channel hydrodynamics is neglected.

[5] A number of observations have been carried out in both channels and tidal flats, shedding light on the main sedimentary mechanisms at play in these environments. Because of the high current velocity, tidal channels are often characterized by a suspended sediment concentration higher than nearby tidal flats during fair weather [Ridderinkhof et al., 2000; Janssen-Stelder, 2000]. During storms, the reduced water depths on tidal flats allow a strong penetration of the wave motion, resulting in higher bottom shear stresses and substrate remobilization [Dyer et al., 2000; Christiansen et al., 2006]. Because of these dynamics, sediments are preferentially transported from channels to tidal flats during fair weather and from tidal flats to channels during storms [Yang et al., 2003; Allen and Duffy, 1998; Boldt et al. 2012].

[6] In order to quantify the sediment exchange between channels and tidal flats, three processes should be taken into account: (1) sediment resuspension in the channels [e.g., Lanzoni and Seminara, 2002], (2) wave resuspension on the tidal flats [e.g., Fagherazzi and Wiberg, 2009], and (3) lateral circulation between the two [e.g., Li and O'Donnell, 1997; Mariotti and Fagherazzi, 2012b]. In addition, sediment exchange with the external coastal environment should be included (e.g., through a closure section that might coincide with a tidal inlet). Indeed, Tambroni and Seminara [2006] concluded that the flux of fine material exiting the Venice Lagoon during stormy weather is a major component of its sediment budget.

[7]Tambroni and Seminara [2012] recently proposed a model for the morphological evolution of a channel-marsh transect that includes many of the abovementioned processes. The model bypasses the problem of imposing an arbitrary sediment concentration at the channel-marsh boundary, includes both currents and wave resuspension, and allows vertical and horizontal marsh migration. However, the model does not include tidal flats adjacent to the channel, which often constitute the bulk of the intertidal area. As a result, this model neglects both the flow partition and the sediment exchange between channels and tidal flats.

[8] The goal of this paper is (1) to explore the morphological coupling between channels and tidal flats driven by the interplay between currents and waves, and (2) to expand the capability to predict how intertidal systems respond to changes in external factors such as sea level rise [Bassett et al., 2005; Sallenger et al., 2012], decrease in sediment supply [Syvitski et al., 2005], or decrease in wind speed [Pryor et al., 2009].

[9] We here develop a “two-point” dynamic model that features a novel procedure to compute flow partition between channels and tidal flats. As with all models, the validity of our results is limited by the number of processes considered and by the simplifications adopted to describe them. In particular, the point approach cannot recreate two- and three-dimensional features typical of the intertidal landscape.

[10] The model is extensively explored as a dynamic system [Phillips, 1992], showing the presence of multiple equilibria and bifurcations associated to changes in the control parameters. The model is then applied to a mudflat in Willapa Bay (Washington State, USA), showing a reasonable agreement in terms of morphological results.

2 Study Site

[11] Willapa Bay, Washington State, USA, is an embayment with mixed-semidiurnal tides having a mean tidal range of 2.72 m [Mariotti and Fagherazzi, 2012b]. The southern part of the bay is characterized by negligible salt marshes and extensive mudflats, incised by tidal channels of a wide range of dimensions [Banas et al., 2004]. The main channel draining the eastern portion of the mudflat has a width of few hundred meters and a depth ranging from few meters to tens of meters; the tidal flat is approximately flat, slightly sloping upward and concave up close to the landward boundaries [Mariotti and Fagherazzi, 2012a, 2012b]. The tidal flat is placed about 2 m below the mean high water level (MHW), and is hence exposed at low tide [Mariotti and Fagherazzi, 2012b].

[12] The mudflat substrate has a high concentration of silt and clay sediments, characterized by flocs with an average size of 0.18 µm [Law et al., 2013]. A common statistic used to characterize storm intensity is the 90th percentile of the wind speed [Pryor et al., 2009]. In Willapa Bay, this reference wind speed is 7.7 m/s, derived from measurements between 2005 and 2011 at the NOAA station Toke Point, placed near the Bay inlet, approximately 20 km north of the tidal flat considered herein.

3 The Channel Tidal-Flat Model

[13] The model considers channels and tidal flats as two sediment reservoirs whose volume varies in time. We consider a schematic basin, with a single rectangular channel bounded by a rectangular tidal flat at the sides, each connected to the rest of the system at a closure section, analogous to the schematization used by Li and O'Donnell [1997] (Figure 1). The total width of the channel, 2bc, the width of a single tidal flat, bf, and their lengths, both equal to L, are kept fixed. The depths of the channel at the closure section, dc, and the depth of the tidal flat, df, are the only dynamic variables of the system. The channel profile is not specified, but it is assumed to only translate vertically during the evolution. As a result, dc and df are univocal indicators of the amount of sediment present in the channel and tidal flat reservoirs. The vertical datum is fixed at MHW and the depths are considered positive below this datum. A single-frequency tide with tidal range r is considered.

Figure 1.

(a) Scheme of the control volume considered. Two sediment reservoirs are selected: channel and adjacent tidal flats. For symmetry, only half domain is considered. The basin length, L, and the channel and tidal flat widths, bc and bf, are constant parameters. The dynamic variables are the average tidal flat depth, df, and the channel depth at the closure section, df. Qc and Qf are the tidally averaged volume-exchange rate [L3T−1] between the sediment reservoirs and the rest of the system through the closure section, Qlat is the tidally averaged volume-exchange rate between the two reservoirs. (b) Image of a channel and tidal flat of Willapa Bay (WA), where the model is applied. The tidal flat depth is ~2 m and the channel depth at the closure section is ~6 m (depth referred to MHW). (c) Schematization of the channel and tidal flat cross section and longitudinal profiles. The tidal volumes for different tidal stages, defined in equations (4) and (5), are highlighted in the cross section.

[14] We assume that the sediment in the basin is composed by very fine material (e.g., D50 = 18 µm), which is prevalently transported in suspension. Because of the elevated mobility of the sediment, we assume that the reference concentrations in the channel, Cc, and on the tidal flat, Cf, are uniform in space.

[15] Because of the system symmetry, we isolate a control volume composed of half channel and one lateral tidal flat (Figure 1a). The channel and tidal flat reservoirs exchange sediment between them and with the rest of the system through the seaward closure section. The concentration at the closure section, Co, is a control parameter. We assume that the tidally averaged transport is dominated by tidal dispersion [Okubo, 1973; Dronkers, 1988; Di Silvio et al., 2010]. Tidal dispersion is driven by the presence of an exchange tidal flow, symmetrical with respect to ebb and flood, and spatial gradients in sediment concentration [Friedrichs, 2012]. In our simplified model, these gradients result from the different sediment concentrations in the channel, on the tidal flat, and at the boundary of our reference volume. Residual transport, i.e., ebb-flood flow asymmetries, is neglected.

[16] We schematize the flow exchange in three parts: flow through the tidal flat closure section, flow through the channel closure section, and flow between tidal flat and channel (Figure 1c). The tidally averaged sediment mass balance for the channel and tidal flat reservoirs reads:

ddcdt=QlatCc−Cf+QcCc−Co/bcLρs+Rddfdt=QlatCf−Cc+QfCf−Co/bfLρs+R(1)

where Qc [L3T−1] is the tidally averaged volume-exchange rate through the channel at the closure section, Qf [L3T−1] is the tidally averaged volume-exchange rate through the tidal flat at the closure section, Qlat [L3T−1] is the tidally averaged volume-exchange rate between the channel and the tidal flat, ρs is the sediment bulk density, and R is the relative sea level rise, which is simulated by keeping the tidal datum fixed and increasing uniformly the channel and tidal flat depths. A uniform sedimentation rate (per unit of area) due to sediment inputs from rivers discharging in the basin can be easily introduced by manipulating the relative sea level rise, but it will not be considered in this work. Here we do not include salt marshes that often border tidal flats, and the basin length and width B = (bc + bf) are kept constant in time. This model is hence suitable for intertidal environments in which extensive salt marshes do not coexist with tidal flats, which is the case of Willapa Bay (Figure 1b) and other estuaries in the Pacific Northwest [Emmett et al., 2000]. In future research, we will relax this hypothesis to study the spatial evolution of the domain, driven for example by the erosion of salt marshes bordering the tidal flats.

[17] We assume that the reference concentrations in the channel, subscript c, and in the tidal flat, subscript f, are such that erosion and deposition are identical. Assuming that the deposition rate is equal to wsCc/f and that the erosion rate is me max[(τc/f − τcr)/τcr, 0] we obtain [Whitehouse et al., 2000]:

Cc/f=maxτc/f−τcrτcr0mews(2)

where ws is the reference settling velocity, me is the erosion coefficient, τc/f is the reference shear stress in the channel or on the tidal flat, and τcr is the critical shear stress for erosion. We assume that waves and currents act independently to determine the reference sediment concentration in both channel and tidal flat:

Cc/f=maxτcur,c/f−τcrτcr0mews+λmaxτw,c/f−τcrτcr0mews(3)

where τcur,c/f is the reference current-induced shear stress in the channel or tidal flat, τw,c/f is the reference wave-induced shear stress in the channel or tidal flat, and λ is a weighting factor that describes the relative frequency of wave resuspension events with respect to currents. Instead of determining λ, here we set this parameter equal to 1, and we select a significant wind speed that is equivalent to the average effect of all wind events. The intensity of such wind event is then the control parameter for the wave processes. An alternative approach would be to fix a wind speed, e.g., 10 m/s, and consider λ as the control parameter.

[18] The model structure is briefly illustrated in Figure 2. The hydrodynamics and sediment transport in our basin are simplified in the following way: (1) the fluxes of water are calculated from the tidal prism; (2) the tidal prism is partitioned between tidal flat and channel as a function of water depth, following a simplification of the momentum equation; (3) from the water fluxes, a characteristic tidal current velocity is obtained; this velocity is used to calculate current bottom shear stresses; (4) wave bottom shear stresses are calculated as a function of wind speed, fetch, and depth; (5) sediment fluxes result from the combination of water fluxes and reference suspended sediment concentrations, calculated as a function of wave and current shear stresses; (6) sediment fluxes determine erosion and deposition at the bottom of the channel and tidal flat and therefore control the vertical morphological evolution.

Figure 2.

Model flowchart. Given the channel and tidal flat geometry, a hydrodynamic model computes the flow partition between channel and tidal flat and related current shear stresses, while a wave model computes wave shear stresses. The reference sediment concentrations are computed from the shear stresses, and they are combined with the water fluxes and the open-sea concentration to obtain the sediment fluxes, which modify the system morphology.

3.1 Hydrodynamic Model for Flow Redistribution

[19] Here we propose a simplified model for the flow in the channel and tidal flat. We assume that the hydrodynamic is tidally dominated, neglecting the presence of wind driven currents, seiches, riverine circulation, and stratification. Restricting the analysis to tidal flats with width, bf, much smaller than their length, L, we assume that the water level is uniform in the cross-channel direction [Mariotti and Fagherazzi, 2012b]. We also consider short tidal flats only (e.g., L << tidal wave length), for which the tide is quasi-static, so that, at the zero-order approximation, the water level is uniform in the along-channel direction [Fagherazzi et al., 2003, Toffolon and Lanzoni, 2010]. We consider two stages: high tide (H), when the tidal flat is submerged, and low tide (L), when the tidal flat is emergent (Figure 1c). Stage H considers water levels between MHW and the tidal flat elevation, while stage L considers water levels between the tidal flat elevation and mean low water level (MLW). The tidal prisms for the channel, PH,c , and tidal flat PH,f , for stage H are as follows:

PH,f=mindfrLbfPH,c=mindfrLbc(4)

[20] According to this equation, the tidal prism changes when the tidal flat elevation is above or below MLW. The total tidal prism for stage H, PH, is the sum of PH,f and PH,c. In addition, the channel tidal prism for stage L is:

PL=minr−minrdf,dc−dfLbc(5)

[21] Because the two-dimensional flow is undetermined at the zero-order approximation, a first-order approximation is required to calculate the flow redistribution between the channel and the tidal flat during stage H. At the first order, the inertial and time derivative terms are neglected [Le Blond, 1978; Friedrichs and Madsen, 1992; Toffolon and Lanzoni, 2010], and the flow in the longitudinal direction at the closure section is in equilibrium with bottom friction, according to the Manning equation:

U=1nh2/3i1/2(6)

where n is the Manning's coefficient, i is the longitudinal barotropic gradient at the closure section, and h is the effective water depth. This depth takes into account that the water level varies during the tidal cycle, and it is set equal to the average between the minimum and maximum water depth in both units (channel and tidal flat):

hc/f=dc/f+max0,dc/f−r/2(7)

[22] If the unit is exposed at low tide, the effective depth is equal to half the high tide depth; if the unit is not exposed, the effective depth is equal to the high tide depth minus the tidal amplitude. Because of the assumption of uniform water level along each cross section, the redistribution of the longitudinal discharge between channel and tidal flat at the closure section is proportional to h5/3 [Mariotti and Fagherazzi, 2012b]. The total volume of water exiting the channel, Vc, and the tidal flat, Vf, in the longitudinal direction becomes

Vc=PL+PHbchc5/3bchc5/3+bfhf5/3Vf=PHbfhf5/3bchc5/3+bfhf5/3(8)

[23] According to this simple model, the partition of the water fluxes between channel and tidal flat increases when the difference in the elevation between the two landforms increases. This partition depends on the ratio between the channel and tidal flat width, but not on their absolute value [Mariotti and Fagherazzi, 2012b]. The volume exchanged between channel and tidal flat is computed by continuity, as the residual between the tidal prism above the tidal flat and the volume exiting the tidal flat along the longitudinal direction [Mariotti and Fagherazzi, 2012b]:

Vlat=PH,f−Vf(9)

[24] The volumes Vc, Vf, and Vlat are divided by the tidal period Tω to obtain the tidally averaged volume-exchange rates Qc, Qf, and Qlat used in equation (1). Noticeably, the proposed equations do not depend on the channel longitudinal profile, but only on its depth at the closure section. This condition is met if the channel depth is lower than MLW at every location, so that the tidal prism can be computed using equation (4).

3.2 Current Shear Stresses

[25] The next step required by the zero-dimensional approximation is the introduction of a reference shear stress, such that its effect, in terms of average sediment resuspension, is equivalent to a spatially distributed shear stress. The choice of a single shear stress for long-term modeling purposes is supported by several studies showing that both channels [Friedrichs, 1995] and tidal flats [Friedrichs and Aubrey, 1996; Pritchard et al., 2002] adjust their profile toward a spatially uniform shear stress.

[26] Here we choose the channel and tidal flat maximum tidal velocity at the closure section as reference for the current shear stress. Such velocity is computed by dividing the volume of water per unit of width flowing in the longitudinal direction by the reference depth and by an effective tidal time Te:

Uc/f=Vc/fbc/fhc/fTe,c/fπ2(10)

[27]Te equals the amount of time during which water is flowing (i.e., the bed is not dry) during ebb or flood in a tidal cycle, i.e.,

Te,c/f=Tω2minrdc/fr(11)

[28] The factor π/2 in equation (10) originates from the ratio between the maximum and average water level temporal derivatives, and it is used to evaluate the maximum rather than the tidally averaged velocity. The current-induced shear stress is then computed from the velocity using the Manning equation:

τcur,c/f=γn2hc/f1/3Uc/f2(12)

where γ is the water specific weight.

3.3 Wave Shear Stresses

[29] Locally generated wind waves are reproduced with a simplified method. Significant wave height, Hs, and peak period, Tp, are computed using the semi-empirical relationships of Young and Veraghen [1996]:

[30] Mariotti and Fagherazzi [2012a] showed that the presence of channels marginally affects the wave regime on tidal flats. Therefore, the tidal flat effective depth, hf, is taken as a reference to compute the wave regime (Hs, Tp) over the entire domain. Fetch and wind speed are considered control parameters. The wave-induced shear stress τw is then computed with the linear wave theory as follows:

τw,c/f=12ρfwπHsTpsinhkhc/f2(14)

where fw is a friction factor, and k is the wave number, derived from the dispersion relationship. Because of the greater water depth, wave shear stresses are lower in the channel than on the tidal flat. These semi-empirical equations [(13) and (14)] produce a non-monotonic relationship between bed shear stress and water depth [Fagherazzi et al., 2006], which is function of both wind speed and fetch [Mariotti and Fagherazzi, 2012a]. For a fixed fetch and wind speed, we refer to the maximum wave shear stress over all the possible depths as τw,max.

4 Comparison With the Results of the Model Delft3D

[31] The simplified model for water fluxes and velocities is tested against the results obtained by solving the depth-averaged shallow water equations with the finite-difference model Delft3D [Lesser et al., 2004], which is equipped with appropriate algorithms to deal with wetting and drying processes [Stelling et al., 1986].

[32] Simulations are performed on the simplified geometry of Figure 1a, by considering a single domain length, L = 10 km, a single total width, B = 1.2 km, and by changing the tidal flat depth, df, channel depth, dc, and the ratio between channel and tidal flat width, bc/B. Two channel profiles are considered: horizontal, with a depth equal to dc, and linearly sloping, with a depth equal to dc at the closure section and equal to the minimum between r and df at the landward edge. A sinusoidal water level is imposed at the closure section, and no flux conditions are imposed on the other boundaries. Friction is computed using the Manning equation, with a coefficient n = 0.016 m1/3/s, which is equivalent to a drag coefficient of ~0.0025, often used in muddy environments [Whitehouse et al., 2000]. Horizontal eddy viscosity is computed with a large eddy simulation technique.

[33] The ebb-flood water volumes exchanged along the longitudinal and transversal directions are computed from the results of the Delft3D model as follows:

where Ux, Uy, and h and h are the velocities and the water depth computed with the model. The maximum tidal velocities are calculated as the average of the maximum ebb and flood velocity, Ux, at the closure section, x = L.

[34] The results for the horizontal and linearly sloping channel are almost identical (Figure 3), suggesting that, while the channel profile affects the velocity distribution along the channel, it does not significantly change the integral flow partition between channels and tidal flats. The results of Delft3D and the simplified model are in good agreement (Figure 3). Both models predict that the maximum lateral exchange generally occurs when the tidal flat elevation is at MLW (df/r = 1 in Figure 3a). When the tidal flat elevation is above MLW, the tidal prism and hence the total amount of water that moves through the system is smaller than the possible maximum value [equation (4)]. When the tidal flat elevation is below MLW, the tidal prism remains constant, equal to the maximum value, but the difference in water depth between channel and tidal flat and hence the flow partition is reduced [equation (8)].

[35] Although the simplified model is not quantitatively accurate, it qualitatively captures the velocities simulated by Delft3D (Figure 3b). In particular, the model reproduces the velocity increase in the channel and the velocity decrease on the tidal flat when the water depth in the tidal flat decreases. The velocity increase in the channel is caused by the tidal flat mass storage effect, which has been extensively investigated [i.e., Dronkers, 1986]. The decrease of velocity on the tidal flat is instead caused by the channel drainage effect. This process is not captured by one-dimensional models of tidal flat evolution [Le Hir et al., 2000; Roberts et al, 2000; Pritchard et al, 2002; Pritchard and Hogg, 2003; Mariotti and Fagherazzi, 2010; Tambroni and Seminara, 2012].

[36] Because here we focus on the gross flow partition between channel and tidal flat rather than on the detailed velocity distribution along the channel (needed instead for the morphological evolution of the channel profile, see Lanzoni and Seminara, [2002]; Seminara et al., [2009]), the good agreement with the Delft3D results supports the use of our simplified schematization.

5 Morphological Results

[37] The simplified morphological model is studied fixing the wave friction coefficient, fw, the Manning's coefficient, n, the erosion parameter, me, the settling velocity, ws, the critical shear stress, τcr, the tidal range, r, and the basin length, L, reflecting the physical settings of the southern Willapa Bay mudflat (Table 1). We therefore explore the model sensitivity to the following parameters: bc/B, Co, R, and Uwind. The model capability to predict the Willapa Bay mudflat morphology will be analyzed in the discussion section.

[38] All morphological equilibria predicted by the model are dynamic, i.e., they are characterized by an instantaneous sediment flux generally different from zero, but vanishing if averaged over a certain time scale, e.g., the tidal cycle or the interval between two formative wind events. Hereafter we will omit the term “dynamic” when referring to the equilibria.

5.1 Case Without Relative Sea Level Rise

[39] We first consider the case with no relative sea level rise. The system admits only one equilibrium solution, Cc = Cf = Co, in which the reference concentration is the same in the channel and tidal flats, and it is equal to the seaward boundary concentration. In order to maintain this equilibrium, the shear stress in both channel and tidal flat is such that C(τeq(w),τeq(cur)) = Co:

λwmaxτeqw−τcr,0+λcurmaxτeqcur−τcr,0+τcr=τeq=Coτcrwsme+τcr(16)

[40] The solution further depends on whether we consider only tidal currents or both currents and waves, as discussed below.

5.1.1 Tidal Currents

[41] When only currents are present, an analytical treatment of the solution is possible. If the channel is deeper than the tidal flat, the current induced bed shear stress and the reference concentration are higher in the channel. As a result, the system admits only trivial equilibria, which correspond to Qlat equal to zero. This condition is met in two simple configurations: tidal flat depth equal to zero or tidal flat depth equal to channel depth. The first scenario is a tidal channel in which the surrounding tidal flats have an elevation equal to MHW and do not contribute to the tidal prism (tidal channel case, Pc). The second scenario is a tidal flat without channels (unchannelized sound case, Ps). In the absence of waves, the equilibrium shear stress is associated with the following equilibrium depth:

dcur=r2+ρgn2rLTω/2π22τeq−13/7(17)

which is valid if dcur > r, i.e., the channel bed is lower than MLW.

[42] The equilibrium point Pc (tidal channel) is stable: an increase in channel depth decreases the bed shear stress and the sediment concentration in the channel promoting infilling, while a decrease in channel depth increases its sediment concentration promoting scouring (Figure 4). The equilibrium point Ps (unchannelized sound) is instead unstable. A channel slightly deeper than the surrounding tidal flat starts draining water, thus increasing the current shear stress in the channel and decreasing it on the tidal flat. This condition creates a gradient of sediment concentration from the channel to the tidal flat, which, combined to the water exchange between the two elements (Qlat), determines a flux of sediment toward the tidal flat. This flux tends to fill the tidal flat, increasing the relative channel depth, thus creating a morphological positive feedback that constitutes the basic mechanism for channel formation [Fagherazzi and Furbish, 2001]. Indeed, this process is the channel spillover mechanism described by Mariotti and Fagherazzi [2012b].

[43] The different stability of the two equilibria results by fixing a priori the dimensions of the channel and tidal flat domains. In fact, if the channel bed in the tidal sound configuration is treated as a tidal flat, it will be unstable with respect to the formation of a channel inside it. The model results hence show that a flat surface is always unstable with respect to the channelization process. However, the model is not able to predict the width of the channel resulting from this instability.

[44] The union of the points with channel depth equal to the tidal flat depth constitutes the stable manifold of the saddle point Ps (Figure 4, 1:1 line). An unstable manifold connects Ps to the stable equilibrium Pc. During the trajectory from Ps to Pc, the tidal flat depth decreases monotonically, while the channel initially deepens, reaches a maximum depth when the tidal elevation is at MLW, and then it shoals (Figure 4). The time needed to move from a neighborhood of Ps to a neighborhood of Pc is on the order of 100 years.

[45] Changing the relative channel width bc/B does not affect the equilibrium depth dcur [Figure 4, equation (17)]. However, the maximum channel depth experienced by the manifold connecting Ps to Pc increases with decreasing values of the relative channel width. Changing the model parameters (L, r, and τcr) affects both the equilibrium points [equation (17)] and the manifold, but does not affect the qualitative bifurcation portrait.

5.1.2 Tidal Currents and Waves

[46] Wind waves induce high shear stresses on tidal flats, triggering complex dynamics. We find that the system experiences four distinct regimes when we increase the wind speed: current dominated, mixed, weakly wave dominated, and strongly wave dominated (Figure 5).

[47] The current dominated regime occurs for low wind speeds, such that τw,max < τeq (Figure 5a). In this regime the equilibria are those due to tidal currents only: the tidal channel configuration, Ps, and the tidal sound configuration, Pc [equation (17)]. In these conditions, the higher sediment concentration in the channel tends to fill the tidal flats, despite the fact that τw on the tidal flats is greater than τcr. Some trajectories are slightly different from the case without waves, but the overall phase space portrait remains unchanged.

[48] The mixed regime occurs for wind speeds such that τw,max > τeq (Figure 5b). In this regime a bifurcation occurs and two non-trivial equilibria appear: the wave stable equilibrium Pw,st and the wave unstable equilibrium Pw,un. The new equilibria are such that τcur,c = τw,f = τeq, i.e., the equilibrium concentrations are determined by currents in the channel and by waves on the tidal flat. The presence of two equilibria Pw,st and Pw,un stems from the fact that there are two depths on the tidal flat capable of generating the same wave shear stress, as shown by Fagherazzi et al. [2006] [equations (14) and (15)].

[49] The stability of the equilibrium Pw,st is determined by the shear stress and sediment concentration modulations associated with variations in channel and tidal flat depth. If the tidal flat becomes shallower than df (Pw,st), the sediment concentration on the tidal flat becomes greater than the sediment concentration on the channel, establishing a sediment flux toward the channel, while the opposite flux establishes if the tidal flat becomes deeper than df (Pw,st). Similar stabilizing fluxes are triggered by variations in channel depths.

[50] The equilibrium Pw,un is unstable with respect to variations in tidal flat depth: if the tidal flat becomes shallower than df(Pw,un), the tidal flat sediment concentration becomes smaller than the concentration on the channel and at the closure section, establishing sediment fluxes toward the tidal flat, while the opposite fluxes establish if the tidal flat becomes deeper than df (Pw,un). Because Pw,un is stable with respect to variations in channel depth, this point is a saddle. One limb of the unstable manifold of Pw,un connects the saddle to Pc, while the other limb connects the saddle to Pw,st (Figure 5b). Because the depth associated with Pw,un is very small (<0.2 m for wind speed between 0 and 15 m/s), the basin of attraction of the stable point Pc is small, and the phase space is dominated by the new equilibrium point Pw,st.

[51] The two wave-dominated regimes occur for wind speeds such that τw(dcur) > τcr, i.e., the wave shear stress is capable of resuspending sediment at the current dominated equilibrium dcur (Figures 5c and 5d). In these regimes the equilibrium point Ps is still characterized by channel and tidal flat with same depth, but now both currents and wave shear stresses contribute to τeq, and hence the depth of Ps is no longer equal to dcur [equation (17)]. This equilibrium becomes stable, being controlled by the wave shear stress, which decreases for increasing depths.

[52] In the weakly wave dominated regime, the stable equilibrium Pw,st persists, and τeq remains determined by currents in the channel and by waves on the tidal flat (Figure 5c). An additional equilibrium appears, Pcw, such that the reference concentration on the tidal flat is determined by a combination of currents and waves. This new equilibrium is a saddle: it is stable with respect to variations of tidal flat depth, but it is unstable with respect to variations of channel depth. Because this regime occurs for a very narrow wind speed window, and because Pw,st remains the point with the greatest basin of attraction, the weakly wave dominated regime is hereafter considered equivalent to the mixed regime.

[53] The strongly wave dominated regime occurs for wind speeds such that the current shear stress on the tidal flat associated with Pw,st is greater than τcr (Figure 5d). In this scenario, it is no longer possible to keep both the wave and the current shear stress on the tidal flat equal or below the equilibrium shear stress, hence, Pw,st and Pcw disappear and Ps becomes the only stable equilibrium point.

[54] In summary, the stable point with the largest basin of attraction is Pc for the current dominated regime, Pw,st for the mixed regime, and Ps for the strongly wave dominated regime (see cartoon in Figure 6c). The bifurcation diagram as a function of wind speed is shown in Figure 6a, reporting the coordinate of the stable point with the largest basin of attraction.

Figure 6.

Stable equilibrium with the largest basin of attraction, (Pc for the current dominated, Pw,st for the mixed, and Ps for the strongly wave dominated regime), as function of wind speed, for (a) different Co and R = 0, and (b) different R and Co = 0.1 kg/m3. In both cases r = 2.72 m, L = 4 km, bc/B = 0.2. The two points on the left panel represent the tidal flat depth (~2 m) and the channel depth at the closure section (~6 m) of the control volume considered for the Willapa Bay mudflat (Figure 1b). (c) Scheme of the channel—tidal flat equilibria under various regimes.

[55] Figure 6a also shows the influence of the boundary condition Co on the model results. If more sediment is available at the closure section, then a higher shear stress is needed to balance the input. As a result, a higher Co decreases the equilibrium depth of Ps, Pw,st and Pc.

5.2 Case With Relative Sea Level Rise

[56] When relative sea level rise is present, the system requires a sediment input in order to maintain a morphological equilibrium. In this scenario the trivial equilibrium Ps persists (Figure 7), and stems from a balance between the longitudinal sediment input-export and the internal gain-loss of sediment: Qc(Cc − Co)/bc = Qf(Cf − Co)/bf = − RLρs. Without waves, the depth of this equilibrium reads as follows:

dcur,RSLR=r2+ρgn2rLTω/2π22−TωρsRr+Coτcrwsme+τcr−13/7(18)

[57] When relative sea level rise is present, the tidal channel equilibrium Pc becomes non-trivial, characterized by a deep channel and a tidal flat placed below MHW (Figures 6b and 7). In this equilibrium, Co > Cc > Cf, and hence sediments are transported from the open sea to the channel and from the channel to the tidal flat to balance sea level rise [Friedrichs and Perry, 2001].

6 Discussion

6.1 Is the Model a Good Description of Intertidal Landscapes?

[58] The proposed model investigates the evolution of a channel—tidal flat complex by using some tools of the dynamic systems analysis, such as equilibrium points, stability, and bifurcations, yet featuring geometries and flow patterns more articulated than those assumed in one-point dynamic models [e.g., Fagherazzi et al., 2006; Fagherazzi et al. 2007; Marani et al., 2007; Marani et al., 2010]. The good agreement with the detailed hydrodynamic model Delft3D suggests that these flow patterns are realistically approximated in our framework.

[59] Our simple model reproduces the flow between tidal flats and channels, which constitutes the main hydrodynamic interaction between these two landforms [Fagherazzi and Furbish, 2001]. As a result, the current shear stresses on the tidal flat and the wave shear stresses in the channel marginally participate to the morphological evolution of the system, and the dynamics are overall driven by currents in the channel and waves on the tidal flat, in accordance with field observations [Allen and Duffy, 1998; Ridderinkhof et al., 2000; Mariotti and Fagherazzi, 2012b].

[60] The choice of only two dynamic variables and the schematization of the tidal dispersion fluxes in three components, Qc, Qf and Qlat, allow the construction of a reduced complexity model which admits some analytical treatments. Clearly many processes are neglected in the model, such as spatial variability in channel and tidal flat hydrodynamics [Lanzoni and Seminara, 2002; Roberts et al., 2000], dead-end channels drainage [Mariotti and Fagherazzi, 2011], ebb-flood asymmetry [Speer and Aubrey, 1985; Friedrichs and Aubrey, 1988], very shallow flow erosion [Fagherazzi and Mariotti, 2012], wind induced currents [Tambroni and Seminara, 2012], settling and scouring lag [Postma, 1961], and biostabilization and bioturbation [Le Hir et al., 2007].

[61] The model does not simulate the dynamics of the tidal flat-marsh boundary, the position of which has been shown to be unstable, prograding or retreating [Mariotti and Fagherazzi, 2010; Tambroni and Seminara, 2012]. As a result, the model should not be applied in landscapes where the marsh boundary is actively migrating. Nonetheless, the model can reproduce the equilibrium between tidal flats and channels only, a configuration eventually achieved when marsh erosion is stopped by geological constraints, i.e., there is no salt marsh left to be eroded.

6.2 Equilibrium Without Waves: A Creek—Marsh Configuration?

[62] The model predicts that, with moderate or no waves, an unchannelized configuration is unstable, thus capturing the autogenic channelization process driven by flow concentration in deeper areas [see also Fagherazzi and Furbish, 2001, equation (8)]. Starting from an unchannelized bed, tidal flat infilling increases the flow partition and hence the velocity in the channel. Once the tidal flat becomes shallower than MLW, the tidal prism decreases, reducing the flow velocity in the channel. During this transition, the channel initially deepens, reaches a maximum when the tidal flat is at MLW, and then shoals converging to its equilibrium depth. Eventually, the tidal flat reaches MHW and does not contribute to the fluxes in the channel. In the presence of relative sea level rise, a similar equilibrium occurs, but with tidal flats slightly below MHW (df <0.5 m for RLSR of 0–10 mm/yr). In this scenario, a net sediment flux from the channel to the tidal flats is established to balance relative sea level rise.

[63] This equilibrium is morphologically similar to a salt marsh, an intertidal environment characterized by a nearly horizontal vegetated platform placed in proximity of MHW and dissected by creeks, although our model does not include vegetation effects, such as organogenic sediment production, decrease in bed turbulence, and sediment stabilization [Fagherazzi et al., 2012]. It is therefore of interest to compare our results to models for creeks-salt marsh dynamics, in which waves are neglected but vegetation is accounted for.

[64] Spatially explicit models for salt marsh evolution without relative sea level rise [D'Alpaos et al., 2006; D'Alpaos et al., 2007] show that the equilibrium configuration is characterized by creeks bounded by salt marshes placed at MHW. Such equilibrium is reached even with a sediment boundary condition greater than zero, in accordance with our predictions. Similar to our model, the creek-salt marsh cross section model proposed by D'Alpaos et al. [2006] predicts a transient channel evolution, caused by the changes in water storage on the shoals. Finally, the equilibrium with the tidal flat just below MHW and dissected by channels qualitatively agrees with the results of Kirwan and Murray [2007], who proposed a 2D model for tidal creeks and marsh platform with relative sea level rise and a boundary sediment concentration greater than zero.

[65] We expect that the introduction of vegetation will accelerate the final stage of the tidal flat accretion, leading to a similar creek-marsh equilibrium. It is important to note that the creek-marsh configuration has no tidal flats, and hence the problem of the marsh boundary migration does not occur. Clearly our model is not able to realistically reproduce the transient phase in which tidal flats and marshes coexist [Mariotti and Fagherazzi, 2010; Tambroni and Seminara, 2012].

6.3 Equilibria With Waves: Predictions and Comparison With Willapa Bay

[66] In the mixed regime, waves trigger the presence of a non-trivial equilibrium, Pw,st, with a deep channel and a tidal flat placed near MLW. This equilibrium occurs for wind speeds such that the maximum wave shear stress is greater than τeq, which is higher than τcr. In other words, the wave regime needs to be strong enough to overcome the channel dominance. In this equilibrium, the reference concentration is maintained by currents in the channel and by waves on the tidal flat. This corresponds to the observation that sediment concentrations are higher in the channel during fair weather [Ridderinkhof et al., 2000; Janssen-Stelder, 2000], and on tidal flats during storms [Dyer et al., 2000; Christiansen et al., 2006].

[67] The proposed model predicts the presence of a stable dynamic equilibrium between a current-dominated channel and a wave-dominated tidal flat. In accordance with models for the evolution of tidal flat profiles [Roberts et al., 2000; Pritchard et al., 2002] and with one-point models [Marani et al., 2007], we predict that the tidal flat deepens with increasing wave intensity and decreasing sediment concentration at the sea boundary. Hence, variations in wave regime and sediment availability can trigger a switch in the tidal flats equilibrium: from intertidal, i.e., above MLW, to subtidal, i.e., below MLW (see Figures 6a and 6b). Such switch has a profound importance, given the number of physical and biological processes affected by tidal flats exposure to air at low water, such as sediment dewatering [Anderson and Howell, 1984], and benthic microorganisms migration and primary production [Consalvey et al., 2004; Spilmont et al., 2007].

[68] The channel-tidal flat equilibrium is similar to the morphological configuration of the inner tidal flat of Willapa Bay. Modeling comparison with the Willapa Bay morphology is performed by considering the control volume in Figure 1b, which features a tidal depth of about ~2 m [Mariotti and Fagherazzi, 2012b] and a channel depth at the closure section of ~6 m (NOAA nautical chart 18504, 1988). Channel and tidal flat depths are reasonably reproduced (Figure 6a) by using a value of Co equal of ~0.2 kg/m3, characteristic of the channel concentration during peak currents and on the tidal flat during moderate waves events [Mariotti and Fagherazzi, 2012b]. It is remarkable to note that the model predicts that the tidal flats of Willapa Bay are indeed intertidal (df < r).

[69] In the strongly wave dominated regime, it is not possible to have equilibrium between a wave dominated tidal flat and a current dominated channel. In this regime, waves are strong enough to suppress the formation and persistence of channels, and the unchannelized sound configuration becomes the only stable equilibrium. When applying this result to real tidal environments, we expect that intense wave conditions will reduce the channel drainage network rather than completely prevent any form of channelization. Such prediction might be investigated by introducing in our model a mechanism governing the evolution of channel width.

6.4 How Strong is the Morphological Coupling Between Tidal Channel and Tidal Flats?

[71] The model predicts that intermediate elevations are not in equilibrium: the platform adjacent to channels can be either close to MHW, resembling a salt marsh, or close to MLW, a typical tidal flat configuration. This result is hence similar to that predicted by one-point models [Fagherazzi et al., 2006; Marani et al., 2007; Marani et al., 2010]. Because the channel quickly adjusts to the sediment boundary condition, its role reduces to a conveyor of sediments entering the system. Channel geometry does affect the fluxes and hence the timescale of tidal flat evolution, but tidal flat equilibrium is solely determined by waves (except for the point Pcw). For example, the switch from the current dominated to the mixed regime occurs when τw,max > τeq, which is independent of channel presence. This result suggests that wave-dominated platform models [Fagherazzi et al., 2006; Marani et al., 2007; Marani et al., 2010] can be decoupled from channels, provided that platform currents are small enough to be negligible in the sediment resuspension process.

[72] On the other hand, channel depth is strongly influenced by the tidal flat elevation, which controls the tidal prism and the flow partition into the channel. Wave regime hence indirectly affects channel depth, through its control on tidal flat morphology. Such effect can be seen in the wide range of channel depths attained in the mixed regime as a function of wind speed (Figures 6a and 6c). The maximum channel depth is found in correspondence of the wind condition such that the tidal flat elevation is at MLW. In this condition, the reference flow velocity in the channel is maximum, because of the optimal combination of tidal prism and flow partition between channel and tidal flat (Figure 3b for df/r = 1). The channel depth decreases for a tidal flat shallower than MLW, i.e., an intertidal flat, and for a tidal flat deeper than MLW, i.e., a subtidal flat. We, hence, suggest that models of tidal channels morphology would benefit from a dynamic coupling with tidal flats.

7 Conclusions

[73] We proposed a novel morphological model for the coupled evolution of channels and tidal flats. The model features few fundamental processes: flow partition between channel and tidal flat and a non monotonic relationship between wave shear stresses and depth. The model uses a novel simplified hydrodynamic approach, validated by numerically solving the shallow water equations applied to a schematic geometry.

[74] According to the model, wind speed is a crucial parameter determining the system dynamics. For low wind speeds, only a trivial stable point is possible, corresponding to a channel without tidal flats. Therefore, tidal flats in equilibrium exist only when waves are present. The unchannelized equilibrium configuration is unstable, suggesting that the model captures the morphological instability leading to channel formation. In the presence of relative sea level rise, the current-dominated equilibrium is characterized by tidal flats placed just below MHW, resembling a channel-marsh configuration.

[75] For intermediate wind speeds, a non-trivial stable equilibrium point characterized by a deep channel and tidal flat placed close to MLW appears. In this scenario, the equilibrium concentration is determined by currents in the channel and by waves on the tidal flats. This equilibrium likely represents real tidal flats, in which both wind waves and tidal currents are relevant erosive processes. For elevated wind speeds, an unchannelized sound is the only stable point. In this scenario, waves suppress the instability that leads to channel incision.

[76] The model suggests that tidal flat elevation dynamics is controlled by waves and can be decoupled from tidal channels, supporting the validity of tidal flat models in which tidal channels are neglected. On the other hand, the model suggests that the tidal flat elevation has a strong influence on the depth of the channel. Indeed, the maximum channel depth is found in correspondence of a tidal flat placed at MLW, and a shallower channel is predicted for both intertidal and subtidal flats. As a consequence, the wave regime becomes an indirect control factor for channel morphology.

[77] Finally, the model predicts that variations in environmental parameters, such as wind regime, sediment availability, and sea level rise can trigger catastrophic morphological shifts, revealing the fragility of intertidal environments. Further studies are required to better understand these transitions between stable states and their dependence on basin geometry, tidal range, and sediment characteristics. Additional dynamic variables and processes, such as tidal flat extension, marsh boundary erosion, organogenic sediment production, and sediment biostabilization can be easily added to the model, broadening its applicability.

Acknowledgments

[78] This research was supported through the VCR-LTER program award DEB 0621014, by the NSF award OCE-0924287, and by the Office of Naval Research award N00014-10-1-0269.