Morphodynamic Equilibria in Double‐Inlet Systems: Existence and Stability

The existence of morphodynamic equilibria of double‐inlet systems is investigated using a cross‐sectionally averaged morphodynamic model. The number of possible equilibria and their stability strongly depend on the forcing conditions and geometry considered. This is illustrated by considering a rectangular double‐inlet system forced by M2 tidal constituents only. Depending on the M2 amplitudes and phases at both entrances, no equilibrium, one equilibrium or multiple morphodynamic equilibria may exist. In case no equilibrium is found, the minimum water depth becomes zero somewhere in the system, reducing the double‐inlet system to two single‐inlet systems. In the other cases, the location of the minimum water depth and the direction of the tidally‐averaged sediment transport, as well as their actual values, depend strongly on the M2 tidal characteristics. Such parameter sensitivity is also observed when including the residual and M4 forcing contributions to the water motion, and when allowing for width variations. This suggests that, when considering a specific system, the number and stability of morphodynamic equilibria, as well as the characteristics of these quantities, can only be assessed by investigating that specific system in detail. As an example, the Marsdiep‐Vlie inlet system in the Dutch Wadden Sea is considered. It is found that, by using parameter values and a geometry characteristic for this system, the water motion and bathymetry in morphodynamic equilibrium are qualitatively reproduced. Also the direction and order of magnitude of the tidally‐averaged suspended sediment transport compare well with those obtained from a high‐complexity numerical model.

DENG ET AL.

10.1029/2021JF006266
3 of 23 One of the great advantages of these semi-empirical models is their computational efficiency, allowing for extensive sensitivity analysis of parameters such as tidal range, phase difference of the tidal forcing at the inlets and basin geometry, factors that strongly influence the morphodynamic stability of the inlet systems. However, these models do not allow the morphodynamic evolution of the back-barrier basin, even though observations and model studies show these regions to be morphodynamically active. For example, Dastgheib et al. (2008) used a high-complexity numerical model to simulate the long-term sediment transport and bottom evolution for a double-inlet system, resulting in the development of channel-shoal systems in the back-barrier basin. After simulating 2,000 years, the system was assumed to be close to a morphodynamic equilibrium with both inlets still open. A detailed analysis of the physical mechanisms resulting in the observed morphodynamic equilibrium is challenging using high-complexity numerical models, as it allows one to gain insight into the sensitivity of the morphodynamic equilibria to variations in parameters and geometric characteristics.
The aim of the present study is therefore to develop a process-based model for double-inlet system on mesotidal barrier coasts (Hayes, 1979) in which both the tidal inlets and back-barrier basin are morphodynamically active. The model will be used to analyze the existence of morphodynamic equilibria, their stability, and sensitivity to parameters and geometry. It should allow for a systematic analysis of the physical mechanisms leading to these morphodynamic equilibria. Motivated by the insights gained from width-averaged equilibrium models for single-inlet systems (Meerman et al., 2019;Schuttelaars & de Swart, 1996;Schuttelaars & de Swart, 2000;ter Brake & Schuttelaars, 2010), this type of models will be extended to allow for two tidal inlets connected to each other by a back-barrier basin of arbitrary length and width. Default parameters used are characteristic for a double-inlet system in the Dutch Wadden Sea, and results will be compared with observations and numerical model results from the Marsdiep-Vlie system. The existence and sensitivity of morphodynamic equilibria to geometry and forcing conditions off the seaward sides of the two inlets, consisting of M 2 and M 4 tidal constituents and a subtidal signal, will be investigated in detail.
In Section 2, the equations governing water motion, transport of sediment and bed evolution are presented. The scaling of the system of equations and the solution method are presented in Section 3. In Section 4, morphodynamic equilibria in double-inlet systems and their linear stability are studied. In Section 5, the results are discussed and in Section 6 conclusions are presented.

Model Description
The geometry under consideration is that of a tidal inlet system that is connected to the open sea by two inlets (for a top view see Figure 2b, which is obtained using the observed geometry in Figure 2a). The x axis in Figure 2b represents the distance along the centerline indicated in Figure 2a, whereas the y axis represents the distances normal to the centerline. The inlet system is aligned with the horizontal x axis, and has a prescribed length L. The basin is connected to the open sea at x = 0 and x = L, respectively. The width is allowed to vary with the longitudinal coordinate, and is denoted by B(x) = B 2 (x) − B 1 (x), with B 2 and B 1 indicating the left and right coastal boundaries, looking from inlet I to inlet II. These coastal boundaries are assumed to be non-erodible. The width of inlet I is denoted by B I = B (0), and the width at inlet II by B II = B(L).
The undisturbed water depth at inlet I is denoted by H I , and at x = L by H II (see Figure 2c for a side view). The free surface is represented by =̂ , with the undisturbed free surface located at z = 0. The erodible bottom, that is assumed to consist of sandy material with a single grain size, is described by the equation =ĥ − , where H I is used as the reference depth and ĥ is the local bed elevation measured from level identified by the reference depth. Hence the instantaneous local water depth is given by −ĥ +̂ .
Since, both the length and width of the inlet system are much larger than the reference water depth, and the width is much smaller than the basin length and the Rossby deformation radius, the water motion can be described by the cross-sectionally averaged equations for a homogeneous fluid (Csanady, 1982): DENG ET AL.
10.1029/2021JF006266 4 of 23 Similarly, u (x, t) denotes the width-averaged horizontal longitudinal velocity and h the width-averaged bed level. Time is denoted by t, and g denotes the gravitational acceleration. Subscripts denote a derivative with respect to that variable. Following Lorentz (1922) and Zimmerman (1992), a linearized formulation of the friction is used with the bottom friction coefficient defined as r * = 8Uc d /3π, with U a characteristic velocity scale and c d a drag coefficient. The characteristic velocity U will be introduced in Section 3.1.

