Sensitivity of tidal characteristics in double inlet systems to momentum dissipation on tidal flats

In a tidal channel with adjacent tidal flats, along–channel momentum is dissipated on the flats during rising tides. This leads to a sink of along–channel momentum. Using a perturbative method, it is shown that the momentum sink slightly reduces the M 2 amplitude of both the sea surface elevation and current velocity and favours flood dominant tides. These changes in tidal characteristics (phase and amplitude of sea surface elevations and currents) are noticeable if widths of tidal flats are at least of the same order as the channel width, and amplitudes and gradients of along–channel velocity are large. The M 2 amplitudes are reduced because stagnant water flows from the flats into the channel, thereby slowing down the current. The M 4 amplitudes and phases change because the momentum sink acts as an advective term during the fall of the tide, such a term generates flood dominant currents. For a prototype embayment that resembles the Marsdiep–Vlie double–inlet system of the Western Wadden Sea, it is found that for both the sea surface elevation and current velocity, including the momentum sink, lead to a decrease of approximately 2% in M 2 amplitudes and an increase of approximately 25% in M 4 amplitudes. As a result, the net import of coarse sediment is increased by approximately 35%, while the transport of fine sediment is hardly influenced by the momentum sink. For the Marsdiep–Vlie system, the M 2 sea surface amplitude obtained from the idealised model is similar to that computed with a realistic three–dimensional numerical model whilst the comparison with regard to M 4 improves if momentum sink is accounted for.


Introduction
In shallow coastal seas, such as the Wadden Sea in Northern Europe, the tidal motion is highly nonlinear, as is evident from the presence of residual and higher harmonics in the primary tide (Buijsman and Ridderinkhof 2007).Wellknown sources of these nonlinear tides are the advection of momentum, quadratic bottom stress, depth dependent friction, divergence of excess mass (i.e.mass stored between mean sea level and actual sea surface) and the hypsometry of the embayment (e.g.Parker 1991;Friedrichs and Aubrey 1994, and references therein).As a result, the sinusoidal shape of the curves of the sea surface height and current velocity at a fixed location is distorted.This tidal asymmetry leads to a net transport of fine and coarse sediment (Groen 1967;Aubrey 1986;Boon and Byrne 1981;van de Kreeke and Robaczewska 1993).
In this study, the focus is on nonlinear tides by hypsometry and in particular by dissipation of momentum on tidal flats.When the water level rises, water flows from the channel onto the flats, carrying with it longitudinal momentum.The mass is temporarily stored on the flats, until it flows back into the channel during the falling tide.The momentum, on the other hand, dissipates rapidly due to high friction on the tidal flats.That is, the flats act like mass storages and momentum sinks (Dronkers and Schönfeld 1959;Dronkers 1964).
In contrast to mass storage, which has been extensively studied (Speer and Aubrey 1985;Friedrichs and Aubrey 1994;Ridderinkhof et al. 2014), the role of the momentum sink on tidal dynamics is less well understood.Numerical model studies (e.g. Brown and Davies 2010) show that tidal asymmetry strongly depends on the distribution of channels and flats, thereby suggesting the importance of processes such as the momentum sink.de Swart et al. (2011) and Alebregtse (2015) 1 studied numerically the momentum sink in a semi-enclosed embayment.They demonstrated that the amount by which the momentum sink deforms tidal curves depends on the ratio of flat-to-channel area, the shape of the tidal flats and the distribution of flats along the channel.However, no details about the underlying processes were given.
The considerations above motivate further investigation of the generation of overtides by a momentum sink, which is conducted here.Similar to Alebregtse (2015), the model of Speer and Aubrey (1985) is extended with a momentum sink term.The model of Speer and Aubrey (1985) has been successfully validated against field data (see also Friedrichs andAubrey 1988, 1994).Here, two main new aspects are introduced.First, the approach is analytic in nature, and approximate solutions of the cross-sectionally averaged shallow water equations are explicitly constructed using perturbation methods that exploit the smallness of the Froude number.This enables the assessment of the impact of the momentum sink and quantification of the difference between solutions that do and do not account for it.Second, a double inlet system (rather than a semi-enclosed bay) with an M 2 tidal wave entering at both ends is considered.In particular, this permits investigating what happens in double inlet systems when the two waves entering at opposite sides differ in amplitude or phase, as typically occurs in nature.Examples of such embayments are the Western Wadden Sea in the Netherlands, with the Marsdiep and Vlie inlets (e.g.Ridderinkhof 1988, also considered here), the Laguna de Términos in México (David and Kjerfve 1998) and the Santa María La Reforma in California (Serrano et al. 2013).
This study has two specific aims.The first is to assess changes in tidal characteristics of sea surface elevation and 1 Chapter 4, submitted.
current velocity due to the momentum sink in double inlet systems with different parameter settings.The parameters that will be varied are the width of the tidal flats and the channel width, the drag coefficient, the slope of the flats, the phase difference between the incoming tidal waves and the length of the embayment.The second aim is to understand the mechanisms behind these changes in tidal characteristics.
The paper is organised as follows.In Section 2, the model is presented.In Section 3, a perturbation method is used to find approximate solutions to the cross-sectionally averaged shallow water equations.The results are presented in Section 4, followed by a discussion in Section 5. Here, mechanisms leading to the model results are explained together with a qualitative comparison with the model output of a complex numerical model.Section 6 contains the conclusions.

