Field Measurements and Modelling of Vessel-Generated Waves and Caused Bank Erosion—A Case Study at the Sabine–Neches Waterway, Texas, USA

Abstract

The Sabine–Neches Waterway (SNWW) is home to the largest commercial port of the United States military and of the refineries that produce 60% of the nation’s commercial jet fuel. The deposited sediments from bank erosion due to wake wash result in frequent dredging to keep the waterway operational. This study investigates vessel-generated waves and their impacts on bank erosion. Surface wave data at Golden Pass and the City of Port Arthur Park dock were measured using a 1 MHz Aquadopp Profiler. Bank properties such as soil strengths were measured and soil samples were collected. Acceptable predictive models for estimating the maximum wave heights caused by vessels sailing through the SNWW were developed and validated with recorded data. Vessel-generated waves are found to produce enough shear forces to mobilize bed sediments and cause bank erosion. The bed erosion rate increases with an increase in wave height or a decrease in water depth. Bank and bank toe erosion occurs at both monitoring locations. Bank stability and toe erosion model (BSTEM) results suggest that potential bank protection options are large woody debris and riprap at Port Arthur. However, other stronger stabilization methods are required at Golden Pass.

1. Introduction

Vessel wake affects the safe and efficient passage of vessels in a navigation channel [1,2,3,4,5,6,7]. Additionally, the vessel-generated waves can seriously erode unprotected banks of navigation channels, rivers, and harbors [8,9,10]. The wake-induced oscillations can also cause damage to moored vessels and pose a safety threat to small boats moving through the wake zone [11,12]. As illustrated in Figure 1, the wake generated by a moving vessel forms a complex wave system, which consists of different wave components that are constantly superimposed. The so-called Kelvin wake pattern introduced by Havelock [13] in the V-shaped wake zone includes sets of diverging waves that move obliquely outward from the vessel’s sailing direction and groups of transverse curved waves with crests across the vessel’s sailing direction (Figure 1). Strong wave interference can be observed along the cusp locus lines, where the transverse waves meet the diverging waves [3,14].

As pointed out by Roo, et al. [15], the primary wave system is the long-period drawdown phenomenon that reflects a significant water level depression along the hull of a moving vessel (Figure 1). The amplification of the drawdown effect has been noticed in restricted shallow water areas such as navigation channels, especially when large displacement ships are in motion [11,16]. The drawdown wave height is generally associated with the maximum vessel-generated wave height. As vessel size and speed increase, the vessel-generated drawdown waves grow significantly, which create a tremendous amount of erosive force on the banks of navigation channels [15]. As a result, bank erosion and sedimentation to the channel bottom caused by large vessel-generated drawdown waves become a major concern in maintaining a navigation channel.

The vessel-generated wave height is difficult to predict because it depends heavily on several parameters: the vessel speed, the vessel geometry, the operating draft of the vessel, the water depth, and the channel cross-sectional area [4,5,6,7]. Under the idealized condition of an unconfined flow domain with a constant water depth and applying the Fourier integral method, Havelock [13,17] presented the analytical solutions of the wave pattern generated by a moving disturbance and the induced wave resistance in deep water. Following a similar approach by formulating the wave elevations as an integral form of a specified moving pressure distribution, the Kelvin wave patterns for cases with a large Froude number and with arbitrary moving surface pressures were investigated analytically by Darmon et al. [18] and Miao and Liu [19], respectively. By solving the Boussinesq equations with the added effects of moving free surface pressure, Torsvik et al. [20] numerically investigated waves generated by moving disturbances in either a rectangular or a trapezoidal channel. Later, Chang and Wang [21] developed a numerical model to simulate the generation and propagation of three-dimensional (3D) fully nonlinear water waves caused by a submerged object moving at speeds varying from a subcritical to supercritical condition.

Numerical models may provide relatively accurate predictions on waves generated by moving vessels. However, the model setup for the inputs of vessel geometry, vessel operating conditions, and 3D grid systems to cover the irregular channel domain and uneven water depths can be tedious, and the numerical calculations are time-consuming and require considerable computational resources. On the other hand, as pointed out by Hartman and Styles [7], the low-order empirical approaches for predicting vessel wake may be appropriate for some applications. These include the use of screening to identify the potential large waves generated by a moving vessel, situations where a detailed wake analysis is not required, or when the resources for a quantitative prediction are not available. A number of simplified vessel-wake estimation models based on field and lab measurements have been developed in the past [3]. In a recent U.S. Army Corps of Engineers (USACE) study, Hartman and Styles [7] provided a table summarizing the empirical models developed over a period from 1965 to 2005, which included models from Balanin and Bykov; USACE; Bhowmik; Bhowmik, Demissie, and Guo; Blaauw, Groot, Knaap, and Pilarczyk (Blaauw et al.); Permanent International Association of Navigation Congresses (PIANC); Sorensen and Weggel; Bhowmik, Soong, Reichelt, and Seddik; and Maynord. Compared to the empirical models described by Sorensen [3], only Maynord’s model, established in 2005, was added to the updated list. This suggests that the development of new empirical models has been very limited. The USACE study [7] also developed a MATLAB based vessel wake prediction tool (VWPT) to estimate vessel-generated maximum wave height and to support navigation studies requiring knowledge of the potential effects of vessel wake, especially on bank erosion. The simplified algorithms in VWPT were based on the empirical models mentioned above as a function of vessel speed and other parameters such as vessel length, beam, and draft.

The erosion due to wake wash can be evaluated by resolving the wave-induced shear forces and the subsequent displacements of the bank, bank slope, and bed material. Soulsby [22] adopted linear wave theory to calculate the bed shear stress for bed erosion study. Limited reports have attempted to estimate the amount of erosion due to wake wash on riverbeds or stream banks [23,24]. The relationships between wave characteristics and the determined shear forces have been applied to study local bed and bank stability [25,26].

A streambank can be divided into toe, bank, and overbank zones. The toe zone is between the ordinary high water and low water levels, which is exposed to strong currents, debris movement and wet-dry cycles; the bank zone is above the toe zone and inundated during periods of moderate flows; the overbank zone is the inland portion and only becomes inundated during periods of high water [27]. Using the same physical process-based routines, the U.S. Department of Agriculture—Agriculture Research Service (USDA-ARS) have developed the bank stability and toe erosion model (BSTEM) [28,29,30] and the conservational channel evolution and pollutant transport system (CONCEPTS) [29,31,32] to simulate bank stability and the fluvial erosion process. The CONCEPTS model incorporates sediment transport processes from beds and banks to evaluate the long-term impact of rehabilitation measures for stabilizing stream systems and reducing sediment yield for a large-scale watershed [31,33]. The widely used BSTEM version 5.4 [34] is a comprehensive and effective tool for evaluating river or stream bank erosion and stability. Although there are limitations on incorporating subaerial processes, the effects of riparian vegetation, spatial, and temporal variability in geotechnical failure, and fluvial erodibility parameters are included [34]. The BSTEM simulates the failure by shearing a soil block with various geometries and the erosion on bank and bank toe material [28,35,36]. The erosion of SNWW due to wake wash is mainly located on the bank and bank toe, and the heavy erosion occurs on the short length of the reaches, so the CONCEPTS model is not applicable and the BSTEM 5.4 is sufficient to determine the erosion control approaches for the SNWW.

