A Baroclinic Fluid Model and Its Application in Investigating the Salinity Transport Process Within the Sediment–Water Interface in an Idealized Estuary

Abstract

Understanding the salinity transport process around the sediment–water interface is important for water resources management in the upper reach of an estuary. In this study, we developed a baroclinic fluid dynamic model for investigating the flow and salt transport characteristics within the sediment–water interface under tidal forcing. The validation showed robust model performance on the salinity transport within the sediment–water interface. The results revealed that the turbulent kinetic energy, dissipation rate, and kinetic energy production rate exhibited periodic variations within the seabed boundary layer. The thickness of the viscous sublayer and the mean flow showed an inverse relationship. Water and salinity exchange within the sediment–water interface occurred predominantly via turbulent diffusion, with extreme turbulent kinetic energy production rates appearing during the tidal reversal, flood, and ebb stages. The sediment acted as a source of salinity release during ebb tides and a sink for salinity absorption during flood tides. As the sediment depth increased, fluctuations in salinity were weakened. These results clearly illustrated that the sediment layer is important in modulating the salinity transport in the upper reach of an estuary. However, such an important process was usually excluded by previous studies. The model developed in this study can be used as a sediment–water interface module that, coupled with other hydrodynamic models, can evaluate the contributions of the sediment layer to the salinity exchange in coastal water.

1. Introduction

Being the interface of river deltas and the open ocean, estuaries are facing significant challenges due to water body salinization and elevated salt content in the soil caused by the combined effect of anthropogenic activity and climate change. In order to promote sustainable water resources development, it is essential to accurately assess the migration of salinity within both the water column and the sediment layer. So far, substantial academic focus has been directed toward understanding the hydrodynamic characteristics of salinity transport processes in the water column, whereas research focusing on the sediment–water interface has been quite limited. Understanding the flow characteristics within the sediment–water interface is crucial as they significantly influence the substance exchange between the sediment layer and the overlying water. Assessing the processes modulating salinity transport within the sediment–water interface could facilitate better water resources management and environment protection.

The sediment–water interface layer encompasses the seabed boundary layer, sediment–water interface, and surface layer of the sediment. The seabed boundary layer consists of the logarithmic layer and the viscous sublayer, while the sediment surface layer extends as far as 0.5 m below the sediment–water interface. Intricate transport mechanisms involve convective diffusion processes that extend beyond molecular diffusion. Substance exchange between sediments and the overlying water, including the exchange of nitrogen, oxygen, nutrients, and metal ions, has been widely studied [1,2,3]. Numerical simulation methods were employed to investigate salinity release at the sediment–water interface, particularly in low-permeability zones. Li et al. [4] conducted field sampling and numerical experiments to calculate hydrogeological parameters and monitor salinity variations in sediment and water. Fu et al. [5] conducted a study to analyze the salt accumulation and release process in shallow underground brine during tidal cycles. After analyzing measured sediment permeability and pore water parameters, they found that hydrodynamic characteristics near the sediment–water interface determine the movement and distribution of bottom sediments, consequently affecting the distribution of nutrients and organic matter within the sediment. Salinity exchange regulates the distribution of salinity in the water and the diffusion of dissolved substances and sediment deposition. Zaghmouri et al. [6] studied the salinity variations within sediment layers and their impact on the microenvironment. The results offered valuable insights into the relationships between salinity changes and their implications for biodiversity, ecological processes, and energy fluxes in estuarine areas. These studies highlight the critical importance of coastal salinity and water treatment solutions. As coastal regions face increasing salinity challenges due to climate change and human activity, developing effective materials and techniques for desalination becomes essential.

The sediment–water interface features a viscous sublayer, which weakens convection and diffusion between the overlying water and sediment. Within this viscous sublayer, the molecular viscosity primarily governs the Reynolds stress associated with the vertical momentum transport. Although the turbulent viscosity is lower than the kinematic viscosity of the water flow, turbulent fluctuations still influence the mass flux. By contrast, within the diffusion sublayer, the resistance to water flow surpasses that of the pores. As a result, the mass flux through the sediment–water interface is predominantly controlled by diffusion, with molecular diffusion dominating turbulent diffusion. Niu et al. [7] investigated sediment dynamics in Western Lake Erie using a sediment model. The results revealed that high-turbidity events during spring and fall were primarily driven by resuspension events, while those in summer were influenced by pulsing river loadings near the mouths. River loading played a crucial role during floods, enlarging high-turbidity areas. Additionally, sediment porosity plays an important role in shaping water flow velocity, diffusion rate, and turbulence intensity. Sediment layers with high porosity display enhanced permeability and larger pore spaces, resulting in a more pronounced salinity gradient in the sediment–water interface layer, faster water flow velocity, higher diffusion rates, and greater turbulence intensity [8,9].

The seabed boundary layer experiences turbulent motion driven by large-scale water flow in the overlying layer. Tidal forcing directly influences the hydrodynamic characteristics of the sediment–water interface in an estuary. Hydrodynamic processes exhibit periodicity under tidal forcing, leading to changes in water flow structure and further influencing the convection and diffusion processes between the overlying water and sediment layer [10]. The flow velocity significantly influences the diffusion flux within the viscous sublayer. As the flow velocity increases, the thickness of the viscous sublayer decreases, increasing the diffusion flux [11]. Conversely, increasing the water depth reduces the bottom shear stress, decreasing the diffusion flux. Periodic water flow under non-steady conditions is closely linked to the substance transport in the interface layer and the dynamic response of the boundary layer. The thickness of the boundary layer exhibits periodic variations with changes in the flow velocity, resulting in vigorous turbulent mixing between the overlying water and sediment pore water at the sediment surface [12].

