Dynamic Response Analysis and Liquefaction Potential Evaluation of Riverbed Induced by Tidal Bore

Abstract

Tidal bores, defined by sudden upstream surges of tidal water in estuaries, exert significant hydrodynamic forces on riverbeds, leading to complex sedimentary responses. This study examines the dynamic response and liquefaction potential of riverbeds subjected to tidal bores in macro-tidal estuaries. An analytical model, developed using the generalized Biot theory and integral transform methods, evaluates the dynamic behavior of riverbed sediments. Key factors such as permeability, saturation, and sediment properties are analyzed for their influence on momentary liquefaction. The results indicate that fine sand reduces liquefaction risk by facilitating pore water discharge, while silt soil increases sediment instability. Additionally, the study reveals that pressure gradients induced by tidal bores can trigger momentary liquefaction, with the maximum liquefaction depth predicted based on horizontal pressure gradients being five times that predicted based on vertical pressure gradients. This research highlights the critical role of sediment characteristics in riverbed stability, providing a comprehensive understanding of the interactions between tidal bores and riverbed dynamics. The findings contribute to the development of predictive models and guidelines for managing the risks of tidal bore-induced liquefaction in coastal and estuarine environments.

1. Introduction

Tidal bores are remarkable natural phenomena in estuaries and bays [1], characterized by a sudden upstream surge of tidal water against the ebb flow [2,3]. This occurs due to the deformation of tidal waves in shallow water, creating a steep front or a series of wave trains that move upstream along rivers or through convergent bays [4,5]. Recent studies have identified tidal bores in 117 rivers across 25 countries [6]. As depicted in Figure 1, notable examples include the Qiantang River in China, the Amazon River in Brazil, Cook Inlet in Alaska, the Hooghly River in India, the Batang Lupar River in Malaysia, the Kampar River in Indonesia, and the Sittang River in Myanmar [7]. During a tidal bore, water levels can rise by over 3 m within seconds, with flow velocities reaching up to 12 m/s as the river rapidly transitions from ebb to flood tide [8,9]. These powerful waves exert increased pressure on the riverbed, affecting pore water pressure within the sediment and potentially leading to liquefaction. The complex interaction between hydrodynamic forces and sedimentary responses highlights the power and complexity of tidal bores, making them a compelling subject for research.

Since the nineteenth century, tidal bore research has advanced from descriptive studies to dynamic mechanism analyses [10]. Early work focused on in situ observation [11], while contemporary research integrates field investigations [12,13,14], laboratory experiments [2,15,16,17], and numerical simulations [17,18,19] to examine the dynamics and turbulent structures of the tidal bore front [20,21,22]. Despite these advancements, limitations in observational and simulation techniques hinder our understanding of the interactions between tidal bores and riverbeds, necessitating further research to improve analytical accuracy and reliability in this field.

Liquefaction, as defined by the Geotechnical Earthquake Engineering Committee of the American Society of Civil Engineers, is the act or process of transforming any substance into a liquid [23]. In cohesionless soils, the shift from a solid to a liquefied state results from increased pore pressure and a decrease in effective stress [24]. Seabed liquefaction under wave loading has been extensively studied, with comprehensive reviews by Jeng [25] and Lin et al. [26]. In natural environments, two mechanisms of wave-induced liquefaction exist [27,28,29]: momentary liquefaction and residual liquefaction. Momentary liquefaction primarily arises from the attenuation and phase lag of oscillatory pore pressure along the seabed depth [30], associated with the momentary volumetric strain of sediments. During this process, pore pressure does not accumulate; however, a significant upward seepage force is generated at wave troughs due to vertical pressure gradients, causing sediment particle movement and overturning [31,32,33]. Residual liquefaction refers to the instability and failure induced by the cumulative rise in pore water pressure and the reduction in effective stress within seabed sediments under wave action [34,35,36], related to the permanent compaction deformation of sediments.

Recent research has also explored tidal flat liquefaction caused by tidal bores from a sedimentological perspective [37]. Tessier and Terwindt [38] investigated significant soft sediment deformations in the tidal channel of Mont-Saint-Michel Bay, characterized by folds and tobacco-pouch structures, attributing these features to the liquefaction of unconsolidated sediments triggered by the passage of a tidal bore. In the Turnagain Arm area of Alaska, Greb and Archer [39] found that tidal bores are a primary cause of soft sediment deformation structures similar to those formed by earthquake-induced sand liquefaction. Fan et al. [40] analyzed sedimentary structures in the northern tidal flat of the Qiantang River estuary, indicating that the rapid increase in water pressure and wave impact induced by tidal bores can cause liquefaction, resulting in features such as convolute bedding and drainage structures. Tidal bore-induced soft sediment deformations, primarily involving liquefaction, are evident in contorted beds and dewatering structures, with geometries ranging from 2D folds to 3D chaotic envelopes. Due to the occurrence of liquefaction, tidal bore-induced sedimentary signatures are more conspicuous than regular tidal facies and are easily recognized within tidal channel infilling successions.

The dynamic response and liquefaction potential of soils in coastal and estuarine environments are strongly influenced by factors such as permeability, saturation, and sediment properties, including grain size and fines content [41]. High permeability facilitates the rapid dissipation of pore water pressure, reducing the risk of liquefaction, whereas low permeability increases susceptibility. Saturation further affects this process, with unsaturated soils more prone to increased pressure gradients [42]. Additionally, sediment characteristics play a critical role: coarse-grained soils generally have a lower liquefaction potential, but the presence of fines can increase susceptibility by restricting drainage [43]. Understanding the interaction of these factors is key to predicting soil behavior under dynamic conditions, such as those caused by tidal bores.