of 23
It should be noted that, when considering the morphodynamic evolution of a tidal basin, the effect of radiation damping on the tidal forcing at the seaward side can only be neglected in the so-called deep sea limit (Roos & Schuttelaars, 2015). However, if one is only interested in morphodynamic equilibria, prescribing the tidal forcing is allowed. One should interpret the morphodynamic equilibria as obtained for that prescribed tidal forcing. For the tidally-averaged hydrodynamic components, we require The constants 2 ( 4 ) and 2 ( 4 ) denote the amplitude of the M 2 (M 4 ) tide at the seaward sides of inlet I and II, respectively. Their corresponding phases are given by 2 ( 4 ) and 2 ( 4 ) . The angular frequency of the M 2 tidal signal is given by σ = 2π/T, with T the M 2 tidal period. Condition (4a) implies that at the first inlet the tidally-averaged water depth is given by H I , whereas the mean water depth at inlet II may deviate with respect to the undisturbed water level due to a mean difference in the mean sea surface elevation on the seaward side of the two inlets. The boundary condition (4b) at inlet II assumes the depth-averaged residual water transport equals Q * . When Q * is positive it represents a residual transport from inlet I to II, a negative value implies a net transport from inlet II to inlet I.
The sediment in the tidal inlet system consists of non-cohesive material with a uniform grain size of 2 ⋅ 10 −4 m. The sediment is mainly transported as suspended load, which is described by the width-averaged and depth-integrated concentration equation (ter Brake & Schuttelaars, 2010: where C is the depth-integrated and width-averaged suspended sediment concentration, k h* is the horizontal diffusivity, k v* is the vertical diffusivity, w s is the settling velocity, α is the erosion parameter (Dyer, 1986) related to sediment properties, and β is the deposition parameter, defined by ter Brake and Schuttelaars (2010) for details.
The tidally-averaged erosion and deposition are assumed to balance each other at seaward entrances, and it is assumed that no diffusive boundary layers develop at these locations in the time-dependent parts of the sediment concentration (Schuttelaars & de Swart, 2000). The resulting boundary conditions read: The width-averaged bed evolution equation is derived from the mass balance in the sediment layer and reads: Here, ρ s is the density of the sediment and p denotes the bed porosity. The first and second term on the right of Equation 8 model the local erosion and deposition of sediment, respectively.
Substituting Equation 5 into Equation 8 results in the following bed evolution equation: 6 of 23 the total depth-integrated and width-averaged sediment transport. This transport consists of three terms: a diffusive contribution related to the diffusion of depth-integrated concentration (F diff ), a diffusive contribution related to topographical variations (F topo ) and an advective contribution (F adv ).
At both boundaries, we prescribe the undisturbed bed level as, Since, erosion and deposition are assumed to balance each other at both entrances (see Equation 7), the undisturbed water depth at the entrances will not change during a morphodynamic experiment. The focus of this paper is on morphodynamic equilibria and not the evolution toward these equilibria. Hence, these boundary conditions are appropriate to obtain morphodynamic equilibria characterized by the specific depths prescribed at the seaward sides of the inlet (see also de Swart (1996, 2000) for a discussion of boundary conditions in a single inlet system). When considering the morphodynamic evolution in time of tidal inlet systems, instead of focusing on morphodynamic equilibria, these boundary conditions are too restrictive and more dynamic ones have to be used (see e.g., Bolla Pittaluga et al., 2014;Lanzoni & Seminara, 2002).

Scaling and Expansion
To make the equations dimensionless, the physical variables are scaled as, where the dimensionless variables are indicated by ⋅ . For the horizontal coordinate x, the length L of the double-inlet system is the appropriate scale as the focus is on basin scale dynamics. The width is made dimensionless using the width at inlet I, B I . Time is made dimensionless using the M 2 angular frequency σ, the surface elevation using the M 2 amplitude at the seaward side of inlet I, denoted by 2 , and the bed level is made dimensionless using the depth H I at inlet I. The typical scale for the velocity is = 2 ∕ . It is obtained by assuming that the amplitude of the sea surface elevation at the first inlet is nonzero, otherwise the amplitude at the seaward side of the second inlet should be used (see details in Text S1 in Supporting Information S1). The suspended sediment concentration is made dimensionless using 2 * ∕ 2 , which follows from assuming an approximate balance between erosion and deposition. In Table 1, characteristic values of the relevant parameters are given for the Marsdiep-Vlie system.
Substituting these dimensionless variables in the equations and suppressing the checks ⋅ , the system of dimensionless equations reads: with the dimensionless deposition parameter β given by,   Ridderinkhof (1988) and Duran-Matute et al. (2014) DENG ET AL.
10.1029/2021JF006266 7 of 23 The parameter = 2 ∕ is the ratio of the M 2 tidal amplitude to the water depth at inlet I. In many tidal inlet systems ϵ ∼ 0.1, indicating that ϵ is a small parameter. Hereafter, it is assumed that ϵζ is always much smaller than the local undisturbed water depth 1 − h, thus neglecting the effects of drying and flooding. Note that ϵ is also equal to the ratio of the tidal excursion length and the length of the tidal inlet system Uσ/L. The parameter λ L = k g L is the product of the frictionless tidal wavenumber = ∕ √ to the length of the inlet system. The dimensionless friction parameter is denoted by r and is defined as r = r * /H I σ. The ratio of the deposition timescale to the tidal period is denoted by = * ∕ 2 , and the sediment Péclet number λ d = H I w s /k v* is the ratio of the typical time it takes for a particle to be mixed through the water column to the typical time it takes to settle in the water column. The parameter δ s = αU 2 /(ρ(1 − p)H I σ) denotes the ratio of tidal period T to the time scale related to suspended load and is typically small. All dimensionless parameters and the assigned values are summarized in Table 2.
The bed evolution Equation 13d shows that at the tidal timescale the bed changes are very small, namely of O (δ s ). This implies that the bed changes significantly at a much larger time scale, the so-called morphodynamic time scale, defined as τ = δ s t. To distinguish between these time scales, a multiple timescale method is used, with t the short tidal time variable and τ the long morphodynamic time variable. Using this multiple timescale approximation, it can be demonstrated that at the leading order of approximation the bed changes are due to tidally-averaged convergences and divergences of sediment (see Krol, 1991;Sanders & Verhulst, 1985). The resulting tidally-averaged bed evolution equation reads, The dimensionless boundary conditions associated to Equations 13 and 15 read: = 2 cos( − 2 ) + 4 cos(2 − 4 ) at = 1, Parameter Definition Value Parameter Definition Value

able 2 The Definition of Dimensionless Parameters and Their Values for the Marsdiep-Vlie System
DENG ET AL. 10.1029/2021JF006266 8 of 23 The parameter γ is the ratio of the amplitude of the M 4 tide and M 2 tide at inlet I, and is the ratio of the amplitude of the M 2 (M 4 ) tide at inlet II and at inlet I. Parameter 2 = 2 − 2 ( 4 = 4 − 2 2 , 4 = 4 − 2 2 ) is the phase difference between the M 2 tide at inlet II (the M 4 tide at inlet I, the M 4 tide at inlet II) and M 2 tide at inlet I. Q = Q * /(B I H I U) is the dimensionless residual water transport.
In most tidal inlet systems, the parameters ϵ and γ are much smaller than 1 (see Table 2), allowing for the introduction of an asymptotic expansion in ϵ and γ of the physical variables Φ ∈ {ζ, u, C}, where the first superscript denotes the order in ϵ while the second denotes the order in γ. Because the water motion is forced with a time-periodical signal, each physical variable can be decomposed as a (infinite) sum of tidal constituents and a residual component, where the subscript "res" denotes the tidally-averaged contribution to the variable Φ, and the contributions that temporally vary as cosines (sines) with frequency k are denoted with the subscript ck (sk).
The system of Equations 13a, 13b, 13c, 15 and the boundary conditions (16) can first be ordered in terms of ϵ and γ, and then in terms of the tidal constituents.
To obtain the tidally-averaged contributions that include both the dominant advective and diffusive transport mechanisms, the system of Equations 13a, 13b, 13c and 15 has to be solved up to the orders ( ) and ( ) . Both the water motion and the suspended sediment concentration consist of a M 2 and a M 4 tidal constituent, and a residual component. This results in tidally-averaged sediment transport contributions due to diffusive and advective processes. The diffusive transport consists of the two contributions 00 diff and 00 topo . The advective transport can be decomposed in two contributions, denoted as 20 adv and 11 adv of order ( 2 ) and ( ) , respectively. The internally-generated advective transport 20 adv is due to the temporal correlation of internally-generated overtides and residual velocities at ( ) with the leading order concentration fields, and the correlation of the leading order velocities with the ( ) concentration. The externally-generated advective transport 11 adv results from the temporal correlation of externally-generated overtides with the leading order concentration fields, and the correlation of the leading order velocities with the ( ) concentration. The resulting tidally-averaged sediment transport contributions are then given by,

Morphodynamic Equilibria and Linear Stability
The resulting system of morphodynamic equations, ordered in terms of the small parameters and expanded in tidal constituents, can be written as, 9 of 23 where Ψ is the vector of physical variables and p the vector of model parameters. The matrix K is a diagonal matrix, with a non-zero element only in the row associated with the bed evolution equation. The nonlinear operator G depends on the parameters p and works on the physical variables Ψ.
To obtain solutions of this system of equations, Equation 23 is discretized using a finite element method with continuous Lagrange elements (Alnaes et al., 2015, see also Text S2.1 in Supporting Information S1). Adopting a weak formulation and using a Galerkin method (Brenner & Scott, 2002), the discretized system of equations reads, in which Ψ is the discretization of Ψ,  is the discretized matrix corresponding to K and  is the discretized nonlinear operator working on Ψ and the parameters p.
This discretized system of equations can be analyzed using two different approaches, the initial value and the bifurcation approach. The latter approach aims at the direct identification of morphodynamic equilibria, whereas with the former approach an initially prescribed bathymetry is integrated in time.
In this paper, the focus is on the bifurcation approach with which we directly solve for equilibrium states Ψ e that satisfy, in which p 0 is a vector of prescribed parameter values. To obtain the equilibrium solutions Ψ e , a Newton-Raphson iterative method is used. An initial guess Ψ i of a morphodynamic equilibrium is updated iteratively until the corrections are small enough (in this study the iteration is stopped when the maximum correction is smaller than 10 −8 ).
For this iterative procedure to converge, the initial guess Ψ i has to be close enough to the morphodynamic equilibrium associated with the parameter vector p 0 . If no information about an approximate morphodynamic equilibrium is available, a good initial guess can be obtained by integrating the discretized system (24) starting from an arbitrary initial bed profile in time (in this study a backward Euler method is used as time-integration method, see the Text S2.3 in Supporting Information S1 for details), until an approximate equilibrium is reached. This approximate equilibrium can then be used as the initial guess Ψ i in the Newton-Raphson iterative procedure. If for a given parameter vector p 0 a morphodynamic equilibrium is found, a so-called continuation method (Seydel, 1994) can be employed to find equilibria for different parameter values p 1 by slowly changing the parameters from p 0 to p 1 . During this continuation process, the previously obtained equilibrium is used as a first guess in the iteration process, resulting in the morphodynamic equilibrium for the new parameter vector p 1 . The continuation method employed in this paper is the arclength method (see Crisfield, 1981;Dijkstra et al., 2014; and the discussion in Text S2.2 in Supporting Information S1).
The linear stability of the morphodynamic equilibria Ψ (corresponding to parameter p e ) can be investigated by substituting Ψ =Ψ + ΔΨexp( ) into Equation 24 and linearizing the resulting equation. The resulting eigenvalue problem reads: in which (Ψ , ) is the Jacobian, ω the eigenvalue and ΔΨ the associated eigenvector. An equilibrium is called linearly stable if all its eigenvalues have a negative real part, and unstable if at least one eigenvalue has a positive real part.

Results
In the numerical experiments presented in this section, the influence of the forcing conditions and geometry on morphodynamic equilibria is presented. Assuming a constant width, that is, a system with a rectangular planform, the morphodynamic equilibria are first obtained for a water motion forced by a M 2 tidal constituent and presented in Section 4.1. Next (Section 4.2), all hydrodynamic forcing contributions are included. Finally (Section 4.3), the DENG ET AL.

10.1029/2021JF006266
10 of 23 influence of width variations is investigated. All results are obtained using parameter values that are representative of the Marsdiep-Vlie inlet system (see Table 1), unless mentioned otherwise.

Constant-Width System Subjected to a M 2 Tidal Forcing
The influence of variations in the prescribed M 2 tidal constituent on the resulting morphodynamic equilibria has been studied by setting the amplitude and phase of the M 2 tide at inlet I equal to 0.62 m and 0° and varying the M 2 amplitude and phase at inlet II. To focus only on the influence of the M 2 forcing, the amplitudes of the externally prescribed overtides ( 4 and 4 ) are assumed to be zero. Furthermore, the undisturbed water depths at both inlets are taken to be 11.7 m and the tidally-averaged water transport Q at the second inlet is set to zero. Sensitivity to inlet depth and tidally-averaged water transport at inlet II, assuming no M 4 tidal forcing, is discussed in the Text S3 in Supporting Information S1. All other parameter values are taken from Table 1.
The resulting morphodynamic equilibrium condition reads, where for the numerical experiment under consideration, the width within the inlet system is constant for all x, that is, B(x) = B I .
To obtain these morphodynamic equilibria, the evolution of a spatially uniform bed is studied (i.e., the system of Equation 24 is integrated in time using a Backward-Euler method). If the divergence of the tidally-averaged sediment transport vanishes, an equilibrium is reached. Two numerical experiments have been carried out with different M 2 tidal forcings applied at the second inlet.
In the first numerical experiment, we take 2 = 0.77 m, 2 = 54 • (default value). The resulting bed evolution is shown in Figure 3a. After approximately 2,500 years, the bed profile reaches its equilibrium with a maximum water depth (WD max ) of 17.8 m, located at approximately 20 km from inlet I.
In the second numerical experiment, we take 2 = 0.77 m and 2 = 0 • . The development of the bed profile for the first 23,200 years is shown in Figure 3b. This figure clearly shows that the water depth decreases over time, and after approximately 23,300 years, it vanishes near km 32. The two inlets are not connected anymore, and a system formed by two single inlets is eventually obtained.
From the previous two numerical experiments, it is evident that a morphodynamic equilibrium in which both inlets are connected does not necessarily exist for all combinations of 2 and 2 . The minimum water depth, denoted by WD min , can vanish at some location, disconnecting the two inlets. To assess the existence of morphodynamic equilibria with both inlets connected for a wider range of forcing conditions, we take 2 = 0.77 m, and vary 2 from −180° to 180°. Figure 3c shows the equilibrium water depth as a function of position along the embayment (horizontal axis) and the M 2 phase at inlet II (vertical axis). The equilibrium water depth is found to strongly depend on the tidal phase at inlet II. By increasing 2 from −180° to −2°, WD min decreases from 11.7 to 2.5 m, and by decreasing 2 from 180° to 13°, WD min varies from 11.7 to 6.0 m. For a phase between −2° and 13°, no double-inlet morphodynamic equilibrium exists.
To quickly characterize these morphodynamic equilibria, the quantities WD min and its location will be used. As an example, using a constant amplitude 2 = 0.77 m, WD min is shown as a function of 2 in Figure 3d. Two distinct branches of morphodynamic equilibrium solutions are found, each consisting of a stable (solid) and unstable (dashed) part. The stable and unstable solutions on each branch are connected by a so-called limit points, denoted by L1 and L2. For phases with values between the two limit points, no equilibria were found with both inlets connected. Note that the stable solutions correspond to the equilibria shown in Figure 3c.
To be more precise, Figure 3d indicates that the number of morphodynamic equilibria and their stability strongly depends on the value of the phase at inlet II. For −1.9 • ≤ 2 ≤ 12.8 • no morphodynamic equilibrium exists, for − 53 • ≤ 2 ≤ −1.9 • and 12.8 • ≤ 2 ≤ 47 • two morphodynamic equilibria (one stable and one unstable) exist, whereas for all other phase values one stable equilibrium exists. Furthermore, the number of morphodynamic equilibria and their stability not only depends on the phase of the M 2 tide at inlet II, but also on its amplitude. Figure 4 shows the number of morphodynamic equilibria as a function of the amplitude 2 and the phase 2 .
DENG ET AL.

10.1029/2021JF006266
11 of 23 To further study these equilibria, the dependency of WD min on the amplitude and phase at inlet II is investigated. For 0.31 m ≤ 2 ≤ 1.24 m and phases | 2 | ≥ 60 • it is found that WD min of stable equilibria occur at the entrances reaching the minimum water depth of 11.7 m. Therefore, in the following we focus on | 2 | < 60 • , while 2 is varied between 0.31 and 1.24 m. In Figure 5a, the depth of the watersheds is shown as a function of the amplitude and phase of the M 2 tide at inlet II. From this figure it follows that by increasing 2 from 0.31 to 1.24 m and/or 2 from −60° to 60° results in a shift of WD min from a location closer to inlet I to a location closer to inlet II. The black line in the figure denotes the location where WD min vanishes and, consequently, the double-inlet system reduces to two single-inlet systems. In two regions in the ( 2 , 2 ) -space, complex bifurcation structures occur, indicating that possibly multiple stable morphodynamic equilibria exist. To assess the number of equilibria, the bifurcation structure is carefully analyzed over two paths in parameter space, with one path characterized by a fixed amplitude and varying phase, and the second path by a fixed phase and varying amplitude.
The values of WD min for the first path in parameter space, that is, 2 = 0.94 m and 2 in the range of 0°-30° is shown in Figure 5b. This figure shows that the number of morphodynamic equilibria depends sensitively on the phase angle. Similarly in Figure 5c, WD min is shown for M 2 amplitudes varied between 0.80 and 1.10 m, with the phase fixed at 15.5°. This figure illustrates that the number of morphodynamic equilibria also depends sensitively on the amplitude. Figures 5b and 5c correspond to the same parameter values ( 2 = 0.94 m, 2 = 15.5 • ). Hence, the equilibrium bed configurations denoted by E1 − E4, are the same in both figures. These bed profiles are shown in Figure 6. Equilibria E1 and E3 are linearly stable, whereas E2 and E4 are linearly unstable. Each equilibrium is characterized by a different balance between the various transport terms, contributing to a different, spatially uniform, and total sediment transport. These contributions, divided in two diffusive contributions and one total advective contribution (see Equation 10), and the total transport are shown in Figures 7a, 7c, and 7e, 7g, while Figures 7b, 7d, 7f, and 7h show the different contributions to the advective contribution. In all cases, the total transport is positive, indicating a net transport from inlet I to inlet II. These plots show that all transport mechanisms result in significant contributions to the total transport. Moreover, Figure 8 indicates that the direction and magnitude of the total transport strongly depends on the forcing conditions at inlet II. Note the presence in the parameter space of a region in which the net sediment transport is negative, indicating a tidally-averaged total transport from inlet II to inlet I. Increasing 2 results in a decrease of net transport from inlet II to inlet I on the left of the thin black lines, that indicate a zero transport, or an increase of net transport from inlet I to inlet II on the right of these lines.

Constant-Width System Subjected to All Forcings
In this section, all hydrodynamic forcings are taken into account. They consist of M 2 and M 4 constituents prescribed at both inlets and a tidally-averaged water transport, prescribed at inlet II. The complete morphodynamic equilibrium condition up to orders O (ϵ 2 , ϵγ) has to be solved. It reads: with the transport terms on the left-hand side given in Equations 20-22. All parameter values are taken from Table 1, except for the M 4 phase at inlet II, 4 , which is varied between −180° and 180°. Only linearly stable morphodynamic equilibria are discussed below.
The sensitivity of morphodynamic equilibria to 4 is illustrated in Figure 9a, where the water depth is shown as a function of distance to inlet I (horizontal axis) and 4 (vertical axis). The total sediment transport of the stable morphodynamic equilibria is shown in Figure 9b as a function of 4 , illustrating the sensitivity of both direction and magnitude of the total transport to 4 . In Figure 9c, the equilibrium water depths are shown for some selected values of 4 . Focusing on −121° (for the other phases, see Text S4 in Supporting Information S1), the total sediment transport, is split into its four main contributions in Figure 9d. The two diffusive contributions, < 00 diff > and < 00 topo > , are only significant near inlet I. The sediment transport related to the internally-generated overtides, < 20 adv > , is mostly directed from inlet II to inlet I, whereas the transport due to external overtides, < 11 adv > , occurs in the largest part of the double-inlet system and is directed from inlet I to inlet II. The advective contributions can be further dissected in different components (see Text S4 in Supporting Information S1). It turns out that only two contributions dominate < 20 adv > , namely the transport of tidally-averaged concentration by the tidally-averaged flow < 10 >< 00 > and the correlation of 00 1 10 1 . Similarly, there are two main contributions to < 11 adv > , that is, 00 1 01 1 and 01 2 00 2 .

Varying-Width System Subjected to All Forcings
To illustrate the influence of width variations on the morphodynamic equilibria, the width is varied as, where c 0 is a parameter that controls the width variation. This width distribution is chosen to be able to systematically vary the width when moving away from the entrances, a geometric feature often observed in tidal Increasing (decreasing) c 0 results in a double-inlet system with a width that increases (decreases) toward the middle of the inlet system as shown in Figure 10a. The associated equilibrium water depths are shown in Figure 10b, while Figure 10c reports the equilibrium water depth as a function of c 0 (vertical axis) and of location in the double-inlet system (horizontal axis). Figure 10d shows the total sediment transport as a function of c 0 . This transport is directed from inlet II to inlet I for a rectangular inlet system (c 0 = 0), and decreases for increasing width variation. If the width in the middle of the inlet system is approximately 50% larger than at the inlets, the total transport vanishes. For even larger width variations, the transport is directed from inlet I to inlet II.
The various contributions to the transport can again be split into its different contributions (see Text S4 in Supporting Information S1). For small enough c 0 (i.e., c 0 < 0.5), the magnitude of < 20 adv > is larger than that of < 11 adv > . Since < 20 adv > is directed from inlet II to inlet I, the total transport is negative. When c 0 is increased, the relative importance of < 11 adv > increases, resulting in a total transport from inlet I to inlet II for large enough c 0 .

Morphodynamic Equilibria
The results presented in Section 4 suggest that, the number and (linear) stability of morphodynamic equilibria in double-inlet systems depends sensitively on the forcing conditions and model geometry. For a rectangular geometry, with the water motion forced by one tidal constituent, the full bifurcation diagram indicates that for most parameter values, either one unique stable equilibrium or no morphodynamic equilibrium may exist in a system with the two inlets connected together. However, in a small part of the parameter space, more than one stable equilibrium was found.
Qualitatively, these observations concerning the number of morphodynamic equilibria and their stability are consistent with model results obtained using the modeling approach employed by van de Kreeke et al. (2008) and Brouwer et al. (2012). In these models, only the tidal inlets are morphodynamically active, while the prescribed bathymetry in the back-barrier basin is assumed uniform and constant in time. In particular, van de Kreeke et al. (2008) investigated the influence of a topographic high on the existence and stability of morphodynamic equilibria. They found that, depending on the height of this topographic high, no equilibrium, one unique stable equilibrium or two stable equilibria could exist. On the other hand, Brouwer et al. (2012) showed that the existence of these equilibria depends sensitively on the amplitude and phase of the M 2 tide at the seaward side of the tidal inlets. Even though these results seem to confirm the findings of Section 4 (i.e., sensitivity of the number and the stability of morphodynamic equilibria to geometry and forcing), it should be realized that the model formulations are quite different and hence that results obtained with the different approaches are difficult to compare.

Comparison With Results of a Complex Numerical Model
In this section, we compare the cross-sectionally averaged water motion obtained with the present model, with the results of a state-of-the-art process-based numerical model, and the modeled equilibrium bathymetry with that observed in the field. The numerical model used to calculate the water motion is the General Estuarine Transport Model (GETM), a three-dimensional model that solves the hydrodynamic equations using a finite difference approach. Effects of drying and flooding are included. For an extensive discussion of the model features and the application site we refer the reader to Duran-Matute et al. (2014). The planform geometry and the bathymetry of the double-inlet tidal system were obtained from field observations and are shown in Figure 2a. In this figure, the Marsdiep inlet is denoted as inlet I and the Vlie inlet as inlet II. The Wadden islands and mainland are colored in white, whereas the water depth is shown in color scale.
The width distribution as a function of distance to inlet I to be used in the present model were obtained by smoothly connecting points with a water depth of 1 m indicated with triangles in Figure 2a. Using these lines, the basin centerline was constructed. This line was parameterized by the spatial coordinate x, which starts at inlet I (x = 0) and ends at inlet II (x = 59 km). The width as a function of x is defined as the length of the lines perpendicular to the basin centerline. Furthermore, the depth is obtained by averaging the observed depths over the width.
The resulting width profile is shown in Figure 11a. The double-inlet system is characterized by a small width near the seaward sides, but also around 31 km. At this location, the width is restricted by the presence of the tidal divide between the Eierlandse gat and the Marsdiep-Vlie system on the northern side, and the mainland in the south. This width profile has then been used to compute the values of the width-averaged depth at both entrances, as well as of amplitudes and phases of the tidal constituents (see Table 1). Figure 11b shows the width-averaged amplitudes of the M 2 and M 4 tidal constituents obtained with the present equilibrium model and computed from the GETM model. The overall trend appears to be well captured. Furthermore, there is a good correspondence with the amplitudes computed with the present model and those presented in Ridderinkhof (1988) and Hepkema et al. (2018), in which the hydrodynamics of the Marsdiep-Vlie system was investigated through cross-sectionally averaged models. The equilibrium water depth is shown in Figure 11c. There is a good qualitative comparison between the observed and modeled equilibrium water depth. In particular, the dramatic depth increase when moving from inlet I a few kilometers into the basin is well-captured. When moving toward inlet II, it is observed that the distribution of the depth variations is well-captured, but the water depth in the model is typically overestimated with respect to observed data.In Figure 11d, the various transport contributions in morphodynamic equilibrium are shown. It is observed that, the total sediment transport is from inlet I to inlet II. In the first few kilometers, the dominant transport contribution is given by topographic diffusion and is directed from inlet I to II, whereas all other contributions are in the opposite direction. The topographic diffusive sediment transport contribution is directly related to the relatively fast increase of the water depth when moving into the basin. When moving further into the basin, the water depth decreases again and topographic diffusion changes sign. However, there is still a net transport toward inlet II due to internally generated advection and diffusion that changed sign as well. Moving even further toward inlet II (x > 24 km), the transport due to the externally-prescribed overtide dominates, resulting in a net transport in the direction of inlet II. The direction of transport at inlet I, the Marsdiep Inlet, is in agreement with  Table 1) are shown in panel (d). Total transport, < 00 diff > , < 00 topo > , < 20 adv > , and < 11 adv > are indicated by the purple, blue, orange, green, and red, respectively.
19 of 23 the transport direction found in Sassi et al. (2015). In our model setup, this directly implies that in the second inlet, the Vlie, sediment should be exported. In Sassi et al. (2015), the median value of the sediment transport at the Vlie is close to zero, with both import and export found, depending on the specific forcing conditions. However, in the model setup of Sassi et al. (2015), the Wadden Sea is modeled as a multiple-inlet system, allowing for sediment transport over watersheds to adjacent basins that are not taken into account in the model schematization used in this paper.
Even though the main characteristics of the water motion (from a refined numerical model), and the main bathymetric features (from observations), are reproduced qualitatively well, there are still some differences. The main differences concerning the hydrodynamic quantities are found in regions with extensive tidal flats and multiple channels. In this study, variations in flow conditions over the cross-sectional area are not explicitly taken into account. It is possible to include these effects by modeling parametrically the effects of mass storage and momentum sinks (Friedrichs & Aubrey, 1994;Hepkema et al., 2018). However, to capture these effects dynamically, the existence and stability of morphodynamic equilibria in a two-dimensional (depth-averaged) model have to be studied. The results obtained with the cross-sectionally averaged model are essential input for models in which the above effects are incorporated explicitly (Boelens et al., 2021;Dijkstra et al., 2014). Furthermore, in reality  Sassi et al., 2015). These transport processes can only be captured by extending the present model to a morphodynamic model for multiple inlet systems (Reef et al., 2020;Roos et al., 2013).
Also, the width-averaged equilibrium bathymetry obtained with the present model exhibits the largest deviations with respect to the observed data in regions with extensive tidal flats. There, the effects of wind and waves can play an important role (de Swart & Zimmerman, 2009;Marciano et al., 2005). Furthermore, the possible import of sediment due to littoral drift along the coast of the sea is neglected in the present study. This transport mechanism was shown to be important by van de Kreeke et al. (2008) and Roos et al. (2013).