It is still challenging to model the hydrodynamic processes around the sediment–water interface layer, particularly with bidirectional flow. Previous studies using unidirectional flow in the seabed layer failed to capture the complete flow dynamics [13,14,15]. Li et al. [16] introduced a numerical method to simulate the seabed flow characteristics under bidirectional tidal flow. As tides undergo continuous fluctuations, they did not obtain a comprehensive understanding of the processes of salt absorption and release in the sediment–water interface layer. In the upper reach of the estuary, the sediment layer is very important in modulating the salinity transport, especially for the changing stratification during flood and ebb tides. However, such an important process was usually excluded in other studies. The purpose of this study is to reveal how the interface layer modulates the salinity absorption and release in the estuarine environment. By incorporating the porosity in a baroclinic fluid model, the non-steady Navier–Stokes equations and the Reynolds stress transport equation are employed to comprehensively examine the hydrodynamic changes in the sediment–water interface layer. The transport and exchange processes of salinity between the water body and the sediment layer will be thoroughly examined. For simplicity, an idealized estuary is configured for model application. Our methods will help improve the prediction and evaluation of salinity transport within the sediment and overlying water and offer insights for water resources management.

The remainder of the paper is organized as follows: Section 2 introduces the numerical model and solving methods. Section 3 presents the validation of the numerical model. Section 4 presents the results and discussion. Conclusions are given in Section 5.

2. Numerical Model Development

2.1. Model Development

There was a lack of existing models addressing the salt absorption and release in the sediment–water interface layer. Therefore, we developed a model to incorporate such processes. The baroclinic fluid dynamic equation was used to develop a fully stratified hydrodynamic model which aimed to investigate the salinity transport process at the sediment–water interface. Unlike the barotropic model in the previous work, the baroclinic model in this study incorporated the salinity transport both in the stratified water column and the sediment layer. By leveraging this modified mathematical framework, we successfully developed a model capable of simulating the flow and salt transport processes in the sediment–water interface layer which have been observed in real partially mixed estuaries under tidal impacts.

2.1.1. The Momentum and Continuity Equations

The continuity equation is

$${\displaystyle \frac{\partial \left(\varphi \rho \right)}{\partial t}}+\left[{\displaystyle \frac{\partial }{\partial x}}\left(\varphi \rho u\right)+{\displaystyle \frac{\partial }{\partial y}}\left(\varphi \rho v\right)+{\displaystyle \frac{\partial }{\partial z}}\left(\varphi \rho w\right)\right]=0$$

(1)

The momentum equations (Navier–Stokes equations) are

$${\displaystyle \frac{\partial \left(\varphi \rho u\right)}{\partial t}}+{\displaystyle \frac{\partial \left(\varphi \rho u^2\right)}{\partial x}}+{\displaystyle \frac{\partial \left(\varphi \rho uv\right)}{\partial y}}+{\displaystyle \frac{\partial \left(\varphi \rho uw\right)}{\partial z}}=-{\displaystyle \frac{\partial \left(\varphi \rho \right)}{\partial x}}+{\displaystyle \frac{{\partial }^2\left(\varphi \mu u\right)}{\partial x^2}}-{\displaystyle \frac{\partial \varphi \rho \overline{u^{\prime }{v}^{\prime }}}{\partial y}}-{\displaystyle \frac{\partial \varphi \rho \overline{u^{\prime }{w}^{\prime }}}{\partial z}}+S_x$$

(2)

$${\displaystyle \frac{\partial \left(\varphi \rho v\right)}{\partial t}}+{\displaystyle \frac{\partial \left(\varphi \rho uv\right)}{\partial x}}+{\displaystyle \frac{\partial \left(\varphi \rho v^2\right)}{\partial y}}+{\displaystyle \frac{\partial \left(\varphi \rho vw\right)}{\partial z}}=-{\displaystyle \frac{\partial \left(\varphi \rho \right)}{\partial y}}+{\displaystyle \frac{{\partial }^2\left(\varphi \mu v\right)}{\partial y^2}}-{\displaystyle \frac{\partial \varphi \rho \overline{v^{\prime }{u}^{\prime }}}{\partial x}}-{\displaystyle \frac{\partial \varphi \rho \overline{v^{\prime }{w}^{\prime }}}{\partial z}}+S_y$$

(3)

$${\displaystyle \frac{\partial \left(\varphi \rho w\right)}{\partial t}}+{\displaystyle \frac{\partial \left(\varphi \rho uw\right)}{\partial x}}+{\displaystyle \frac{\partial \left(\varphi \rho vw\right)}{\partial y}}+{\displaystyle \frac{\partial \left(\varphi \rho w^2\right)}{\partial z}}=-{\displaystyle \frac{\partial \left(\varphi \rho \right)}{\partial z}}+{\displaystyle \frac{{\partial }^2\left(\varphi \mu w\right)}{\partial z^2}}-{\displaystyle \frac{\partial \varphi \rho \overline{w^{\prime }{u}^{\prime }}}{\partial x}}-{\displaystyle \frac{\partial \varphi \rho \overline{w^{\prime }{v}^{\prime }}}{\partial y}}+S_z$$