In this study, the tasks of measurements and model development for the predictions of the maximum wave height generated by large vessels sailing through the SNWW were performed. Estimations of the erosion rates and eroded bank profiles caused by the vessel-generated waves were also investigated. However, it should be noted that the effects of channel flow and tidal oscillation were not included in the model development and bank erosion analysis. The field study included the deployment of instruments at two selected locations to collect wave data, extract soil samples, and survey the bank conditions. The measured wave data were analyzed to determine the maximum wave heights at the time of vessels sailing through the wave data sampling sites. The predictive models used to establish the relationship between the maximum wave height and vessel characteristics were developed through the calibration and validation procedures. The eroded bank profiles caused by the vessel-generated maximum waves were calculated using the BSTEM model. Potential erosion control approaches to reducing the amount of erosion and consequent dredging were also analyzed and discussed.

2. Methods

2.1. Study Area

The Sabine–Neches Waterway (SNWW) is of strategic importance to the economy of not only southeast Texas but the entire nation; it is home to the largest commercial port of the United States military and refineries that produce 60% of the nation’s commercial jet fuel [37]. In 2017, the waterway was set to undergo dredging to increase the depth of the channel from 12.2 m (40 ft) to 14.6 m (48 ft) to accommodate larger cargo vessels [37]. The deepening project is currently underway by the U.S. Army Corps of Engineers (USACE). The industrial entities along the waterway have applied various erosion protection measures around their properties. This leaves the rest of the banks to experience more severe erosion. The banks of the SNWW are mostly subject to wave action caused by vessel passages. At the Sabine Pass station, the mean and diurnal tidal ranges are 0.41 m (1.33 ft) and 0.59 m (1.93 ft) [38], respectively. They decrease, respectively, to 0.26 m (0.84 ft) and 0.32 m (1.04 ft) at the Port Arthur station [38]. The effect of tides with the strong diurnal component is small when compared to the vessel-generated large local waves. The freshwater inflows are mostly from the Sabine and Neches Rivers, which contribute an average flow rate of 415.16 m³/s (14,650 ft³/s) [39]. On average, 125 million tons of cargo are transported each year along the waterway, with thousands of vessels passing through each month [40,41]. Significant bank erosion due to vessel-generated waves (also called “wake wash”) has been observed to cause loss of land and damage to existing structures along the SNWW. The eroded sediment also settles at the bottom of the channel, necessitating further dredging. To keep the waterway available for commercial waterborne traffic, the USACE Galveston District signed a $3.6 million dredging contract for maintenance dredging of the SNWW in 2015. This kind of dredging contract occurs approximately once every two years [42]. This recurring expenditure to counteract the bank erosion caused by wake wash makes the maintenance of the SNWW to be very costly.

2.2. Field Study

Two wave-monitoring locations (Figure 2, left panel) were chosen based on the waterway channel and bank conditions, instrument installation requirements, and vessel transportation characteristics. Location A was a concrete dock (Figure 2, right bottom panel) at the Golden Pass LNG terminal with heavy erosion upstream on both the left and right banks facing downstream. Location B was also a concrete dock (Figure 2, right middle panel) on the right bank of the SNWW by the City of Port Arthur Park, where heavy erosion was observed on the left bank while the right bank has hurricane walls to protect it from erosion. The instrument was mounted on the pier of the dock to avoid bank effects on the wave measurements. The channel width at Location A is much larger than at Location B.

The bank and channel geometries were surveyed at a site (29°46′29″ N, 93°56′28″ W) about 200 m upstream of Location A on 3 November 2017 and at a site (29°52′13″ N, 93°55′44″ W) on the opposite bank of Location B on 23 March 2016. The channel flow moved in the downstream direction at the time the field survey was conducted. The variation of discharge was small. The survey site for Location A (Golden Pass) was the closest upstream point to Location A experiencing severe erosion (Figure 3, left panel) where bank erosion control measures have been implemented at the Golden Pass LNG terminal property. Severe erosion was also observed at the site close to Location B, as shown in Figure 3 (right panel), due to sheet piles constructed upstream. Soil strengths were measured and soil samples were collected at three different layers of the bank, as marked in Figure 3.

2.3. Wave and Vessel Data Collection and Significant Wave Height Calculation