Physical domain and geometry
The tidal embayment considered in this study consists of a channel of constant width b * c with adjacent tidal flats.The embayment is connected to an open sea on both sides, has length L * and geometry symmetric with respect to the centerline of the channel (see Fig. 1).
The along-channel coordinate is x * , the lateral coordinate y * , the vertical coordinate z * and the time coordinate t * .The origin is placed in the middle of the embayment at the height of the undisturbed sea level.The undisturbed water depth on the flats is d * f and in the channel it is h * c .The width of the channel or of the flats henceforth refers to the width of the wetted part of the channel or flats.The total width b * of the embayment is dependent on the sea surface elevation ζ * , and hence also on t * and x * .The geometry of the cross section is based on that of Speer and Aubrey (1985).It is assumed that the lowest part of the flats are located at the level of low water.The total width reads where, b * max is the width of the embayment including the dry parts of the flats and ϕ the angle of inclination of the flats as in Fig. 1.It follows from the figure that partly horizontal flats with a kink are possible if This motivates the introduction of a parameter q ≥ 1 defined as The total width of the embayment now reads The parameter q ≥ 1 determines flat steepness; q = 1 means linear flats and larger q values imply steeper flats (in Fig. 1, q = 5).

Governing equations
The hydrodynamics in the embayment are modelled by the cross-sectionally averaged shallow water equations.They ( 5 ) Here, u * (ms −1 ) is the current velocity averaged over the cross-section of the channel (excluding the flats) in the along-channel (x * ) direction, ζ * (m) is the sea surface elevation, g * (m s −2 ) is the gravitational acceleration, r * (ms −1 ) is a friction parameter and H is the step or Heaviside function, which is one (zero) when its argument is positive (negative).Further, M * is the momentum sink term which is derived and put into context of previous studies in Appendix A.
In Eqs.3-5, it is assumed that the channel width b * c is constant, the water in the channel has a constant density, there is no effect of wind or waves and the stresses that result from averaging nonlinear momentum fluxes over the cross section can be ignored.On the tidal flats, the velocity in the along-channel direction is assumed to be zero.Instead of a quadratic bottom stress τ * b = ρ * c d |u * |u * , with c d a drag coefficient and ρ * the water density, a linear stress τ * b = ρ * r * u * is used to allow for analytical approximations.The friction parameter r * is determined by demanding that the energy dissipated in one tidal cycle (averaged over the spatial domain) is to be the same in both formulations (e.g.Lorentz 1922;Terra et al. 2005).Thus, in which • denotes the average over one tidal cycle.As, a priori, u * is not known, r * is determined iteratively.
The system of Eqs.1-5 will be considered in the weakly nonlinear regime, i.e. magnitudes of nonlinear terms are small compared to those of linear terms.In that case, as will be shown in Section 3, solutions consist of waves that propagate in both the positive and negative x-direction.The (Riemann) boundary conditions are that at both open sides of the embayment, a prescribed T *periodic tidal wave (where T * = 2π/σ * , with σ * the M 2 tidal frequency) propagates into the embayment and that waves may freely exit the domain.This choice implies that boundary conditions are independent of the dynamics inside the domain.
Furthermore, it is imposed that the average volume of water over one tidal period is constant and that there is no net transport of water through the embayment.The former poses a condition for the residual sea surface elevation and the latter a condition for the residual current velocity.

Model analysis
Equations 3-5 are nonlinear and cannot be solved analytically, necessitating the derivation of approximate solutions, either numerical or perturbative.The former have the advantage that they remain valid over extended regions in parameter space, as they do not require nonlinear terms to be much smaller than the linear ones.The benefit of perturbative methods is that they result in closed form solutions that yield information about the effect of different physical forcing terms on solutions.As mentioned in Section 1, here a perturbation analysis is made.For this, dimensionless equations containing a small parameter are needed in which all variables are of order 1.Let be dimensionless variables, where The scales U and Z are typical values for the current velocity and the sea surface elevation amplitude, respectively.The scales L and T are typical values for a frictionless tidal wave length and period of the principal constituent (both divided by 2π ), quantifying the length and time scale of the dominant dynamics.The scaled versions of Eqs.3-4 are The Froude number is denoted by ε: In many tidal embayments (e.g., the Wadden Sea), ε is around 10 −1 and here assumed to be a small parameter.The dimensionless length of the channel is The parameter β is the ratio of the tidal flat width to the main channel width, the parameter α controls the steepness of the sides of the flats and r is the dimensionless friction parameter.
It is assumed that the area of the flats is of the same order as the channel area, that the slope of the tidal flats is steep and that friction is moderate.The consequence of these choices are that the magnitude of the momentum sink and mass storage are comparable with that of other nonlinear terms in the equations of motion, such as advection.Note that Friedrichs and Aubrey (1994) and Speer and Aubrey (1985) assume β to be small.However, inlet systems like those in the Wadden Sea are characterised by β ≈ 2 (Dronkers 2005).Thus, it is assumed that where the symbol O(1) denotes that they are of order 1.Moreover, O(ε) parts of the friction parameter r and the O(ε 2 ) parts of the last term in Eq. 10 are neglected.All dimensionless variables are assumed to be of order 1, so the magnitude of the terms in the dimensionless equations is determined by the parameters they involve.The Lorenz linearisation condition (6) in dimensionless variables reads with the additional assumption that