(4)

where $u$, $v$, and $w$ are the velocity component in the $x$, $y$, and $z$ directions, respectively; $\varphi$ is the porosity; $\rho$ is a function of density concerning temperature and salinity; $-\rho \overline{u\prime v\prime }$, $-\rho \overline{v\prime w\prime }$, and $-\rho \overline{w\prime u\prime }$ are the Reynolds stress; $t$ is the time; $\mu$ is the dynamic viscosity coefficient; and $S_x$, $S_y$, and $S_z$ are the resistance source term of the momentum equation in the $x$, $y$, and $z$ directions.

In order to account for the momentum losses incurred by the movement of water within the sediment layer, a resistance source term was added to the momentum equation. This was performed to quantify the momentum loss of the water flow. The momentum equations were modified by including a resistance source term for the flow resistance and movement loss within the sediment layer, reflecting the pressure gradient generated within each porous medium cell. This results in a pressure drop proportional to the velocity (or velocity squared) and includes viscous and inertial loss terms. The source term in the momentum equation [17] is expressed as follows:

$$S_X=-\left({\varphi }^2{D}_X\mu \left(u+v+w\right)+{\varphi }^3{C}_X{\displaystyle \frac{1}{2}}\rho \left(\left|u\right|u+\left|v\right|v+\left|w\right|w\right)\right)$$

(5)

where $X$ represents the $x$, $y$, and $z$ directions, respectively; $\left|u\right|, \left|v\right|, \left|w\right|$ represent the magnitude of the velocity; and $D_X$ and $C_X$ represent the coefficients of viscous resistance and inertial resistance, respectively.

In the calculation domain, the actual area simulated by the calculations is primarily composed of silty clay, which is why the porosity is 0.5. The horizontal permeability is 1.4 × 10⁻¹⁶ m², and the vertical permeability is 1.1 × 10⁻¹⁶ m² [18]. These parameters were used to calculate the values of the viscous loss and inertial loss coefficients using empirical formulas [19], as shown in

Table 1. Values of the viscous loss and inertial loss coefficients.

Direction	${\mathit{D}}_{\mathit{X}}$ (1/m2)	${\mathit{C}}_{\mathit{X}}$ (1/m2)
$x$	7 × 1015	280,000
$y$	7 × 1015	280,000
$z$	9 × 1015	280,000

below.

This study used the Reynolds stress model (RSM) to represent the anisotropy of turbulent viscosity. The Reynolds stress transport and turbulent dissipation rate transport equations are expressed as follows:

$${\displaystyle \frac{\partial }{\partial t}}\left(\varphi \rho \overline{{u^{\prime }}_i{u^{\prime }}_j}\right)+{\displaystyle \frac{\partial }{\partial x_r}}\left(\varphi \rho u_r\overline{{u^{\prime }}_i{u^{\prime }}_j}\right)=D_{T,ij}+D_{L,ij}+P_{ij}+G_{ij}+{\varphi }_{ij}+{\epsilon }_{ij}+F_{ij}+S_{RS}$$

(6)

In the context of turbulence modeling, $i$ and $j$ typically represent the $x$, $y$, and $z$ coordinates of velocity components. $D_{T,ij}$ represents the turbulent diffusion term, $D_{L,ij}$ is the molecular diffusion term, $P_{ij}$ corresponds to the stress generation term, $G_{ij}$ represents the buoyancy production term, ${\varphi }_{ij}$ is the stress–strain term, ${\epsilon }_{ij}$ denotes the dissipation term, $F_{ij}$ denotes the stress generation owing to system rotation, $S_{RS}$ represents the Reynolds stress source term, and ${\mu }_t$ represents turbulent viscosity [20,21].

2.1.2. Salt Transport Equation

Considering the density variation due to salinity, the salt transport equation and equation of state were modeled using the baroclinic model. The sediment layer is treated as a rigid porous medium, and the salinity exchange between the sediment layer and water column is represented by turbulent diffusion and convective diffusion processes. The continuity between the overlying water and the pore water within the porous medium is maintained, enabling the investigation of the salinity exchange mechanism at the sediment–water interface.

The salt transport equation is as follows:

$$\begin{array}{l}{\displaystyle \frac{\partial \left(\varphi \rho S\right)}{\partial t}}+\left[{\displaystyle \frac{\partial }{\partial x}}\left(\varphi \rho uS\right)+{\displaystyle \frac{\partial }{\partial y}}\left(\varphi \rho vS\right)+{\displaystyle \frac{\partial }{\partial z}}\left(\varphi \rho wS\right)\right]\\ =\varphi \left(\rho \lambda +{\displaystyle \frac{{\mu }_t}{S{c}_t}}\right)\left[{\displaystyle \frac{{\partial }^2}{\partial x^2}}\left(\varphi \rho uS\right)+{\displaystyle \frac{{\partial }^2}{\partial y^2}}\left(\varphi \rho vS\right)+{\displaystyle \frac{{\partial }^2}{\partial z^2}}\left(\varphi \rho wS\right)\right]+S_S\end{array}$$