Conclusions
An idealized model has been developed to systematically investigate the existence of cross-sectionally averaged morphodynamic equilibria in double-inlet systems. A morphodynamic equilibrium is defined by the condition for which the bottom does not evolve anymore, requiring the tidally-averaged sediment transport to be spatially uniform. The model consists of the cross-sectionally averaged shallow water equations, a width-averaged and depth-integrated advection-diffusion equation for the suspended sediment dynamics and a width-averaged bed evolution equation.
The morphodynamic equilibria have been computed by using a so-called bifurcation method. This approach allows to obtain the equilibrium bed profiles (and associated water surface elevation, velocity and concentration fields) directly for given values of the relevant parameters, of the tidal forcing constituents, of the water depth at the seaward boundaries and the system geometry, without resorting to time-integration techniques.
Considering a double-inlet system with a constant width and forced by a M 2 tidal constituent at the seaward boundaries, it was shown that, depending on the prescribed M 2 amplitudes and phases, no equilibrium, one equilibrium or more than one morphodynamic equilibrium configurations exist. In the absence of a morphodynamic equilibrium, the water depth vanishes somewhere within the tidal basin and the double-inlet system reduces to two uncoupled single-inlet systems. In case of multiple morphodynamic equilibria, most often two equilibrium solutions are found, one linearly stable and the other one linearly unstable. However, for very specific regions of the parameter space, more than one stable equilibrium can be found, with typically one equilibrium having a much smaller minimum water depth than the other. When only one morphodynamic equilibrium was found, the equilibrium was always linearly stable.
The location of the watershed, the local depth at the watershed, and the total sediment transport are used to characterize the morphodynamic equilibria, since these properties are sensitive to changes in tidal forcing. Typically, the watershed tends to get closer to the inlet with a larger tidal amplitude. The total transport is usually directed from the inlet with the largest tidal amplitude to the one with the smallest tidal amplitude, with a weak dependency on the phase difference between the tidal forcing at the two inlets.
The model results have been compared with observations by forcing the system with the first two tidal constituents and considering a tidally-averaged water velocity at one of the inlets. Furthermore, the observed width distribution along the basin has to be prescribed in the model. Taking parameter values representative of the Marsdiep-Vlie system one stable morphodynamic equilibrium was found. To assess the robustness of this equilibrium, the phase of the M 4 tide at one of the inlets and the width profile of the double-inlet system were systematically varied. The resulting bathymetry and direction of sediment transport were found to strongly depend on these parameters. This finding suggests that information concerning morphodynamic equilibria for double-inlet systems with characteristics different from those of the Marsdiep-Vlie system can only be obtained by dedicated numerical experiments to explore the parameter space.
Taking the large-scale width variations observed in the Marsdiep-Vlie system into account, the main characteristics of the observed width-averaged bathymetry in this double-inlet system were qualitatively reproduced. The large water depth near the Texel Island was found in the model results, as well as two shallower regions, divided by a deeper part, when moving toward the Vlie inlet. Even though the residual water transport is directed from the Vlie inlet to the Marsdiep inlet, the model predicts a sediment transport in the opposite direction. This transport from the Marsdiep to the Vlie system is also found in much more refined numerical models. Furthermore, the order of magnitude of this transport, predicted by the numerical model, is well-reproduced by the model presented in this paper.