Substituting the perturbation series
in Eqs.9-10 and collecting terms of order 1, one obtains the so-called O(1) problem.Likewise, collecting terms of order ε yields the O(ε) problem.The O(1) and O(ε) problem will be treated separately in the following sections.

The O (1) problem
The O(1) problem is linear and reads hence, it admits nontransient solutions of the form with ζ0 (x) and û0 (x) complex amplitudes and Re {•} denoting the real part of a complex number.Equation 16describes the primary tidal wave in the system.By eliminating u 0 from Eqs. 14-15 and substituting ζ 0 from Eq. 16, an ordinary differential equation (ODE) for ζ0 is found: with complex wave number k satisfying where A and B are complex-valued integration constants.When Eq. 18 is substituted into Eq.16, and choosing the root with positive real part (i.e, Re {k} > 0), it follows that the first term on the right hand side of Eq. 18 describes the spatial structure of a right-propagating wave, while the second term describes the structure of a left-propagating wave.Equations 14-15 allow for solutions that consist of waves propagating in opposite directions.Following the considerations in Section 2, conditions at both the left (x = −L/2) and right (x = L/2) boundary are imposed on the waves that enter the domain.At the left boundary, the incoming wave progressing to the right has a sea surface elevation amplitude z l .At the right boundary a leftpropagating tidal wave enters with a sea surface elevation amplitude of z r .The phase difference between the incoming tidal waves at the boundary is θ .Hence, the sea surface variations corresponding to right-propagating wave at x = −L/2 and the left-propagating wave at x = L/2 are, respectively, Thus, high water related to the incoming wave at the left boundary occurs at time t = −Re {k} L/2 and at the right boundary it occurs a time θ later.These boundary conditions are met by choosing with Im {•} denoting the imaginary part of a complex number.Note that, when θ = 0 and z l = z r , the two incoming waves have equal phase and amplitude.In that symmetric case, the velocity u 0 vanishes in the middle of the channel (at x = 0) and water motion in half the embayment (from x = −L/2 to x = 0) behaves as if in a semi-enclosed embayment with the left-propagating wave representing the wave reflected at the closed boundary.Once the sea surface elevation ζ 0 is known, the along-channel current velocity u 0 follows from either Eq. 14 or Eq. 15 using Eq.16, As mentioned in Section 2.2, the friction parameter r is determined from the Lorentz linearisation condition (13).
First, u 0 is calculated with Eq. 20-16, using an initial guess for r.Next, r is updated to satisfy Eq. 13, with u = u 0 and the integrals are calculated numerically.A new value for u 0 is calculated using this updated value for r, and the process is iterated until the relative change in r between successive iterations is less than 0.0001%.This r-value is also used in the O(ε) problem below.

The O (ε) problem
The O(ε) problem reads These equations have the same form as those of the O( 1) problem, but with the known O(1) solutions forcing the unknowns u 1 and ζ 1 through the inhomogeneous terms on the right hand side of Eqs.21-22.Every term in the inhomogeneity originates from nonlinear terms in Eqs.9-10.
In Table 1, the physical source of every nonlinear term is presented.
Since the inhomogeneities I, II, III and IV in the O(ε) problem consist of products of the O(1) solutions, which themselves are M 2 signals, they are composed of M 0 and M 4 harmonics.Term V is different from terms I-IV in the sense that it is not quadratic.Because of the presence of the discontinuous step function H , that term consists of an infinite number of harmonics, where the complex coefficients p m are given in Appendix B. Note that p m depends in a nonlinear way on β through ζ 0 and u 0 .The nontransient solution of the O(ε) problem has the form (24) In the remainder of this section, ζ 1 , u 1 , ζ1,m and û1,m are calculated.The residual current velocity u 1 is determined by averaging Eq.21 over a tidal cycle using the condition that there is no net transport of water, u(1 + εζ ) = 0. So, the residual current only compensates for Stokes transport created by the tidal wave.This yields The residual sea surface elevation ζ 1 results from the tidal average of Eq. 22 and imposing the condition that the average volume of water over one tidal period remains the same, In Eq. 26, the variable s is used as a dummy variable for x.It follows from Eq. 27 that advection of momentum (III), depth dependent bottom friction (IV) and momentum sink (V) lead to residual sea surface elevation.
In order to obtain solutions for ζ1,m , Eq. 24 is substituted into Eqs.21-22.Subsequently, u 1 is eliminated, which results in an ODE for every ζ1,m , with and To study the overtides generated inside the domain by incoming M 2 waves, z l,m and z r,m are chosen to be zero, except when the comparison is made with a complex numerical model in Section 5.3.The complex amplitudes û1,m of the O(ε) current velocity u 1 are found by substituting Eq. 24 in Eq. 22.This yields It follows from these expressions that advection of momentum (III), depth dependent bottom friction (IV) create M 4 velocities and the momentum sink (V) generates M 2 , M 4 , M 6 , M 10 , . . .velocities.