(7)

The state equation [22] is as follows:

$$\begin{array}{cc}\hfill \rho \left(T,S\right)& =a_1+a_{2}T+a_3{T}^2+a_4{T}^3+a_5{T}^4+a_6{T}^5\hfill \\ & +S\left(a_7+a_{8}T+a_9{T}^2+a_{10}{T}^3+a_{11}{T}^4\right)\hfill \\ & +S^{2/3}\left(a_{12}+a_{13}T+a_{14}{T}^2\right)+a_{15}{S}^2\hfill \end{array}$$

(8)

where $\varphi$ is the porosity; $\rho$ is a function of density concerning temperature and salinity; $t$ is the time; $S_S$ is the source term of the salt transport equation; $S$ is the salinity; $\lambda$ is the molecular diffusion coefficient; ${\mu }_t$ is the turbulent viscosity; $S_{C_t}$ is the turbulent Schmidt number, value 0.7; $T$ is the seawater temperature, where the temperature change in the calculation process can be neglected and the fixed value taken; and $a_1~a_{15}$ is the empirical coefficient.

The salinity diffusion flux at the sediment–water interface is determined by the salinity gradient across the interface and can be calculated based on Fick’s first law [23]:

$$F=\varphi D_s{\left({\displaystyle \frac{\partial C}{\partial Z}}\right)}_{Z=0}$$

(9)

where $F$ represents the molecular diffusion flux; $D_s$ is the diffusion coefficient; and $\varphi$ is the porosity of the sediment layer when $\varphi \le 0.7$, $D_s=\varphi D_0$, and $\varphi >0.7$ are $D_s={\varphi }^2{D}_0$. The calculation area was based on the actual porosity of 0.5, and the sedimentary layer was defined as a rigid region, so the porosity remained constant. $D_0$ is the diffusion coefficient of the molecule in seawater, which is 19.8 × 10⁻⁶ cm²/s, and ${\left({\displaystyle \frac{\partial C}{\partial Z}}\right)}_{Z=0}$ represents the concentration gradient at the sediment–water interface.

2.2. Solution Methods

Two commonly used macroscopic approaches were employed in the numerical computation (Figure 1). The first approach involves a two-domain computational method that couples the transport equations of the upper free water flow with the time-averaged transport equations in a porous medium, incorporating appropriate boundary conditions [24]. However, this method neglects the water releasing from the sediment layer to the overlying water column, resulting in one-way coupling. The second approach, known as the full-water stratification method, dispenses with the Brinkman correction term in the general form of the Brinkman equation. It represents the flow resistance in a porous medium using the Darcy term, which accounts for viscous and inertial resistance. This method introduces a Darcy–Forchheimer modified resistance source term to the momentum equation [25]. The second method allows for a more comprehensive two-way coupling between the upper water body and the sediment layer. This study adopted the second method to simulate the variations in hydrodynamics and salt transport in the sediment–water interface layer. This approach enabled a more accurate representation of the interactions between water and sediment, contributing to a deeper understanding of the dynamic processes occurring in the interface layer.

The finite volume method was used to present a comprehensive numerical simulation of two-way periodic water flow. The dynamic behavior of water flow was accurately captured by solving the Reynolds-averaged Navier–Stokes and Reynolds stress transport equations. The previous work successfully validated the numerical calculation method by comparison with measured data, confirming its reliability for computing two-way unsteady flow processes in an upper water body. In this study, we adopted the concept of full-water stratification by considering the sediment layer as a porous medium. This innovative approach allows for a more sophisticated representation of the complex interactions between water and sediment in upper water bodies. As a result, the model can effectively investigate intricate processes occurring in the sediment–water interface layer.

The finite volume method was employed to discretize the governing equations, enabling accurate numerical simulations of the hydrodynamic phenomena. The values at the control volume grid boundaries and their derivatives were obtained by interpolating them at the nodes, ensuring a precise representation of the flow characteristics. The Pressure-Implicit with Splitting of Operators (PISO) method was used to calculate the pressure and velocity terms. The diffusion term was discretized using a central differencing scheme to guarantee the stability and accuracy of the diffusion processes. To ensure the model converged quickly and maintained its stability, the Courant number was set to 0.5. The convection term was discretized using a second-order upwind scheme that effectively captured the dominant flow direction for enhanced accuracy. The remaining terms were discretized using a first-order backward differencing scheme to achieve computational efficiency. This approach balances accuracy with computation time, making it well suited to simulating complex hydrodynamic scenarios.

2.3. Numerical Model Application

2.3.1. Grid Partitioning

The model was used to simulate the salinity transport at the sediment–water interface of an idealized estuary (Figure 2). The open ocean boundary was located at the mouth of the estuary ($x=1000 \mathrm{m}$). Structured grids were refined at the free surface and sediment–water interface to accurately capture the intricate hydrodynamic processes. The Courant number in the model was set to 0.5. The grid was 2 m in the along-channel direction and 1 m in the cross-channel direction. In the vertical direction, the grid size ranged from 0.002 m to 0.2 m.