Analytical solutions for the response of seabeds to wave loading have garnered significant attention since the 1940s [25], leading to the development of various computational models. These models include uncoupled (drained) models, consolidation (quasi-static) models, dynamic models, and poro-elastoplastic models [41]. The core of seabed response analysis under wave loading involves examining the coupled effects of soil stress and pore fluid pressure induced by wave trains [44]. Building on Biot’s pioneering work on consolidation theory [45,46,47], analytical solutions for seabed response under wave loading have been proposed for both two-dimensional [31,48] and three-dimensional spaces [42,49]. A set of analytical solutions for wave-induced seabed response under cnoidal waves was developed to evaluate the impact of high-order nonlinearity in shallow water [44,50], with a parametric study revealing that this nonlinearity significantly influences shallow water waves and amplifies the effects of soil characteristics and liquefaction potential. All the aforementioned analytical studies were formulated based on the assumption of simple harmonic waves. However, the transient response of the seabed under the abrupt application of wave loads has not been accounted for in these analyses.

Wave action can lead to localized liquefaction and erosion of the seabed, which is particularly evident in transient responses. Selecting appropriate integral transforms facilitates the conversion of differential and integral equations into solvable algebraic equations. Consequently, several analytical or semi-analytical solutions for seabed transient response under wave action have utilized integral transform methods [51,52]. For layered seabeds, the state-space method can be extended to improve computational efficiency [53,54]. In macro-tidal estuary regions, a tidal bore is a typical transient wave. It is characterized by a sudden rise in water level at the front of the tidal wave. When a breaking bore occurs, the front becomes nearly vertical and propagates upstream as a form of moving load. However, previous studies have often been limited by inadequate consideration of the complex interactions between tidal bores and the riverbed, particularly in terms of transient response and the liquefaction potential. These limitations highlight the need for more comprehensive models that can accurately capture the dynamic behavior of the riverbed under such extreme conditions. The purpose of this study is to address these gaps by developing a more detailed analysis of the transient response of the riverbed and to better understand the mechanisms leading to momentary liquefaction induced by pressure gradient.

The structure of this paper is as follows: Section 2 provides a detailed description of the theoretical model and methodologies used in this study, including the governing equations and the integral transform method. Section 3 presents the results of the study, with subsections focusing on model verification, riverbed response analysis, seepage characteristics, and liquefaction analysis. The discussion and limitations of these findings are addressed in Section 4. Finally, Section 5 concludes the study, summarizing the key contributions and outlining directions for future research. The findings of this study are expected to provide new insights into the stability of riverbeds under extreme hydrodynamic conditions, contributing to improved predictive models and risk assessment frameworks for coastal engineering and management practices.

2. Materials and Methods

2.1. Theoretical Model

The generalized Biot theory provides a comprehensive framework for describing the interaction between skeleton deformation and pore fluid flow in multiphase porous media under transient loading conditions [55]. This study employs a theoretical model to analyze the dynamic response of a riverbed subjected to tidal forces, based on the generalized Biot theory. The model operates under the following assumptions:

• The riverbed is considered horizontal, homogeneous, highly saturated, isotropic in permeability, and of finite thickness.
• The compressibility coefficients of the soil skeleton and pore water are constant, while the soil particles are assumed to be incompressible.
• The stress–strain relationship of the soil skeleton adheres to Hooke’s law.
• The seepage of pore water complies with Darcy’s law.
• The energy loss during tidal bore propagation is disregarded.
• The water pressure exerted on the riverbed surface is equivalent to the water pressure on the surface of an impermeable, rigid horizontal riverbed at the same depth.

2.1.1. Governing Equations

In a two-dimensional Cartesian coordinate system (x, z), we analyze the propagation of an incident tidal bore at a water depth d over a submerged porous riverbed with a thickness h, as illustrated in the side view sketch in Figure 2. Following Zienkiewicz et al. [55], the equilibrium equations for a unit total volume can be expressed as

$${\sigma }_{ij,j}+{\rho }_0{g}_i={\rho }_0{\ddot{u}}_i+{\rho }_w{\ddot{w}}_i,$$

(1)

where σ_ij is the total stress component, u_i is the displacement of the solid matrix, w_i is the average relative displacement of the fluid to the solid, ρ₀ represents the combined density, ρ₀ = (1 − n) ρ_s + n ρ_w, ρ_w, and ρ_s denotes the fluid and solid density, respectively, n is the porosity of the solid phase, and g_i refers to the body force acceleration, typically defined by the relevant components of gravity.

The equilibrium equations for the fluid phase, incorporating viscous resistance as defined by Darcy’s law and assuming isotropic permeability, are established as follows:

$$-p{,}_i+{\rho }_w{g}_i={\rho }_w{\ddot{u}}_i+{\displaystyle \frac{{\rho }_w}{n}}{\ddot{w}}_i+{\displaystyle \frac{{\rho }_{w}g}{k_z}}{\dot{w}}_i,$$

(2)

where p is the pore pressure in the riverbed; k_z is the permeability coefficient in the z direction, considering a seabed with isotropic permeability k_z = k_x; and k_x is the permeability coefficient in the x directions.

The continuity equation for mass conservation in the fluid is as follows:

$$n\beta \hspace{0.33em}\dot{p}+{\dot{u}}_{i,i}+{\dot{w}}_{i,i}=0,$$

(3)

where β is the volumetric compressibility coefficient of the fluid, which can be expressed as follows for a highly saturated riverbed [56]

$$\beta ={\displaystyle \frac{1}{K_w}}+{\displaystyle \frac{1-S_r}{P_{w0}}},$$

(4)

where K_w is the true modulus of elasticity of water, P_w0 is the absolute water pressure, and S_r is the degree of saturation.