The wave data were collected using an Aquadopp Profiler, which is a small acoustic doppler current profiler (ADCP) primarily used in shallow coastal water and estuaries (https://www.nortekgroup.com/products/aquadopp-profiler-1-mhz, accessed on 1 January 2016). Aluminum rigid frames with lengths of 29.0 m and 7.3 m for housing the Aquadopp Profiler were constructed (Figure 2, right top panel) and applied successfully for wave measurements at Locations A and B, respectively. The ADCP was mounted at a depth of 15 m below the dock at Location A (Golden Pass, right bottom panel of Figure 2) and 3 m below the dock at Location B (Port Arthur, middle panel of Figure 2). The wave data for Location B were collected from 4 May 2016 14:21:01 to 11 May 2016 13:11:01 and from 12 May 2016 14:31:01 to 26 May 2016 10:01:01, located 107 m (350 ft) from the waterway centerline. Location A was 298 m (978 ft) from the waterway centerline and the wave data were collected from 11 January 2017 15:01:01 to 25 January 2017 13:41:01 and from 30 January 2017 15:52:01 to 15 February 2017 10:02:01. The raw data consisted of pressure signals sampled at 2 hertz that were grouped into 10-min blocks of data.

The QuickWave software developed by Nortek (https://www.nortekgroup.com/software, accessed on 1 January 2016) was used to convert the submerged pressure measurements into water-surface elevations using the linear wave theory. This localized linear approach of converting pressure data into wave data has been commonly adopted for the development of software for wave elevation sensors such as the Aquadopp Profiler used in the present study. This approach has the advantages of convenience and simplicity in the calculation of wave elevation and wave properties such as the mean wave period, the peak wave period, and the directional spread. Additionally, the integrated frequency domain analysis in QuickWave superimposes all wave components to calculate the wave energy spectrum. For very large waves, the use of nonlinear analysis [43] may improve the conversion of pressure measurements into wave elevations and individual wave characteristics. However, as pointed out by Sobey and Hughes [43], nonlinear analysis is computationally intensive and is probably not suited for routine analysis and storing data for instruments with a lengthy deployment. In this study, since the instruments were deployed in a confined shallow waterway and the linked software QuickWave was used, it was considered reasonable to neglect the nonlinear influence on the conversion of each individual wave measurement collected. Additionally, test runs were carried out to compare the reconstructed wave elevations using both the linear formula and Bonneton and Lannes’s [44] nonlinear nonlocal implicit formula, or equivalently Oliveras et al.’s [45] formula (See also in Marino et al., [46]). Both formulas produced nearly identical wave elevations and shapes.

The data were filtered by setting a frequency range from 0.001 Hz to 0.99 Hz to remove unrealistic data and low frequency tidal waves. Therefore, the tidal wave effect was not included in the wave statistics analysis. However, the long period drawdown waves were included in the data analysis and the maximum wave height modeling study. The significant wave height, H_s, was calculated using wave data collected at each 10-min interval as a representation of the wave height for waves occurring within that 10-min wave sample. Two methods were used to calculate the significant wave height. One was the direct wave-count method according to the definition of the significant wave height. After ranking the wave heights from the highest to the lowest values using all waves obtained from the 10-min wave data, the significant wave height was calculated as the average wave height of the highest one-third of the waves. The other method was the energy spectrum method based on the direct correlation between significant wave height and the sum of the energy spectrum curve, which is formulated as [47]

$$H_s=4 \sqrt{\underset{0}{\overset{\infty }{{\displaystyle \int }}}\mathrm{S}(f)df}$$

(1)

where $S\left(f\right)$, as a function of wave frequency f, represents the wave energy spectrum of the waves within a specified time interval, which is 10 min in this study. Once the variation of wave energy as a function of wave frequency for all waves is determined, a representative wave-energy spectrum can be defined. In general, its distribution can be approximated as the Rayleigh distribution. As a result, the significant wave height can be correlated with the wave energy spectrum as shown in Equation (1).

Vessel transit information was provided by the U.S. Coast Guard, which contained automatic identification system (AIS) data and the types of moving vessels. The AIS data gave the temporal and spatial information for each vessel at about a 1-min interval and the vessel type included the vessel geometry. The maritime mobile service identity (MMSI) number was used to uniquely identify each vessel passing through the SNWW for the vessel geometry information. The vessel data and AIS data were merged into one table based on each vessel’s MMSI number, including vessel’s MMSI, name, geometry (length and beam), draft, heading, course over ground, speed over ground, latitude and longitude, and timestamp of record. The data were imported into a PostgreSQL database. A shape file was exported based on the latitude and longitude (reported from AIS) to verify that the recordings were from the same locations as the wave monitoring locations. The vessel passage times were matched with corresponding wave measurement periods.

2.4. Empirical Models for Predicting Vessel-Generated Maximum Wave Height

Similar to the general vessel wake patterns described in the Introduction section, in confined and depth-limited waterways such as the SNWW, the vessel-generated waves include the short-period Kelvin wave pattern (secondary waves) and the long-period drawdown (primary waves). The Kelvin wave characteristics, in terms of wave celerity and wavelength, have been found to be heavily dependent on the vessel speed and water depth, or on the combined effect of vessel Froude number [13,48], $Fr$, which is defined as

$$Fr=\frac{V}{\sqrt{gh}}$$

(2)

where V is the vessel speed, h is the water depth, and $g$ is the gravitational acceleration. Therefore, the celerity of Kelvin wave wake systems, i.e., the divergent waves and the transverse waves, can be determined in a more straightforward way. It has been shown that the celerity (C) of the transverse waves in principle is equal to the vessel speed (V). However, the celerity of the diverging waves can be estimated from an empirical formula, $C=Vcos\left({35.27}^{°}\left(1-e^{12\left(Fr-1\right)}\right)\right)$, proposed by Weggel and Sorensen [48], where 0 <$Fr$ < 1. Based on the determined wave celerity, the wavelength L of either the transverse or diverging waves can be calculated according to the linear wave-based dispersion relationship:

$$C^2=\frac{gL}{2\pi }tanh\left(\frac{2\pi h}{L}\right)$$

(3)

and the wave period T can be determined as T = L/C.

However, for the estimation of vessel-generated maximum wave height, the primary wake system of the long-period drawdown waves was used to develop the predictive models, as the drawdown wave height generally represents the vessel-generated maximum wave height. An empirical approach using collected maximum wave height ($H_{max}$) data from field observations or laboratory measurements has been adopted with regression analysis to establish model equations. In Sorensen [3], various maximum wave height predictive models based on the wave conditions, vessel data, and channel geometries have been tested in a study of the upper Mississippi River system.

A series of test runs using the VWPT [7] were conducted with the field wave data measurements from Locations A and B in the SNWW. Using the model equations listed in Hartman and Styles’ [7] report, it was found that the USACE model for the upper Mississippi River performed better than other listed empirical models for determining maximum wave height. The value of the coefficient of determination ${\widehat{R}}^2$ from the USACE model was −0.427. However, the ${\widehat{R}}^2$ values for Blaauw et al., PIANC, and Sorensen and Weggel models were −1.116, −1.397, and −6.801, respectively. It is noted that the negative ${\widehat{R}}^2$ value was not the result of regression but rather an application of an empirical model that did not follow the trend of the data. The original USACE model, though producing a negative ${\widehat{R}}^2$ value with underestimated results, showed a better trend than the other three possible model selections. Therefore, the USACE model equation was adopted for the further development of the SNWW models through calibration and verification.

To extend the application of the USACE model to The Sabine–Neches Waterway, its empirical formulation is generalized as

$$H_{max}=\mathsf{\alpha } V^{\mathsf{\beta }}{\left(\frac{D}{L_v}\right)}^{\mathsf{\gamma }}{\left(\frac{S_c}{S_c-1}\right)}^{\mathsf{\delta }}$$

(4)

where $H_{max}$ is the maximum wave height, $\mathsf{\alpha }$, $\mathsf{\beta }$, $\mathsf{\gamma }$, and δ are calibrated coefficients and in meters, S_c is the channel section coefficient, which is the channel cross-sectional area divided by the wetted cross-sectional area of a vessel at midship, L_v is the vessel length, V is the vessel speed, and D is the vessel draft. The original USACE model coefficients are $\mathsf{\alpha }=0.0448$, $\mathsf{\beta }=2$, $\mathsf{\gamma }=0.5$, and $\mathsf{\delta }=2.5$. In the present study, a regression approach was followed by using the combined wave and vessel data collected at Location B to obtain calibrated coefficients that best predict the vessel-generated maximum wave height in the SNWW navigation system.

2.5. Stream Bed Erosion and Bank Stability

Erosion rates (E) at the stream bed due to wave induced bed shear stress are a function of the critical shear stress of the soil (τ_c), the induced bed shear stress (τ_b), and an erosion coefficient (k) as defined in Equation (5) [49,50,51]:

$$E=k\left({\tau }_b-{\tau }_c\right)$$

(5)

The critical shear stress is a function of plasticity, grain distribution, and percent clay and can be estimated with the shields diagram for initiation of particle motion [52,53]. The erosion coefficient (k) can be estimated as $k=0.2{\tau }_c{}^{-0.5}$ [54], or $k=1.62{\tau }_c{}^{-0.838}$ for cohesive materials [30]. The wave induced bed shear stress in Equation (6) can be calculated with water density (ρ), the rough bed wave friction factor (f_wr) and the peak value of orbital velocity (U_orb$)$ [22]:

$${\tau }_b=\frac{1}{2}\rho f_{wr}{U}_{orb}^2$$

(6)

where $U_{orb}=\frac{\pi H}{T·sin\mathrm{h}\left(\frac{2\pi h}{L}\right)}$, which is calculated using H = significant wave height, T = wave period, h = water depth, and L = wavelength [55]. The wave friction factor $(f_{wr})$ is estimated as $f_{wr}=1.39{\left(\frac{A}{z_o}\right)}^{-0.52}$, where A = length scale (m), calculated as $A=\frac{U_{orb}T}{2\pi }=\frac{2H}{sinh\left(\frac{2\pi h}{L}\right)}$, and $z_o$= roughness length (m). The $z_o$ values of 0.2 mm for cohesive sediment, 0.7 mm for sand and clay, and 0.4 mm for sand have been estimated in Soulsby [22].

Stream bank failures in fluvial systems includes hydraulic failures, geotechnical failure, or a combination of both [56]. Hydraulic failure occurs when flowing water exerts a tractive force that exceeds the critical shear stress of the bank material, while geotechnical failure is the result of bank moisture when the driving force (stress) exceeds the resisting force (strength). However, a combination of both failures is the most common due to surface erosion from the additional cantilever failures of undercut bank [57]. Slope failure modes can be generally categorized as fall, topple, slide, spread, and flow [58]. Failure modes for the two locations on this study fall more closely into the combinations of “slide/rotation failure” and “flow” modes due to flow-induced bottom shear stress.

To model bank stability, BSTEM calculates a factor of safety (FS) using horizontal layers, vertical slices, and cantilever shear failure models [28,59]. It accounts for up to five user-input horizontal soil layers and determines the shear stress of saturated soil with unique geotechnical properties and driving forces controlled by bank height and slope [28,59]. Along vertical slices, the model examines the flow-induced normal and shear forces active in slices of the failure blocks. The factor of safety is calculated as the ratio between the resisting forces and the driving forces along a potential failure plane. The resisting forces determined by shear strength of soil $({\tau }_f)$ are as follows [28,33]:

$${\tau }_f=c^{\prime }+\left(\sigma -{\mu }_w\right)tan{\varphi }^{\prime }\left(\mathrm{saturated} \mathrm{soil}\right)$$

(7)

$${\tau }_f=c^{\prime }+\sigma tan{\varphi }^{\prime }+\psi tan{\varphi }^b\left(\mathrm{unsaturated} \mathrm{soil}\right)$$

(8)

where c’ is the effective cohesion, σ is the normal stress, μ_w is the pore-water pressure, ϕ’ is the effective internal angle of friction in degrees, ψ is the matric suction, and ϕ^b is an angle that describes the relationship between shear strength and matric suction. Driving forces for streambank instability are controlled by bank height and slope, the unit weight of the soil and the mass of water within it, and the surcharge imposed by any objects on the bank top [60,61]. An FS value greater than one indicates stability, while an FS value less than one indicates instability. Various combinations of failure plane angle and shear emergence elevation can be estimated in BSTEM 5.4 [34] using an iterative procedure to automatically determine when the FS approaches unity. The shear emergence elevation is defined as the elevation on the bank face where shear surface emergence is located [34]. The BSTEM computes an average erosion rate for each node by utilizing an excess-shear stress approach [49], which is then integrated with respect to time to yield an average erosion distance, E_L [60]:

$$E_L=k\Delta t({\tau }_0-{\tau }_c)$$

(9)

where k is the erodibility coefficient, estimated as $2\times {10}^{-7}{\tau }_c^{-0.5}$ for silt-clays [54], ∆t is time step, τ_c is the critical shear stress, and ${\tau }_0$ is an average shear stress induced by fluid flow. To correct the effect of channel curvature, the BSTEM accounts for secondary currents that can cause the bottom of the bank to possibly experience shear stress many times that of the first order approximation [62]. It separately calculates shear stress induced by the stream flow on the bank using “no-lag kinematic model” to account for secondary currents as [62]

$${\tau }_0=\frac{{\gamma }_w{n}^2{(u+U)}^2}{R^{1/3}}$$

(10)

where n is Manning’s roughness coefficient, ${\gamma }_w$ is the specific weight of water, R is the hydraulic radius, u is the reach-averaged water velocity, and U is the increase in the near-bank velocity due to superelevation caused by channel curvature [62]. One of the main limitations of BSTEM is that it can only assign the same soil properties for the same horizontal soil layer without considering its moisture, dispersion, and flow conditions [33]. Therefore, the model was developed using the soil type of the different bank layers identified during the field survey to account for the soil’s non-homogeneity and with the consideration of soil to be saturated as the worst-case scenario.

The model was developed to investigate potential erosion caused by waterway flow with the supplement of the analysis of wave and soil data. The BSTEM evaluates the bank stability using Equations (7) and (8) with the inputs of bank and channel profiles, bank and bank-toe material data, water depth, and water table depth. The model also estimates the erosion rates using Equation (9) based on the shear stress calculated from Equation (10) by accounting for the calculations of u and U with the inputs of channel curvature, reach length, Manning’s roughness coefficient, and bed slope. Additionally, it predicts the eroded area at bank, bank toe, and channel bed at different times. To account for the wake wash, the input water depth was determined by adding the maximum wave height and the measured field water depth. To determine possible erosion protection approaches, the critical shear stress caused by the flow and wake wash maximum wave on the bank material can be determined and compared with the protection options in the BSTEM. For example, if the riprap protection approach is chosen, the critical shear stress of the bank would be estimated as the value of 204 N/m² obtained from experimental studies [63].

3. Results

3.1. Field Survey and Soil Characteristics

The bank profile and the water elevation below the bank was measured at both locations to provide the input data for the BSTEM base model. The water elevation measured in the field was 2.23 m below the bank at Location A on 3 November 2017 and 1.86 m below the bank at Location B on March 23 2016. The soil samples for Location A were taken at the water’s edge (Location A.3), from the road below the shoulder of Highway 87 (Location A.1), and approximately 3 m from Location A.1 at Location A.2 as indicated in Figure 3. The soil samples for Location B were obtained from dry land (Location B.1), at the shoreline (Location B.2), and at the water line (Location B.3) as shown in Figure 3, respectively. The lab analysis on specific gravity, medium size, liquid limit and plastic limit, porosity, hydraulic conductivity, and soil type based on AASHTO and Unified Soil Classification Systems (USCS) are summarized in

Table 1. Analyzed soil characteristics from soil samples collected at study sites.

Location	Specific Gravity	Medium Size (D50)	Liquid Limit	Plastic Limit	Porosity	Hydraulic Conductivity	Soil Type (AASHTO/USCS)
Location A.1	2.60	0.19 mm	21.0	18.7	0.31	0.0006 cm/s	A-2-4/poorly graded sand-silt mixture
Location A.2	2.54	2.5 mm	23	17.9	0.39	1.4 cm/s	A-1-a/poorly graded sand with gravel, silty and clay mixture
Location A.3	2.43	0.7 mm	25	20.7	0.41	0.004 cm/s	A-1-b/silty sand with gravel
Location B.1	2.70	0.018 mm	63	19	0.51	0.000003 cm/s	A-7-6/Fat Clay
Location B.2	2.58	0.97 mm	28.5	17.3	0.38	0.005 cm/s	A-2-6/clayey sand
Location B.3	2.63	0.42 mm	28	20	0.42	0.003 cm/s	A-1-b/silt sand

. The unusually small D₅₀ and hydraulic conductivity found at Location B.1 in Table 1 is because the soil was backfilled from dredge waste. The bank experiences severe erosion and soil shear strength plays an important role in the stability analysis of the embankment. Field tests including the pocket penetrometer and vane shear tests were performed. The shear strengths were approximately 239.4, 143.6, and 47.9 kN/m² at Locations A.1, A.2, A.3, respectively, and 33.5, 143.6, and 47.8 kN/m² at Locations B.1, B.2, and B.3, respectively. Since the soil strength strongly relates to the moisture content, lab tests including Standard Proctor tests with corresponding unconfined compression tests were also performed. These tests investigated the possible range of soil shear strength under different moisture contents and dry densities for soil samples retrieved from each location. As expected, the strength values obtained from the lab and field tests vary greatly, ranging from 240 kN/m² to near zero strength depending on the moisture content and dispersion conditions of the soils. The shear strengths used in the analyses were estimated based on the best engineering judgement available through past experience related to the embankment stability analyses of SNWW.

3.2. Wave Data Analysis and Significant Wave Height

In this study, all vessel-generated waves, including short-period Kelvin waves (secondary waves, Figure 1) and long-period drawdown (primary waves, Figure 1), were collected and analyzed. A typical example showing the time variations of the significant wave height, which represents a wave energy related average wave height for all waves (including short- and long-period waves) collected at a defined time interval, is presented in this section.

After applying the filter to the pressure data (See Section 2.3), the data set that includes only the waves with frequencies from 0.001 Hz to 0.99 Hz was used for the data analysis. An example plot showing the filtered wave related pressure measurements is presented in Figure 4. An example plot showing the time variations of significant wave height, calculated using the direct-wave count method from the recorded data at Locations A and B in the Sabine–Neches Waterway, is presented in Figure 5. As shown in Figure 5, the maximum significant wave height is about 0.34 m at Location A and about 0.46 m at Location B. The Froude number as defined in Equation (2) for Locations A and B were in the subcritical category (Fr < 0.7), indicating that the moving vessels generated the mixed Kelvin wake pattern and drawdown [15,16,64].

The measurements from the SNWW suggest that the wave period of the drawdown associated maximum wave-height waves range from 85 s to 155 s. For subcritical vessel speeds in a confined channel with shallow water depth, as in the SNWW, the drawdown waves of greatest height and period can cause the most bank erosion [65]. According to the recorded waves generated by vessels, the estimated surge waves had an average maximum wave height of 0.3 m at Location A and 0.33 m at Location B. In addition, the vessel-generated maximum wave height is about 0.45 m at Location A and about 0.55 m at Location B. As shown in Figure 6, the surge wave transitioned into a long-period wave after a vessel left. The green dots indicate the time history of a moving vessel and the free-surface elevation is referenced to the mean water surface level.

3.3. Predictive Models for Long-Period Drawdown Induced Maximum Wave Height

In this study, predictive models for estimating the vessel-generated maximum wave height in the SNWW were developed by calibrating the proposed empirical formula of Equation (4) and validating the calibrated formula. The maximum wave height $\left(H_{max}\right)$ recorded for each passing vessel are summarized in

Table 2. Measured maximum wave heights and their corresponding passing vessel information.

Vessel Name	Vessel Type	Speed (knots)	Draft (m)	Length (m)	Beam (m)	Hmax (m)
SILVIO	LNG Tanker	6.5	8.80	220	36.6	0.202
PARAMOUNT HATTERAS	Tanker	7.7	8.50	250	44.0	0.204
NS LAGUNA	Tanker	7	9.00	250	43.8	0.207
LUBARA	LNG Tanker	6.7	12.05	255	36.6	0.212
LONE STAR STATE	Tanker	6.7	11.43	186	32.2	0.217
GULF MIST	Tanker	7.8	10.65	184	27.0	0.218
GUARDIANSHIP	Cargo	6	12.19	190	32.3	0.221
ALTHEA	Tanker	7.1	10.10	248	43.0	0.223
ORIENT TRADER	Cargo	7.7	10.10	180	30.0	0.223
LEFKARA	Tanker	5.3	11.40	183	32.2	0.230
LONE STAR STATE	Tanker	7.5	11.43	186	32.2	0.234
EAGLE TORRANCE	Tanker	8.4	8.50	236	42.1	0.236
GENER8 ORION	Tanker	6.9	10.60	274	48.0	0.238
NS ARCTIC	Tanker	6	10.50	250	44.0	0.241
EMERALD SPIRIT	Tanker	6.9	11.70	243	42.0	0.247
NORDBAY	Tanker	6.3	8.79	249	44.0	0.249
HOUSTON	Tanker	8	11.13	187	27.4	0.250
MARAN SAGITTA	Tanker	6.4	8.60	244	42.0	0.251
MARAN ATLAS	Tanker	6.5	8.61	244	42.0	0.252
SILVIO	LNG Tanker	7.1	12.00	220	36.6	0.256
EAGLE SIBU	Tanker	6.7	8.80	244	42.0	0.260
GULF MIST	Tanker	6.4	10.65	184	27.0	0.268
BOW CARDINAL	Tanker	7.2	8.10	183	32.2	0.277
MARAN ATLAS	Tanker	7.4	8.61	244	42.0	0.279
ALPINE ETERNITY	Tanker	6.5	10.80	183	32.2	0.282
TELLEVIKEN	Tanker	6.9	8.80	250	43.8	0.289
HELLAS SPARTA	LNG Tanker	7.1	12.00	225	36.0	0.289
BRITISH ROBIN	Tanker	8.4	11.10	252	43.8	0.305
BOW TONE	Tanker	8.1	9.50	170	26.6	0.306
CORRIDO	Tanker	6.9	10.90	183	32.0	0.309
ENERGY PATRIOT	Tanker	6.8	11.10	183	32.2	0.315
SCF SAMOTLOR	Tanker	6.9	9.90	274	48.0	0.319
EAGLE SIBU	Tanker	6.2	12.19	244	42.0	0.339
PELICAN STATE	Tanker	8.4	9.12	183	32.2	0.346
TEXAS	Tanker	7.7	10.80	183	32.2	0.352
PAZIFIK	LNG Tanker	6.2	11.10	205	32.0	0.357
NS ARCTIC	Tanker	7.8	10.50	250	44.0	0.383
OVERSEAS ANACORTES	Tanker	8.4	11.41	183	32.2	0.395
STAR KVARVEN	Bulk Carrier	9	9.38	209	32.3	0.396
AURORA VAR	LNG Tanker	8.2	11.60	225	36.6	0.397
EMERALD SPIRIT	Tanker	8.5	11.70	243	42.0	0.404
EAGLE KANGAR	Tanker	7.2	9.00	244	42.0	0.428
BBC KIMBERLEY	Cargo	9.4	5.45	139	20.0	0.431
HOUSTON	Tanker	6.7	11.13	187	27.4	0.440
BW LEO	LNG Tanker	9.1	11.70	225	36.0	0.445
NS CHAMPION	Tanker	8.8	10.20	244	42.0	0.449
EAGLE SIBU	Tanker	8.8	8.80	244	42.0	0.460
AGATHONISSOS	Tanker	8.8	8.30	244	42.0	0.465
ELKA BENE	Tanker	8.4	11.50	189	32.2	0.473
EAGLE FORD	Tanker	8.4	11.38	270	42.0	0.474
ORIENT TRADER	Cargo	8.5	10.10	180	30.0	0.527
SEABULK CHALLENGE	Tanker	9.2	11.51	187	32.3	0.529
OMODOS	Tanker	9.3	7.50	183	32.2	0.552

. It was observed that the maximum wave height was generally proportional to the vessel speed and vessel draft, while inversely proportional to the vessel length. The most common vessel types are tanker, LNG tanker, and cargo. With the calculated channel cross-sectional area, the parameter S_c can be determined with reference to a given vessel. The vessel speed and other relevant parameters can be found in Table 2.

Part of the data listed in Table 2 were used for calibration and the remaining data were used for model validation. The coefficients of the wave height predicting model (Equation (4)) $\mathsf{\alpha }$, $\mathsf{\beta }$, $\mathsf{\gamma }$, and δ were calibrated in four scenarios using the nonlinear least squares method from MATLAB. The first scenario (SNWW-1 model) assumed that the parameters $\mathsf{\beta }$, $\mathsf{\gamma }$, and δ are the same as those proposed in USACE model, i.e., $\mathsf{\beta }$= 2, $\mathsf{\gamma }$ = 0.5, and δ = 2.5. The calibrated coefficient $\mathsf{\alpha }$ was 0.063, which is slightly larger than the value used in the original USACE model of 0.0448. The second scenario is a two-parameter model (SNWW-2 model), which fixes $\mathsf{\gamma }$ = 0.5 and δ = 2.5 and allows $\mathsf{\alpha }$ and $\mathsf{\beta }$ to be adjusted based on the measured data. The calibrated coefficients are $\mathsf{\alpha }$ = 0.055 compared to USACE model’s value of 0.0448 and $\mathsf{\beta }$ = 2.11 compared to USACE model’s value of 2.0. The third scenario is a 3-parameter model (SNWW-3 model), where $\mathsf{\alpha }$, $\mathsf{\beta }$, and $\mathsf{\gamma }$ are fully calibrated while δ is set as 2.5. The obtained best-fit values are $\mathsf{\alpha }$= 0.0076, $\mathsf{\beta }$ = 2.19, and $\mathsf{\gamma }$ = −0.1. The fourth scenario is a 2-parameter model (SNWW-4 model), where $\mathsf{\beta }$ is set as 2.0 (with the assumption that the maximum wave height is proportional to velocity squared) and $\mathsf{\gamma }$ is set as 0.5. The unknown parameters are $\mathsf{\alpha }$ and δ, which were calibrated as 0.076 and 1.67, respectively. To measure the goodness of fit between the model results and calibration data, the coefficient of determination ${\widehat{R}}^2$ was calculated for each scenario. The values of ${\widehat{R}}^2$ for SNWW-1, SMWW-2, SNWW-3, and SNWW-4 models are 0.664, 0.726, 0.763, and 0.769, respectively. With the obtained ${\widehat{R}}^2$ values, the correlation between measured data and predicted values from the four models developed is visible. Although the SNWW-4 model was the best model in wave height prediction based on the ${\widehat{R}}^2$ value, all four models were considered acceptable for predicting the vessel-generated maximum wave height in the SNWW.

Table 3. Calibrated model coefficients for the proposed maximum wave height predictive models.

Maximum Wave Height Predictive Models	Model Coefficients				${\widehat{\mathit{R}}}^2$
	α	β	γ	δ
SNWW-1 model (α is a calibrated value.)	0.063	2	0.5	2.5	0.664
SNWW-2 model (α and β are calibrated values.)	0.055	2.11	0.5	2.5	0.726
SNWW-3 model (α, β, and γ are calibrated values.)	0.0076	2.19	−0.1	2.5	0.763
SNWW-4 model (α and δ are calibrated values.)	0.076	2	0.5	1.67	0.769

summarizes the calibrated model coefficients and the corresponding ${\widehat{R}}^2$ values for the proposed maximum wave height predictive models for the SNWW.

Figure 7 presents the comparison between the original USACE model predictions and the measured maximum wave height data (including the data for calibration and validation). A perfect-fit 45° line is also plotted in Figure 7 to reflect the level of agreement between predicted values and measured ones. The USACE model underestimated nearly all of the maximum wave height measurements recorded at the study area in the SNWW. The results from the calibrated models (SNWW-1, SNWW-2, SNWW-3, and SNWW-4) proposed in this study were also plotted to compare them to the measured maximum wave heights. The comparison plots using the results from SNWW-1, SNWW-2, SNWW-3, and SNWW-4 models are presented, respectively, in Figure 8, Figure 9, Figure 10 and Figure 11. The data used for model calibration and validation were all included in the comparison plots. The predicted maximum wave heights from the four proposed models generally follow the ideally fitted 45° line nicely, with the results from the SNWW-4 model fitting the best.

Typically, two statistical values are calculated to evaluate the performance of predictive models in validating stage. One is the root mean squared error (RMSE) and the other is the Nash–Sutcliffe model efficiency coefficient (NSE). In this study, the NSE value was calculated as a measure of the performance of the developed predictive models when compared to the validation data of maximum wave height. By using the event data, the NSE formula is given as

$$\mathrm{NSE}=1-\frac{{{\displaystyle \sum }}_{i=1}^N{H_{p_i}-H_{m_i}}^2}{{{\displaystyle \sum }}_{i=1}^N{\left(H_{m_i}-H_{m-ave}\right)}^2}$$

(11)

where N is the total number of data points, $H_{m-ave}={{\displaystyle \sum }}_{i=1}^N{H}_{m_i}/N$, and H_p and H_m are predicted and measured wave heights, respectively. When NSE = 1, the model is suggested to have a perfect prediction.

To evaluate the accuracy of the model performance, the RMSE and NSE values for each of the four models developed in this study were calculated by comparing the values calculated from each model with the data used for model validation. For reference, the RMSE and NSE values of the original USACE model were also calculated. The computed RMSE and NSE values for the USACE model, SNWW-1 model, SNWW-2 model, SNWW-3 model, and SNWW-4 model are summarized in

Table 4. Comparisons of the values of the root-mean-squared error (RMSE) and Nash–Sutcliffe model efficiency coefficient (NSE) obtained from the USACE model, SNWW-1 model, SNWW-2 model, SNWW-3 model, and SNWW-4 model.

Model	Nash–Sutcliffe Model Efficiency Coefficient (NSE)	RMSE (m)
USACE Model	−0.574	0.127
SNWW-1 Model	0.481	0.073
SNWW-2 Model	0.505	0.072
SNWW-3 Model	0.483	0.073
SNWW-4 Model	0.535	0.069

. The original USACE model, with a negative NSE value, completely fails to predict the maximum wave heights occurring in the SNWW. All of the calibrated models show a better performance than the original USACE model, with the NSE values for the proposed SNWW models (SNWW-1 to SNWW-4) varying from 0.481 to 0.535. The SNWW-4 model has the highest NSE value of 0.535. Additionally, the SNWW-4 model has the smallest RMSE. This demonstrates that the SNWW-4 model has the best performance in estimating the vessel-generated maximum wave height in The Sabine–Neches Waterway. The effect of the S_c value on the predicted maximum wave height was also examined. Using the SNWW-4 model, the predicted maximum wave height versus velocity squared (V²) for various S_c values is presented in Figure 12. As expected, the predicted maximum wave height increases with an increase in V². However, with an increase in S_c value, the predicted maximum wave height decreases. This suggests that when the channel cross-sectional area increases, the vessel-generated maximum wave height decreases. Alternatively, vessels with a larger wetted cross-section generate a greater maximum wave height. It is interesting to note from Table 4 that the SNWW-3 model was not the best model when compared to the validation data. This 3-parameter model has a calibrated negative value of γ, which reflected a positively proportional trend of the maximum wave height versus the vessel length rather than a negatively proportional tendency as indicated mostly by the data. This suggested that this 3-parameter model was over-fitting the calibration data. The 4-parameter model was not included in the list of predictive models since it was not as good as a 2-parameter model and some of the determined parameter values, because of over-fitting data, did not follow the suggested physical trends.

3.4. Channel Bed Erosion due to Waves

For the primary sandy soil in the field, the bed-roughness length z_o was determined to be 0.4 mm and the critical shear stress was estimated as 0.30 N/m² [22]. Therefore, the erosional coefficient of the soil can be found to be $k=0.2\times {0.3}^{-0.5}=0.365 {\mathrm{cm}}^3/\mathrm{N}.\mathrm{s}$ [54]. With the wave period and wavelength calculated based on the selected wave characteristics, the bed shear stress can be estimated using Equation (6) for different water depths near the bank. Under the short-period wave condition, the total number of waves in a wave package is generally considered to be around 10–15 waves [23] and the wave height is typically small comparing to the drawdown surge waves as observed in SNWW. We considered the duration of wave action is about 50% of fourteen waves, i.e., the first seven waves. The remaining seven waves are expected to have even smaller wave height. Considering the worst scenario in the calculation, we assumed that seven waves with the largest wave height from a short-period wave package were encountering the bank. For the drawdown surge waves, the duration was assumed to be two times the wave period. Due to the effects of larger wave height and action duration, the surge waves have a tendency to cause much greater bed erosion than an average short-period wave. For example, at 1 m of water depth, a surge wave at Location A is predicted to mobilize 195–480 g/m² of sediment, while an average short-period wave is predicted to mobilize 42 g/m² of sediment. The predicted erosion rate in g/m² versus water depth for the four selected surge waves at Locations A and B are presented in Figure 13. The erosion rate on the channel bed decreases with an increase in water depth but increases with an increase in wave height. It is shown in the left subplot of Figure 13 that the wave period has little effect on the erosion rate when similar wave height was used in the calculation. In those cases, the long shallow-water waves propagate with a similar phase speed of $\sqrt{gh}$. Since the drawdown surge wave plays an important role in mobilizing the stream-bed sediment, it becomes necessary to take into account its effect on sediment erosion when evaluating bank protection measures.

3.5. Bank Erosion due to Waterway Flow

The reach widths were measured as 440 m at the survey site of Location A and 274 m at Location B. A slope of 0.0004 m/m was selected as the bed slope of the Lower Sabine River [66]. Due to the sediment accumulation in the waterway, one must note that the slope of the waterway is likely less than that of the Lower Sabine River. The waterway radius of curvature was estimated from an aerial map to be 365 m for the 30 m reach length at both locations, and the combined Manning’s roughness coefficient for the bank material data was estimated to be 0.025 [64]. In this analysis, the bank level is treated as the zero elevation. The water-surface elevations were estimated with the field high water marks as 1.37 m below the bank at Location A and 1.31 m below the bank at Location B. The shear surface angle was determined to be 39.37° using the BSTEM based on the bank height and angle, bank toe length, and bank toe angle. Field survey on the bank geometry was used to generate the initial bank profile. The measured maximum wave height of 0.552 m in Table 2 was considered the worst-case scenario in the evaluation of the bank and toe erosion using the BSTEM. Therefore, the maximum water-surface elevation was calculated as a sum of the measured water-surface elevation and the maximum wave height, resulting in an elevation of −1.68 m (−2.23 + 0.552= −1.68) at Location A and an elevation of −1.31 m (−1.86 + 0.552 = −1.31) at Location B. The bank material was determined based on the soil characteristics data shown in Table 1. At Location A, the bank material at layer 1 was classified as erodible silt because its D₅₀ was 0.19 mm. Layers 2 and 3 were classified as fine rounded sand due to their D₅₀ being 2.5 mm and 0.7 mm, respectively. At Location B, the bank material at layer 1 was classified as erodible soft clay because its D₅₀ was 0.018 mm. Layers 2 and 3 were classified as erodible silt sand as their D₅₀ was 0.97 mm and 0.45 mm, respectively. Based on the field survey observation, the soil types erodible clay and resistant stiff clay were chosen as layers 4 and 5, respectively, for both locations. Since an unconfined aquifer is in the study area, the water table elevation was set to be the same as the water elevation in BSTEM.

The total erosion area changing with time was determined by adding the amount of bank and bank toe erosion in meters squared (m²) estimated by BSTEM under different time durations. As indicated in Figure 14, the erosion at Location A is more severe than at Location B. The erosion area increases at a fast pace during the first 100 days at Location A before gradually approaching a constant area. The total erosion area reaches 189 m² in 10 days, 268 m² in 50 days, and 364 m² in one year. At Location B, the variation of bank erosion also shows a rapid increase in erosion area in the first 50 days before transitioning to a nearly constant area. The total erosion area reaches 30 m² in 10 days, 66 m² in 50 days, and 110 m² in one year. The model outputs shown in Figure 15 indicate the initial profile (black line) and the eroded profile (red line) at 10, 50, and 365 days. The eroded profile of Location A presented in Figure 15 (upper panel) suggests that erosion occurs very quickly above layer 4 and can reach a near vertical profile in 10 days. A stable profile is observed above layer 4 after 50 days and 365 days. After 365 days, the erosion has continued until a near vertical profile is formed due to the resistance of the stiff clay in layer 5. At Location B, as illustrated in Figure 15 (lower panel), the bank profile does not show much change in the first 10 days. The significant erosion was observed in layer 4 in 50 days and the bank reaches a stable profile in layer 5 after 365 days.

In addition, when the water-surface elevation decreases from −1.68 m to −3.68 m at Location A, the total erosion area decreases from 189 m² to 77 m². Similarly, when the water-surface elevation decreases from −1.31 m to −2.20 m at Location B, the total erosion area decreases from 30 m² to 23 m². This suggests that the total erosion area increases as the water-surface elevation increases.

3.6. Bank and Toe Protection Measures

The average boundary shear stresses at the toe were estimated by BSTEM to be 90 N/m² at Location A and 32 N/m² at Location B. The permissible shear stresses of geotextiles, large woody debris, and riprap (D₅₀ = 0.256 m) are 144 N/m², 192 N/m² and 204 N/m², respectively. These erosion protection measures were chosen because their permissible shear stresses increase the bank and bank toe critical stress to be greater than the average boundary shear stress of 90 N/m² at Location A. According to Equation (9), the erosion rate becomes negative, or effectively no erosion. These erosion protection measures, applied either on the bank only, on the bank toe only, or on both the bank and bank toe, were evaluated by calculating the erosion area for 10 days at both locations. The erosion area of Locations A and B with each of the protection measures is summarized in

Table 5. Erosion area (m²) with the use of erosion protections of either geotextiles, large woody debris, or riprap.

Erosion Protections Applied On	Erosion Location	Erosion Area (m2)
		Location A			Location B
		Geotextile	Large Woody Debris	Riprap	Geotextile	Large Woody Debris	Riprap
Bank only	Bank	0	0	0	0	0	0
	Toe	87.651	87.779	84.114	15.866	12.859	12.846
Toe only	Bank	0.231	0.231	0.231	0.525	0.525	0.525
	Toe	5.396	2.338	1.657	1.203	0.921	0.092
Bank and Toe	Bank	0	0	0	0	0	0
	Toe	4.936	1.879	1.198	0.127	0	0

. It suggests that none of the protection measures are effective at Location A as erosion can occur at the bank toe no matter how the protection measures are applied. However, at Location B, no erosion is observed when either the large woody debris or the riprap (D₅₀ = 0.256 m) is applied on both the bank and bank toe. The results also indicate that the shear stress at some nodes of the bank toe exceed the critical shear stress of the protection measures.

4. Discussion

4.1. Predictive Models for Estimating Vessel-Generated Maximum Wave Height

As seen in Section 3.3, empirical models for predicting the potential vessel-generated maximum wave height occurring in the SNWW were developed for this study. The adopted model equation is based on the USACE model, as it generally outperformed the other empirical models described in Section 2.4 after a series of test runs were conducted and compared to measured wave data. However, the original USACE model without calibration and validation was found to fail in predicting the maximum wave heights occurring in the SNWW. Four models, named SNWW-1, SNWW-2, SNWW-3, and SNWW-4, were developed through a regression process by calibrating the model coefficients with part of the measured wave data. The RMSE and NSE value for each model was calculated based on the remaining measured wave data, and the SNWW-4 model was identified as the best model for estimating the vessel-generated maximum wave height. The proposed SNWW-4 model equation is given as

$$H_{max}=0.076 V^2{\left(\frac{D}{L_v}\right)}^{0.5}{\left(\frac{S_c}{S_c-1}\right)}^{1.67}$$

(12)

It was noticed that the vessel-generated maximum wave height is proportional to the term of V², indicating the maximum wave height is strongly influenced by kinetic energy from the vessel. The maximum wave height is also affected by the vessel length, the vessel draft, the channel cross-sectional area, and the wetted cross-sectional area of a vessel at midship. The developed model for calculating vessel-generated maximum wave height, as shown in Equation (9), can be a practical prediction tool to provide routine screening and, by identifying potential large waves that could be generated by a moving vessel in the SNWW, it has certain safety applications. The model’s predictions on maximum wave height can also assist in modeling and analysis related to the erosion of the bank and channel bottom, which is important to the design selection of bank protection measures and decision-making regarding the dredging of The Sabine–Neches Waterway.

4.2. Wake Wash Erosion and Protection

The results from the field survey, field wave measurements, and model study demonstrate that this part of the SNWW waterway bank is unstable and the erosion can be increased by wake wash. After evaluating the predicted erosion results at Locations A and B, it appears that soil type and the angle of the erosion profile significantly affect the erosion of the bank and bank toe. The total erosion area is larger at Location A because the erodible layer is thicker than at Location B. The erosion causes the slope of the bank toe to increase from an acute angle to 90°, after which the soil at the bottom of the profile continues to erode due to flowing water and eventually causes the soil above to collapse. The slide or rotation failure mode changes the bank slope with the action from the flowing water and surge waves as observed in the field surveys. The bank and bank toe erosion increases with an increase in water-surface elevation and water table level. Therefore, the protection measures are needed. According to the BSTEM results, the viable materials for erosion protection at Location B are large woody debris and riprap, or a combination of the two. However, none of the BSTEM recommended erosion protection measures work at Location A. Instead, stronger stabilization measures, such as concrete revetments or a sheet pile wall, can be considered, with proper designs following the USACE design manual [67] for the design of coastal revetments, seawalls, and bulkheads. The sheet pile wall was selected by the Texas Department of Transportation at Location A in 2019 as part of improvements to State Highway 73 along the SNWW. At Location B, the riprap was constructed in 2019. As shown, the BSTEM is a useful tool to evaluate measures for bank erosion protection.

5. Conclusions

This study focuses on the investigation of vessel-generated waves, their associated wake wash impacts on the streambank erosion, and the evaluation of erosion control and protection measures for reducing the bank erosion and consequent dredging. Vessel-generated waves at the City of Port Arthur Park dock (Location B) and the Golden Pass LNG terminal (Location A) were collected over a period of three weeks to a month using the Aquadopp Profiler. Based on the recorded wave data, it was found that the vessel-generated maximum wave height in the SNWW was dominated by the drawdown waves. At Location A, the vessel-generated maximum wave height was about 0.45 m and, at Location B, the value reached up to 0.55 m. Four models (named as SNWW-1, SNWW-2, SNWW-3, and SNWW-4) with calibrated model coefficients were developed and validated using the collected maximum wave-height data. The NSE value for each model in validation varied from 0.481 to 0.535. Based on both the calibrated and validated results, the SNWW-4 model was identified as the best predictive model for estimating the vessel-generated maximum wave height in the SNWW. The maximum wave height was found to strongly correlate with the square of the vessel speed, a kinetic energy related variable. It also appeared to be affected by the vessel geometry, vessel draft, and the channel cross-sectional area. The proposed SNWW-4 model is anticipated to provide reasonable estimation of the vessel-generated maximum wave heights in the SNWW. The SNWW-4 model can be applied to assist in general screening and has safety applications by identifying the potential large waves generated by a moving vessel in the SNWW. It can also provide wave information for a bank and channel-bottom erosion-related modeling study, which is important to the evaluation of effective bank protection measures and the resulting frequency of channel dredging. Vessel-generated waves at both monitoring locations in the SNWW were found to produce enough shear force to mobilize bed sediment and cause erosion. The surge waves can mobilize between 125–480 g/m² of bed sediment at a water depth of 1 m at the two data sampling sites. The bed erosion rate generally increases as the wave height increases and as the water depth decreases. Erosion protection measures such as large wood debris and riprap are suitable erosion control solutions at the Port Arthur site according to the results obtained by applying the bank stability toe erosion model (BSTEM). However, stronger stabilization measures, such as using the concrete revetments, a sheet pile wall, or an articulating block, are required to resist the erosive forces along the bank of Golden Pass site, where a frontal attack by vessel-generated larger waves is expected. On a final note, this paper presents a specific case study at The Sabine–Neches Waterway, Texas, USA. The wave measurements, vessel information, and soil samples are unique SNWW data sources used for the investigation of vessel-generated waves and consequent bank erosions. It is anticipated that the methodologies and modeling techniques established in this study may be applicable to the calculation of vessel-generated maximum wave height and bank erosion rate at other navigation channels if similar data can be collected there.