This research focuses on the sediment–water interface and its surrounding area, which is referred to as the sediment–water interface layer. The sediment–water interface layer consists of the benthic boundary layer, the sediment–water interface, and the sediment surface layer. Here, the benthic boundary layer is defined as the region from the layer where the velocity reaches 99% of the average flow velocity of the overlying water to the seabed, encompassing the logarithm law layer and the viscous sublayer. The sediment–water interface refers to the interface between the seabed and the seawater, while the sediment surface layer is the region below the sediment–water interface to a depth of 0.5 m. Monitoring station A, located in the center of the domain ($x=500 \mathrm{m}$, $y=25 \mathrm{m}$), was chosen to display the salinity transport process within the sediment—water interface.

2.3.2. Boundary Conditions

Although an idealized estuary was configured in this study, its 3D features represented most common estuaries in the world. Since the target of the numerical simulation was to investigate the processes of salinity absorption and release at the sediment–water interface influenced by tides, the periodic forcing in the tidal cycle needed to be prescribed at the open boundaries. In order to obtain reasonable tidal information, the tidal current and tidal water level from in situ observations can be adopted for reconstructing open boundary conditions. The data used in this study are shown in Figure 3. Similar observations can be obtained for any local estuary. For the inflow, these data were prescribed at the open boundary. For the outflow, the free radiation open boundary condition was used. For the vertical velocity, it deviates from the commonly assumed logarithmic velocity distribution. Through nonlinear fitting of the coefficients in the double-logarithmic equation, the results conform to a double logarithmic distribution. Both upstream and downstream boundaries of the sediment layer were set as rigid wall boundaries. The boundary conditions used in this study were based on reasonable physical assumptions. These settings were justified by their minimal influence on the computational structure at the calculation points [26] and had been validated by Zhao et al. [27].

3. Validation of Numerical Simulation Results

Our previous numerical results for the water level, vertical mean flow velocity, and vertical velocity distribution are in excellent agreement with the measured data. The simulation and validation of these aspects had been conducted. Additional validation on the variations of water level, flow velocity, saltwater intrusion distance, turbulence kinetic energy, and turbulence dissipation rate had been conducted by Zhao et al. [27] and Li et al. [28]. They all indicate that the model can show robust performance in reproducing the hydrodynamic processes in the water column and sediment. Therefore, here, we avoid duplicating the relevant validation analyses in this study. These simulations effectively captured the salinity transport mechanisms at the salt–freshwater interface.

The operating conditions used in the model validation are same as the physical model conditions established by Zhang et al. [29]. The primary focus of this study was to investigate the hydrodynamics and salinity transport within the sediment–water interface. To ensure the model robustness, the computed salinity was compared and validated against physical flume experiments conducted by Zhang et al. [29]. The resulting salinity time series for the entrance of the estuary (10 m away from the downstream boundary) are shown in Figure 4. This station was selected as a validation site because the salinity in this area undergoes rapid and substantial variations in the vicinity of the estuary entrance. A remarkably consistency in the salinity variation for both the bottom layer and the middle layer can be observed from our numerical simulation and physical experiments. Slight deviations still can be discerned. Owing to the environmental factors affecting the flume data of the physical model, such as solid wall influences on the water flow velocity and tidal wave reflections, these discrepancies were acceptable. The relative errors remained within an acceptable range, confirming the reasonable performance of our numerical model. The salinity variations observed in this study are consistent with the patterns shown by Zhang et al. [30] and Zhu et al. [31]. Notably, we observed periodic fluctuations in salinity with varying amplitudes within the two tidal cycles, similar to the trend in the water level changes. This comparison confirmed the model’s reliability in simulating salinity exchange at the sediment–water interface.

4. Results and Discussion

4.1. Variations of Salinity at the Sediment–Water Interface Layer

To investigate salinity variations at the sediment–water interface layer, simulations were conducted in the idealized estuary from 13 August, 19:00, to 14 August, 06:00, encompassing one tidal cycle of salinity fluctuation. To address the potential issue of model warm-up, we applied appropriate initial conditions that were representative of typical conditions in August. Specifically, the initial salinity distribution and flow conditions were set by the model output from a repeated model run. Thus, the tidal cycle effect could be reasonably included. For short-term salinity migration processes, we focused on the hydrodynamic response of the sediment–water interface to the tidal cycle. Thereby, using the model output from two flood–ebb cycles was reasonable to represent the dynamic response of the entire sediment and overlying water column.

Figure 5 presents the vertical transect of salinity along the channel of the estuary (y = 25 m). The flooding tide presented from 18:00 to 00:00 and the ebb tide lasted from 00:00 to 06:00. The salinity transects at 20:00 and 02:00 were selected to show the pattern at a typical flood tide and ebb tide, respectively. During the flood tide, seawater flows into the estuary, creating a saltwater wedge that influences both the flow and salinity patterns in the region. During the flood tide, as seawater enters the estuary, a saltwater wedge forms at the estuarine boundary. Along this boundary, various parameters such as salinity and flow velocity exhibit changes with increasing distance. The interface between seawater and freshwater gradually moves upstream during this process, leading to shifts in hydraulic characteristics at the interface. As seawater flows into the freshwater region, it accelerates, generating strong impacts and inducing intense turbulence and eddies due to differences in salinity and density. These variations foster the mixing of seawater and freshwater, increasing salinity changes at the interface. Furthermore, the rise in tidal flow causes changes in water depth, affecting flow velocity. Typically, flow velocity decreases with increasing water depth. As a result, in the estuarine area, the flood tide phase involves significant changes in water depth, leading to more pronounced variations in flow velocity. This dynamic hydrodynamic variation intensifies the mixing of saltwater and freshwater, enhancing the fluctuations in salinity. At this stage, the sediment layer begins to absorb the salt as the salinity variations in the sediment–water interface layer evolve.

