A semianalytical three-dimensional model is set up to dynamically calculate the coupled water motion and salinity for idealized well-mixed estuaries and prognostically investigate the influence of each physical mechanism on the residual salt transport. As a study case, a schematized estuary with an exponentially converging width and a channel–shoal structure is considered. The temporal correlation between horizontal tidal velocities and tidal salinities is the dominant process for the landward residual salt transport. The residual salt transport induced by residual circulation is locally significant, but the induced salt transport integrated over the cross section is small. The impacts of the estuarine geometry, Coriolis force, and bathymetry on the salt dynamics are studied using three dedicated experiments, in which the impact of each of these factors is studied separately. To assess the impact of width convergence, a convergent estuary without bathymetric variations or Coriolis force is considered. In this experiment, the temporal correlation between tidal velocities and salinities is the only landward salt transport process. In the second experiment, Coriolis effects are included. This results in a significant residual salt transport cell due to the advection of the tidally averaged salinity by residual circulation, with salt imported into the estuary from the left side and exported on the right (looking seaward). In the last experiment, a lateral channel–shoal structure is included while the Coriolis effects are excluded. This results in a significant landward salt transport through the deeper channel and a seaward salt transport over the shoals due to the advection of the tidally averaged salinity by residual circulation.
Estuaries are dynamically complex systems that are important from both an economical and ecological point of view. To be able to address issues related to these aspects, a good understanding of the physical processes resulting in the estuarine exchange flow is fundamental, as this flow drives long-term transport of sediment and other suspended particle matter in estuaries. Since the exchange flow is significantly affected by the salt dynamics, it is key to improve our understanding of the essential three-dimensional processes that result in the spatially and temporally varying salinity fields.
Salt transport processes can vary significantly from estuary to estuary, since estuaries are subject to different forcing conditions and have greatly varying geometries and bathymetries (Fischer 1972). The classical subtidal salt balance, in which the river-induced seaward transport is balanced by the landward transport due to the longitudinal–vertical gravitational circulation and dispersion, is extensively described in literature (see, e.g., Hansen and Rattray 1965; MacCready 2004; MacCready and Geyer 2010). In this way, both lateral processes and processes that vary on the tidal time scale are parameterized. This classical balance is challenged in partially to well-mixed estuaries where gravitational circulation is relatively weak (Jay and Smith 1990) and other processes are likely to dominate. In well-mixed estuaries, for instance, the tidal oscillatory transport due to temporal correlations on the tidal time scale between the longitudinal velocity and salinity can be the predominant salt import mechanism (McCarthy 1993; Wei et al. 2016). The tidal oscillatory salt transport can also be significant in partially stratified estuaries such as the Hudson estuary during neap tides (Chen et al. 2012; Wang et al. 2015). Even in salt wedge estuaries such as the Merrimack River, the landward salt transport is mainly associated with the tidal oscillatory processes rather than with the steady sheared gravitational circulation (Ralston et al. 2010). This implies the potential importance of considering residual salt transport mechanism due to tidal processes (Hughes and Rattray 1980; Geyer and MacCready 2014).
Apart from the importance of processes on the tidal time scale, the transverse circulation (also known as secondary estuarine circulation) can be of great importance for the landward salt transport (Fischer 1972; Dyer 1974; Smith 1977). Significant transverse circulations are observed in well-mixed, partially mixed, and stratified estuaries (Nunes and Simpson 1985; Becherer et al. 2015; Valle-Levinson et al. 2000; Lacy and Monismith 2001). These circulations are usually associated with lateral variations in bathymetry (Nunes and Simpson 1985; Valle-Levinson et al. 2000), Coriolis force (Lerczak and Geyer 2004; Valle-Levinson 2010), lateral density gradients (Huijts et al. 2006, 2011), and curvature-induced centrifugal forces (Lacy and Monismith 2001; Buijsman and Ridderinkhof 2008). Transverse circulations are also affected by the estuarine width-to-depth ratio (Schulz et al. 2015) and temporal/lateral variations of mixing (tidal and lateral straining; see Burchard et al. 2011; Becherer et al. 2015).
To investigate the influence of transverse processes on estuarine salt transport, the salt flux is usually decomposed into different contributions associated with transverse/vertical steady or oscillatory shear, using velocities and salinities from measurements (Fischer 1972; Dyer 1974; Hughes and Rattray 1980; West et al. 1990; Guymer and West 1992) or numerical models (Oey et al. 1985; Ralston and Stacey 2005; Ralston et al. 2010). This method, however, strongly depends on the quality/quantity of the measurements and the resolution and accuracy of numerical models. Moreover, it is difficult to relate the shear-induced salt transport to a specific physical process since the shearing itself can be a result of the interplay of many different processes.
To unravel the dominant salt transport mechanisms, idealized models are appropriate tools, as they are specifically developed to investigate the influence of each physical mechanism in isolation. This type of model has been successfully used to gain insights into the influence of lateral processes on water motion (Nunes and Simpson 1985; Li and O’Donnell 1997; Wong 1994; Valle-Levinson et al. 2003; Kumar et al. 2016) and lateral sediment trapping (Huijts et al. 2006, 2009, 2011). Existent idealized models aiming at understanding the salt dynamics (explicitly including processes at the tidal time scale), however, are limited to width-averaged models (McCarthy 1993; Wei et al. 2016). In these models, the coupled width-averaged water motion and salinity are dynamically resolved while the lateral processes are parameterized as diffusion. Since these models do not resolve the lateral dynamics, the relative importance of the transverse circulation on the residual salt balance are not investigated.
In this paper, we therefore develop a three-dimensional idealized model by extending the model of Wei et al. (2016). This model is specifically aiming at dynamically resolving the temporal and spatial (both longitudinal and lateral) variability of the water motion and salinity (i.e., prognostic in salinity, including lateral processes). With this model, the influence of each physical process on the residual salt balance can be systematically investigated. As a first step, the estuaries are assumed to be well mixed, that is, the top-to-bottom salinity difference is one order of magnitude smaller than the bottom salinity, and the effects of tidal straining are assumed to be small. The nonlinear system of the three-dimensional water motion and salinity is solved at the tidal time scale, employing a perturbation method to analytically obtain the vertical structures of the physical quantities and a finite-element method to calculate their horizontal variations. An iterative approach is used to calculate the coupled gravitational circulation and salinity field. Using this model, the influence of the longitudinal/transverse processes on the residual salt transport can be investigated, and their sensitivity to parameters/forcing conditions can be systematically studied.
The paper is structured as follows: The equations governing the coupled water motion and salinity (see section 2a) and the solution method (see section 2b) are introduced in section 2. In section 3, the salt dynamics for the default experiment (representative for the Delaware estuary; see section 3a) and the influences of estuarine convergence, bathymetry, and Coriolis force on salt transport mechanisms (see section 3b) are investigated. In section 4, the salt transport induced by the tidal advective diffusion (see section 4a) and gravitational circulation (see section 4b) are analyzed considering the influence of the estuarine convergence, Coriolis, and bathymetry, respectively. Conclusions are drawn in section 5.
2. Model description
a. Governing equations
To investigate the salt dynamics and transport mechanisms in tidal estuaries, an idealized three-dimensional model is developed for estuaries that are assumed to be tidally dominated and well mixed. The shape of the estuary (i.e., its geometry) is arbitrary (see Fig. 1). The seaward boundary ∂SΩ is connected with the open sea and forced by tides. A weir is assumed to be located at the landward boundary denoted by ∂RΩ, where a river flows into the estuary. The closed boundaries ∂CΩ are impermeable. The undisturbed water level is located at z = 0, and the free-surface elevation is denoted by z = η. The estuarine bed is located at z = −H(x, y), with H an arbitrary function of the horizontal coordinates (x, y).
The water motion is described by the three-dimensional, shallow-water equations, assuming hydrostatic equilibrium and using the Boussinesq approximation (Cushman-Roisin and Beckers 2011):
Here, t denotes time and U = (u, υ, w) is the velocity vector, with u, υ, and w as the velocity components in x, y, and z directions, respectively. The acceleration of gravity is denoted by g. The quantities Aυ and Ah are, respectively, the vertical and horizontal eddy viscosity coefficients. The Coriolis parameter is denoted by f, and the estuarine water density is denoted by ρ. The density is assumed to depend only on the salinity S as ρ = ρc(1 + βsS), with βs = 7.6 × 10−4 psu−1 and ρc as the background density, taken to be 1000 kg m−3. To dynamically calculate the density, the salinity equation has to be solved:
Here, Kh and Kυ are the horizontal and vertical eddy diffusivity coefficients, respectively. In this paper, the vertical eddy diffusivity Kυ is assumed to be equal to the vertical eddy viscosity Aυ.
At the seaward boundaries, the water motion is forced by a prescribed sea surface elevation that consists of a semidiurnal constituent M2 and a residual water level M0:
where and φ(x, y) are the prescribed amplitude and phase of the semidiurnal tidal elevation at the seaward boundary ∂SΩ, with σ ~ 1.4 × 10−4 s−1 (the M2 tidal frequency). The term is the prescribed residual water level at ∂SΩ, where averaged over the seaward boundary vanishes. Furthermore, the tidally averaged salinity at the seaward boundary is prescribed (denoted by Se):
where the overbar denotes a tidal average.
At the free surface, the kinematic and no stress boundary conditions are prescribed:
At the bottom, the normal velocity is required to vanish and a partial-slip condition is applied. Here, we use a linear relationship between the bed shear stress and the velocity at the top of the bottom boundary layer (Schramkowski and De Swart 2002):
The slip parameter s depends on the bed roughness and can vary from 0 in frictionless cases (free slip) to large values in strongly frictional cases (no slip). At the free surface and the bottom, the salt flux is required to vanish:
Since the main aim of the model is to analyze the salt dynamics at the estuarine scales, the boundary layers where viscous effects are dominant are not resolved [see Winant (2008) for a detailed discussion]. Instead of prescribing a velocity at each location of the closed boundaries and the river boundaries, the normal component of the water flux and the tidally averaged salt flux integrated over the depth are prescribed. At the closed boundaries, the tidally averaged water/salt flux integrated over the water depth is required to vanish due to impermeability. At the river boundary ∂RΩ, a river discharge R is prescribed, and the tidal discharge is assumed to vanish; the tidally averaged, depth-integrated salt flux through the river boundary is also required to vanish. The resulting conditions read
Here, nh is the horizontal unit normal vector pointing outwards. The depth-integrated river flux per unit width Ri is assumed to be proportional to the local water depth:
No boundary conditions for the temporally varying salinity transport are prescribed at ∂SΩ, ∂RΩ, or ∂CΩ, as the horizontal and vertical structure of the time-dependent salinity will be directly related to the tidally averaged salinity field (a detailed discussion will follow in section 2b).
b. Solution method
1) Perturbation method
The system of Eqs. (1)–(14) is solved using an asymptotic expansion of the physical variables in a small parameter ε, which is the ratio of the mean M2 tidal amplitude and the mean water depth , obtained by averaging and H over the seaward boundary ∂SΩ. To make a consistent expansion, as a first step, the governing equations are made dimensionless, and a scaling analysis is performed to order various terms in the equations with respect to this small parameter ε. As a next step, the physical variables are written as asymptotic expansions in ε:
where Γ represents any of the physical variables (η, u, υ, w, and S). The subscript denotes the order in the expansion. Finally, all terms at the same order of ε are collected and required to balance, resulting in a system of equations at each order of ε.
Using this perturbation approach, the salinity equation at leading order becomes (see details in appendix A)
Since there is no horizontal transport of salt at this order, the boundary conditions given in Eq. (14) at ∂CΩ or ∂RΩ are automatically satisfied. The remaining boundary conditions read
It is found that the leading-order salinity field S0 is time independent and vertically homogeneous (McCarthy 1993; Wei et al. 2016), that is, the tidally averaged salinity is an unknown function of the horizontal coordinates (x, y):
To get a closed formula for S0, the first-order salinity equation and the depth-integrated, tidally averaged, second-order salinity equation have to be solved. The first-order salinity equation reads
which reveals that the temporal and vertical variations of the first-order salinity S1 are forced by the advection of the leading-order salinity S0 by the leading-order horizontal velocities u0 and υ0. The leading-order velocities are driven by the M2 sea surface elevation that is prescribed at the seaward boundary. This tidal motion is independent of salinity and can be calculated explicitly by numerically solving a two-dimensional wave equation for the leading-order, free-surface elevation η0 (Kumar et al. 2016). Then, the horizontal tidal velocities can be expressed as
where C1 and C2 are analytical expressions that describe the vertical structure of the velocity fields [for details, see appendix B and Kumar et al. (2016)]. The boundary conditions for salinity equation at this order read:
Since the leading-order water motion only consists of a M2 tidal constituent, and the leading-order salinity field is time independent, the first-order salinity field consists only of an M2 tidal constituent as well. Hence, the tidally averaged salt flux at O(ε) vanishes, and boundary conditions Eq. (22a) and Eq. (22b) are automatically satisfied. By substituting Eq. (21) into Eq. (20) and applying boundary condition Eq. (22c), an analytical expression for the first-order salinity S1 is obtained:
Equation (23) shows that the horizontal structure of S1 is directly related to the gradients of the leading-order water elevation and salinity, while its vertical structure, denoted by and , can be explicitly calculated (for details see appendix C).
Next, the tidally averaged, depth-integrated salinity equation at O(ε2) is derived:
that is, the divergence of the depth-integrated residual salt transport due to all processes has to vanish at each location. The tidally averaged salt transport at this order is a result of advection of the first-order salinity S1 by the leading-order tidal flow u0, υ0 (first term), the leading-order salinity S0 by the first-order residual flow u1, υ1 together with transport due to the time-varying sea surface elevation (second term), and diffusive processes (third term).
At the closed boundaries ∂CΩ and the river boundaries ∂RΩ, the normal component of the depth-integrated residual salt flux at O(ε2) has to vanish:
while at the seaward boundary, the tidally averaged salinity S0 is prescribed [see Eq. (18a)]:
is the effective total diffusivity tensor, and F is the residual water transport vector. The entries in consist of the prescribed horizontal eddy diffusivity Kh and the tidal advective diffusion coefficients , , , and , which are induced by the temporal correlation between the leading-order horizontal tidal velocities and the first-order salinity (Wei et al. 2016). The tidal advective diffusion coefficients can be explicitly evaluated using the leading-order water motion:
Here, indicates that the real part is used. The complex amplitude of the M2 tidal surface elevation is denoted by , and the asterisk (*) denotes the conjugate of a complex variable. Here, and measure the longitudinal and lateral diffusion, respectively, due to the temporal correlation between u0 and S1; and measure the longitudinal and lateral diffusion due to the temporal correlation between υ0 and S1. Inclusion of these advective diffusivities in the effective diffusivity tensor shows that in a tidally averaged model, diffusion is generally nonisotropic. The tidal advective diffusion coefficients , , and are associated with lateral processes/variations in the leading-order water motion due to geometrical features such as estuarine convergence, Coriolis deflection, and lateral bathymetric variations. In case of no lateral variations of the tidal surface elevation and no Coriolis forcing, that is, and f = 0, these three diffusion coefficients vanish, and becomes identical to the width-averaged tidal advective diffusivity as used in Wei et al. (2016). The term Kh is the prescribed horizontal diffusivity, which parameterizes unresolved processes such as mixing induced by river bends and tidal straining.
The residual water transport in the longitudinal and lateral directions is given by
The first-order residual velocity , consists of various contributions including tidal rectification of the leading-order M2 tide (denoted by AD), density-driven gravitational circulation (GC), the stress-free surface condition (NS), river discharge (RD), and a return flow (TRF). Hence, F can be decomposed into various components:
All contributions to F, except for FGC, can be calculated explicitly (without information of the salinity field) by numerically solving a two-dimensional equation for that specific contribution to the residual free-surface elevation [see section 2b(2)], while the vertical structure is obtained analytically [see details in appendix D, section a, and Kumar et al. (2017)]. The residual contribution FGC due to the gravitational circulation,
itself depends on salinity S0 (here, and are 2 × 2 matrices, whose elements are known functions of the horizontal coordinates; see details in appendix D, section b). Hence, to obtain S0, the coupled system of equations, given by
has to be solved simultaneously; here, is the barotropic residual water transport independent of salinity. Equation (36) follows from substituting Eqs. (34)–(35) into Eq. (27), and Eq. (37) follows from requiring the divergence of the horizontal residual water transport by gravitational circulation to vanish (see appendix D, section b for details). Note that since the leading- and first-order water motion (except for gravitational circulation) can be calculated explicitly, the only unknown variables in this system of equations are the leading-order salinity field S0 and .
2) Numerical method
All physical variables of water motion and salinity will be calculated using a standard finite-element method (Gockenbach 2006). As a first step, the weak form of the equations is derived. These equations are all two-dimensional (as the vertical structure of all physical variables is derived analytically). Hence, the domain is discretized in the horizontal directions using triangles, and the unknown variables are approximated using polynomial basis functions. Then, these approximations are substituted into their weak form. The numerical solution method for the leading-order tidal motion and the first-order residual flow components (except for gravitational circulation) are discussed in detail in Kumar et al. (2016, 2017). However, as the gravitational circulation and the tidally averaged salinity field are coupled, the solution procedures as discussed in Kumar et al. (2016, 2017) cannot be used, and an iterative approach has to be taken. To eliminate any accuracy loss that results from taking derivatives of S0 and , a special iteration scheme is adopted, in which no information of the gradients of S0 or from the previous iteration step is required. This iterative scheme reads
where the superscripts k and k + 1 denote the physical variables at the kth and (k + 1)th iteration step, respectively.
The weak form of Eqs. (38) and (39) is derived by multiplying these equations with a test function V, integrating them over the domain Ω, and applying the no normal salt transport boundary conditions and integration by parts
Next, the unknown sea surface elevation related to the gravitational circulation and salinity are approximated using Lagrangian basis functions ϕj:
where the ϕj are chosen to be quadratic polynomials, which are equal to 1 at node j and 0 at other nodes, and Nf is the total number of nodes. Substituting Eq. (42) into Eq. (40) and Eq. (41), and taking V = ϕi, the following system of equations (for i = 1, 2, …, Nf) is found:
This system of equations is linear in the unknown coefficients and using a known . The iterative procedure is stopped when the relative difference between the solutions from the (k + 1)th and kth iteration step
becomes smaller than 5 × 10−4, with . To start the iteration, a prescribed salinity field [e.g., ] is used. When solving this coupled system of equations, special care has to be taken of the seaward boundary, since numerical instabilities can occur when using an arbitrary prescribed water level at this boundary (Chen and Sanford 2009); see detailed procedures in appendix D, section c.
In this paper, schematized estuaries with a simplified geometry and bathymetry are studied. The estuarine width B is assumed to decrease exponentially in the landward direction:
with B0 as the estuarine width at the seaward entrance and Lb as the estuarine convergence length. The seaward boundary is located at x = 0, and the landward boundary is located at x = L. The midaxis of the idealized estuary is located at y = 0, and the lateral boundaries are located at y = ±B/2. In this geometry, the contribution of each physical process to the along-channel residual salt transport can be obtained by integrating Eq. (24) from the left bank of the channel y = y1 to the right bank y = y2 and subsequently integrating from the weir located at x = L (where the cross-sectionally integrated longitudinal salt transport vanishes) to any longitudinal location x:
This expression shows that the cross-sectionally integrated tidally averaged salt transport due to tidal rectification (AD), gravitational circulation (GC), river discharge (RD), tidal return flow including Stokes’ drift (TRF), the stress-free surface condition (NS), tidal advective salt transport (TASF), and diffusion (DIF) balance each other, and their relative importance at each longitudinal location can be assessed individually. Here, the tidal advective salt transport represents the residual salt transport caused by the advection of tidal salinity by the tidal velocity, which is called tidal oscillatory salt flux in Lerczak et al. (2006) and Wang et al. (2015).
a. Salt dynamics for the default experiment
For the default experiment, parameters representative for the Delaware estuary are used (L = 215 km, B = 39 km; see values of all default parameters in Table 1). In this default experiment, the water motion is forced at the mouth by an M2 tide with tidal amplitude (0.75 m). A river discharge R = 72 m3 s−1 (low discharge condition) is prescribed at the landward end of the estuary (at x = L), where a weir is located (Wei et al. 2016). The bathymetry of the Delaware estuary is characterized by a deep channel in the center and shallow shoals on the sides, with the lateral depth difference between the channel and the shoals decreasing in the landward direction (Aristizábal and Chant 2013). To qualitatively capture this spatial dependency, the bathymetry is described by
with Hmax (15 m) and Hmin (3.6 m) as the maximum and minimum water depths at the mouth. The width-averaged water depth Hm is longitudinally constant (Hm = 8 m). By varying the parameter Cf, the width of the tidal flats is varied, with larger values of Cf corresponding to wider shoals and a narrower channel and smaller values to narrower shoals and a wider channel.
The vertical eddy viscosity Aυ is assumed to be linearly proportional to the local water depth (Friedrichs and Hamrick 1996), and the same is assumed for the partial-slip parameter s:
with m2 s−1 and sm = 0.039 m s−1 as the width-averaged friction parameters, which are the same as those used in Wei et al. (2016). Using these parameter settings, the leading-order tidal surface elevation obtained from the present 3D model is nearly the same as that obtained in the width-averaged model of Wei et al. (2016); the tidal surface gradients and velocities, however, change significantly by including the three-dimensional effects. The horizontal eddy diffusivity has to be prescribed and is assumed to decrease from the mouth to the landward side, indicating that more unresolved processes have to be parameterized near the mouth (such as rapid depth variations, multiple channels, complex sea–estuary interactions, and tidally varying mixing effects). These unresolved processes are assumed to be proportional to the estuarine width as
with 10 m2 s−1 and = 40 m2 s−1. An overview of all default parameters is shown in Table 1.
1) Tidally averaged salinity field and longitudinal salt transport
The tidally averaged salinity is prescribed to be 31 psu at the seaward side, and it decreases to 2 psu at ~120 km inside the estuary, as shown by the background colors in Fig. 2, where the white color indicates a salinity of ~31 psu and black indicates a vanishing salinity. Hence, the salt intrusion length Ls, which in this paper is defined as the distance between the estuarine mouth to the location where the cross-sectionally averaged salinity is 2 psu, is 120 km. This salt intrusion length agrees well with the observed salinity profile under low river discharge (Kuijper and Van Rijn 2011; Wei et al. 2016), confirming that the longitudinal dependency of the horizontal diffusivity Kh, as given in Eq. (51), is appropriate for the Delaware estuary. A lateral variation of S0 (of ~1 psu) is found, with larger tidally averaged salinities on the left shoals than those on the right, looking seaward. This lateral distribution qualitatively agrees with the observation of Wong (1994) near the mouth of the Delaware Bay. Furthermore, the top-to-bottom salinity difference (below 0.5 psu) is much smaller than the bottom salinity during the tidal cycle (not shown), justifying the well-mixed assumption.
To visualize the residual salt transport pattern, streamlines of the depth-integrated salt transport are used. Since the depth-integrated residual salt transport is divergence free [see Eq. (27)], the amount of salt transported through any area connecting the same streamlines is constant. The streamlines of the residual salt transport are shown in Fig. 2a, where solid lines indicate counterclockwise salt transport, and dashed–dotted lines indicate clockwise salt transport (looking into the paper). It is found that, in the idealized Delaware estuary, the salt is mostly transported counterclockwise, meaning the salt is transported into the estuary from the left side of the estuary and out of the estuary from the right side of the deeper channel looking seaward.
Since the salt transport between the same streamlines is constant, large distances between these streamlines indicate small (normal) salt transports, and small distances indicate large transports. As shown in Fig. 2a, the distance between the two neighboring streamlines first increases from the left bank to the midaxis (at y = 0) and then decreases sharply toward the right bank. This suggests an interesting steady-state salt transport pattern: the salt is slowly transported landwards from the left side of the mouth toward the midaxis and gets quickly transported out of the estuary on the right of the deeper channel. To evaluate the relative importance of various longitudinal salt transport mechanisms, the along-channel salt transport at each cross section is calculated using Eq. (48). As shown in Fig. 2b, the seaward salt transport due to river discharge (line with upside-down triangles) is mainly balanced by the landward salt transport due to the tidal advective diffusion (line with triangles) and prescribed diffusion (line with asterisks). Gravitational circulation (solid line) and tidal rectification (dashed–dotted line) transport salt landward in the central region of salt intrusion. The stress-free surface condition (line with crosses) and tidal return flow (dashed line) transport salt landward only near the mouth. The total amount of salt transport is zero (dotted line in 2b) by definition, as the salt content in the estuary is constant in steady state.
2) Salt transport due to each process
In equilibrium, the divergence of the total (depth integrated) salt transport is zero at every location. However, the divergence of salt transport due to each physical transport process separately is not necessarily zero, indicating a local change of the salt content due to each process. If the divergence of the salt transport due to one process has locally a negative (positive) value, it implies this process contributes to an increase (a decrease) of salt content at that location.
Figures 3a–e show the divergence of the salt transport due to diffusive processes (including tidal advective diffusion and prescribed diffusion), gravitational circulation, tidal rectification, stress-free surface condition, and tidal return flow, respectively. Local convergences of salt are shown in dark gray, and divergences of salt are shown in white. The lines in Figs. 3a–e show the direction of salt transport due to each process. The solid (dashed–dotted) lines in Figs. 3b–e indicate that the salt is transported in a counterclockwise (clockwise) direction. Note that these lines are not the streamlines for the residual salt transport but for the residual water transport, since the divergence of the residual salt transport for each process individually does not vanish.
The diffusive processes transport salt from the shoals toward the deeper channel and subsequently toward the upstream (Fig. 3a). The diffusive processes result in an increase of salinity on the right side of the estuary and a decrease on the left. The convergences and divergences of the diffusive salt transport are significant, especially in the central region near the deeper channel, with a salt accumulation of up to 1 × 10−4 psu m s−1. Moreover, the salt accumulation on the right exceeds the salt decrease on the left side, resulting in an increase of salt content integrated over a cross section. This increase in salt content agrees with the decreasing longitudinal salt transport by diffusive processes, as shown in Fig. 2b.
The two-cell streamlines of gravitational circulation show that the density-driven flow transports salt landward into the estuary through the deeper channel, reverses direction in the central region, and then transports salt out of the estuary over the shoals (see Fig. 3b). The salt transport induced by gravitational circulation is asymmetric with respect to the midaxis due to the Coriolis deflection; the landward salt transport is deflected toward the left bank, and the seaward transport is deflected toward the right bank. Even though the longitudinal salt transport (integrated over a cross section) induced by gravitational circulation is much smaller than that induced by the diffusive processes, the gravitational circulation makes a significant contribution to local changes in salinity. The gravitational circulation induces a significant decrease of salt on the shoals, and a large increase of salt in the deeper channel, which reaches up to 5 × 10−5 psu m s−1 in the central region of salt intrusion.
The salt transport induced by residual flow due to tidal rectification, stress-free surface condition, and tidal return flow all show a similar pattern, with salt transported into the estuary over the shoals and out through the deeper channel (see Figs. 3c–e). The tidal rectification results in an accumulation of salt on the left side of the deeper channel, accompanied with a strong decrease of salt on the right. The stress-free surface condition and tidal return flow, however, contribute to a decrease of salt in the deeper channel and an accumulation of salt on the shoals.
b. Influence of estuarine geometry, Coriolis, and bathymetry
The salt transport processes in the default experiment are affected by the interplay between estuarine geometry, Coriolis deflection, and bathymetry. To investigate the influence of estuarine convergence on salt dynamics separately, an exponentially convergent channel with a horizontal bottom (no Coriolis) is studied in experiment I. In experiment II, the Coriolis force is added and hence the influence of the Coriolis force on salt dynamics can be found by comparing the results obtained in experiment II with those found in experiment I. In experiment III, the default bathymetry profile with a lateral channel–shoal structure [see Eq. (49)] is used, but now the Coriolis force is neglected. Hence, the influence of bathymetric variations on salt dynamics can be found by comparing the results from experiments III and I. The parameters used to define estuarine bathymetry and the Coriolis parameter are summarized in Table 2.
1) Influence of estuarine convergence (experiment I)
In a convergent channel without bathymetric variations or Coriolis force, the salt intrusion length is only 105 km, which is much smaller than that in the default experiment (120 km). The tidally averaged salinity is symmetric with respect to the midaxis, with the salinity near the banks slightly larger than that in the middle (see Fig. 4a). Since no circulation cells of salt transport are formed, the depth-integrated residual salt transport is zero at each location, and no streamlines are obtained.
The residual salt balance is identical to that obtained from the width-averaged model (Wei et al. 2016), where the seaward salt transport induced by river discharge is balanced by the landward transport by tidal advective diffusion and prescribed diffusion (Fig. 4b). The residual flow induced by gravitational circulation, tidal rectification, stress-free surface, and the tidal return flow, however, makes no contributions to the residual salt balance (integrated over the cross section).
The divergence of the residual salt transport induced by the diffusive processes, as shown in Fig. 5, indicates that diffusion results in a convergence (accumulation) of salt in the central region of the estuary. The salt transport paths induced by the diffusive processes are fully adapted to the estuarine geometry, with no horizontal circulations formed. It is important to note that even though no salt transport circulation cells are induced by including estuarine convergence, the salt dynamics can be significantly influenced by this factor. For instance, the estuarine convergence has a strong impact on the tidal water motion (Friedrichs and Aubrey 1994; Lanzoni and Seminara 1998) and river-induced flow due to changes of the cross-sectional area, thus affecting the salt transport contributions due to tidal advective diffusion (Wei et al. 2016) and river-induced flow (not shown).
2) Influence of Coriolis deflection (experiment II)
Including the Coriolis force in the convergent channel slightly increases the salt intrusion length to 108 km, with the salinity on the left larger than that on the right (see Fig. 6a). The salt transport pattern is significantly changed by the Coriolis deflection, with salt transported into the estuary from the left side of the channel and out from the right. Comparing the streamlines of the total salt transport in experiments I–II, it is found that the geometry-induced salt transport pattern as found in experiment I is replaced by a Coriolis-induced pattern in experiment II by including Coriolis force. It suggests that the influence of the Coriolis deflection on the total salt transport is much stronger than that of the estuarine convergence. Because of the Coriolis deflection, the seaward salt transport induced by river discharge is mainly balanced by the landward salt transport induced by diffusive processes, while the residual flow due to the stress-free surface condition and the tidal return flow makes a small but nonnegligible contribution (see Fig. 6b). The gravitational circulation, however, does not contribute to the residual salt balance.
Because of the Coriolis force, the landward salt transport by diffusive processes is deflected toward the left side of the estuary (see Fig. 7a), decreasing the cross-sectionally integrated longitudinal salt transport due to this contribution. It results in a significant divergence of salt on the left side and a convergence of salt on the right. The magnitude of local change in salinity induced by diffusive processes reaches up to 2 × 10−5 psu m s−1, much larger than that in experiment I where no Coriolis is included. The residual flow due to the stress-free surface condition and the tidal return flow tends to transport salt landward from the left side of the estuary and seaward from the right (see solid lines in Figs. 7b,c). It results in a remarkable accumulation of salt on the left side of the estuary and a divergence of salt on the right. The local change of salinity induced by the stress-free surface condition and the tidal return flow are up to 1 × 10−5 and 2 × 10−6 psu m s−1, respectively (see gray scales in Figs. 7b,c). The salt transport induced by tidal rectification and gravitational circulation are negligible (not shown).
3) Influence of estuarine bathymetry (experiment III)
By including a lateral channel–shoal structure, the salt intrusion length is significantly increased to 117 km (see Fig. 8a), compared to experiment I. Furthermore, the tidally averaged salinities become higher in the deeper channel than those on the shoals. The residual salt transport patterns are also completely changed: salt is transported landward over the shoals and seaward through the deeper channel. This clearly shows the influence of bathymetric variations is larger than that of the estuarine convergence. By including the bathymetric variations, the landward salt transport contributions due to tidal advective diffusion, prescribed diffusion, and gravitational circulation are increased (see Fig. 8b). On the other hand, the residual flow induced by the stress-free surface condition and tidal return flow tends to transport salt seaward, in contrast to their landward salt transport contributions in experiment II. The salt transport due to tidal rectification is very small compared to the other processes.
Because of the lateral channel–shoal structure, the diffusive processes transport salt from the shoals to the deeper channel and then transport salt landward (see Fig. 9a). This is accompanied with an accumulation of salinity in the deeper channel and a decrease of salinity on the shoals. The gravitational circulation presents a landward flow in the deeper channel and a seaward flow on the shoals, as found by Wong (1994), transporting salt into the estuary through the deeper channel and out over the shoals (see Fig. 9b). It results in a strong convergence of salt in the deeper channel and a small divergence of salt on the shoals. The residual flow induced by tidal rectification, stress-free surface condition, and tidal return flow is characterized by a net landward flow over the shoals and balanced by a return (seaward) flow in the deeper channel (see Figs. 9c–e). The pattern of these three residual components is consistent with the topography-induced exchange flow, as found by Li and O’Donnell (1997) (without considering density), which tends to transport salt landward over the shoals and seaward through the deeper channel. This results in a divergence of salt in the deeper channel and an accumulation of salt on the shoals. The amount of local salt loss induced by the tidal return flow or stress-free surface condition is up to 5 × 10−5 psu m s−1. The salt change induced by the tidal rectification is very small (less than 1 × 10−5 psu m s−1) compared to the default experiment. This confirms that effects of earth rotation are very important in generating the tidally rectified flow, as reported in Huijts et al. (2009).
The importance of various physical processes on salt intrusion length is evaluated by excluding each process individually and comparing corresponding changes in the salt intrusion length. It is found that the tidal advective diffusion is the predominant physical process for longitudinal salt intrusion, which is hardly influenced by residual flows (except for the river-induced flow). Since tidal advective diffusion is strongly influenced by many factors such as bottom friction, vertical mixing, external tidal forcing, estuarine bathymetry, and geometry (Wei et al. 2016), it is deduced that the salt intrusion length changes significantly with these factors. Though the residual circulations make no big contributions to salt intrusion lengths, they are crucial for the transverse structure of the tidally averaged salinity, which is essential for the transverse gravitational circulation and thence lateral sediment trapping in tidal estuaries (Huijts et al. 2006).
In section 4a, the salt transport due to the tidal advective diffusion is discussed: in section 4a(1), the salt transport induced by tidal advective diffusion in the idealized Delaware estuary is explained by comparing results of the experiments I–III with those of the default experiment; in section 4a(2), the resulting tidal advective diffusivity coefficients are discussed in relation to estuarine geometry, bathymetry, and Coriolis. The salt transport induced by the gravitational circulation, the strongest residual circulation process in the idealized Delaware estuary, is discussed in section 4b.
a. Tidal advective diffusion
1) Tidal advective salt flux
The salt flux induced by tidal advective diffusion (i.e., tidal advective salt flux) is a result of the temporal correlation between the tidal salinity S1 and the three-dimensional tidal currents u0, υ0, and w0. The longitudinal tidal advective salt flux can be expressed in terms of the amplitudes of u0 and S1 and their phase difference (Wei et al. 2016)
where Φs, Φu, |S1|, and |u0| are the phases and amplitudes of S1 and u0, respectively. Similarly, the transverse tidal advective salt flux (, ) is related to the amplitudes of the transverse tidal velocities, the tidal salinity, and their phase differences. The cross-sectional distributions of the longitudinal and transverse tidal advective salt fluxes are shown in Fig. 10 by gray scales and by arrows, respectively. All cross sections are located at x ~ 40 km.
In experiment I, the longitudinal tidal advective salt flux is similar to that from a width-averaged model by Wei et al. (2016). The longitudinal tidal advective salt fluxes are nearly laterally uniform, with landward fluxes near the surface and seaward fluxes near the bottom (see the contours in Fig. 10a). Moreover, the estuarine convergence results in a lateral tidal current that interacts with the tidal salinity and transports salt toward the midaxis near the free surface and toward the banks near the bottom. These geometry-induced transverse tidal advective salt fluxes are smaller than 0.005 psu m s−1 (see arrows in Fig. 10a). When including the Coriolis effects (experiment II), |u0| and |S1| become larger on the left than the right (see Figs. 11a,b), while their phase differences remain almost unchanged (see Fig. 11c). This results in larger longitudinal tidal advective salt fluxes on the left than on the right (see Fig. 10b). In the transverse direction, the lateral tidal current induced by the Coriolis effects interacts with the tidal salinity and yields salt fluxes toward the left near the surface and toward the right near the bottom. This Coriolis-induced circulation, together with the geometry-induced circulation, results in a slightly enhanced counterclockwise salt transport on the right side of the estuary and a slightly weakened clockwise transport on the left (see arrows in Fig. 10b).
In experiment III, tides propagate faster in the deeper channel than on the shoal, resulting in a larger |u0| and |S1| in the deeper channel than on the shoals (see Figs. 11d,e). Meanwhile, u0 in the deeper channel lags that on the shoals, and S1 on the shoals lags that in the deeper channel. Thus, in the deeper channel, the phase difference between u0 and S1 is smaller than 90° near the surface and slightly larger than 90° near the bottom; on the shoals, it is larger than 90° near the bottom, and smaller than 90° near the surface (see Fig. 11f). This results in enhanced landward tidal advective salt fluxes near the surface (relatively small seaward fluxes near the bottom) in the deeper channel. On the shoals, since |u0| and |S1| are much smaller than those in the deeper channel, the induced salt fluxes are much smaller (see gray scales in Fig. 10c). Integrated over the depth, the landward tidal advective salt transport is much stronger in the deeper channel. Moreover, the interaction between the tidal salinity and the lateral tidal current induced by the lateral depth variations results in axial-convergent tidal advective salt fluxes (up to 0.01 psu m s−1) toward the midaxis at all depths (see arrows in Fig. 10c).
Comparing the tidal advective salt fluxes in experiments I–III with the default experiment suggests that the lateral bathymetric variations dominate the longitudinal tidal advective salt fluxes in the idealized Delaware estuary. The lateral tidal advective salt fluxes, however, are influenced by estuarine convergence, Coriolis, and bathymetry. Because of the overall effects, the tidal advective salt fluxes toward the left on the top right of the channel as well as the fluxes toward the right on the bottom left are augmented (see Fig. 10d).
2) Tidal advective diffusivity
The strength of the (depth integrated) salt transport due to tidal advective diffusion is measured by the tidal advective diffusion coefficients, which result from the temporal correlation between the longitudinal/lateral tidal velocities and tidal salinities [see Eqs. (29)–(32)]. Since more processes are resolved in a three-dimensional model compared to a width-averaged approach, it is to be expected that the prescribed diffusion coefficient (including unresolved processes) has to be smaller compared to the values that have to be used in width-averaged models to qualitatively reproduce the observed salt intrusion in the Delaware estuary. Figure 12 shows the tidal advective diffusion coefficients for all experiments. Here, only , , and are shown because ( m2 s−1) is much smaller than the other diffusion coefficients.
The tidal advective diffusivity is positive in all experiments, confirming that the temporal correlation between u0 and S1 always results in a landward salt transport integrated over the depth (Wei et al. 2016). In experiment I, is almost the same as obtained from the width-averaged model in Wei et al. (2016), which is laterally uniform (see Fig. 12a). This implies that the estuarine convergence has little influence on . However, different from the width-averaged model, and become nonnegligible (see Figs. 12b,c).
Coriolis effects significantly change the spatial pattern of compared to experiment I, with larger magnitudes on the left than on the right (see Fig. 12d), resulting in a larger along-channel tidal advective salt transport on the left. The terms and show slightly larger magnitudes near the entrance compared to experiment I (see Figs. 12e,f) and are almost unchanged in the landward region. In experiment III, is strongly increased in the deeper channel and decreased on the shoals (see Fig. 12g). This lateral difference of contributes to a stronger along-channel tidal advective salt transport in the deeper channel than on the shoals. The lateral bathymetric variations also significantly increase the magnitudes of and (see Figs. 12h–i). Comparing experiments I–III with the default experiment suggests that lateral bathymetric variations dominate the magnitudes of the tidal advective diffusion coefficients. The Coriolis force hardly changes (see Fig. 12j); however, it slightly increases the magnitudes of and on the left and decreases them on the right (see Figs. 12k,l). The magnitudes of the tidal advective diffusion coefficients are much larger than those of the prescribed horizontal diffusivity Kh. It means that by resolving the lateral processes in the three-dimensional model, the significance of the tidal advective diffusion is generally magnified, resulting in fewer processes that need to be parameterized by a prescribed horizontal diffusion coefficient than in the width-averaged model of Wei et al. (2016).
b. Gravitational circulation
As found by Wei et al. (2016), the gravitational circulation does not influence the residual salt transport in convergent estuaries with a horizontal bed (with negligible lateral processes), since the near-bottom landward salt transport is canceled by the near-surface seaward transport (see gray scales in Figs. 13a,b). It suggests in estuaries with no depth variations, the gravitational circulation and salinity are hardly coupled, and the gravitational circulation can be understood straightforward using the obtained salinity.
The estuarine convergence induces a lateral gradient of the tidally averaged salinity with S0 near the banks larger than that in the center. To balance this gradient, a barotropic pressure gradient is induced with a larger in the center than that near the banks (see upper panel of Fig. 13a). These gradients result in a small but nonnegligible transverse gravitational circulation, which transports salt (smaller than 0.02 psu m s−1) toward the banks near the surface and toward the midaxis near the bottom (see arrows in Fig. 13a). Under the influence of Coriolis effects, however, S0 becomes larger on the left than on the right, and becomes larger on the right than on the left in order to balance the salinity gradient (see upper panel of Fig. 13b). This results in a transverse gravitational circulation that transports salt toward the left near the surface and toward the right close to the bottom (see arrows in Fig. 13b). The salt fluxes induced by the Coriolis-induced transverse gravitational circulation (up to 0.1 psu m s−1) are much stronger than those induced by the estuarine convergence.
Because of the lateral bathymetric variations, the gravitational circulation transports salt landward through the deeper channel and seaward over the shoals (see gray scales in Figs. 13c,d). The longitudinal gravitational circulation is qualitatively consistent with Lerczak and Geyer (2004), with the strongest landward fluxes near the bottom of the deeper channel and the strongest seaward fluxes near the surface of the shoals. Nevertheless, there is a discrepancy: the gravitational circulation is generally landward near the bottom and seaward near the surface in Lerczak and Geyer (2004); in the present work, however, the landward flow is located at all depths of the deeper channel with the seaward flow on the shoals. This discrepancy is probably related to different vertical mixing coefficients used in two models; Lerczak and Geyer (2004) used a relatively small Aυ (0.002 m2 s−1) for their well-mixed runs, while in the present model Aυ is up to 0.005 m2 s−1 in the deeper channel. It is important to note that in the presence of the channel–shoal structure, the gravitational circulation and the salinity distribution become strongly coupled, and hence they cannot be interpreted straightforwardly, as in experiments I–II.
In experiment III, the diffusive processes and the residual flows act together, resulting in a larger (smaller) S0 () in the deeper channel than on the shoals (see upper panel of Fig. 13c). This induces a two-cell transverse circulation with near-bottom salt fluxes toward the banks and near-surface return fluxes (up to 0.06 psu m s−1) toward the midaxis (see arrows in Fig. 13c). This pattern is consistent with the bathymetry-controlled secondary circulation reported by Nunes and Simpson (1985). Including the Coriolis force hardly influences the longitudinal gravitational circulation; however, it significantly enhances the transverse gravitational circulation (see arrows in Fig. 13d). Because of the combined effects of the estuarine convergence, lateral bathymetric variations, and Coriolis force, S0 () becomes larger (smaller) on the left than on the right, inducing near-surface salt fluxes (up to 0.25 psu m s−1) toward the left and near-bottom fluxes toward the right (see Fig. 13d). The transverse salt fluxes induced by gravitational circulation have almost the same magnitude as the longitudinal fluxes, implying the transverse gravitational circulation is as significant as the longitudinal circulation.
5. Response to river discharge
To systematically explore the response of salt intrusion to river discharge, the model is applied to experiment III using three different river discharge values: 72, 288, and 864 m3 s−1, respectively. It is found that the simulated, tidally averaged, salt intrusion length varies significantly from 120 km (for Q = 72 m3 s−1) to 80 km (Q = 288 m3 s−1) and 60 km (Q = 864 m3 s−1). This strong dependency of salt intrusion length on river discharge is not consistent with observations of Garvine et al. (1992), who found the salt intrusion length only weakly responds to the increase of river discharge in the Delaware estuary. Earlier studies have shown that increasing river discharge can increase stratification and lead to weaker vertical mixing, consequently enhancing the residual circulation (Monismith et al. 2002; Ralston et al. 2008). This enhanced residual circulation is found to play an important role in transporting salt landward and weakening the dependency of salt intrusion on the variations of river discharge (Ralston et al. 2008; Lerczak et al. 2009). Therefore, to understand the overestimated sensitivity of the simulated salt intrusion to river discharge, the stratification (measured by the top-to-bottom salinity difference) at each location is calculated for different river discharges (not shown). It is found that the stratification in the central region of salt intrusion increases with increasing river discharge, consistent with model studies (Ralston et al. 2008; Lerczak et al. 2009) and observations (McSweeney et al. 2016). The top-to-bottom salinity differences (less than 1.5 psu), however, are much smaller than the observed values during large river discharges (McSweeney et al. 2016). This implies that the relatively strong stratification in case of large river discharges is not resolved in this model, which is probably the reason for the overestimated impact of river discharge on salt intrusion. Therefore, to reproduce the observed salt intrusion for large river discharges using the model presented in this paper, a larger horizontal eddy diffusivity coefficient needs to be prescribed. This is consistent with Gay and O’Donnell (2009), who captured the weak response of the salt intrusion length to river discharge in the Delaware estuary by allowing the along-channel dispersion coefficient to exponentially increase with river discharge.
A semianalytical, three-dimensional model is developed to investigate the three-dimensional salt dynamics in idealized well-mixed estuaries. The adopted perturbation method allows for a systematic decomposition of the total salt transport into contributions due to various physical processes. Each salt transport process and its relative importance can be investigated in isolation.
First, the salt dynamics for a schematized estuary is investigated. The schematized estuary is exponentially converging, and the bathymetry shows a channel–shoal structure. In this experiment, the tidal advective diffusion, induced by the temporal correlation between horizontal tidal velocities and salinity (including the lateral processes), is the dominant landward salt transport process. The cross-sectionally integrated salt transport induced by residual circulations due to gravitational circulation, tidal rectification, stress-free surface condition, and tidal return flow is small compared to that induced by the tidal advective diffusion. The salt transport due to these residual circulations, however, is locally significant. The gravitational circulation transports salt landward through the deeper channel and seaward over the shoals. The residual circulations due to other mechanisms are smaller and transport salt in the opposite direction compared to the transport pattern of gravitational circulation.
The residual salt balance is influenced by the estuarine convergence, bathymetry, and Coriolis. Three experiments are conducted in order to investigate these three factors on the estuarine salt transport in separation: in experiment I, an exponentially convergent estuary is used with a horizontal bed, no Coriolis; in experiment II, the estuarine geometry and bathymetry are the same as in experiment I, but the Coriolis force is included; and in experiment III, a channel–shoal structure is included while the Coriolis effects are excluded. The tidal advective diffusion is found to be the most dominant landward salt transport process for all experiments. The salt intrusion length slightly increases, compared to the results of experiment I, when Coriolis is included (experiment II) and increases significantly by including the channel–shoal structure (experiment III). The contribution of tidal advective diffusion to estuarine salt transport can be measured by the tidal advective diffusion coefficients, which can be explicitly calculated and are strongly dependent on the Coriolis effects and lateral bathymetric variations.
The residual circulations hardly influence the salt intrusion length, but are crucial for the local salt transport. In experiment I, the depth-integrated residual circulations vanish, with the river-induced seaward salt transport balanced by diffusive processes only. No significant lateral salt transport patterns are induced by the width convergence of the estuary. In experiment II (including Coriolis deflection), the residual circulations due to tidal rectification and the stress-free surface condition transport salt into the estuary from the left side of the estuary (looking seaward) and out on the right side, resulting in a significant contribution to the landward salt transport. In experiment III, salt transport due to all processes becomes stronger due to the channel–shoal structure. The gravitational circulation, as the most significant residual circulation process, is characterized by a symmetric two-cell circulation pattern. Residual circulations due to other processes are relatively weak. The residual salt transport due to residual circulations is essential for the lateral salinity structure, which can drive a strong transverse gravitational circulation within the salt intrusion region.
The authors appreciate the suggestion from Dr. George P. Schramkowski (Flanders Hydraulics Research) about visualizing salt transport patterns using streamlines. This project is mainly funded by the China Scholarship Council (CSC, File 201206710049). This work was sponsored by NWO Exacte Wetenschappen (Physical Sciences) for the use of supercomputer facilities, with financial support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek [Netherlands Organization for Scientific Research (NWO)]. The work is part of the research programme NWO–ALW project 8843.10.005, which is financed by the NWO and the Chinese Organisation for Scientific Research (NSFC).
A perturbation method is used to analytically solve the water motion and salinity. First of all, variables are scaled with their typical scales (see Table A1), denoted by a tilde (~). This model does not resolve the boundary layer, and processes in typical estuarine length scales are focused. The variations of water motion and salinity in the longitudinal and lateral directions are assumed to scale with the typical estuarine length scale Lb. For estuaries with an exponentially decreasing width, Lb is equivalent to the estuarine convergence length [see Eq. (46)]. The variations in the vertical are assumed to scale with the estuarine water depth.
The scale for velocities in the horizontal directions (denoted by U) is obtained by approximately balancing the terms in the depth-averaged continuity equation:
which yields . The scale for the vertical velocity (denoted by W) is obtained by assuming the terms in the continuity equation [Eq. (A2)] to be at the same order; hence, we found W = UH0/Lb. The river discharge per unit width at the landward boundary is scaled by the width-averaged river discharge: , with BL as the estuarine width at the landward boundary.
Substituting these scales into the shallow-water equations and the salinity equation and replacing the density in terms of salinity in the momentum equations, we derived the dimensionless governing equations for the water motion and salinity, which read
The corresponding dimensionless boundary conditions at the free surface and the bottom are given by
At the seaward boundary, we have
At the closed boundaries and the landward boundary, the dimensionless boundary conditions are
As a next step, the magnitudes of the scaling parameters are compared with a small parameter ; the ratio of the width-averaged M2 tidal amplitude to the mean water depth at the mouth. This results in an ordering of scaling parameters of each term in the dimensionless equations, as summarized in Table A2.
Then, we substitute the magnitudes of these scaling parameters into Eqs. (A2)–(A9), which yields the dimensionless governing equations and boundary conditions for water motion and salinity in terms of ε. Substituting the expansions of each dimensionless variable in terms of ε and collecting the terms at the same order of ε, we derived the dimensional governing equations and boundary conditions at different order of ε, which will be used to calculate the water motion and salinity semianalytically.
Leading-Order Water Motion
The leading-order water motion is forced by the semidiurnal tide; thus, it consists of only the M2 tidal constituent.
As shown by Kumar et al. (2016), u0 and υ0 can be written in terms with the gradients of η0:
In the case of no Coriolis force (f = 0), we have C2 = 0, yielding
First-Order Salinity Equations
Since the leading-order salinity is independent of depth, the salinity equation at the first order reduces to
As discussed before, the boundary condition of S1 at the top and the bottom is required:
where and are known functions given by
In the case of f = 0, we find , and becomes the same as Sz, which is derived by Wei et al. (2016).
First-Order Residual Flow
a. Governing equations
The first-order residual flow follows from the first-order tidally averaged shallow-water equations
with boundary conditions
Here, the underbrace denotes various mechanisms that force the residual flow. Equations (D1)–(D2) show that the residual flow is forced by advective contributions of the leading-order M2 tide (tidal rectification, denoted by AD), density-driven gravitational circulation (GC), the stress-free surface condition (NS), river discharge (RD), and a return flow (TRF). Since the first-order water motion Eqs. (D1) and (D2) are linear, the residual flow components due to these forcing mechanisms can be solved separately. Hence, the solution of the residual water motion can be written as
with the solution vector . All residual contributions can be calculated explicitly without information of the salinity field, except gravitational circulation, which is dependent on salinity itself.
b. Gravitational circulation
To solve the gravitational circulation, the following system of equations has to be solved:
with boundary conditions
Introducing two rotating flow variables and (Kumar et al. 2017),
Substitute Eq. (D7) into Eq. (D4a), integrate it over depth, and apply corresponding boundary conditions. And then, we obtained an elliptic partial differential equation for the GC-induced residual sea surface elevation:
c. The seaward boundary condition
When solving this coupled system of equations, special care has to be taken of the seaward boundary, since numerical instabilities can occur when using an arbitrary prescribed water level at this boundary (Chen and Sanford 2009). This is because an inconsistent prescribed sea surface elevation at the mouth can result in a rapid adjustment to a consistent solution in the interior domain. This rapid adjustment will result in large unphysical gradients in the water motion close to the boundary and thus cause numerically inaccurate solutions.
To avoid that, the computational domain is extended by 30 km in the seaward direction. At the open boundary of the computational domain, a constant M2 tidal water level is prescribed, and the tidal water motion is calculated. This results in a physically consistent water motion within the domain of interest. Then, the resulting M2 tidal water elevation at the seaward boundary of the physical domain is rescaled, such that the width-averaged tidal amplitude at the original open boundary (now located at x = 30 km) is the same as the tidal amplitude that is prescribed at the open boundary of the physical domain. This will result in a consistent M2 water motion in the domain of interest. The extended domain is then used to calculate the residual flow using the calculated tidal water motion, except for the residual flow generated by the baroclinic pressure. To satisfy the physical boundary, a constant residual water level is used at the open boundary of the computational domain such that the width-averaged residual water level is zero at the open boundary of the physical domain. After that the calculated water motion in the physical domain is used to obtain the salinity field and gravitational circulation.