In plane strain, the constitutive relationship is expressed as follows:

$${\sigma }_{ij}=\lambda {\delta }_{ij}\theta +G(u_{i,j}+u_{j,i})-{\delta }_{ij}p,$$

(5)

where λ and G are the Lamé constants, δ_ij is the Kronecker delta denotation, and θ is the volume strain, θ = u_i,i.

2.1.2. Integral Transform Method

Omitting the static gravity terms and focusing only on excess pressures and stresses above the static state, apply the Laplace transform to generalized Biot’s Equations (1) through (3) and constitutive Equation (5), respectively [41]

$$\widehat{f}(x,z,s)={\displaystyle {\int }_{\hspace{0.17em}0}^{\hspace{0.17em}+\infty }f(x,z,t)e^{-st}dt},$$

(6)

where f(x, z, t) is the original function in the time domain; $\widehat{f}(x,z,s)$ is the Laplace transformed function; and s is the complex frequency variable in the Laplace domain. After simplification, the generalized Biot’s equations can be written in the following form:

$$G{\nabla }^2{\widehat{u}}_x+(\lambda +G){\displaystyle \frac{\partial \widehat{\theta }}{\partial x}}-(1-\vartheta ){\displaystyle \frac{\partial \widehat{p}}{\partial x}}-s^2({\rho }_0-{\rho }_w\vartheta ){\widehat{u}}_x=0,$$

(7)

$$G{\nabla }^2{\widehat{u}}_z+(\lambda +G){\displaystyle \frac{\partial \widehat{\theta }}{\partial z}}-(1-\vartheta ){\displaystyle \frac{\partial \widehat{p}}{\partial z}}-s^2({\rho }_0-{\rho }_w\vartheta ){\widehat{u}}_z=0,$$

(8)

$${\nabla }^2\widehat{p}-{\displaystyle \frac{n\beta {\rho }_w{s}^2}{\vartheta }}\widehat{p}-{\displaystyle \frac{1-\vartheta }{\vartheta }}{\rho }_w{s}^2\widehat{\theta }=0,$$

(9)

where

$$\vartheta ={\displaystyle \frac{{\rho }_w{s}^{2}n}{{\rho }_w{s}^2+nbs}},$$

(10)

$$b={\displaystyle \frac{{\rho }_{w}g}{k_z}}.$$

(11)

From Equations (7) through (9), the following equations can be derived

$${\nabla }^4\widehat{p}+{\beta }_1{\nabla }^2\widehat{p}+{\beta }_2\widehat{p}=0,$$

(12)

where M = λ + 2G,

$${\beta }_1=-{\displaystyle \frac{s^2(Mn\beta {\rho }_w+{\rho }_w-2{\rho }_w\vartheta +\vartheta {\rho }_0)}{M\vartheta }},$$

(13)

$${\beta }_2={\displaystyle \frac{({\rho }_0-{\rho }_w\vartheta )n\beta s^4{\rho }_w}{M\vartheta }}.$$

(14)

Apply the following Fourier transforms to Equations (7), (8), and (12):

$$\tilde{f}(\xi ,z,s)={\displaystyle {\int }_{\hspace{0.17em}-\infty }^{\hspace{0.17em}+\infty }\widehat{f}(x,z,s)e^{i\xi x}}dx,$$

(15)

where $\widehat{f}(x,z,s)$ is the original function, $\tilde{f}(\xi ,z,s)$ is the Fourier transformed function, ξ is the frequency variable, and i is the imaginary unit.

After converting the partial differential equation into an ordinary differential equation, the operator method is used to obtain the following general solution:

$$Y_1(\xi ,z,s)=T_1(\xi ,z,s)X_0(\xi ,z,s),$$

(16)

where

$$Y_1={[\tilde{p}(\xi ,z,s)\begin{array}{c} \end{array}\tilde{\theta }(\xi ,z,s)\begin{array}{c} \end{array}\hspace{0.17em}{\tilde{u}}_z(\xi ,z,s)\begin{array}{c} \end{array}\xi (i{\tilde{u}}_x)(\xi ,z,s)]}^T,$$

(17)

$$X_0={[A_1{e}^{{\gamma }_{1}z}\begin{array}{c} \end{array}{A}_2{e}^{-{\gamma }_{1}z}\begin{array}{c} \end{array}{A}_3{e}^{{\gamma }_{2}z}\begin{array}{c} \end{array}{A}_4{e}^{-{\gamma }_{2}z}\begin{array}{c} \end{array}{A}_5{e}^{{\gamma }_{3}z}\begin{array}{c} \end{array}{A}_6{e}^{-{\gamma }_{3}z}]}^T,$$

(18)

$$T_1=\left[\begin{array}{cccccc}1& 1& 1& 1& 0& 0\\ -G{\kappa }_1& -G{\kappa }_1& -G{\kappa }_2& -G{\kappa }_2& 0& 0\\ {\gamma }_1{\alpha }_1& -{\gamma }_1{\alpha }_1& {\gamma }_2{\alpha }_2& -{\gamma }_2{\alpha }_2& 1& 1\\ {\chi }_1& -{\chi }_1& {\chi }_2& -{\chi }_2& {\gamma }_3& -{\gamma }_3\end{array}\right],$$

(19)

$${\gamma }_j=\sqrt{{\xi }^2-L_j^2},$$

(20)

$$L_1^2={\displaystyle \frac{1}{2}}({\beta }_1+\sqrt{{\beta }_1^2-4{\beta }_2}),$$

(21)

$$L_2^2={\displaystyle \frac{1}{2}}({\beta }_1-\sqrt{{\beta }_1^2-4{\beta }_2}),$$