As the ebb tide began, seawater started to flow out, reversing the hydrodynamic conditions and leading to a different pattern of mixing and salinity fluctuations. During the ebb tide, seawater started to flow out of the estuarine region, causing a reversal in the flow direction and influencing the hydrodynamic characteristics. First, as the seawater flowed out, the water depth gradually became shallower, increasing the water flow velocity. Although the water flow velocity and discharge during the ebb tide were smaller than those during the flood tide, the rate of increase in flow velocity was faster. Second, during the ebb tide, the saltwater wedge moved downstream. In this process, mixing between seawater and freshwater was reduced, resulting in a relatively small variation in salinity at the interface. However, a salinity difference between seawater and freshwater still existed, although the magnitude of the change was not as drastic as that during the flood tide. Reduced mixing between seawater and freshwater leads to smaller variations in salinity. Salt had accumulated to some extent in the sediment layer during this period, causing the salinity variation to be larger than during the flood tide stage.

Figure 6 shows the salinity variations at different depths within the sediment–water interface at monitoring station A. Between 19:00 and 00:00, the flood tide phase occurred and was characterized by the intrusion of saline water that vertically mixed with freshwater, leading to a significant salinity exchange between the sediment layer and the water body. This exchange created a distinct salinity gradient between the sediment surface and overlying water, resulting in a salinity discontinuity. Subsequently, from 01:00 to 06:00, the ebb tide phase ensued, during which saline water and freshwater exhibited distinct vertical stratifications. As the saline water receded, the salinity in the sediment layer became higher than that in the overlying water, prompting salt diffusion into the water body. Moreover, the salinity gradient between the sediment layer and the overlying water was comparatively smaller during the ebb tide than during the flood tide, leading to less pronounced salinity discontinuities.

The spatial-temporal variations of salinity at the sediment–water interface mirrored the patterns observed in the overlying water, fluctuating in tandem with the variations in the overlying water salinity. Because of the salinity gradient between the sediment and overlying water, the salinity at the sediment surface gradually increased but with a noticeable time lag. At depth within the sediment layer, salinity experiences relatively minor variations. However, toward the end of the ebb tide, the salinity in the deep layer can be higher than that of the upper layers. These highlight the dynamic interactions between the sediment and overlying water during different tidal phases. Tidal phases modulate the increase or decrease in salinity and their relative change in different layers. The response of salinity in the sediment layer shows an obvious delay in comparison with those in the overlying water.

4.2. Salinity Flux at the Sediment–Water Interface

The variation in salinity flux over time is shown in

Table 2. Time variation of salinity flux. Negative (positive) values represent the salinity absorption (release) of the sediment layer. The flooding tide was from 18:00 to 00:00 and the ebb tide was from 00:00 to 06:00.

Time (Flood Tide)	19:00	20:00	21:00	22:00	23:00	00:00
F $(\times {10}^{-2}\cdot {\mathrm{m}}^{-2}\cdot {\mathrm{h}}^{-1}$)	0	−2.673	−6.237	−6.459	−6.548	−6.593
Time (Ebb Tide)	01:00	02:00	03:00	04:00	05:00	06:00
F $(\times {10}^{-2}\cdot {\mathrm{m}}^{-2}\cdot {\mathrm{h}}^{-1}$)	−6.637	−1.782	0.7128	1.381	0.668	0.0267

. The flux is calculated based on the Fick’s first law. Negative values indicate the salinity absorption of the sediment layer, whereas positive values indicate the salinity releasing from the sediment layer. The flooding tide was from 18:00 to 00:00 and the ebb tide was from 00:00 to 06:00. The results clearly indicate that the salinity absorption increased continuously from the beginning of the flood tide to the first 2 h of the ebb tide. The absorption flux reached its maximum at 01:00, then decreased rapidly. The salinity release was relatively slow and the magnitude was smaller than the salt absorption.

Figure 7 illustrates the salinity distribution near the sediment–water interface during the flood and ebb tide cycles. The interaction between the sediment layer and the overlying water is governed by the tidal cycle, with the salinity absorption and release by the sediment being highly influenced by the evolving salinity gradient. According to the table and figure, during the initial 1–2 h of the flood tide phase, the salinity absorption rate of the sediment layer increased rapidly and stabilized later. Similarly, during the early stage of the ebb tide, as the saltwater did not completely recede, a salinity concentration gradient persisted, and the sediment layer continued to absorb salinity. However, after the saltwater receded, the sediment layer maintained a relatively stable salinity release. Under these conditions, the salinity gradient was the main factor influencing the salt flux. However, there remained a gap between the absorption and release of salinity by the sediment layer. The impact of the salinity of the bottom sediment layer on the overlying water is a long-term process. In this model, the sediment layer started with zero salinity, and there was an accumulation of salinity in the sediment layer, consistent with the patterns observed in the physical model by Wang [32]. The hydrodynamic characteristics and salinity changes between the sediment layer and the overlying water were mainly influenced by the tidal currents.

