The Otzumer Balje embayment ( Figure 1 ), which is located in the mesotidal southern German Bight (North Sea), is used as an example. There are three major sources of river runoff, the Ems to the west and the Weser and Elbe to the east. Moreover, terrestrial fresh water enters the Wadden Sea from the south via four floodgates at a rate of 80×10⁶ m³ a^-1 (Rupert et al., 2004) at the surface and via groundwater flows in an unknown amount. The largest freshwater contributor to the Otzumer Balje system is the flood gate Neuharlinger Siel located on the coast opposite the island of Spiekeroog (Kölsch et al., 2003). Therefore, a horizontal salinity gradient stretches along the tidal channel toward the open ocean in the north.

Acoustic Doppler current profiler (ADCP) measurements in a cross-section of Otzumer Balje (Valle-Levinson et al. (2018), see their figure 7 are used to validate the numerical model that studies the mechanism responsible for the secondary circulation. The measurements were collected in 2011 along a ~5 km-long transect between Janssand and Spiekeroog (Otzumer Balje, see Figure 1 ) during neap tide on May 11 and 12 (three tidal cycles), and during spring tide between May 17 and May 19. In total, R/V Otzum traversed the inlet 245 times (136 times during neap tide and 209 times during spring tide), maneuvering at speeds between 2 and 2.5 m s^-1 while a downward-pointing 1200 kHz Teledyne RD Instruments ADCP recorded a velocity profile every 0.8 s. Ensemble averages were taken every 9 profiles, yielding a mean horizontal resolution of ~ m; with a bin size of 0.25 m (Valle-Levinson et al., 2018).

Hydrodynamic simulations are carried out for the area of the German Bight ( Figure 1 ) from January 1 to July 11, 2011, using the Semi-implicit Cross-scale Hydroscience Integrated System Model (SCHISM) (Zhang et al., 2016), which is a derivative product of the Semi-implicit Eulerian-Lagrangian Finite Element (SELFE) model (Zhang and Baptista, 2008) that solves the Reynolds-averaged Navier–Stokes equations on an unstructured grid under the application of the hydrostatic and Boussinesq approximations. The semi-implicit time step allows numerical stability and efficiency. Moreover, the use of a higher-order Eulerian-Lagrangian method (ELM) to compute the momentum advection reduces constraints on numerical stability. Tracer transport, e.g., salinity, is solved with a 2nd-order total Variation diminishing (TVD) scheme.

The turbulence closure model follows the generic length scale (GLS) formulation of Umlauf and Burchard (2003) as expressed in the k–ϵ parameterization. The parametrization was chosen due to its widespread use and because it was among the most tested closure models. A deeper analysis of the simulated quantities such as turbulent kinetic energy (TKE), dissipation and velocity shears (not shown here) illustrates a reasonable performance. This is further supported indirectly by a suspended sediment modeling study (Stanev et al., 2019), which relies on the performance of the closure model. An important consideration for this study involving the shallow Wadden Sea is the natural treatment of wetting and drying areas. Since most of the Wadden Sea area (such as the tidal basin of the studied channel cross-section) is subject to periodic inundation and drying throughout the tidal cycle, capturing these processes is particularly important for the replication of nonlinear tidal transformations. This was demonstrated in detail by Stanev et al. (2003) for the East Frisian Wadden Sea.

The model setup for the German Bight is based on the model that was presented and validated in Stanev et al. (2019). The mesh size is scaled up to 200 m near the shore of the inlet and the resolution increases up to 50 m in the area inside Otzumer Balje. The horizontal mesh consists of 438k nodes connected within approximately 1 million triangles, while the vertical dimension is resolved with 21 terrain-following sigma coordinates. To facilitate a comparison with the ADCP observations in Otzumer Balje, the model results are mapped from the nodes of the triangular grid cells onto the corresponding transect using inverse distance-weighted horizontal interpolation.

The initial forcing and hydrodynamic open boundary forcing of the German Bight model are interpolated from the 3.6 km Hereon operational Geesthacht Coupled Coastal Model System 35 (GCOAST35)-Nucleus for European Modeling of the Ocean (NEMO) model setup for the northwest European Shelf [for the application and validation of GCOAST, see, e.g., Bonaduce et al. (2020) and Staneva et al. (2021)]. The forcing encompasses sea surface height, three-dimensional velocities, and temperature and salinity at temporal resolutions of one, three, and six hours, respectively.

The atmospheric forcing is derived from the German Weather Service (DWD) Consortium for Small-scale Modeling (COSMO)-Europe (EU) model with an hourly resolution for the atmospheric pressure, 10 m wind, 2 m air temperature, specific humidity and solar radiation (https://opendata.dwd.de/). River discharge is applied for the Ems, Weser and Elbe Rivers from daily observations provided by the German Waterways and Navigation Administration (WSV) (https://www.kuestendaten.de/DE/Startseite/StartseiteKuestendatennode.html). The smaller Eider River is introduced with its annual mean discharge.

To quantify the dominant forcing of the secondary circulation at subtidal timescales, the tidally averaged momentum equation is analyzed in the lateral direction, as has been done in other studies (Geyer, 1993; Lacy and Monismith, 2001; Nidzieko et al., 2008; Zhu et al., 2017; Cui et al., 2018; Chen et al., 2019). In this study, we computed each forcing term for the transect across the tidal inlet (the red line in Figure 1 ) with the expression given in sigma coordinates (Cheng, 2014):

Here, the relative depth is σ=(z-η)/D in which D = H + η is the total depth, H is the absolute value of local tidal mean depth and η is the free water surface that deviates from the tidal mean. The variables u and v are the longitudinal (x) and lateral (y) velocity components, respectively, in the σ-coordinates ( Figure 1 ). The positive x-axis is directed in the landward flood direction, the positive y-axis is directed to Spiekeroog, the positive σ-axis is directed upward with the origin at the water surface. The vertical velocity component w(m s^-1) is defined as w = Ddσ/dt with t denoting time. The overbar ( · ¯ ) and prime ( · ′ ) denote the tidal mean and tidally varying part of a variable, respectively. Note that Z ¯ = 1 / H 2 and Z' = 2η/H ³ are obtained by estimating 1/D ² to the first-order, using Taylor expansion. The Coriolis parameter is f = 10^-4 s^-1, and the gravitational constant of acceleration is g = 9.8 m s^-2. The water density is ρ and has a reference value of ρ ₀ = 1000 kg m^-3. In baroclinicity term, σ ˜ is a dummy variable. The radius of curvature, R was a magnitude of approximately 4 km, with the sign changing from positive (when the bend turns to the positive transverse direction) during flood tide to negative during ebb tide because of the S-shape of the tidal inlet Becherer et al. (2015). The coefficient A_v denotes the vertical eddy viscosity, which is obtained from the k–ϵ model output.

On the left-hand side of Eq. 1, the first term is the friction term. This term is the response of the fluid to the local derivative of inertia, the Coriolis acceleration and other accelerations represented by the remaining terms on the left-hand side of the equation: the centrifugal acceleration (CFA), the baroclinic pressure gradient, nonlinear advection, and eddy viscosity-shear covariance (ESCO). Note that the tidally averaged local acceleration remains in Eq. 1 due to asymmetric tidal currents of flood and ebb. The right-hand side of Eq. 1 represents the barotropic (surface slope-driven) pressure gradient, which balances the terms on the left-hand side of the equation.

The next step is to derive the equations that govern residual currents induced by individual forces. The method follows that of Cheng (2014) and Chen and de Swart (2018), but with a focus on the across-channel direction. The tidally averaged conservation of mass is

The barotropic residual flow is denoted as v ¯ BT , while the non-barotropic is decomposed into four constituents: the flow due to Coriolis deflection of longitudinal residual current ( v ¯ C ), the baroclinic (density-driven) flow ( v ¯ D ) the ESCO-induced flow ( v ¯ ESCO ) and the advection-induced flow ( v ¯ ADV ). The inertia term is usually negligible (which will be clear in the next section), and thus not included in the non-barotropic constituents. Therefore, the residual currents and the corresponding residual water elevations read

Substitution of Eq. 3 into Eq. 1 and Eq. 2 yields the equations for the across-channel residual current due to individual forces:

Here, the subscript i represents an individual component of the residual flow. Integrating Eq. 4 vertically twice and applying zero stress at the free surface and no-slip bottom condition yields

with

In the across-channel direction, it is considered that the net water transports are barotropic, i.e., F _BT=0. Thus, integrating Eq. 2 across the section from the side boundary (–B) yields the barotropic transport, which can be computed numerically from the model output:

where v ¯ * is the tidally averaged across-channel velocity prescribed from the model at each grid point. Substitution of Eq. 5 into Eq. 7 yields the lateral barotropic water-level slope

Considering that the ratio of width to length of the tidal inlet ( Figure 1 ) is small, the continuity equation regarding flow components due to non-barotropic forces is simplified as (Cheng, 2014)

Hence, applying boundary conditions at the surface and the bottom to the corresponding reduced momentum equations (Eq. 4) and the continuity equations (Eq. 9), the expressions p_i corresponding to Coriolis (“C”), density gradient (“D”), curvature (“CFA”), ESCO mechanism (“ESCO”) and the advection (“ADV”) are

where "< · >" denotes taking vertical variation of the forces.

The model has been extensively validated in Stanev et al. (2019) at the scale of the German Bight. Hence for the model validation, the observation transect (see Figure 1 for the location) is taken to compare the modeled residual currents with observations. As shown in Figure 2 , the spatial patterns of the observed streamwise residual current are well represented by the model. In both neap tide and spring tide, the residual current is inward (from the North Sea to the Wadden Sea) with a maximum of 6 cm s^-1 in the middle of the transect and outward near the edge of the western and eastern shoals. The direction change occurs at approximately 0.4 km and 0.9 km from the western boundary. However, differences exist between measurements and the simulation. The maximum inflow is observed inside the water column near the deep trough of the channel, at ~0.7 m from the western boundary, while the model predicts the maximum current at the water surface near ~0.5 m. Moreover, the model overestimates the outflow, especially over the western shoal from 0 to 0.3 km.

For the lateral component, the model is able to represent major features of the observations of the transect. In neap tide, a two-layer current is observed in the main channel of the transect, i.e., the lateral distance between 0.6 km and 1.0 km ( Figure 2 ). The lateral flow is mainly to the left (looking seaward) at 0.6 km, showing a maximum in the middle of the water column. Leftward flow converges with rightward flow at approximately 0.3 km ( Figure 2 , upper left), where the longitudinal current reverses direction. In spring tide, convergence of lateral flow is found at the same place. However, east of this location, from 0.3 to 1.0 km, the surface lateral flow reverses to rightward. These features are captured by the model. Regarding model-data differences, in neap tide model shows a larger leftward current near the surface of the right side of the transect (0.6~1.4). The observed leftward currents are 2~3 m s^-1 while the model flows are 4 ∼ 6 cm s^-1. During spring tide, the modeled rightward current near the water surface is larger than the in-situ data and extends deeper in the water column. The measurement shows a rightward current with a magnitude of 1 - 2 cm s^-1 hear the bed on the right shoal (0.8 ~ 0 km), while the modeled current is leftward. Moreover, the currents over the slope near the right-hand boundary (1.0 ~ 0.4 km) are overestimated by the model.

Considering the subtidal secondary circulation is on a scale of 2~3 m s^-1 validating the model performance in the cross-channel plain is a challenge. Minor errors in, e.g., bathymetry and measured velocities, may cause differences. Despite the model and observations differences, the basic patterns in both neap and spring tides remain similar. This comparison demonstrates the reliability of the model to study lateral dynamics over the “analyzed transect” (the red transect in Figure 1 ), which maximized the longitudinal flows following the principal axis (perpendicular to the “analyzed transect”).

To quantify the degree of stratification, the potential energy anomaly (Θ) is computed across the analyzed transect:

where ρ ^ denotes the vertical mean water density. Figure 3 suggested that destratification is stronger during spring tide than during neap tide. Stratification is strongest at the end of neap ebb tide. During spring tide, the water column is mostly well-mixed. A short period of stratification is observed (with Θ being approximately 1~5 J m^-3 toward the end of flood, when the cross-channel current is close to zero at high water. The late flood straining observed in the model is consistent with that reported by Becherer et al. (2015), who applied a time-dependent dynamic equation to the potential energy anomaly (Burchard and Hofmeister, 2008). The overall stratification at the Otzumer Balje inlet transect is weak, with a maximum Θ of 4 J m^-3 during neap tide and a maximum of 2 J m^-3 during spring tide.

Figure 4 shows the along- and cross-channel transect velocities on the analyzed transect at four phases of a tidal cycle during neap and spring tides. The largest along-channel tidal currents (the peaks of both flood tide and ebb tide) occur at the surface on the right side (looking seaward) of the thalweg. Phase differences are also observed between the currents at shoals and those in the deep channel (e.g., end of the ebb and flood during neap tide), as well as at the surface and bottom (end of the flood tide during spring tide).

Different spatial patterns arise in the cross-channel flow at different tidal stages from neap to spring tides. For example, at the end of ebb during neap tide, a two-layer secondary flow develops in the deep channel, with currents directed toward the left shoal near the surface and a return flow at the bottom ( Figure 4A ). At peak flood, the tidal current diverges at approximately 0.45 km ( Figure 4B ). On the right shoal the maximum cross-channel current is located near the bottom and decreases toward the surface. At the end of the flood tide, the current near the surface of the right shoal reverses, while the current in the lower layers remain toward the right ( Figure 4C ). A anti-clockwise rotation develops akin to that at the end of ebb tide. Later, in the ebb phase, the tidal current flows from the left side (looking seaward) and converges with the leftward current from the right shoal at 0.45 km. The current structure at the ebb peak reverses compared to that at the flood peak ( Figure 4D ).

The cross-channel current pattern during spring tide is similar to that during neap tide at the flood and ebb peaks ( Figures 4G, I , comparing with b and d), but the velocities are larger during spring tide due to the stronger external tidal forcing. Spring-tide flood and ebb ends, however, show distinct structures compared to that during neap tide. At the end of ebb tide, a leftward current is observed across the transect as the velocity increases toward the bottom ( Figure 4F ). At the end of the flood tide, a two-layer flow appears with the opposite sign to that observed during neap tide ( Figure 4H ).

During both tidal periods, the along-channel residual current consists of inflow along the left side of the channel and outflow on the right side ( Figures 4E, J ). In the cross-channel direction, the residual currents show a two-layer structure both in the neap and spring tides, but with opposite signs. Considering that the inlet is predominately well-mixed, the lateral density gradient generated by differential advection is weak. The density difference over the entire cross-section has a magnitude of ~10^-2 g m³ in both neap and spring tide. Further analysis reveals that in both periods, the leftward current is much larger than the rightward current. The difference is that during neap tide, the strong leftward current is near the surface while in spring tide it is close to the bottom. Integrating lateral residual currents over the entire water column yields a net water transport from the right shoal to the left shoal for both neap and spring tides (not shown).

The individual terms of the residual lateral momentum budget (see Eq. 1) are computed on the transect analyzed (see Figure 1 ). Figure 5 compares the depth-averaged acceleration terms with the barotropic pressure gradient. In both neap and spring tides, the sum of all accelerations (left-hand side of Eq. 1) approximately equals the barotropic pressure gradient (right-hand side of Eq. 1) along the transect. Advection is dominant on the vertically averaged subtidal momentum budget, implying a small contribution to the net transverse transport due to the other accelerations. This result is consistent with the assumptions made for the decomposition method in the previous section. Differences between the sum and the barotropic pressure gradient are caused mainly by numerical errors from interpolating the model output to the grids on the transect analyzed.

Individual acceleration terms of the momentum budget as a function of depth and position, are shown in Figure 6 . During neap tide, the local derivative (inertia) is the smallest term among all contributions to the residual lateral momentum budget. The baroclinic pressure gradient induces clockwise accelerations (looking seaward) on the slopes and anti-clockwise accelerations over the main channel ( Figure 6B ). The contribution of baroclinicity reaches a value of ±1 × 10⁻⁵ m s^-2, which is one order of magnitude larger than the Coriolis accelerations ( Figure 6C ). The ESCO mechanism shows a three-layer structure, with leftward accelerations both near the surface and bottom, and a rightward acceleration in between. The ESCO magnitude is as large as that from the baroclinic pressure gradient ( Figure 6D ). The CFA is clockwise on the left slope, while it is anti-clockwise in the thalweg and the right slope ( Figure 6E ). In comparison, lateral advection dominates the contributions to the depth-dependent lateral momentum budget ( Figure 6F ).

During spring tide, the local derivative remains small ( Figure 6I ). The lateral baroclinic pressure gradients drive clockwise-rotating accelerations over the main channel and the left slope, and a relatively small anti-clockwise rotation on the right slope ( Figure 6J ). The Coriolis acceleration is comparable to the baroclinicity but shows a single rotating cell directed toward the left near the surface and toward the right near the bottom ( Figure 6K ). The ESCO mechanism causes a single rotating cell ( Figure 6L ), in which the acceleration is from left to right near the surface and to the right underneath. Its contribution to the momentum budget is much larger than that of baroclinicity ( Figure 6J ). This result reveals that the ESCO mechanism is more relevant than the baroclinicity in predominantly well-mixed coastal waters. This finding is similar to that of Cheng et al. (2013), who focused on the along-channel direction. Remarkably, compared to the neap tide, CFA reverses in the thalweg and the right slope ( Figure 6M ). It is smaller than the ESCO but has a similar spatial structure. In spring tide, advection is similar to that of the neap tide, regarding both the spatial pattern and the magnitude, except for the near-surface layers of the deep channel.

The sum of the individual acceleration terms ( Figures 6G, O ) is equivalent to the ‘Friction’ (with a minus sign) ( Figures 6H, P ). With a perfect model performance and momentum budget decomposition, the ‘Sum’ and the negative ‘Friction’ terms should be identical (in balance). Overall, the differences in this application are small.

Figure 7 further shows the spatial structure of the cross-channel residual flow induced by the individual accelerations displayed in Figure 6 . The barotropic flow v _BT. ( Figures 7A, I ) and the flow induced by the non-barotropic part of advection, v _ADV ( Figures 7F, N ), are dominant residual components. They reach more than ~4 m s^-1. The barotropic flow diverges at 0.7 km in neap tide and at 1 km in spring tide. It also shows a convergence at the left shoal during spring tide. The baroclinic flow and the Coriolis flow are mostly less than 0.3 cm s^-1 except for v _D in the main channel, where it reaches 1 cm s^-1 The difference between the two tidal periods is that the ESCO flow v _ESCO changes the spatial structure and becomes more important from neap to spring tide. Over the left side of the deep channel (0.4 km), the ESCO flow direction reverses ( Figures 7D, K ). The maximum velocity exceeds 3 m s^-1 but the surface in spring tide. In neap tide, the CFA-induced flow is similar to that driven by advection, but with a small magnitude. The CFA-induced flow reverses its direction over the thalweg and the right slope during the spring tide. The flow is rightward near the surface and leftward near the bottom, which is similar to the two-layer structure of the ESCO flow, but with a magnitude of approximately 50% lower. Summing all individual residual flow components yields flow with spatial patterns ( Figures 7G, O ) similar to that obtained by averaging the modeled lateral flow in tidal periods ( Figures 7H, P ), demonstrating a valid performance of the decomposition method.

The lateral residual flow due to advection, i.e., the summation of the barotropic part v _BT ( Figures 7A, I ) and the non-barotropic part v _ADV ( Figures 7F, N ), is shown in Figure 8 . It has a comparable structure and magnitude as the v _m in the neap tide, revealing its dominance in subtidal secondary circulation. In spring tide, the flow driven by advection presents a similar structure as that during the neap tide, e.g., currents leftward near 0.1 km and rightward near 1.1 km with a velocity of approximately 4 cm s^-1. Moreover, between 0.3 km and 0.9 km the flow is leftward near the surface and rightward near the bottom, despite the difference in flow intensity (4 cm s^-1 n neap tide and 2 cm s^-1 in spring tide). This implies that the change in advective process from neap to spring tide is small.

Comparing the advective flow with the mean flow (v _m) clearly demonstrates the contribution of the ESCO mechanism. The clockwise circulation induced by ESCO acts against advection and reverses the lateral flow near the surface ( Figure 7O ). Meanwhile, the ESCO flow below the surface enhances the leftward flow. The baroclinic flow (v _D) as a pattern similar to that of the ESCO flow but with a magnitude 2~3 times smaller ( Figure 7I ). The CFA also plays a role in subtidal secondary circulation, especially during spring tide, which together with ESCO, contributes to the reversal of the subtidal secondary circulation. The analysis reveals that the reversal of the flow due to the ESCO mechanism ( Figures 7D, K ) plays the most important role in the change of the subtidal secondary circulation pattern between neap and spring tide ( Figures 7H, P ).

The decomposition analysis is also applied to the observation transect (the black line in Figure 6 ), where the barotropic transport is mainly driven by the advection and could be isolated from the other processes. On the depth-dependent momentum budget, the advection also plays a dominant role for both neap and spring tides. However, similar to the findings on the transect analyzed, the ESCO mechanism explains the clockwise circulation observed at 0.5 km of the observation transect during spring tide (see Figure 2 ). Figures illustrating structures of the decomposed flow are provided in the Supplementary Material .

To understand the changes in the ESCO-induced acceleration, as well as the flow, from neap tide to spring tide, the ESCO mechanism is further scrutinized. The ESCO mechanism (see Eq. 1) consists of four covariance components, namely, those between the tidal velocity (v'), tidal elevation ( Z ¯ ) and the tidally varying component of the eddy viscosity ( A v ′ ). Note that A v ′ changes with respect to the tidal mean, and is negative when turbulent mixing is relatively weak in the tidal cycle. To quantify the contributions of different covariances to the ESCO mechanism, each term is computed separately. The ESCO acceleration mainly results from the first term (not shown), i.e., the covariance between the tidally varying part of the eddy viscosity ( A v ′ ) and the vertical shear of the lateral tidal velocity ( ∂ v / ∂ σ ). The contributions of the other terms are one order of magnitude smaller.

Given the complexity of the correlations in ESCO, the structures of ∂ v / ∂ σ and A v ' it different phases of a tidal cycle are analyzed for both neap and spring tides. Figure 9 shows ∂v/∂σ, A v ' and their product A v ' ∂ v ∂ σ at four phases of the neap tide: end of the ebb, peak flood, end of flood and peak ebb. Because the tidal mean of the eddy viscosity is removed, negative values of A v ′ indicate a smaller eddy viscosity value than the tidal mean. This occurs from late ebb to early flood and from late flood to early ebb. At the end of the ebb tide, the velocity shear exhibits a structure with values of opposite signs with respect to the deep channel, which, when multiplied with a negative A v ' for the entire cross-section, yields A v ′ ∂ v ∂ σ with a similar structure as the shear. At flood peak, the ∂ v / ∂ σ values change signs over the left shoal and the upper layer of the right shoal, meanwhile A v ' becomes positive. Their product again yields a cross-sectional pattern similar to that of end of ebb. At the end of flood, the vertical shear has two-layer structures but displays opposite signs over the two shoals. A similar feature is observed in A v ' with respect to the deep channel, with values of different signs. This leads to a three-layer structure of A v ' ∂ v ∂ σ over the two shoals, despite a relatively weaker product compared to the other phases. At maximum ebb, A v ' has a two layer structure on the right shoal and the deep channel, with positive values in deep water columns and negative at the surface. Over the left shoal, A v ' is mostly negative (except for the section close to the left side boundary). The spatial structure of the vertical shear at this phase is negative over the right shoal and positive over the left shoal. The product of A v ' ∂ v ∂ σ reflects similar features as ∂ v / ∂ σ . The residual A v ' ∂ v ∂ σ bottom plot of Figure 9 ) has a two-layer pattern on the right shoal and the channel, while it is almost three layers over the right shoal.

A similar analysis is made for spring tide ( Figure 10 ). Both the cross-sectional distribution and the magnitude of A v ' it different tidal phases are similar to those at neap tide. The only difference occurs at flood slack, when A v ' is negative along the entire cross-section. However, A v ' ∂ v ∂ σ during spring tide shows a pattern rather different compared to that during neap tide. It is negative over almost the entire tidal cycle. This is because the velocity shear has patterns that are different from those at neap tide. For example, ∂ v / ∂ σ is mostly positive at the end of ebb while it becomes negative throughout the section near the flood peak. This pattern suggests that the maximum lateral velocity is near the bottom and decreases toward the surface. A similar situation occurs near peak ebb. The velocity shear is negative from the thalweg to the right shoal, where A v ' is positive. On the left shoal, the eddy viscosity is still smaller than the tidal mean and the vertical shear remains positive. As a result, the product retains the same sign. Taking an average over the tidal cycle yields a negative value of the product.