(22)

$$L_3^2={\displaystyle \frac{({\rho }_w-{\rho }_0)\vartheta s^2}{G}},$$

(23)

$${\kappa }_j={\displaystyle \frac{\vartheta L_j^2+n\beta {\rho }_w{s}^2}{G(1-\vartheta ){\rho }_w{s}^2}},$$

(24)

$${\alpha }_j={\displaystyle \frac{{\kappa }_j(\lambda +G)+(1-\vartheta )/G}{L_3^2-L_j^2}},$$

(25)

$${\chi }_j=G{\kappa }_j+{\gamma }_j^2{\alpha }_j.$$

(26)

Using Equation (16) and the constitutive Equation (5), the total stress components can be expressed in the Laplace–Fourier transform domain as follows:

$$Y_2(\xi ,z,s)=T_2(\xi ,z,s)X_0(\xi ,z,s),$$

(27)

where

$$Y_2={[{\tilde{\sigma }}_x(\xi ,z,s)\begin{array}{c} \end{array}{\tilde{\sigma }}_z(\xi ,z,s)\begin{array}{c} \end{array}i\xi {\tilde{\tau }}_{xz}(\xi ,z,s)]}^T,$$

(28)

$$T_2=\left[\begin{array}{cccccc}{c}_1& c_1& c_2& c_2& -2G{\gamma }_3& 2G{\gamma }_3\\ d_1& d_1& d_2& d_2& 2G{\gamma }_3& -2G{\gamma }_3\\ g_1& -g_1& g_2& -g_2& {\gamma }_3^2+{\xi }^2& {\gamma }_3^2+{\xi }^2\end{array}\right],$$

(29)

where

$$c_j=-(\lambda {\kappa }_j+2{\chi }_j)G,$$

(30)

$$d_j=-(M{\kappa }_j-2{\chi }_j)G,$$

(31)

$$g_j={\gamma }_j({\xi }^2{\alpha }_j+{\chi }_j)G.$$

(32)

At this stage, the analytical solution of the dynamic response control equation for the riverbed in the Laplace–Fourier transform domain is obtained. In Equations (16) and (27), A_j (j = 1, ..., 6) in $X_0(\xi ,z,s)$ are undetermined constants, which will be solved based on the specific problem in the next subsection.

2.2. Solving the Boundary Value Problem

According to [57], the tidal bore-induced dynamic water pressure at the water–riverbed interface can be expressed as [52]

$$\hspace{0.17em}{p}_b=\rho gH\left\{1-\tanh b(x-ct)\right\}/2,$$

(33)

where

$$b={(3H/4d^3)}^{1/2}.$$

(34)

The boundary conditions for the riverbed response problem under the action of tidal bore can be described as follows:

$$z=0,\hspace{0.17em}\hspace{0.17em}{\sigma }_z+p=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\tau }_{xz}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}p=p_b,$$

(35)

$$z=-h,\hspace{0.17em}\hspace{0.17em}{u}_x=u_z=0,\hspace{0.17em}\hspace{0.17em}{\displaystyle \frac{\partial p}{\partial z}}=0.$$

(36)

Applying the Laplace transform followed by the Fourier transform to Equations (35) and (36) yields

$${\tilde{\sigma }}_z(\xi ,0,s)=-{\tilde{p}}_b\begin{array}{c},\end{array}\tilde{p}(\xi ,0,s)={\tilde{p}}_b\begin{array}{c},\end{array}{\tilde{\tau }}_{xz}(\xi ,0,s)=0,$$

(37)

$${\tilde{u}}_x(\xi ,-h,s)=0\begin{array}{c},\end{array}{\tilde{u}}_z(\xi ,-h,s)=0\begin{array}{c},\end{array}{\displaystyle \frac{\partial \tilde{p}}{\partial z}}(\xi ,-h,s)=0,$$

(38)

where

$${\tilde{p}}_b=2\pi \rho gH\left\{{\displaystyle \frac{1}{s}}\delta (\xi )+{\displaystyle \sum _{k=1}^{\infty }\frac{{(-1)}^k}{s+2kbc}}\delta (\xi -2kbi)\right\}.$$

(39)

The detailed derivation process of Equation (39) can be found in Appendix A.

Substitute Equations (37) and (38) into Equations (16) and (27), solve the system of equations, and obtain the expressions for the coefficients A_j (j = 1, ..., 6)

$$A_j={\tilde{p}}_b{B}_j\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(j=1, \cdots , 6),$$

(40)

where coefficients B_j are functions of ξ and s; the specific forms can be found in Appendix B.

The solution for the riverbed response under tidal bore in the Laplace–Fourier transform domain is given as follows:

$$\tilde{f}(\xi ,z,s)=2\pi \rho gH\left\{{\displaystyle \frac{1}{s}}\delta (\xi )+{\displaystyle \sum _{k=1}^{\infty }\frac{{(-1)}^k}{s+2kbc}}\delta (\xi -2kbi)\right\}\Gamma (\xi ,z,s),$$

(41)

where $\Gamma (\xi ,z,s)$ is a rational function of ξ, s and z, and the specific expression varies depending on the response variable.

By applying the inverse Fourier transform to Equation (41) and utilizing the sifting property of the δ function, we obtain the following result:

$$\widehat{f}(x,z,s)=\rho gH\left\{{\displaystyle \frac{\Gamma (0,z,s)}{s}}+{\displaystyle \sum _{k=1}^{\infty }\frac{{(-1)}^k\Gamma (2kbi,z,s)}{s+2kbc}}{e}^{2kbx}\right\}.$$

(42)