4.3. Characteristics of Turbulent Kinetic Energy Within the Sediment–Water Interface Layer

The dynamics of the seabed boundary layer are governed by various factors, with tidal currents playing a central role in modulating its thickness and the associated turbulent parameters. The thickness of the seabed boundary layer varies depending on specific conditions, typically forming a relatively thin layer of water, ranging in thickness from millimeters to meters. Several factors influence this thickness, including water velocity, viscosity, bottom topography, and marine environmental conditions. Real-time thickness was calculated, and vertical averages of the observed physical quantities were computed as necessary. Figure 8 depicts the variations in the mean flow velocity, turbulent kinetic energy, turbulent dissipation rate, and viscous sublayer thickness within the seabed boundary layer over the two tidal cycles. Notably, the variation in turbulent kinetic energy closely mirrored the turbulent dissipation rate, with secondary peaks evident within each tidal cycle. However, it is noteworthy that the peaks in the turbulent kinetic energy occurred approximately 1–2 h later than the peaks in the mean flow velocity.

Tidal currents play a significant role in shaping the turbulent parameters and thickness of the seabed boundary layer. The empirical formula ${\delta }_v=11.1v/u_{\ast }$ proposed by Chriss and Caldwell [33] was used to determine the viscous sublayer thickness (Figure 8d). Our numerical simulation results demonstrated that the velocity distribution within the seabed boundary layer followed a double-logarithmic profile [34]. During the acceleration phase, the boundary layer underwent compression, resulting in a thinner sublayer that expanded and thickened during deceleration. These findings underscore the direct influence of tidal variations on the hydrodynamic conditions in the seabed boundary layer, which contributes to the observed tidal periodicity in the turbulent kinetic energy, turbulent dissipation rate, and viscous sublayer thickness.

The shear force of the overlying water causes turbulent effects at the sediment–water interface, resulting in the pulsating motion of turbulent eddies that cascade into smaller eddies. The transfer rate of the mean flow energy to the turbulent kinetic energy is referred to as the turbulent kinetic energy production rate (P_TKE) [35]. It is calculated as follows:

$$P=-\rho \left(\overline{u^{\prime }{w}^{\prime }}{\displaystyle \frac{\partial u}{\partial z}}+\overline{\overline{v^{\prime }{w}^{\prime }}{\displaystyle \frac{\partial v}{\partial z}}}\right)$$

(10)

Figure 9 shows the temporal variation in turbulent kinetic energy production rate (P_TKE). The distribution of turbulent kinetic energy (P_TKE) in the sediment–water interface layer reflected complex interactions influenced by tidal phases, with flow velocity and water level playing key roles in energy transfer. Within the water layer situated at a distance of 0.3 m from the sediment–water interface, P_TKE ranged from 0.008 to 26.842 Wm⁻³ over one tidal cycle (08–06 T 20:00–08–07 T 07:00, 2018). Notably, P_TKE reached its maximum value when the tidal water level was at its highest or lowest. In contrast, it reached its minimum value during the tidal transition phase. The vertical distribution of P_TKE closely followed half of the tidal cycle. In the sediment layer below the sediment–water interface, within a depth of 0.3 m, P_TKE values spanned from 1.54 × 10⁻¹⁴ to 5.92 × 10⁻¹⁰ Wm⁻³ over one tidal cycle. P_TKE values in the sediment layer were significantly lower than those near the sediment–water interface. This discrepancy can be attributed to the presence of pores in the sediment, which affect the strength of turbulence, resulting in less pronounced turbulent activity in small-scale pores. During the flood tide phase, the maximum value of P_TKE occurred at the time of the highest tidal water level, whereas the minimum value was observed during the transition phase. Conversely, during the ebb tide phase, P_TKE values at the surface of the sediment layer decreased as the flow velocity decreased from the highest tidal water level to the transition phase. Subsequently, from the transition phase to the beginning of the ebb jet, P_TKE values at the surface of the sediment layer initially increased and then decreased with increasing flow velocity. At the lowest tidal water level, P_TKE values at the surface of the sediment layer increased as the flow velocity decreased. Remarkably, the maximum and minimum values of P_TKE occurred during the transition, extreme water level, and ebb jet periods. These findings highlight the significance of the flow velocity and water level as key factors influencing kinetic energy transfer between the overlying water and the surface of the sediment layer.

The turbulent dynamics at the sediment–water interface are influenced by the varying impedance of the sediment, leading to distinct patterns of turbulence during different tidal phases. The turbulence within the interface layer is constrained by the impedance, reducing the turbulent kinetic energy (k) and turbulent dissipation rate (ɛ). Specifically, during the flood tide phase, the turbulent dissipation rate in the interface layer was 2.15 × 10⁻⁶ m²/s at the transition phase, which was smaller than the turbulent dissipation rate in the overlying water. Conversely, during the flood slack phase, the turbulent dissipation rate increased to 5.80 × 10⁻⁴ m²/s, surpassing the overlying water. As for the ebb tide phase, the turbulent dissipation rate decreased to 5.20 × 10⁻⁶ m²/s at the transition phase, again being smaller than that in the overlying water. Subsequently, during the ebb slack phase, the turbulent dissipation rate rose to 2.91 × 10⁻³ m²/s.