Results
In order to address the two objectives set out in Section 1, parameter values will be chosen that are representative of a prototype double inlet system, viz. the Marsdiep-Vlie system in the Western Wadden Sea.The embayment is roughly 60 km long, 15 km wide and the channel has an averaged depth of 10 m.The primary tide is the semidiurnal lunar tide M 2 and its first overtide M 4 is clearly present (Buijsman and Ridderinkhof 2007).The total tidal range is around 1.4 m at the Marsdiep and 1.8 m at the Vlie and the moment of high water at the Marsdiep differs by approximately 90 min from that at the Vlie (Royal Netherlands Navy 2016).The default parameter values used are listed in Table 2.

Impact of the momentum sink on tides in the default setting
In Fig. 2, results are shown of the primary tide M 2 and its first overtide M 4 of the sea surface elevation and current velocity in the regime where all forcings in Eqs.21-22 are taken into account and for the case where the momentum sink is turned off (term V = 0).Panels (a) and (c) of this figure reveal that including the momentum sink lowers the M 2 amplitude and increases the M 4 amplitude of both the sea surface elevation and the current velocity of the order of centimetres, respectively, centimetres per second.
For the M 2 amplitudes, that change is approximately 2%, for the M 4 amplitudes, approximately 25%.The increase in absolute value of the residual sea surface elevation due to the momentum sink is of the order of 0.5 cm, corresponding to a relative change of 30%.Note that in Fig. 2a, c, the M 2 amplitude of both the sea surface elevation and the current velocity at the left boundary differs from that at the right boundary, even though the amplitudes of the incoming tidal waves at the boundaries are equal.
Figure 2b, d reveal that the M 2 phase of both the sea surface elevation and the current velocity is hardly influenced.However, the phase of the M 4 component of the sea surface elevation is lowered by approximately 15 • .The phase of the M 4 harmonic of the velocity is lowered by approximately 10 • .Figure 2b illustrates that the M 2 sea surface elevation has the character of a propagating wave near the boundaries, but it turns into a standing wave in the middle of the channel as the phase becomes constant there.The M 4 component of the sea surface elevation has the character of a standing wave with a jump in phase around x * = − 6 km.Both the M 2 and M 4 components of the current velocity have the character of standing waves with jumps in phase around x * = 15 km.

Impact of the momentum sink on sediment transport
Changes in tidal harmonics as described above may lead to a difference in sediment transport.The net sediment transport due to the primary tide and its first overtide is calculated as in de Swart and Zimmerman (2009) (after Groen 1967), with ν * (s m −1 ) an erosion parameter and γ * (s −1 ) a settling parameter (1/γ * is the timescale in which sediment settles to the bed), a = σ * /γ * , τ = arg(ε û1,2 ) − 2 arg( û0 + ε û1,1 ) the relative phase of the M 2 and M 4 harmonic, and U * M 2 = U| û0 +ε û1,1 | and U * M 4 = U|ε û1,2 | the M 2 and M 4 amplitude of the current velocity.Coarse sediment is characterised by a small a ≈ 0.01, while a ≈ 5 is representative for fine sediment.Figure 3 shows the tidally averaged volumetric sediment transport per unit mass and width S * , divided by ν * /4γ * , both for when the momentum sink is included and when it is not.The purple lines shows the situation with coarse sediment and the red ones the situation with fine sediment.The arrows indicate the direction of the sediment transport.The figure reveals that including the momentum sink increases the absolute value of the transport of coarse sediment by approximately 35% and has little impact on the transport of fine sediment.On the left (right) side of the channel, the flood current is to the right (left).Hence, the momentum sink increases transport of coarse sediment in the flood direction, that is, it leads to an import of coarse sediment.Furthermore, Fig. 3 shows that fine sediment is transported through the embayment, from the right boundary to the left boundary.Coarse sediment accumulates around x * = 5 km and suggests the formation of a tidal watershed at this location.