By applying the inverse Laplace transform, using the numerical method for inverse Laplace transform proposed by Durbin [58], a corresponding program was developed to perform the calculations and obtain the results.

Meanwhile, let

$${\widehat{n}}_1(x,z,s)={\displaystyle \frac{\Gamma (0,z,s)}{s}},$$

(43)

$${\widehat{n}}_2(x,z,s)={\displaystyle \frac{{(-1)}^k\Gamma (2kbi,z,s)}{s+2kbc}}{e}^{2kbx}.$$

(44)

Considering Equation (43), it can be written as a reduced rational function

$${\widehat{n}}_1(x,z,s)={\displaystyle \frac{W_{11}(x,z,s)}{W_{12}(x,z,s)}}.$$

(45)

where coefficients W₁₁(x, z, s) and W₁₂(x, z, s) are irreducible proper rational functions.

Assuming the denominator has N zeros, s₁, s₂, ... s_N, with s = 0 clearly being one of them, let us assume s₁ = 0. Then, according to Heaviside’s first expansion theorem for the inverse Laplace transform [59], the following result can be obtained:

$$n_1(x,z,t)=\Gamma (0,z,0)+{\displaystyle \sum _{n=2}^N\frac{W_{11}(x,z,s_n)}{{W^{\prime }}_{12}(x,z,s_n)}}\exp [s_{n}t].$$

(46)

where

$${W^{\prime }}_{12}(x,z,s)={\displaystyle \frac{d}{ds}}\left[W_{12}(x,z,s)\right].$$

(47)

Similarly, from Equation (44), we obtain the following result:

$$n_2(x,z,t)={\displaystyle \sum _{k=1}^{\infty }{(-1)}^k\Gamma (2kbi,z,-2kbc)e^{2kb(x-ct)}}+{\displaystyle \sum _{k=1}^{\infty }{\displaystyle \sum _{n=2}^N\frac{W_{21}(x,z,s_n)}{{W^{\prime }}_{22}(x,z,s_n)}}\exp [s_{n}t]}.$$

(48)

where coefficient W₂₁(x, z, s) and W₂₂(x, z, s) are irreducible proper rational functions.

The second term on the right-hand side of Equations (46) and (48) decays rapidly in the time domain, i.e.,

$${\left.{\displaystyle \sum _{n=2}^N\frac{W_{11}(x,z,s_n)}{{W^{\prime }}_{12}(x,z,s_n)}}\exp [s_{n}t]\right|}_{t\to \infty }\to 0,$$

(49)

$${\left.{\displaystyle \sum _{k=1}^{\infty }{\displaystyle \sum _{n=2}^N\frac{W_{21}(x,z,s_n)}{{W^{\prime }}_{22}(x,z,s_n)}}\exp [s_{n}t]}\right|}_{t\to \infty }\to 0.$$

(50)

Therefore, the steady-state solution for the dynamic response of the riverbed under tidal bore action is given by

$$f^{\ast }(x,z,t)=\rho gH\left\{\Gamma (0,z,0)+{\displaystyle \sum _{k=1}^{\infty }{(-1)}^k\Gamma (2kbi,z,-2kbc)e^{2kb(x-ct)}}\right\}.$$

(51)

2.3. Liquefaction Triggered Criterion

The propagation of tidal bores leads to sudden changes in water depth and rapid fluctuations in physical variables such as velocity and pressure fields. When a tidal bore arrives, the flow quickly shifts from ebb to flood in a very short time. For instance, in the Qiantang River in China, a tidal bore can cause a sudden rise of 2–4 m within just 5 s. The resulting pressure gradient is strong enough to cause the momentary failure of riverbed sediments, potentially leading to the liquefaction of the bed material. This study examines the momentary liquefaction of riverbed triggered by the vertical and horizontal pressure gradients.

When fluid flows through a porous medium, momentary liquefaction can occur if the vertical or horizontal pressure gradient exceeds a certain threshold. Bear [60] proposed that when the vertical resultant force acting on the particles composing a porous medium becomes zero, localized momentary liquefaction can occur

$$-{\displaystyle \frac{1}{{\rho }_{w}g}}\left({\displaystyle \frac{\partial p}{\partial z}}\right)={\displaystyle \frac{s_z}{k_z}}>{\displaystyle \frac{{{\gamma }^{\prime }}_s}{{\gamma }_w}},$$

(52)

where S_z represents the vertical seepage velocity, γ_s′ stands for the effective unit weight of the soil particles, and γ_w is the unit weight of water.

Effective unit weight is the unit weight of the soil when buoyancy from water is considered. It is calculated by subtracting the buoyant force from the gravitational force acting on the soil particles within a unit volume of soil. For most silty soil,

$${\displaystyle \frac{{{\gamma }^{\prime }}_s}{{\gamma }_w}}\approx 0.7.$$

(53)

Meanwhile, Madsen [61] suggested that when the horizontal pressure gradient exceeds the intergranular stress, localized momentary liquefaction can also occur

$$-{\displaystyle \frac{1}{{\rho }_{w}g}}{\left({\displaystyle \frac{\partial p}{\partial x}}\right)}_{crit}={\left({\displaystyle \frac{s_x}{k_x}}\right)}_{crit}={\displaystyle \frac{{{\gamma }^{\prime }}_s}{{\gamma }_w}}\tan \varphi ,$$

(54)

where S_x represents the horizontal seepage velocity and φ is the internal friction angle.

The internal friction angle of soil represents the soil’s inherent frictional properties. It is generally considered to comprise two components: the surface friction between soil particles and the interlocking force generated by the embedding and interlocking of particles. Thus, it is not only related to the material properties of the soil particles but also to the stress conditions. It is typically, for clay, 0 ≤ φ ≤ 20°; for saturated silty soil, it is 0 ≤ φ ≤ 35°. For the silt in the Qiantang River, a particular value of φ = 35° was chosen for the calculations.