In the layer where the sediment met the water, the rates of turbulent kinetic energy production were as follows: 8.70 × 10⁻⁴ Wm⁻³ during the transition phase of the flood tide, 0.25 Wm⁻³ during the slack phase of the flood tide, 9.77 × 10⁻⁴ Wm⁻³ during the transition phase of the ebb tide, and 1.29 Wm⁻³ during the slack phase of the ebb tide. Remarkably, the turbulent kinetic energy production rates in the sediment–water interface layer were two to three orders of magnitude greater than the turbulent kinetic energy dissipation rates, clearly indicating the dominance of turbulent production over turbulent dissipation within this layer.

Figure 10a presents the variation of convective mass flux per unit area, calculated using the formula $F=\rho w$. Figure 10b presents the variation of diffusive flux per unit area at the interface, calculated using the formula $D=\left(\mu {\displaystyle \frac{\partial w}{\partial z}}-\rho \overline{w^{\prime }{w}^{\prime }}\right)/{\displaystyle \frac{\partial w}{\partial z}}/\partial z$. The comparison between convective and diffusive fluxes revealed the overwhelming dominance of turbulent diffusion in the water exchange process at the sediment–water interface. Figure 10 reveals a striking contrast between the magnitudes of the convective and diffusive fluxes, with the convective being significantly smaller by several orders of magnitude than the diffusive. It is important to note that this study does not account for the impact of terrain variation. The presence of pores and the hindering influence of the viscous sublayer substantially restrict the magnitude of molecular diffusion in the diffusive process, making it much smaller than turbulent diffusion. Consequently, turbulent diffusion emerges as the dominant form during the water exchange process, while the convective effect remains relatively weak.

In comparison with the previous study, our case only focuses on the short-term variations. Actually, environment problems associated with sediment absorption and release usually cause long-term effects through accumulative diffusion. For example, Liu et al. [36] revealed that tidal cycling influenced the metal release from estuarine sediment. The combined effect of sediment and hydrodynamic processes controlled the long-term accumulation. From the perspective of short-term effect, water and salinity exchange within the sediment–water interface occurred predominantly via turbulent diffusion, with extreme turbulent kinetic energy production rates appearing during the tidal reversal, flood, and ebb stages. This information is valuable for understanding the complex processes within the sediment–water interface.

4.4. Broad Implications

This study provides insights towards understanding the material exchange around the sediment–water interface. Firstly, our results clearly illustrated that the sediment layer is important in modulating the salinity transport in the upper reach of an estuary, which was usually excluded by previous studies. Secondly, we introduced a numerical model that incorporates bidirectional flow in the interface layer. The resulting dynamic impact on the salinity diffusion might be underestimated by previous models. Thirdly, by including sediment porosity, this research addresses a gap, as conventional models fail to account for how porosity affects salinity diffusion in sediment and water. It should be noted that, for practical application, the transport processes simulated here are not only specific to the salinity absorption and release, but also applicable to dissolvable pollutants and nutrients. Thereby, this research offers practical applications for substance exchange around the sediment–water interface. The method presented here can be used to explore broad interfacial exchange in other coastal areas, which will help to solve problems in the coastal environment.

5. Conclusions

This study presented a comprehensive baroclinic fluid dynamic model for investigating the flow and salt transport characteristics around the sediment–water interface under tidal forcing. The model treated the sediment layer as a fluid layer and incorporated the Darcy–Forchheimer correction term. In an idealized estuary, this model showed good performance in revealing the dynamic process that controls the water and salinity exchange within the sediment–water interface.

The results indicated that the turbulent kinetic energy and turbulent dissipation rate showed semidiurnal tidal variations in the upper water layer. Specifically, during tidal reversal and ebb phases, these parameters decreased with increasing water depth, while, during tidal flood and neap phases, they increased near the sediment–water interface. Additionally, the seabed boundary layer exhibited lagged extreme values compared to the mean flow velocity. Turbulent dissipation at the sediment–water interface exhibited phase-dependent variations, being stronger during the tidal flood phase and weakening during the tidal reversal and ebb phases. These variations were attributed to the constraining impact of pore spaces and the inhibiting effect of the viscous sublayer. Turbulent diffusion emerged as the dominant water exchange mechanism, with molecular diffusion playing a smaller role. Salinity changes in the sediment layer closely followed those in the overlying water during the tidal flood phase, driven by convective diffusion with larger salinity gradients. In contrast, during the tidal ebb phase, the sediment layer exhibited higher salinity than the overlying water, leading to weaker diffusion effects due to a smaller salinity gradient.

These results provided insights into the hydrodynamic and salinity transport characteristics of the sediment–water interface layer influenced by tides in the upper estuary where freshwater and saltwater meet. Future research should consider incorporating additional factors such as wind fields, waves, and sediment types in order to obtain a more comprehensive understanding of the exchange occurring in the sediment–water interface layer. This will contribute to a broader understanding of estuarine and nearshore environments.