Sensitivity study
In this section, it is investigated how sensitive the results of the preceding section are to changes in model parameters.
The difference in amplitude due to the momentum sink is henceforth denoted by ΔAmpl.
Figure 4 shows ΔAmpl in sea surface elevation for various harmonics as a function of x * and for different values of flat-to-channel area ratio β, slope parameter α and drag coefficient c d .Red colours denote an increase and blue colours a decrease in amplitude.Figure 5 is similar to Fig. 4, but for the current velocity instead of the sea surface elevation.The difference in phase due to the momentum sink is not shown since it remains close to zero for the M 2 harmonic and close to 10 • -15 • for the M 4 harmonic when the parameters β, α and c d are varied.The M 0 component of the current velocity is also left out since Eq. 25 already shows that Ampl = 0.
What stands out in Figs.4-5 is that ΔAmpl in both the sea surface elevation and current velocity is most sensitive to changes in flat-to-channel area ratio β, and less to changes in slope parameter α and drag coefficient c d .Figure 4 reveals that the momentum sink lowers the M 2 amplitude and increases the M 4 one for almost all β values considered.Furthermore, it illustrates that ΔAmpl depends nonlinearly on β.In particular, the sensitivity of ΔAmpl to β decreases, as β increases.Figure 5 shows that these conclusions also hold for the current velocity.
Furthermore, Fig. 4 reveals that, although ΔAmpl in M 0 and M 2 are rather insensitive to changes in α, the response of the M 4 harmonic of the sea surface elevation on changing α values is more interesting.On the right side of the embayment, for small α, the momentum sink increases (red colours) the M 4 sea surface elevation amplitude while for larger values of α, the M 4 amplitude is decreased (blue colours).From the bottom panels of Figs.4-5 it follows that,   2 for smaller values of the drag coefficient, the momentum sink lead to a larger difference in sea surface elevation and current velocity amplitude.
The sensitivity of ΔAmpl in sea surface elevation and current velocity to the embayment length L * and the phase difference θ between the incoming tidal waves is presented in Figs. 6 and 7.For increasing θ values, the spatial pattern of ΔAmpl shifts towards the right, and the M 2 and M 4 sea surface elevation amplitude and M 4 current velocity amplitude decreases.Particularly noticeable about the sensitivity of ΔAmpl on L * is that the M 2 sea surface elevation amplitude has a local maximum at L * ≈ 45 km.

Physical interpretation of the results
In this section, some results obtained in the preceding section are interpreted by analysing the expressions from Section 3. To that end, it is helpful to rewrite term V in the O(ε) problem ( 22) using the continuity Eq. 14.That term becomes: This shows that the momentum sink is similar to the advection term when the tide falls and, in particular, that the momentum sink term is always smaller than the advection term, by virtue of β/(1 The impact of the momentum sink on the residual sea surface elevation becomes apparent by this similarity.In fact, since H (−∂ζ 0 /∂t) = 1/2, integration over one tidal period of the O(ε) momentum Eq. 22 using Eq.35 results in It follows from this expression that the momentum sink increases the absolute value of the gradient of residual sea surface elevation.This clarifies the results in Fig. 2a.
In Section 4, it was found that the momentum sink reduces the amplitude of the M 2 harmonic for most parameter values.To understand why, note that term M * in Eq. 5 (which arises as term V in the O(ε) momentum equation) represents a force per mass unit, as it appears on the right hand side in the dimensional momentum Eq. 4.During the falling tide, this force acts against the direction of the current (as shown in Eq. 35, term V and u 0 have opposite signs when ∂ζ 0 /∂t < 0), whilst during the rising of the tide it is inactive (term V is zero).In other words, when the water level drops, stagnant water enters the channel from the flats, thereby slowing the channel current.During the rising of the tide, water flows from the channel onto the Fig. 7 As Fig. 6, but for different embayment lengths L * .The bold value on the vertical axis represent the default setting L * = 57 km flats, which leaves the channel current unaffected.All in all, this implies that term V reduces the current of the principal tidal component.The continuity Eq. 14 implies, in turn, that the same holds for the sea surface elevation.
In a system resembling the Marsdiep-Vlie system, the momentum sink increases the M 4 amplitudes (Fig. 2).This can also be explained using Eq.35.This equation expresses that the momentum sink term is the product of a constant, the advection term and the step function.The advection term consists of an M 0 and M 4 harmonic and the step function consists of an M 0 and odd harmonics.The M 4 component of the momentum sink is therefore the product of the M 4 harmonic of the advection term and the M 0 harmonic of the step function.That is, the momentum sink generates an M 4 harmonic with the same phase as the M 4 generated by advection (see Appendix B for details).Ridderinkhof et al. (2014) showed that advection leads to stronger and shorter flood currents (flood dominance).The momentum sink therefore also favours flood dominance.Since Fig. 3 shows that the transport of coarse sediment is in the flood direction, it implies that the current is flood dominant.Thus, the M 4 generated by the momentum sink increases the total M 4 current velocity amplitude.
Figure 4 showed that, when the slopes of the flats are small (large α), the M 4 amplitude of the sea surface elevation is decreased by the momentum sink while it is increased for the default value of α.The reason for this is similar to the above.When α is large, the mass storage is the dominant mechanism producing the M 4 overtide in the sea surface elevation.Since this term generates M 4 that is out of phase with the M 4 generated by the momentum sink, including the momentum sink reduces the total M 4 amplitude of the sea surface elevation.
Figure 2 revealed that the amplitudes of the sea surface elevation (and the current velocity) are different at the two boundaries.In these experiments, the amplitudes of the incoming tidal waves are equal, so even then, the resulting amplitudes at the boundaries differ.This is already included in the linear dynamics.The O(1) solution ( 18)-( 20) of the linearised problem yields, if z l = z r , that Hence, for equal amplitudes of incoming waves, the amplitude at the boundaries differ when θ and Re {k(r, β)} L are not a multiple of π, with the difference depending on the phase difference θ , the friction parameter r, the length of the channel L and the flat-to-channel area ratio β.Furthermore, it follows from Eq. 37 that, in a longer embayment, the difference in amplitude due to the phase difference of the incoming tidal waves, becomes smaller.Physically, this makes sense, since in a longer embayment, a wave needs more time to propagate from one side to the other.All this time it is subject to friction.Hence, for large L, ζ0 (±L/2) is mainly determined by the incoming wave.Figure 8 1) sea surface elevation amplitude at the left and right boundary for different flat-to-channel area ratios β and phase differences θ between the incoming tidal waves.The difference varies between + 10 and − 50 cm and is largest for small β and large θ values (blue area).