$${\displaystyle \frac{{{\gamma }^{\prime }}_s}{{\gamma }_w}}\tan \varphi \approx 0.5.$$

(55)

3. Results

3.1. Verification

Assuming that the soil skeleton is incompressible, meaning a rigid soil bed, the dynamic water pressure on the bed surface and the seepage velocity under the influence of a tidal bore can be approximately expressed as [57]

$$p_b/\rho gd=\epsilon \left\{1-\tanh bX\right\}/2,$$

(56)

$${\displaystyle \frac{s_x}{k_x}}={\displaystyle \frac{\epsilon \sqrt{3\epsilon }}{4}}{\sec \mathrm{h}}^{2}bX,$$

(57)

$${\displaystyle \frac{s_z}{k_z}}={\displaystyle \frac{3h{\epsilon }^2}{4d}}\tanh bX{\sec \mathrm{h}}^{2}bX,$$

(58)

where X = x − ct, ε = H/d.

The results presented in Equation (51) can easily degenerate to the case of a rigid soil bed. A comparison between the degenerate solution and the analytical solution by Packwood and Peregrine [57] is shown in Figure 3. The calculation parameters are set as follows: tidal bore height H = 1.2 m, water depth d = 3.0 m, and soil bed thickness h = 3.0 m. For the rigid soil bed, G→∞, β→0. As illustrated in the Figure 3, the two solutions show good agreement.

3.2. Riverbed Responses

The degradation solutions of pore water pressure and seepage velocity caused by tidal bores have been validated using the analytical solution provided by Packwood and Peregrine [57]. Here, we provide more general and broader results regarding the riverbed response induced by tidal bores. In the example of riverbed response to be presented, the following common values are assumed: water depth d = 3 m, tidal bore height H = 3 m, riverbed thickness h = 9 m, permeability coefficient k_x = k_z = 10⁻³ m/s, Poisson’s ratio μ = 1/3, porosity n = 0.3, and shear modulus G = 10⁷ N/m². Given Poisson’s ratio μ and the shear modulus G, the Lame constant λ can be calculated using the following formula:

$$\lambda ={\displaystyle \frac{2\mathrm{G}\mu }{1-2\mu }}.$$

(59)

The response of riverbed pore water pressure and effective stresses to tidal bore events is significantly influenced by the degree of saturation. Figure 4 presents the contours of excess pore water pressure within the riverbed for two saturation conditions: fully saturated (S_r = 1) and nearly fully saturated (S_r = 0.99). The results indicate that in the fully saturated scenario, the pore spaces are entirely filled with water, leading to higher excess pore water pressures near the riverbed surface. This uniform pressure distribution is attributed to the incompressibility of water. Conversely, in the nearly saturated case, the presence of a minimal amount of air in the pores reduces the overall pore pressure and results in more dispersed pressure contours. This suggests a more gradual dissipation of energy within the soil matrix. Unlike waves, the pore water pressure induced by tidal bores is always positive; it decreases with increasing depth. In an unsaturated riverbed, the gradients of pore pressure change both horizontally and vertically are greater than those in a saturated riverbed.

Figure 5 illustrates the contours of effective stress for two saturation conditions: fully saturated and nearly fully saturated. In the fully saturated scenario, the stress distribution is more localized, with concentrated regions of effective stresses occurring directly beneath the front of the tidal bore. In contrast, the nearly saturated condition results in a more dispersed stress pattern, with higher magnitudes of both compressive and tensile stresses extending across a broader area of the riverbed. The presence of air in the pores reduces pore water pressure, leading to an increase in soil compressibility and a wider distribution of stresses. When the riverbed is fully saturated, the effective stress distribution induced by the tidal bore is generally symmetrical about the front at X = 0. The impact of saturation on effective normal stress is greater than that on shear stress.

3.3. Seepage Characteristics

Figure 6 illustrates the time histories of tidal bore-induced horizontal and vertical seepage velocities within the fully saturated riverbed at three elevations: z = 0, z = −0.25 h, and z = −0.5 h. The dynamic water pressure (Figure 6a) is uniformly affected across all elevations. However, the horizontal and vertical seepage velocity (Figure 6b,c) shows a pronounced decrease with depth, being highest at the seabed surface and minimal at z = −0.5 h.

Figure 7 illustrates the contours of both vertical and horizontal seepage velocities induced by a tidal bore within the riverbed for two different degrees of saturation: fully saturated and nearly fully saturated. In the fully saturated case (Figure 7a,c), both vertical and horizontal seepage velocities are concentrated near the surface, showing higher peak velocities and sharper distributions, with the vertical seepage displaying alternating upward and downward flow zones and the horizontal seepage exhibiting intense and localized flow. In contrast, when the saturation is slightly reduced (Figure 7b,d), the seepage velocity distributions become broader and more diffuse, with lower peak velocities extending deeper into the riverbed. This difference highlights the impact of saturation on seepage dynamics, where a fully saturated riverbed exhibits more intense and localized seepage, while reduced saturation leads to a more widespread but weaker seepage response.

3.4. Liquefaction Analysis

To examine the influence of tidal bore and soil parameters on the momentary liquefaction of the riverbed, the following five sets of parameters are considered:

• Tidal bore height: 0.6 m to 4 m;
• Water depth: 2 m to 5 m;
• Bed thickness: 0.6 m to 10 m;
• Degree of saturation: Gβ from 0 to 0.6;
• Permeability coefficient: k_x = k_z, from 10⁻⁶ m/s to 10⁻² m/s.

The other parameters remain constant, including Poisson’s ratio μ = 1/3, porosity n = 0.3, and shear modulus G = 10⁷ N/m². Figure 8 illustrates the boundary between the momentary liquefaction zone and the stable zone within the riverbed, as predicted by two different models: Equation (53) by Bear [60] and Equation (55) by Madsen [61].

According to Bear’s criterion, momentary liquefaction occurs when the vertical resultant force on the particles becomes zero, resulting in a shallow liquefaction zone with a depth of z_l = 0.1 d. In contrast, Madsen’s criterion suggests that liquefaction occurs when the horizontal pressure gradient exceeds the intergranular stress, leading to a deeper liquefaction zone extending to z_l = 0.5 d. The maximum liquefaction depth predicted based on horizontal pressure gradients is five times that predicted based on vertical pressure gradients. This difference in liquefaction depth highlights the varying influences of vertical and horizontal seepage forces, with Bear’s criterion predicting shallow liquefaction under rapidly changing vertical seepage forces, while Madsen’s criterion indicates deeper liquefaction under significant horizontal pressure gradients. For conservative reasons, the liquefaction criterion proposed by Madsen [61] is adopted in the subsequent parameter analysis.

Figure 9 illustrates the relationship between maximum liquefaction depth (z_l/d) and tidal bore height (H/d) for a partially saturated riverbed composed of isotropic silt soil and fine sand. The results indicate a distinction in liquefaction behavior between the two soil types. Silt soil exhibits a greater liquefaction depth across all tidal bore heights compared to fine sand, highlighting its higher susceptibility to liquefaction. This behavior can be attributed to the finer particle size and lower permeability of silt, which result in higher pore pressure gradient under loading conditions such as those induced by tidal bores. When the water depth is 2 m and the tidal bore height is less than 1.4 m, no liquefaction occurs in the fine sand. However, when the tidal bore height exceeds 1.4 m, the liquefaction depth increases sharply. This indicates that when the relative tidal bore height is low, it is insufficient to impact the fine sandy seabed.

Figure 10 illustrates the relationship between the maximum liquefaction depth (z_l/d) and the water depth (d/h) before the arrival of the tidal bore for a partially saturated riverbed composed of isotropic silt soil and fine sand, with a constant tidal bore height. The data reveal a pronounced decrease in liquefaction depth as the initial water depth increases relative to the tidal bore height for both soil types. Silt soil, due to its finer particle size and lower permeability, initially leads to a higher liquefaction depth. However, as the water depth increases, the relative impact of the tidal bore on the silt soil’s stability lessens, leading to reduced liquefaction. As the water depth increases, the liquefaction depth in fine sand decreases more rapidly than in silt, indicating that fine sand exhibits relatively better stability under deeper water conditions.

Figure 11 illustrates the relationship between the liquefaction depth z_l/d and the relative riverbed thickness h/d, comparing isotropic silt soil and fine sand. As h/d increases, the liquefaction depth for both soil types also increases, indicating that thicker riverbeds are more prone to deeper liquefaction zones. For silt riverbeds, when h/d is less than 1, the liquefaction depth increases significantly with the increase in riverbed thickness; when h/d exceeds 1, the impact of riverbed thickness on liquefaction depth diminishes, tending towards a stable depth. The pattern is similar for fine sand riverbeds, but the critical value is around h/d = 2.5. This suggests that there is a most unstable thickness in the riverbed, resulting from the interaction between incident tidal bore and reflected stress waves from the riverbed bottom. Therefore, in engineering site selection, these sensitive areas should be avoided, and structural foundations should be placed in relatively safer locations.

Figure 12 illustrates the relationship between the maximum liquefaction depth z_l/d and the soil stiffness parameter Gβ, comparing the liquefaction behavior of a riverbed composed of isotropic silt soil and fine sand. The results indicate that as Gβ increases, the liquefaction depth for both types of soil also increases, though the rate of increase differs. The liquefaction depth for silt soil rises significantly with increasing Gβ, suggesting that the liquefaction potential of silt soil is highly sensitive to changes in soil stiffness. In contrast, the liquefaction depth for fine sand increases more gradually, indicating that its liquefaction potential is less affected by variations in soil stiffness. As the saturation degree increases, the Gβ value decreases and the soil stiffness increases, leading to a reduction in liquefaction depth, meaning that the more saturated the riverbed, the more stable it becomes. However, when the saturation degree S_r→1, the soil stiffness Gβ→0 and the maximum liquefaction depth of fine sand riverbeds does not vary significantly with changes in saturation.

Figure 13 illustrates the relationship between the maximum liquefaction depth z_l/d and the soil permeability k_x, examining the liquefaction behavior of partially saturated isotropic riverbeds (k_z = k_x). The results show that the liquefaction depth initially increases with increasing permeability, reaches a peak, and then decreases, indicating that there is a specific permeability range where the liquefaction depth is maximized. Specifically, at lower permeability values typical of silt soil (k_z < 1.6 × 10⁻⁴ m/s), the liquefaction depth increases as permeability increases due to the limited drainage capacity, which leads to a higher pore water pressure gradient under tidal bore, thereby enhancing liquefaction potential. However, as it continues to increase to values typical of fine sand, the liquefaction depth decreases because higher permeability facilitates drainage, reducing the pore water pressure gradient and thus lowering the risk of liquefaction. When the permeability coefficient is greater than 7 × 10⁻³ m/s, no liquefaction occurs in the riverbed because the pore water is easily discharged.

4. Discussion

4.1. Liquefaction Potential Evaluation