Model limitations
In this study, the bottom stress is linearised to obtain analytically tractable (linear) O(1) and O(ε) problems.When the M 2 tide is dominant, nonlinear bottom stress is known to especially generate odd overtides M 6 , M 10 ,. . .(Parker Fig. 8 The difference between the O(1) sea surface elevation amplitude at the left (x = −L/2) and right (x = L/2) boundary versus flat-to-channel area ratio β and phase difference θ between the tidal waves entering the embayment at the left and right boundary.The remaining parameters have the same value as in Table 2.In particular, the amplitudes z l and z r of the incoming tidal waves are equal.The red dot denotes the default value of β and θ 1984).It is possible to include an approximation of the M 6 generated by the nonlinear bottom stress.However, as stated in Friedrichs and Aubrey (1994), one-dimensional numerical models, even with nonlinear bottom stress, simulate observed M 6 variations rather poorly.The authors remark that this might be due to the invalid assumption of a constant drag coefficient c d .Therefore, as in Friedrichs and Aubrey (1994), the focus of this study is on the primary tide and its first overtide.
In natural tidal embayments, along-channel velocity u is strongly reduced on the flats, but does not vanish completely.The choice of vanishing u on the flats is made to consider the extreme case.As is shown in Section 4, there are (physically relevant) parameter regimes in which the momentum sink generates noticeably overtides.It would be interesting to relax this assumption in a future study.
The trapezoidal cross section is chosen to approximate concave-up tidal flats which typically occur in nature (Friedrichs 2011).The approximation is such that the partly linear flat is tangential to the convex bed at y * = ±b * max /2 and y * = ±b * c /2 (see Fig. 9).In this study, steep slopes and wide flats are chosen to obtain a system in which all nonlinear terms, and in particular both the momentum sink and the mass storage, are O(ε), i.e., β = O(1) and q = O(1/ε).Other choices are possible.For example, Speer and Aubrey (1985) considered moderate slopes and narrow flats, i.e. β = O(ε) and q = O(1).In that case, the momentum sink term is 2 ) and the mass storage is O(ε).
Furthermore, for simplicity, an embayment with a constant depth and width is chosen.Also, effects of wind, waves, density differences and radiation damping (as in Roos and Schuttelaars 2015) are ignored.These are possible extensions of the model.

Comparison with results of a complex numerical model
In this section, the results of Section 4 are compared with hydrodynamics simulations of the Dutch Wadden Sea performed using the General Estuarine Transport Model Fig. 9 The trapezoidal cross section considered in the idealised model (solid line), which is an approximation of the concave-up flats (dashed line) A transect across the Western Wadden Sea is chosen (see Fig. 10a) across which the tidal channels connect the Marsdiep and Vlie inlet (similar as in Ridderinkhof 1988).In Fig. 10b, the M 2 and M 4 amplitude of the sea surface elevation in April 2009 (a calm month in terms of wind) along this transect are depicted together with the ones from the idealised model (as in Fig. 2).The magnitude and spatial distribution of the M 2 harmonic of the sea surface elevation as calculated by the idealised model and that of GETM roughly agree.The M 4 amplitude calculated by GETM is approximately four times larger as that calculated by the idealised model.However, when the momentum sink is accounted for, the error becomes smaller (reduction of approximately 7%).A possible explanation for this difference is that, in the idealised model, there is no M 4 harmonic in the tidal wave entering the domain.In GETM, the M 4 amplitude of the sea surface elevation 14 km northwest of Texel is approximately 0.082 m.The dasheddotted line in Fig. 10b represents the M 4 harmonic of the sea surface elevation when such an M 4 tidal wave entering at the boundaries is imposed.That is, the sea surface elevation of the right-propagating wave at x = −L/2 and the left-propagating wave at x = L/2 are chosen to be, respectively, z l cos(Re {k} x − t) + εz l,2 cos(Re {k} x − t) and z r cos(Re {k} x + t − θ) + εz r,2 cos(Re {k} x + t − θ m ), with z l = z r = z l,2 = z r,2 = 0.7 and θ 2 = 0.In that case, the difference between the M 4 amplitude of the idealised model and that of GETM decreases at the boundaries, but increases in the middle of the channel.
Figure 2 suggests that the amplitude of the sea surface elevation at the Vlie will be higher than at the Marsdiep.This agrees with data presented in the Royal Netherlands Navy (2016) and the results of GETM also demonstrate this behaviour (Fig. 10b).In Section 5.1, it was shown that this amplitude difference is also possible without a difference in the external forcing.