The results of this study reveal the critical impact of tidal bores on the liquefaction potential of riverbeds, particularly in macro-tidal estuaries, where the sudden influx of water creates significant hydrodynamic forces that challenge the structural integrity of the sediment layers. By applying the generalized Biot theory and integral transform methods, this study not only captures the transient response of riverbed sediments under these dynamic pressures but also offers a quantitative framework for predicting the conditions under which liquefaction is likely to occur. The study demonstrates that momentary liquefaction is not merely a theoretical possibility but a realistic outcome influenced by key factors such as the permeability coefficient and the degree of saturation of the riverbed. This understanding bears similarity to the seabed liquefaction phenomenon induced by wave loading [29,52].

Fine sand facilitates pore water discharge, significantly reducing liquefaction risk, whereas silt soil increases sediment instability. This finding emphasizes the critical need to consider local soil characteristics when evaluating the vulnerability of riverbeds to tidal bore events. In riverbeds with low permeability, the retention of pore water can result in pore pressure gradients reaching critical levels, thereby increasing the risk of liquefaction. Conversely, in riverbeds with high permeability, the efficient dissipation of pore water pressure can significantly mitigate this risk, identifying a potential area for targeted engineering interventions aimed at enhancing riverbed stability.

Additionally, the study’s insights extend beyond the specific context of tidal bores, offering broader implications for the management of coastal and riverine environments subjected to extreme hydrodynamic events. The consistency of these findings with previous research [25,27] reinforces the critical role of soil properties in determining sediment stability under rapid hydrodynamic changes. At the same time, the study advances the field by providing a more detailed analytical approach that considers the specific interactions between tidal bores and riverbed sediments. This approach not only enhances the understanding of sediment behavior in extreme conditions but also lays the groundwork for future studies that might explore more complex scenarios, including the effects of sediment layering, heterogeneity, and multi-bore interactions. Such advancements will be crucial in developing more effective strategies for mitigating the impacts of tidal bores and preserving the integrity of vulnerable estuarine and riverine ecosystems.

4.2. Addressing Limitations

Although this study offers a comprehensive analysis, several limitations indicate the need for further research to fully understand the complex dynamics of tidal bore-induced liquefaction. One key limitation is the assumption of homogeneity and isotropy in riverbed properties. This oversimplifies real-world conditions, where riverbeds typically display significant stratification and variability in material properties, such as grain size, permeability, and sedimentary structure. As highlighted by Hazirbaba [62], Zuo and Baudet [63], and Porcino et al. [64], such heterogeneities can lead to diverse responses to dynamic forces. While this study focused on clean sand and silty soils, the role of fines content in sand–silt mixtures was not explicitly considered. Another limitation lies in the use of a constant shear modulus G. Liquefaction involves large deformations, and G typically varies with shear strain and the number of applied cycles. While this study employed a constant G, we have acknowledged that strain-dependent variations in G are crucial for accurately modeling soil behavior under cyclic loading. Future work will incorporate these variations to enhance the precision of liquefaction predictions, especially in layered sediments where mechanical behavior may differ significantly from homogeneous assumptions.

Additionally, this study does not account for the effects of turbulent flow and sediment transport, which are critical in natural tidal bore environments [7]. Turbulent flow can generate complex stress patterns, leading to localized high-pressure zones that may trigger unexpected liquefaction. Sediment transport processes, such as erosion and deposition, also alter riverbed morphology over time, affecting stability in subsequent tidal events. The exclusion of these processes limits the model’s applicability in actual scenarios. To address this, we have expanded our discussion on the assumptions made and emphasized the importance of validating the model through field studies and in situ experiments. These additions strengthen the conclusions and highlight the need for empirical validation to support theoretical findings.

Moreover, the study primarily focuses on the dynamic response to a single tidal bore event, overlooking the long-term effects of repeated events. In many estuarine environments, riverbeds are subjected to daily tidal bores, leading to cumulative effects that may weaken sediment structures over time. The current model does not fully capture the progressive weakening of the riverbed under sustained cyclic loading. Future research should explore these long-term effects to better predict the potential for failure in riverbed stability.

Finally, we have addressed the limitations of the analytical solution, particularly when coefficients change over time, and have introduced numerical methods, such as the Gauss elimination method and the Newmark time integration scheme. These approaches offer alternatives when analytical solutions are impractical, providing a pathway for future refinements of the model. A numerical analysis method that combines stochastic processes with the finite element method (FEM) can also be employed to assess the quantitative impact of the random distribution of soil material properties [65].

5. Conclusions

This study presents a comprehensive analysis of the dynamic response and liquefaction potential of riverbeds induced by tidal bores, with a focus on macro-tidal estuaries. The analysis highlights the importance of permeability and saturation in determining the stability of riverbed sediments under the dynamic pressures generated by tidal bores. It also provides valuable insights into the role of different soil types in influencing liquefaction risk in estuarine environments.

The key findings demonstrate that fine sand reduces liquefaction risk by facilitating pore water discharge, while silt soil increases the likelihood of sediment instability. Additionally, the degree of saturation was found to significantly affect the distribution of effective stresses and seepage characteristics within the riverbed. Partially saturated riverbeds exhibit higher pore pressure gradients, which raise the risk of liquefaction compared to fully saturated conditions.

The assumptions of homogeneity and isotropy in riverbed properties, along with the exclusion of turbulent flow, sediment transport, and repeated tidal bores, limit the applicability of the study in complex natural environments. Future research should develop models that account for the layered structure of riverbeds, turbulent flow, sediment transport, and long-term effects of tidal bores to improve liquefaction risk predictions. Field validation in estuarine environments with stratified riverbeds and frequent tidal bores is also crucial. Gathering high-resolution data on sediment properties, flow dynamics, and riverbed behavior will further refine these models and provide more accurate risk assessments.