Conclusions
The first objective of this study was to assess changes in tidal characteristics of sea surface elevation and current velocity due to the momentum sink in double inlet systems.It is found that the momentum sink decreases the M 2 amplitude of both the current velocity and the sea surface elevation and favours flood dominant tides.In a prototype system, viz. the Marsdiep-Vlie system in the Western Wadden Sea, it implies a decrease of M 2 amplitudes of approximately 2% and an increase of both M 4 amplitudes by 25%.The absolute value of the residual sea surface elevation increased by 30% due to the momentum sink, whereas the residual current velocity was unaffected.The phases of the M 2 tidal harmonic of both sea surface elevation and current velocity were not influenced by the momentum sink; only the M 4 phases were lowered by approximately 10 • -15 • .In total, this amounts to an increase of net import of coarse sediment by approximately 35%.Accumulation of coarse sediment was found inside the domain, indicating the formation of a tidal watershed.The transport of fine sediment was hardly influenced.
It is found that the changes in tidal characteristics are most sensitive to the flat-to-channel area ratio.The change in amplitude depends nonlinearly on that ratio; the sensitivity is smaller for larger values of the ratio.Furthermore, it is found that the momentum sinks affects the tidal harmonics of the sea surface elevation and current velocity more strongly when the drag coefficient is small and when the slope of the tidal flats decreases, at the right side of the channel (0 < x < L/2), the sea surface elevation M 4 amplitude does the same.A phase difference between the incoming tidal waves leads to spatial shifts of the curves that show amplitudes and phases versus alongchanel distance and to a small decrease of the M 2 and M 4 amplitudes.Finally, it is found that the difference in M 2 sea surface elevation amplitude due to the momentum sink has a local maximum for a length of the embayment at approximately 45 km.
The second objective was to understand the mechanims behind the changes in tidal characteristics due to the momentum sink.The reduction of the M 2 harmonic amplitude of the sea surface elevation and current velocity by the momentum sink is attributable to that, during the fall of the tide, still water enters the channel and slows down the current.The increase in M 4 amplitude of the sea surface elevation and current velocity is explained by noting that the momentum sink acts as an advective term during the fall of the tide, but with a smaller amplitude.Since advection is known to favour flood currents (e.g.Friedrichs and Aubrey 1988;Ridderinkhof et al. 2014), the momentum sink does so as well and increases the M 4 amplitude in a flood dominant system.
A comparison with the complex numerical model GETM showed that the two models produce similar M 2 amplitudes of the sea surface elevation.When the momentum sink is accounted for and an M 4 harmonic is added to the incoming tidal waves, the amplitude of the M 4 harmonic of the sea surface elevation shows a closer resemblance to the one modelled by GETM.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/),which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Appendix A: Derivation of term M * in Eq. 4 In this section, the expression for the momentum sink term M * in Eq. 4 is derived.The derivation only concerns dimensional equations and the asterisks are omitted to keep notation as simple as possible.Equations 3-4 follow from integrating the depth averaged shallow water equations over the channel width.The latter equations read where j 1 = 0 when the water level rises and otherwise depends on the velocities on the flats.Speer and Aubrey (1985) and Friedrichs and Aubrey (1988) assumed the term (b − b c )/b c to be small, and M was therefore neglected.Speer (1984) also considered the case with j 1 = 1.In that case, the tidal flats act as momentum storage regions.The choice made in Alebregtse (2015), de Swart et al. (2011) and the current study is j 1 = H − ∂ζ ∂t .

Appendix B: Fourier coefficients of momentum sink term
In the main text term V is written as In order to see that the phase of the M 4 harmonic of term V equals that of the advection term, note that term V can be written as a of a constant, the advection term and the step function as in Eq. 35.From Eqs. 16-20, it follows that the advection term consists of an M 0 and M 4 harmonic, where A M 0 is a real number and A M 4 a complex number.Since the even Fourier coefficients of the step function c m are zero, multiplying ( 46) and (45) yields that the M 4 harmonic of the product equals Hence, the phase of the M 4 harmonic of term V equals arg(A M 4 ) and thus the phase of the M 4 harmonic of the advection term.

Fig. 1
Fig. 1 Geometry of the open embayment (a sideview at any longitudinal position, b topview).Here, b * max is the width of the embayment including the dry parts of the flats, b * (t * ) the width of the wetted part of the embayment and b * c the width of the channel.Furthermore, h * is the undisturbed depth (being h * c in the channel and d * f on the flats), ϕ the angle of inclination of the flats, L * the length of the embayment and ζ * the sea surface elevation.The origin is at the intersection of the dotted lines.In (a), the x * axis points into the paper.The blue area in (b) is the wetted part ) where | û0 | and φ û0 are the amplitude and phase of the O(1) current velocity û0 .Likewise, | ζ0 | and φ ζ0 are the amplitude and phase of the O(1) sea surface elevation ζ0 .

Fig. 2 Fig. 3
Fig. 2 a, b Amplitude and phase of the sea surface elevation ζ * versus distance x * , for different harmonics M 0 , M 2 and M 4 and for the default parameter setting (Table2).The solid line represents the solution accounting for all the nonlinear terms (I-V), including the momentum sink (term V).The dashed line represents the solution with the momentum sink neglected.The shaded area marks the difference between them.c, d As (a) and (b), but for the current velocity u *

Fig. 4
Fig. 4 Contours of the momentum sink-induced amplitude difference ΔAmpl in sea surface elevation versus flat-to-channel area ratio β (top row), slope parameter α (middle row), drag coefficient c d (bottom row) and along-channel variable x * for the M 0 (left column), M 2 (middle column) and M 4 (right column) harmonics.The bold values on the vertical axis represent the default setting, β = 2, α = 1 and c d = 0.0025.The remaining parameters have the same value as in Table 2

Fig. 5
Fig. 5 Contours of the momentum sink-induced amplitude difference ΔAmpl in current velocity versus flat-tochannel area ratio β (top row), slope parameter α (middle row), drag coefficient c d (bottom row) and along-channel variable x * for the M 2 (left column) and M 4 (right column) harmonics.The bold values on the vertical axis represent the default setting, β = 2, α = 1 and c d = 0.0025.The remaining parameters have the same value as in Table 2

)Fig. 6 a
Fig.6a Contours of the momentum sink-induced amplitude difference ΔAmpl, of the sea surface elevation versus the phase differences between the incoming tidal waves θ and along-channel variable x * for the M 0 (left panel), M 2 (middle panel) and M 4 (right panel)

Fig. 10 a
Fig. 10 a Transect across the Marsdiep-Vlie inlet systems located in the Western Wadden Sea.Colours indicate the bathymetry used in GETM.b The M 2 and M 4 amplitudes of the sea surface elevation calculated by GETM on this transect (dotted lines).The solid lines are the (ms −1 ) is the depth averaged longitudinal velocity, v (ms −1 ) the depth averaged lateral velocity and τb (Nm −2 ) the bottom stress.Note that ζ does not depend on y.Integration of Eq. 39 over the constant width of the channel [−b c /2, b c /2], using the continuity Eq. 3 and neglecting the stresses ( ū−u)( ū−u) arising from averaging over the width, yields b c (ζ + h c u (without the bar) represents the cross-sectionally averaged velocity.Term A arises from integration of the first term in Eq. 39 over the width of the channel and term B from integrating the second term in Eq. 39 and subsequently substituting the continuity Eq. 3. Term C represents the lateral exchange of longitudinal momentum between the channel and the flats, and it equalsC = (ζ + h) v ū| y=b c /2 − (ζ + h) v ū| y=−b c /2 .(41)This expression is subsequently rewritten in terms of crosssectionally averaged velocities.Integration of the continuity Eq. 38 over the right and left flats yields the volume transport (per unit length) through the boundary of the main channel,(ζ + h) v| y=±b c /2 = ± b of thewater level, v is positive at y = b c /2 and along-channel momentum is transferred towards the flats.When the water level falls, water moves from the flats to the main channel.Since the momentum is dissipated on the flats, no along-channel momentum returns to the channel.This motivates the choice ū| y=±b c /2 = uH ∂ζ ∂t .(43) Substitution of Eqs.42-43 in (41) yields C = (b − b c )uH ∂ζ ∂t ∂ζ ∂t .(44) The momentum sink term M as in the main text is now obtained by combining terms A, B and C, using 1 − H (x) = H (−x) and dividing by b c (ζ + h c ), which yields sink is neglected ( ū| y=±b c /2 = u), then A + B = −C and hence, M = 0.To frame Eqs.3-5 in context of previous studies, note thatDronkers (1964) suggested that M || ζ0 | sin(φ û0 − φ ζ0 ),where (in this appendix) • * denotes the complex conjugate and

Table 1
Source of the forcing terms

)
It follows from Eqs. 29-31 that advection of momentum (III), depth dependent bottom friction (IV) and divergence of excess mass (I) create an M 4 harmonic in the sea surface elevation.The mass storage on the tidal flats (II) generates M 2 and M 4 .The momentum sink (V) generates M 2 , M 4 , M 6 , M 10 , M 14 , M 18 , . . . .The momentum sinks skips generation of M 8 , M 12 , M 16 , . . ., because p m = 0 if m > 2 and even.Solutions of the ODEs in Eq. 28 are obtained using variation of parameters, ζ1,m = C m e ik m x + D m e −ik m x m and D m are complex-valued integration constants.The boundary conditions for the O(ε) problem (21)-(22) are similar as those for the O(1) problem.At x = −L/2, it is imposed that the M 2m tidal wave that enters the domain has sea surface amplitude z l,m .Similarly, at the right boundary (x = L/2) a left-propagating M 2m tidal wave enters with a given sea surface amplitude z r,m and a phase shift θ m with respect to the incoming wave at the left boundary.These conditions are met when