수리역학 및 알 이동 모델을 활용한 외래종 잉어 생존에 대한 생태적
추론

Ruichen Xu, Duane C. Chapman, Caroline M. Elliott, Bruce C. Call, Robert B. Jacobson, Binbin Wang

Abstract

Bighead carp (Hypophthalmichthys nobilis), silver carp (H. molitrix), black carp (Mylopharyngodon piceus), and grass carp (Ctenopharyngodon idella), are invasive species in North America. However, they hold significant economic importance as food sources in China. The drifting stage of carp eggs has received great attention because egg survival rate is strongly affected by river hydrodynamics. In this study, we explored egg-drift dynamics using computational fluid dynamics (CFD) models to infer potential egg settling zones based on mechanistic criteria from simulated turbulence in the Lower Missouri River. Using an 8-km reach, we simulated flow characteristics with four different discharges, representing 45–3% daily flow exceedance. The CFD results elucidate the highly heterogeneous spatial distribution of flow velocity, flow depth, turbulence kinetic energy (TKE), and the dissipation rate of TKE. The river hydrodynamics were used to determine potential egg settling zones using criteria based on shear velocity, vertical turbulence intensity, and Rouse number. Importantly, we examined the difference between hydrodynamic-inferred settling zones and settling zones predicted using an egg-drift transport model. The results indicate that hydrodynamic inference is useful in determining the ‘potential’ of egg settling, however, egg drifting paths should be taken into account to improve prediction. Our simulation results also indicate that the river turbulence does not surpass the laboratory-identified threshold to pose a threat to carp eggs.

Introduction

Bighead carp (Hypophthalmichthys nobilis), silver carp (H. molitrix), black carp (Mylopharyngodon piceus), and grass carp (Ctenopharyngodon idella), are considered invasive in North America. These species were imported into North America in the 1970’s to support aquaculture and escaped into the wild where they alter aquatic environments and food webs, resulting in undesirable ecological consequences^1,2,3. On the other hand, these carp species are important food sources in China, yet their populations in their native environment have been declining due to over-fishing and the negative effects on fish habitats resulting from dam construction^4,5. As either native or invasive species, it is of great importance to understand their life cycles in order to identify potential intervention strategies to control their populations⁶.

These rheophilic, broadcast-spawning carps exhibit prolific reproduction, with a single female carp capable of producing between 100,000 and one million eggs annually⁷. Carps typically engage in spawning during the spring and summer months when the temperature is within a range favorable for successful reproduction (peaking at roughly 20–24 ∘C) and during periods of high flows^8,9. They select specific locations for spawning characterized by high turbulence, including rocky rapids, riffles, islands, river confluences, and bends. This choice helps prevent the settling of eggs onto the riverbed, as sediment burial causes high mortality¹⁰. Within 3–5 h after spawning, eggs absorb a large amount of water in a process known as water hardening, leading to an increase in egg size and decrease in egg density. The water-hardening process leads to a decrease in settling velocity by approximately 70%, making eggs more likely to suspension in the water column^10,11.

After spawning and fertilization, the drift stage of carp eggs begins, a critical early-life stage in carp recruitment. Eggs hatch in approximately 30 h at optimal temperatures^10,12. During the drift stage before hatching, eggs are susceptible to predation, relying entirely on river currents and turbulence to remain suspended until hatch. After hatching, larval carp remain in the drift for a period, but they can behaviorally avoid settling^10,12. Because hydrodynamics plays a critical role in the suspension, dispersion, and transport of carp eggs across various scales in rivers, numerous studies have been conducted to explore river hydraulics and turbulence in relation to suitable carp spawning grounds, survival potential, and hatch locations^13,14,15,16. A key survival condition is the necessity for eggs to remain suspended in the water column throughout the entire egg drift stage, or at the very least, to avoid settling and being buried by sediment. Consequently, assessing whether river hydrodynamics can support this condition is a fundamental step in gauging recruitment success.

Flow velocity has been used as a simple indicator for assessing the suspension of eggs in rivers. For instance, Kocovsky et al.¹⁷ used a threshold velocity of 0.7 m/s as suitable for the spawn-to-hatch environment. Selection of 0.7 m/s is based on early literature with limited mechanistic studies^9,18. Lower critical flow velocities were also reported in the literature. Tang et al.¹⁹ suggested a value of 0.25 m/s based on a flume experiment, which agreed with some early field observations in the Yangtze River. Murphy and Jackson²⁰ found that mean velocities of 0.15–0.25 m/s allowed for egg suspension in four tributary rivers of the Great Lakes. Guo et al.²¹ suggested a critical flow velocity of 0.3 m/s in a flume experiment. Because rivers are largely non-uniform and vary in size and morphology, selecting a specific flow velocity as the sole empirical indicator for assessing suitability of carp recruitment is rather challenging.

While using flow velocity as an indicator for examining egg suspension or settling might be practical, it does not fully represent the underlying physics, especially in areas where turbulence is not well correlated with mean flow velocity. To account for the mechanism of egg suspension, Garcia et al.²² proposed three different criteria involving the ratio of shear velocity and egg settling velocity, the ratio of vertical turbulence intensity and egg settling velocity, and the Rouse number to predict the suspension and settling of carp eggs. In their laboratory experiment, they observed that 65% of eggs remained in suspension with a mean flow velocity of 0.07 m/s, corresponding to a Rouse number of 1.32 and shear velocity of 0.004 m/s. At higher flow velocities of 0.2 and 0.4 m/s, with Rouse numbers of 0.57 and 0.58 and shear velocity of 0.008 and 0.016 m/s, respectively, all eggs were in suspension. These observations agree well with the empirical values of Rouse number classification for sediment transport for bedload, partial suspension, full suspension, and washload²³. Therefore, using these parameters is better supported by the mechanism of particle suspension compared to velocity alone.

Given the above simple criteria of using shear velocity or Rouse number, hydraulic models or measurements can be used to infer whether a stream or a river reach can support a favorable environment for egg suspension in the egg-drift stage¹⁷. In addition, three dimensional hydrodynamic models can provide additional insights into the spatial distributions of potential egg settling zones, given the strong spatial heterogeneity of river turbulence^24,25,26. In this paper, we use an 8-km reach in the Lower Missouri River as representative of channelized segments of the Upper Mississippi River basin where carps are established. We used computational fluid dynamics (CFD) modeling to explore the overall suitability for egg drift and to infer potential egg settling zones, with an emphasis on understanding the spatial distributions of hydrodynamics associated with in-stream hydraulic structures, river morphology, and strong topographic gradients on the riverbed. Specifically, we examine the criteria of egg suspension and evaluate the locations where the hydrodynamics are unfavorable for suspending eggs. Our objective is to evaluate whether the potential egg settling zones based on hydrodynamic inference would agree with entrapment locations that can be estimated using drift models. We additionally evaluate whether turbulence conditions indicated in the model approach criteria for turbulence-induced damage to carp eggs as determined in laboratory studies.

Methods

Study site

The study site is a selected reach in the Lower Missouri River near Lexington, Missouri (Fig. 1). The reach is approximately 8 km long with a sinuosity index of 1.12. The mean bankful width is 331.4 m. The bed is mostly covered by medium and coarse sand (D50 = 0.55 mm) with fine muddy materials (< 0.125 mm) near the banks and close to the dike fields^27,28. The mean annual discharge is approximately 1700 m3/s measured at a U.S. Geological Survey (USGS) gaging station approximately 24 km downstream (station no. 06895500, Waverly, Missouri, USGS). The reach is representative of rivers that have been highly engineered to support navigation and bank stability, with complex hydraulic conditions where water flows around and over the rock channel-training structures^29,30. This reach has been used as the main site for model development stage of SDrift^31,32, an egg drift model used in this study. The previous studies have accumulated substantial data for the bathymetric-topographic digital elevation model (DEM), water surface elevations, and cross-channel velocity profiles³³, which have been used for calibration and validation of our CFD model.

Hydrodynamic model

The flow was simulated using FLOW-3D HYDRO with a Reynolds-averaged Navier-Stokes (RANS) solver and a Re-Normalization Group (RNG) modified k−ε turbulence sub-model. The model was set up for solving the steady-state flows under four discharge conditions ( Q = 1342, 2282, 3060 and 4219 m3/s, referred to as Q1 to Q4 conditions), which correspond to approximately 45–3% daily flow exceedance during spawning season. A Cartesian mesh with a final size of 4×4×0.4 m in the east-north-up coordinate system was used after a mesh independence study to evaluate optimal mesh dimensions³¹.

The upstream and downstream boundary conditions were set to the measured flow discharge and calculated hydrostatic pressure from the measured water-surface elevation, respectively. The model was calibrated by adjusting the roughness coefficient until the simulated water-surface elevations agree with the measured data, where the water-surface elevations were measured using a ship-mounted, real-time corrected kinematic global navigation satellite system (RTK-GNSS). The measured cross-channel velocities at 22 locations at two flow conditions (Q=2282 and 3060 m3/s) were used to evaluate model performance, where the velocities were measured using a ship-board acoustic Doppler current profiler (ADCP, Workhorse Rio Grande, Teledyne, Inc) at each cross section with four repeated transects. The ADCP had a vertical resolution of 0.5 m and horizontal resolution of 1 m. The velocities within 1 m below the water surface and within 1 m above the river bed were not measured due to instrument blanking distance and measurement noise. Additional details on model calibration and evaluation are in Li et al.³¹.

Egg drift model

The egg drift model SDrift was used for egg transport modeling in this study³¹. This model uses Lagrangian particle tracking to simulate the transport of carp eggs, where turbulent fluctuations are modeled using an explicit solution for the Langevin equation, i.e., the Markov-chain continuous random walk (CRW) algorithm^34,35,36. The density and diameter of carp eggs were determined as a function of post-fertilization time and water temperature based on the regression equation to the laboratory measured data¹¹. The details of regression can be found in³¹. The time-varying characteristics of eggs result in evolving egg settling velocity in the water, which is determined based on the drag law for spherical particles³⁷.

SDrift was incorporated with the CFD model outputs to predict transport of silver carp eggs in the selected reach. A broad surface-spawning event across the entire cross section at an upstream location in the model (x= 427,130 m, near River Mile 314) was simulated by releasing 6600 model eggs on the water surface at 33 locations³¹. All eggs were tracked until they were transported outside the downstream boundary or ‘entrapped’ in the model domain determined by the model criterion.

Criterion of egg entrapment from the egg drift model

SDrift allows the simulated eggs to be ‘entrapped’ if they are stationary for a pre-defined duration. The entrapment would occur if a simulated egg is transported into a low velocity zone and eventually loses its momentum. From the model evaluation, entrapment primarily occurs in the region with high topographic gradients, e.g., near the bank and hydraulic structures. A duration of 30 s was used here to determine the entrapment, i.e., if a simulated egg does not move for 30 s, it would be considered entrapped and would no longer be tracked. Although the entrapment does not necessarily provide a certain prediction of egg settling, it offers insight into locations where the eggs may be stopped and eventually buried by bed sediment. The selection of a 30-s duration is somewhat arbitrary. From a physics standpoint, this duration should ideally exceed the largest turbulent time scale. However, due to the extensive spatial scale of the modeled reach and the river-training structures, the turbulent time scale varies significantly across space. Furthermore, both the spatial resolution in the CFD simulation and the temporal resolution in particle tracking have the potential to influence particle movements and their entrapment. Therefore, determining the optimal duration requires further investigation in future studies.

Criterion of egg suspension and settling from the hydrodynamic model

Suspension of carp eggs depends on whether the flow can provide adequate upward motions that overcome their settling. Analogous to sediment suspension and transport³⁸, several means have been used to quantify the settling and suspension of carp eggs in turbulent flows. Here we analyze three parameters following Garcia et al.²²: the ratio between shear velocity and settling velocity, the ratio between vertical turbulence intensity and settling velocity, and the Rouse number.

Shear velocity

Shear velocity (u∗) is a velocity scale defined from the bed shear stress. The ratio of shear velocity and particle terminal velocity (wt), a so-called movability number (M∗=u∗/wt), has been used to classify sediment transport³⁹. Different critical values have been proposed to define particle suspension^38,39. Here, the critical value of 1.0 is used following the studies of carp eggs^20,22: locations with u∗/wt<1 are the potential settling zones of carp eggs, where particle terminal velocity is the egg settling velocity (wt=Vegg).

Because shear velocity only represents the bed shear but does not provide the vertical variability in the water column, we applied a scaling method so that potential egg suspension and settling can be evaluated in the entire water column. Using the relationship between bed shear and turbulence kinetic energy (TKE)^40,41, i.e., τb=C1ρk with C1=0.19⁴⁰, the movability number can be estimated at every grid point using the TKE determined from the CFD simulation:

The potential egg settling zones were then determined based on M∗<1.

Vertical turbulence intensity

The vertical turbulence intensity (wrms′) is a direct parameter to quantify the turbulent velocity scale in the vertical direction, which can be used to define the initiation of particle suspension³⁸. Therefore, we also calculated the ratio between wrms′ and V{egg} as the second indicator for egg settling: locations with wrms′/Vegg<1 are the potential settling zones of carp eggs. Here, we estimated w′ based on anisotropy of turbulent fluctuations in open channel flows:

with Du=2.30, Dv=1.27, and Dw=1.63⁴². This gives wrms′/Vegg=0.75TKE/Vegg where TKE was obtained from the CFD simulations.

Rouse number

In sediment transport, the Rouse number has been used to describe the suspended load³⁸. The Rouse number is defined as Ro=wt/(βκu∗) with wt=Vegg for carp eggs, where κ is von Kárman constant and β is a coefficient related to diffusion of particles^22,23:

The Rouse number (Ro, also used as Z or P in the literature), can be used to classify the sediment transport similar to the movability number. Hearn²³ suggested that sediment particles are in 100% suspension or wash load when Ro<1.2; particles are partially suspended when 1.2<Ro<2.5; particles are predominantly transported by bedload if Ro>2.5. Here, we use 1.2 as the criterion, such that the potential egg settling zones were determined based on Ro>1.2.

Results and discussion

Model calibration and evaluation

The model calibration results for water-surface elevation are shown in Fig. 2 for four flow conditions³¹. The elevation of river bed in the main channel is also plotted for reference. The root-mean-square-error (RMSE) in the water surface elevation between the measurement and modeling is 0.07, 0.03, 0.04, and 0.03 m, for Q1 to Q4, respectively. The RMSE is considered to be small compared to the length of the reach and the water depths.

The measurement-modeling comparison of double-averaged velocities over the flow depth and the cross section in both streamwise (Us) and transverse (Ut) directions is given in Fig. 3 for two measured conditions (Q2 and Q3). The RMSE of Us and Ut is 0.055 and 0.028 m/s, much smaller than the mean flow of 1.29 and 1.38 m/s in the measured cross sections for Q2 and Q3, respectively. The direct measurement-modeling comparison in all 22 cross sections is given in the supplementary file (Figs. S1 and S2).

Mean flow characteristics

The CFD simulated flow depth and depth-averaged flow velocity for two out of four conditions are shown in Figures 4 and 5. Greater depths are located downstream from the dikes (i.e., in scour holes) and near the right bank at the upstream bend (i.e., Easting 431,000–432,000 m, downstream of river mile 311). Shallower depths are located upstream from the dikes and along the left bank in the downstream bend (i.e., Easting 432,000–433,500 m, in the vicinity of river mile 310).

Flow velocities are greater at a higher discharge, and are strongly related to the in-stream hydraulic structures: high velocities are located within the main channel and low velocities are located close to the dike areas and both sides of the bank. For Q1, the L-head dikes on the left bank around Easting 430,500–431,000 m (upstream of river mile 311) block the flow into the left bank, resulting in channel narrowing and an area of localized higher velocity. Relatively faster velocities are also located close to the right bank from Easting 432,000–433,500 m (in the vicinity of river mile 310) and then shaped by the L-head dike at Easting 433,500–434,500 m (between river miles 309 and 310). When water enters the L-head dike area at Easting 430,500–431,000 m (between river miles 311 and 312) in high discharge conditions (e.g., Q4), the localized fast flow is not observed.

Turbulence quantities

Two turbulence quantities were selected to elucidate the turbulence in the reach: the depth-averaged TKE (Fig. 6) and the depth-averaged dissipation rate of TKE (Fig. 7). For Q1, TKE shows a similar spatial pattern as the flow velocity, indicating that the high TKE is usually associated with high velocities. For Q4, additional high TKE regions are located within the low velocity zones near the dikes. These high turbulence regions are caused by the interaction of flow with the hydraulic structures. For instance, enhanced turbulence may occur within wakes downstream from the flows over the dikes. Strong shear-induced turbulence may also occur at the water surface near the edge of the dikes close to the main channel. Similar to TKE, the locations of high TKE dissipation rate are coincident with high velocity in the main channel and near the dikes where strong flow-structure interactions occur.

To examine the correlation between turbulence and the mean flow in the reach, Fig. 8 elucidates the ratio between TKE and the mean kinetic energy (MKE) where MKE is defined based on mean velocity values, MKE = 0.5(U2+V2+W2). The data show that the TKE/MKE ratio is much smaller than 1 in the main channel, a typical open-channel feature. However, near the river bank and in the dike fields, greater TKE than MKE is common, with the spatial distribution of TKE/MKE>1 being dependent on discharge. This result documents strong interactions between water flow and the solid boundaries, which generate substantial turbulence comparing to the reduced mean velocity in these regions. Within these regions, particles would be expected to have longer residence times³².

Egg suspension and settling

The CFD modeling results allow for analysis of potential egg settling zones based on the criteria of particle suspension outlined in section “Criterion of egg suspension and settling from the hydrodynamic model”. In Fig. 9, the potential egg settling locations are plotted based on the Rouse number criterion for all four discharge conditions. The plot shows that potential settling zones are located near the river banks, in dike fields, and even in the channel at locations with strong gradients in the bed morphology. We note that the criterion was applied to all data points simulated in the CFD. Therefore, the settling zones represent the x–y locations where turbulence is inadequate to suspend eggs. Not surprisingly, the estimated potential settling zones become smaller with increasing discharge. Results using shear velocity and vertical turbulence intensity criteria show similar results, which are plotted in the supplementary file (Figs. S3 and S4).

Figure 10 shows the predicted locations of entrapped eggs using the egg drift model, SDrift³¹. Comparing Fig. 10 with Fig. 9, we found that both hydrodynamic-inferred potential settling locations and drift-model predicted locations include the regions near the dike fields and the sparse areas in the channel where strong topographic gradients are present. However, careful examination of the wing dike areas (Fig. 11 under Q1 condition and Fig. 12 under Q4 condition), shows that the predicted egg settling zones using two methods are located in different regions near the dike areas. SDrift results indicate that egg entrapment is mainly located adjacent to the dikes, whereas the hydrodynamic inference indicates strong egg settling potential downstream from the dikes under low-flow conditions, such as the discharge condition Q1 (Fig. 11). The potential egg settling zones are substantially decreased by increasing discharge (Fig. 12). SDrift results indicate that egg entrapment is primarily due to interception of egg movement due to strong topographic gradients near the dikes while being tracked in the model under these hydrodynamic conditions. Although this does not directly imply that the eggs would settle in these areas, higher probability of egg-dike interaction would occur that could potentially affect egg survival. In contrast, the hydrodynamic inference only suggests hydrodynamic conditions that are favorable for egg settling, which differs from the drift models.

In addition, the drift model predicts substantial egg entrapment near the left bank upstream of the bend located around x=43,100 m (upstream of river mile 311), where these regions were not inferred from hydrodynamic data. The differences indicate that eggs can be entrapped within locations where hydrodynamics would indicate suspension. The potential entrapment in the drift model is likely due to the reduction in egg-drift speed close to the left bank, which increases the probability of egg settling. In curved rivers reaches, the unevenly distributed flow in the cross section and secondary flow may push eggs towards the outer side of the channel, which can increase the probability of the particle-bank interaction.

The drift trajectories of 200 simulated eggs released near the left bank for discharge Q1 and Q4 can be used to visualize drift dynamics simulated in SDrift (Fig. 13). The modeling results show that, under Q1, there is minimal egg drift into the low-flow region between the L-head dikes and the left bank in Area 1, as well as into the high-riverbed region close to the left bank in Area 2. This restriction occurs because the elevation of the dikes in Area 1 are higher than the water surface elevation during low-flow conditions, preventing eggs from entering these areas. As a result, the drift model predicts minimal entrapment of eggs in these areas. However, the hydrodynamic inference only takes into account favorable conditions for egg settling, implying significant settling in these regions even when trajectories would fail to transport eggs into the areas. Nevertheless, under higher-flow conditions that permit eggs to enter these areas (see Fig. 13b), particularly in Area 1, entrapment of eggs can occur (see Fig. 10), even though the hydrodynamic inference does not indicate significant settling compared to other low-velocity areas.

Vertical distribution of potential egg settling zones

To examine the likelihood of egg settling based on vertical position in the water column, the number of cells were counted that satisfy the criterion of egg settling based on hydrodynamic inference at the same vertical height above the riverbed (z) under the four simulated discharges. Figure 14 illustrates an example based on Rouse number criterion. The results show that the flow condition of Q1 has substantially more counts (about one order of magnitude) due to weaker turbulence compared to the other three flow conditions (Fig. 14a and b). We interpret this large change between Q1 and higher discharges as a threshold resulting when flows begin to overtop the wing dikes. Overtopping flows substantially decrease low-turbulence areas downstream and landward of wing dikes.

The modeling data also indicate that egg settling is more likely to occur in the lower part of water column but not near the riverbed. Taking Q1 as an example, the peak of the number of counts are located about 2 m above the riverbed, with the number of counts decreasing both towards surface and towards the riverbed (Fig. 14a). In the normalized water column profile (Fig. 14b), substantial counts are located within the bottom 20% of the water column. We note that various water depths occur across the river reach, and hence the number of counts on the x-axis of the plots (Fig. 14a and b) are different before and after the water column normalization.

Examining the probability distribution function (PDF), we found that four discharge conditions show similar vertical profiles: egg settling has more than 10% probability within approximately the bottom 5 m (Fig. 14c), corresponding to approximately the bottom 20% of water depth (Fig. 14d). This result suggests that when eggs are transported to the bottom 20% layer, the hydrodynamic condition is less favorable for them to be re-suspended compared to higher-up in the water column. Similar results of profiles were found for the criterion using shear velocity and the vertical turbulence intensity, albeit the number of counts and the PDF values are different due to different criteria (see supplementary file, Figs. S5 and S6).

Discussion on the egg survival

Examining river hydrodynamics in three dimensions through well-calibrated models yields valuable insights into the spatial distribution of flow velocity, water depth, and associated turbulence. These parameters can be used to identify potential locations where carp eggs may settle. However, using and interpreting results based on hydrodynamic criteria must be exercised with careful consideration. For instance, the Rouse number classification for particle suspension involves a broad range of values. In this study, we adopted Ro>1.2 as an indicator of egg settling, with Ro=1.2 representing the lower Rouse number bound for partial suspension. Conservatively, a critical value of Ro=2.5²³ is recommended for assessing predominantly bedload particle transport, indicating minimal to no suspension in the water column. Hence, at Rouse numbers between 1.2 and 2.5, partial suspension would be expected. In addition, the analysis using three-dimensional drift model results indicates that carp eggs would not drift into the egg settling zones within the L-head dikes and left bank (Area 1 in Fig. 13), for example, which would have predicted settling using hydrodynamic inference under the low-flow condition. This is because the actual egg drift pathway is governed by various parameters including egg spawning locations, streamlines of water flows, and interactions of flow and hydraulic structures. Consequently, predictions relevant to invasive carp management would improve when using the hydrodynamic-inferred egg settling zones if these additional parameters were taken into account.

Although egg settling zones based on hydrodynamic inference may not represent the actual conditions for egg settling, those predictions provide valuable information about the local hydrodynamics and suitability for egg settling at lower computational cost compared to drift modeling (for example SDrift). Therefore, this information could be useful for managers in determining the desirability of implementing hydraulic controls for egg settling. For example, if flow patterns can be adjusted to guide eggs into low-turbulence zones with adequate residence time, the hydrodynamics would facilitate the desired settling of eggs, aligning with management objectives for controlling aquatic invasive species. However we noted that solely using hydrodynamic inference may be misleading in invasive carp management without knowledge of drift pathways.

While high turbulence zones are the necessary environment for carp eggs to be suspended, eggs can be damaged or killed if turbulence exceeds a certain threshold. Prada et al.⁴³ found an increased mortality in drifting grass carp eggs when exposed to turbulence with TKE greater than 2 m2/s2 for 1 minute in a grid-stirred turbulence tank. When TKE reaches 2.7 m2/s2, the mortality rate increased by nearly 30%. The corresponding maximal shear stresses were found to be 20 and 30 N/m2 near the grid for these two TKE values respectively. From our hydrodynamic model, mean TKE in the simulated reach under discharges Q1 to Q4 ranges from 0.01 to 0.02 m2/s2, with maximal depth-averaged TKE ranging from 0.16 to 0.21 m2/s2. The maximal TKE in the water column is found within 0.31–0.38 m2/s2 under four discharge conditions. These values are much smaller than the reported values that are harmful for carp eggs. Therefore, in a typical egg drift process, it is unlikely for eggs to experience persistent, extreme turbulence that could cause direct damage or mortality.

However, strong turbulence often generates high suspension and transport of sediment in the river. The abrasion between carp eggs and the suspended sediment may affect the egg survival rate. In the laboratory experiment conducted by Prada et al.¹⁵, carp eggs were found to drift within the lower 75% of the water column with lower flow velocity in the flume (0.08 m/s). When the flow velocity was increased to 0.22 m/s, the egg distribution in the water column was uniform, indicating a well-suspended condition for carp eggs. With further increasing flow velocity, Prada et al.¹⁵ observed that eggs were drifting more towards the bottom where they collided with the sediment particles. This indicates that the suspension of sediment could affect the vertical distribution of suspended eggs. They also observed reduced survival rate in medium and high flows compared to the control, while the survival rate was almost the same in low flow compared to the control. They also observed different larvae behaviors in different flow velocities, which may also contribute to the survival of carps. In our simulated Missouri River reach, the river turbulence may not pose a threat to carp eggs, but the suspended sediment could have negative effects. There has been limited study on the quantitative effects of sediment abrasion on egg mortality, indicating a fruitful subject for future studies.

Conclusions

In this study, we analyzed the simulated hydrodynamics of an 8-km reach in the Lower Missouri River, a site characterized by extensive channelization and river training. Four discharges representing 45–3% daily flow exceedance were examined. Calibration and validation of the simulations were conducted based on field observations. Flow depth, mean flow velocity, and turbulence quantities were investigated through computational fluid dynamics modeling. Simulated results show highly varied spatial distributions of mean flow and turbulence characteristics, primarily attributed to the curvature of the channel, variation in bed morphology, and the presence of river-training hydraulic structures, including wing dikes and L-head dikes.

To investigate the use of hydrodynamics for inferring the settling and suspension of carp eggs, we applied three criteria established in previous carp egg studies to analyze the spatial distribution of potential settling zones. The simulation results enabled the identification of low turbulence zones where insufficient suspension may hinder carp egg development. When comparing these hydrodynamic-inferred egg settling zones with the entrapment predicted by a Lagrangian egg-drift model, we observed that egg drift paths significantly influenced the locations where eggs may settle or be intercepted by in-stream hydraulic structures. Therefore, it is crucial to consider additional factors, such as spawning locations and drift paths, when using hydrodynamic inference to identify potential egg settling zones and larval nursery locations for invasive carp management.

Lastly, river turbulence may also influence carp egg survival through shear stresses and interactions with suspended sediment. Our data indicate that turbulence kinetic energy in the river does not surpass the laboratory-identified threshold associated with direct egg damage. However, abrasion from suspended sediment and the complex interactions between eggs and hydraulic structures, riverbed, and banks, accentuated by high morphological variations as demonstrated in the entrapment areas in the egg drift model, could affect the overall survival rate of carp eggs.

Data availibility

The data of field measurements and modeling are available in the online repository doi:10.5066/P9X5M3WH³³.

References

1. Irons, K. S., Sass, G. G., Mcclelland, M. A. & Stafford, J. D. Reduced condition factor of two native fish species coincident with invasion of non-native Asian carps in the Illinois River, USA—is this evidence for competition and reduced fitness?. J. Fish Biol. 71, 258–273. https://doi.org/10.1111/j.1095-8649.2007.01670.x (2007).
2. Cudmore, B., Mandrak, N. E., Dettmers, J. M., Chapman, D. C., & Kolar, C. S. Binational ecological risk assessment of bigheaded carps (Hypophthalmichthys spp.) for the Great Lakes Basin. Technical Report 2011/114 (2012).
3. Chick, J. H., Gibson-Reinemer, D. K., Soeken-Gittinger, L. & Casper, A. F. Invasive silver carp is empirically linked to declines of native sport fish in the Upper Mississippi River system. Biol. Invasions 22(2), 723–734. https://doi.org/10.1007/s10530-019-02124-4 (2020).
4. Chapman, D. C. et al. Bigheaded carps of the Yangtze and Mississippi Rivers. In Fishery Resources, Environment, and Conservation in the Mississippi and Yangtze (Changjiang) River Basins (eds. Chen, Y. et al.) 113–127 (American Fisheries Society, 2016). https://doi.org/10.47886/9781934874448.
5. Tang, C., Yan, Q., Li, W., Yang, X. & Zhang, S. Impact of dam construction on the spawning grounds of the four major Chinese carps in the three gorges reservoir. J. Hydrol. 609, 127694. https://doi.org/10.1016/j.jhydrol.2022.127694 (2022).
6. Chapman, D. C. et al. Status of the major aquaculture carps of China in the Laurentian Great Lakes Basin. J. Great Lakes Res. 47(1), 3–13. https://doi.org/10.1016/j.jglr.2020.07.018 (2021).
7. Garcia, T., Jackson, P. R., Murphy, E. A., Valocchi, A. J. & Garcia, M. H. Development of a fluvial egg drift simulator to evaluate the transport and dispersion of Asian carp eggs in rivers. Ecol. Model. 263, 211–222. https://doi.org/10.1016/j.ecolmodel.2013.05.005 (2013).
8. Deters, J. E., Chapman, D. C. & McElroy, B. Location and timing of Asian carp spawning in the Lower Missouri River. Environ. Biol. Fishes 96(5), 617–629. https://doi.org/10.1007/s10641-012-0052-z (2013).
9. Yih, P. & Liang, T. Natural conditions of the spawning grounds of the domestic fishes in Yangtze River and essential external factor for spawning. Act Hydrobiol. Sin. 5, 1–15 (1964).
10. George, A. E. & Chapman, D. C. Embryonic and larval development and early behavior in grass carp, Ctenopharyngodon idella: Implications for recruitment in rivers. PLoS ONE 10, 3. https://doi.org/10.1371/journal.pone.0119023 (2015).
11. George, A. E., Garcia, T. & Chapman, D. C. Comparison of size, terminal fall velocity, and density of bighead carp, silver carp, and grass carp eggs for use in drift modeling. Trans. Am. Fish. Soc. 146(5), 834–843. https://doi.org/10.1080/00028487.2017.1310136 (2017).
12. George, A. E. & Chapman, D. C. Aspects of embryonic and larval development in bighead carp Hypophthalmichthys nobilis and silver carp Hypophthalmichthys molitrix. PLoS ONE 8(8), 1932–6203. https://doi.org/10.1371/journal.pone.0073829 (2013).
13. Kasprak, A., Jackson, P. R., Lindroth, E. M., Lund, J. W. & Ziegeweid, J. R. The role of hydraulic and geomorphic complexity in predicting invasive carp spawning potential: St Croix River, Minnesota and Wisconsin, United States. PLoS ONE 17(2), e0263052. https://doi.org/10.1371/journal.pone.0263052 (2022).
14. Zhu, Z. et al. Using reverse-time egg transport analysis for predicting Asian carp spawning grounds in the Illinois River. Ecol. Model. 384, 53–62. https://doi.org/10.1016/j.ecolmodel.2018.06.003 (2018).
15. Prada, A. F., George, A. E., Stahlschmidt, B. H., Chapman, D. C. & Tinoco, R. O. Survival and drifting patterns of grass carp eggs and larvae in response to interactions with flow and sediment in a laboratory flume. PLoS ONE 13(12), 1932–6203. https://doi.org/10.1371/journal.pone.0208326 (2018).
16. Heer, T., Wells, M. G. & Mandrak, N. E. Asian carp spawning success: Predictions from a 3-d hydrodynamic model for a Laurentian Great Lake tributary. J. Great Lakes Res. 47(1), 37–47. https://doi.org/10.1016/j.jglr.2020.07.007 (2021).
17. Kocovsky, P. M., Chapman, D. C. & McKenna, J. E. Thermal and hydrologic suitability of Lake Erie and its major tributaries for spawning of Asian carps. J. Great Lakes Res. 38(1), 159–166. https://doi.org/10.1016/j.jglr.2011.11.015 (2012).
18. Abdusamadov, A. S. Biology of white amur, Ctenopharyngodon idella, silver carp, Hypophthalmichthys molitrix, and bighead, Aristichthys nobilis, acclimatized in the Terek region of the Caspian basin. Vopr. Ikhtiologii 3, 425–433 (1987).
19. Tang, M., Huang, D., Huang, L., Xiang, F. & Yin, W. Preliminary forecast of hydraulic characteristic test of grass, green, silver carp, bighead carp egg incubation conditions in the three gorges reservoir area. Reserv. Fisher. 4, 26–30 (1989).
20. Murphy, E. A. & Jackson, P. R. Hydraulic and water-quality data collection for the investigation of Great Lakes Tributaries for Asian carp spawning and egg-transport suitability. Report 2013-5106 (2013).
21. Guo, H. et al. Settling and transport properties of grass carp and silver carp eggs in the water-hardened phase: Implications for resource protection and invasion control during early life period. Ecol. Ind. 148, 110064. https://doi.org/10.1016/j.ecolind.2023.110064 (2023).
22. Garcia, T. et al. A laboratory investigation of the suspension, transport, and settling of silver carp eggs using synthetic surrogates. PLoS ONE 10(12), 1–19. https://doi.org/10.1371/journal.pone.0145775 (2016).
23. Hearn, C. J. The Dynamics of Coastal Models (Cambridge University Press, 2008).
24. Sukhodolov, A. N. et al. Turbulent flow structure at a discordant river confluence: Asymmetric jet dynamics with implications for channel morphology. J. Geophys. Res. Earth Surf. 122(6), 1278–1293. https://doi.org/10.1002/2016JF004126 (2017).
25. Le, T. B. et al. Large-eddy simulation of the Mississippi River under base-flow condition: Hydrodynamics of a natural diffluence-confluence region. J. Hydraul. Res. 57(6), 836–851. https://doi.org/10.1080/00221686.2018.1534282 (2019).
26. Li, G. et al. Turbulence near a sandbar island in the lower Missouri River. River Res. Appl. 39(9), 1857–1874. https://doi.org/10.1002/rra.4180 (2023).
27. Gaeuman, D. & Jacobson, R. B. Acoustic bed velocity and bed load dynamics in a large sand bed river. J. Geophys. Res.-Earth Surface 111, 2. https://doi.org/10.1029/2005jf000411 (2006)
28. .Poulton, B. C. & Allert, A. L. An evaluation of the relative quality of dike pools for benthic macroinvertebrates in the Lower Missouri River, USA. River Res. Appl. 28(10), 1658–1679. https://doi.org/10.1002/rra.1558 (2012).
29. Galat, D. L. & Lipkin, R. Restoring ecological integrity of great rivers: Historical hydrographs aid in defining reference conditions for the Missouri River. Hydrobiologia 422, 29–48. https://doi.org/10.1023/A:1017052319056 (2000).
30. Jacobson, R. B. & Galat, D. L. Flow and form in rehabilitation of large-river ecosystems: An example from the Lower Missouri River. Geomorphology 77(3–4), 249–269. https://doi.org/10.1016/j.geomorph.2006.01.014 (2006).
31. Li, G. et al. A three-dimensional Lagrangian particle tracking model for predicting transport of eggs of rheophilic-spawning carps in turbulent rivers. Ecol. Model. 470, 110035 (2022).
32. Li, G. et al. Evaluations of lagrangian egg drift models: From a laboratory flume to large channelized rivers. Ecol. Model. 475, 110200 (2023).
33. Elliott, C. M., Call, B. C., Li, G., & Wang, B. Field data and models of the Missouri River at Sheepnose Bend, near Lexington, Missouri, 2019–2021. In U.S. Geological Survey data release. https://doi.org/10.5066/P9X5M3WH (2022).
34. Bocksell, T. L. & Loth, E. Random walk models for particle diffusion in free-shear flows. AIAA J. 39(6), 1086–1096. https://doi.org/10.2514/2.1421 (2001).
35. Wang, B., Huijie, W. & Wan, X.-F. Transport and fate of human expiratory droplets-a modeling approach. Phys. Fluids 32(8), 083307. https://doi.org/10.1063/5.0021280 (2020).
36. Wang, B., Sullivan, L. L. & Wood, J. D. Modeling wind-driven seed dispersal using a coupled Lagrangian particle tracking and 1-D k-ε turbulence model. Ecol. Model.486, 110503. https://doi.org/10.1016/j.ecolmodel.2023.110503 (2023).
37. Goossens-Walter, R. A. Review of the empirical correlations for the drag coefficient of rigid spheres. Powder Technol. 352, 350–359. https://doi.org/10.1016/j.powtec.2019.04.075 (2019).
38. van Rijn, L. C. Sediment transport, part ii: Suspended load transport. J. Hydraul. Eng. 110(11), 1613–1641. https://doi.org/10.1061/(ASCE)0733-9429(1984)110:11(1613) (1984).
39. Dey, S. A., Zeeshan, S. K. & Padhi, E. Terminal fall velocity: the legacy of stokes from the perspective of fluvial hydraulics. Proc. R. Soc. A: Math. Phys. Eng. Sci. 475(2228), 20190277. https://doi.org/10.1098/rspa.2019.0277 (2019).
40. Kim, S.-C., Friedrichs, C. T., Maa, J.P.-Y. & Wright, L. D. Estimating bottom stress in tidal boundary layer from acoustic doppler velocimeter data. J. Hydraul. Eng. 126(6), 399–406. https://doi.org/10.1061/(ASCE)0733-9429(2000)126:6(399) (2000).
41. Liao, Q., Wang, B. & Wang, P.-F. In situ measurement of sediment resuspension caused by propeller wash with an underwater particle image velocimetry and an acoustic doppler velocimeter. Flow Meas. Instrum. 41, 1–9. https://doi.org/10.1016/j.flowmeasinst.2014.10.008 (2015).
42. Nezu, I. & Nakagawa, H. Turbulence in Open Channel Flows (Routledge, 1993).
43. Prada, A. F. et al. Influence of turbulence and in-stream structures on the transport and survival of grass carp eggs and larvae at various developmental stages. Aquat. Sci. 82, 1. https://doi.org/10.1007/s00027-019-0689-1 (2020).

단일 파동 하에서 투과성 포장지가 있는 현수 파이프라인의 동적 성능

Youkou Dong, Enjin Zhao, Lan Cui, Yizhe Li, Yang Wang

Abstract

Submarine pipelines are widely adopted around the world for transporting oil and gas from offshore fields. They tend to be severely ruined by the extreme waves induced by the natural disaster, such as hurricanes and tsunamis. To maintain the safety and function integrity of the pipelines, porous media have been used to wrap them from the external loads by the submarine environment. The functions of the porous wrappers under the hydrodynamic impact remain to be uncovered before they are widely accepted by the industry. In this study, a numerical wave tank is established with the immersed boundary method as one of the computational fluid dynamics. The submarine pipelines and their porous wrappers are two-way-coupled in terms of displacement and pressure at their interfaces. The impact from the solitary waves, which approximately represent the extreme waves in the reality, on the pipelines with different configurations of the porous wrapper is investigated. The results present significant protective functions of the wrappers on the internal pipelines, transferring the impact forces from the pipelines to the wrappers. The protective effects tend to be enhanced by the porosity and thickness of the wrappers. The influence of the pipeline configurations and the marine environment are then analysed. As for the front pipeline, an increase in the gap leads to a slight increase in the horizontal forces on both the wrapper and the pipeline, but a significant increase in the vertical forces. As for the rear pipeline, because of the shield function of the front pipeline, the velocity within the gap space and the forces on the pipes are decreased with the decrease in the gap size. The complex flow fields around the pipelines with wrappers are also illuminated, implying that the protection function of the wrapper is enhanced by the wave height reduction.

Keywords

extreme wave; submarine pipeline; external wrapper; coupling analysis; computational fluid dynamics

1. Introduction

Pipelines that are laid on or below the seabed and continuously transport large amounts of oil (or gas) are collectively referred to as submarine pipelines. They constitute the main transporting structures and currently they are the most economical and reliable selections in the design of transportation tools. Pipelines are usually installed within the seabed sediments under the protection of rock berms [1]. However, the sediments around the pipelines may be scoured by contour currents and internal waves, which expose the pipelines to the threat of complex marine environments [2]. The scour mechanism and its evolution process around the in-position pipelines were investigated by many scholars, such as Reference [3]. Occasionally, segments of a pipeline may be suspended between high points through continental slopes due to an uneven seabed profile. For example, suspended pipelines were widely used in the Ormen Lange projects, with massive depressions and landslide blocks scattered along the 120-km-long route [4].
Natural disaster, such as hurricanes and tsunamis, may induce extreme waves that generate enormous impact loads on the pipelines and may cause serious ruins to the whole production and transportation system [5,6,7]. Tsunamis, one of the major marine disasters caused by earthquakes and submarine landslides [8,9], send surges of water with extremely long waves that are not especially steep [10]. The tsunami triggered by a 9.0-Mw earthquake in 2011 extensively destroyed 70% of the total 200,000 structures along the Miyagi coastline, including submarine pipelines, seawalls, and coastal bridges. A tsunami is typically composed of several transient waves with varying amplitudes, wave-lengths, and wave periods during propagation. Solitary waves were proposed to simulate the tsunami waves by decomposing them into N-waves through the Korteweg-de Vries equation [11,12,13,14]. Since then, the run-up process of the tsunami waves along the shoreline was investigated with the depth-averaged smooth particle hydrodynamics method [15,16]. References [17,18] quantified the impact loads over cylinders from a tsunami wave.
To protect the marine structures from potential damages due to extreme marine conditions, engineers have developed outer protections in terms of wrappers made of porous media. A porous medium enhances the buffering performance of the structures and dissipates part of the incoming wave energy [19]. For example, the turbulent intensities on a permeable breakwater were significantly attenuated in the numerical analysis by References [20,21,22]. Naturally, porous media are expected to be protective to submarine pipelines under extreme marine conditions, although thermal insulation and erosion prevention were mainly considered in designing pipeline coatings in the industry [23,24]. Reference [25] quantified the wave forces on pipelines buried in an impermeable bed with coverings of porous media. References [26,27] evaluated the protective performance of a porous polymer coating on subsea pipelines under sudden impacts. The drag reduction function of the porous coatings over cylinders were then quantified by Reference [28]. Two factors were considered to influence the stabilization effect of the porous coatings on pipelines: the production of an entrainment layer through the coating and the triggering of turbulent transition of the detaching shear layers. In engineering practice, applications of porous coatings on submarine pipelines are limited. Concrete wrappers, mainly designed to counteract the buoyancy forces of pipelines, can be considered as one kind of porous wrapper with medium permeability. In addition, porous wrappers made with woven carbon-fiber materials or polyurethane foam may be designed in future for pipeline protection.
The above literature review revealed that few studies were performed to examine the protective effect by the porous media on submarine pipelines, which is the main aim of this study. The porous wrapper and the submarine pipeline modules are simulated in a numerical wave tank (NWT) with the immersed boundary (IB) method. The numerical methods and equations will be provided in Section 2. Verification of the numerical model is provided in Section 3. The parametric simulations are in Section 4, in which the effects of different waves on various pipelines with porous wrappers are analysed. The conclusions are given in Section 5.

2. Numerical Methods

For simulating the interactions between pipelines and waves, the finite volume methods have been widely used. In this study, the commercial finite volume package FLOW-3D® (version 11.1.0; 2014; https://www.flow3d.com (accessed on 10 December 2022); Flow Science, Inc., Santa Fe, NM, USA). Flow-3D aims to solve the transient response of fluids under interactions with structures, internal and external loads and multi-physical processes. It features some advantages in terms of a high level of accuracy in solving the Navier-Stokes equation with the volume of fluid (VOF) method, efficient meshing techniques for complex geometries, and high efficiency level for large-scale problems. Also, Flow-3D provides the flexibility and utility for flowing through porous media. A two-dimensional numerical wave tank was constructed by using the immersed boundary (IB) method and an in-house subroutine termed as IFS_IB. A submarine pipeline and porous medium were two-way coupled at the interface described by the individual volume fractions [29]. The pipeline was wrapped with a layer of a porous medium. A solitary wave was generated at the inlet boundary of the tank to simulate an approaching tsunami. Non-slip wall conditions were assigned at the bottom of the tank and the pipe surface, which was also specified with a roughness coefficient. The top boundary was defined as a free boundary and configured with the atmospheric pressure. A Neumann-type absorbing boundary condition, a stable, local, and absorbing numerical boundary condition for discretized transport equations [30], was imposed on the outlet boundary to attenuate the reflections of the outgoing waves. A transition zone is set within a certain range from the boundary to reduce the horizontal gradient force of the elements near the boundary and suppress the calculation wave caused by this boundary condition. Through the relaxation coefficient, the predicted value on the inner boundary of the transition zone and the initial value on the outer boundary are continuously transitioned to achieve the purpose of reducing the reflection of propagating waves. The CUSTOMIZATION function of the software FLOW-3D was utilised to impose the Neumann-type absorbing boundary condition. The FLOW-3D distribution includes a variety of FORTRAN source subroutines that allow the user to customize FLOW-3D to meet their requirements. The FORTRAN subroutines provided allow the user to customize boundary conditions, include their own material property correlations, specify custom fluid forces (i.e., electromagnetic forces), add physical models to FLOW-3D, and have additional benefits. Several “dummy” variables have been provided in the input file namelists that users may use for custom options. A user definable namelist has also been provided for customization. Makefiles are provided for Linux and Windows distributions and Visual Studio solution files are provided for Windows distributions to allow users to recompile the FLOW-3D code with their customizations.

2.1. Governing Equations

The governing equations involved include the continuity equations and Reynolds-averaged Navier-Stokes equations. The mass and momentum are conserved in a two-dimensional zone [31]:

where U is the velocity vector, X is the Cartesian position vector, g denotes the gravitational acceleration vector, and ρ represents the weighted averaged density. The term μ is the viscosity. σκ∇α identifies the surface tension effects with σ as the surface tension and α as the fluid volume fraction. Each cell in the fluid domain has a water volume fraction (α) ranging between 0 and 1, where 1 represents cells that are fully occupied with water, while 0 represents cells that fully occupied with air. Values between 1 and 0 represent free surface between air and water. The free surface elevation is defined by using the volume of fluid (VOF) function:

where V_F is the volume of fluid fraction, F_SOR is the source function, F_DIF is the diffusion function; A_x, A_y, and A_z represent the fractional areas; and u, v, and w are the velocity components in the x, y, and z directions.

2.2. Porous Media Module

In FLOW-3D, the porous medium’s flow resistance is modelled by the inclusion of a drag term in the momentum equations (Equation (2)). Coarse granular material is used in most coastal engineering applications, in which case the Forchheimer model is suitable. Using this model, a drag term F_du_i is added to the righthand-side of Equation (2):

where |U| is the norm of the velocity vector, n the porosity, and a and b are the factors.

2.3. Solitary Wave Boundary

The solitary wave is generated in terms of variations of the surface elevation η and velocities u and v by following McCowan’s theory [32]:

where h is the still water depth; Q is the reference value

where X = x − c₀t; 𝑐0=√𝑔𝐻+ℎ; H is the wave height; and t is the elapsed time.

3. Validation

3.1. Propagation over a Porous Breakwater

An experimental test on the propagation process of a solitary wave over a permeable breakwater was performed by Reference [20], which was simulated in this study to validate the adopted two-way coupling model (Figure 1a). The length, width, and depth of the flume tank were 25, 0.5, and 0.6 m, respectively. A permeable breakwater was mounted at the bottom of the flume, which had dimensions of 13 cm and 6.5 cm in the length and height, respectively. The porous breakwater with an average porosity of 0.52 was configured by glass beads with a constant diameter of 1.5 cm. Two wave gauges were fixed before (WG1) and behind (WG2) the breakwater, respectively. The initial still water depth h was assumed to be 10.6 cm. Height of the solitary wave H was considered to be 4.77 cm. In the numerical model, the calculation zone had dimensions of 5 m in length and 0.25 m in height. The second order quadrilateral mesh elements were adopted. The grid around the breakwater was the finest of 0.001 m. The adopted time step size was 0.05 s. The numerical predictions of the water elevations at the locations WG1 and WG2 by the adopted numerical tool FLOW-3D are close to both the experimental measurements and the numerical predictions from another CFD FLUENT version 14.0.1 [33] (Figure 1). Figure 1b,c show the comparison of monitored water levels at the two water level monitoring points in Figure 1a. It can be seen that the experimental results of the two monitoring points are consistent with the numerical simulation results, indicating that the propagating solitary wave energy is basically completely dissipated and then flows out. If the propagating wave energy is not dissipated, the phenomenon of wave reflection will occur. The waves monitored at the two monitoring points will appear superposition of propagating waves and reflected waves. The numerical simulation results do not agree with the physical experiment results. The fluctuations of the water surface elevation after the bypass of the incoming wave are due to its residual reflection at the right absorbing boundary condition, which arrives at WG2 at an earlier time than WG1. Evolution of the wave surfaces was also compared between the experimental and the numerical models (Figure 2), which demonstrates that the numerical tool is sufficiently reliable. The velocity of the wave is reduced by the porous medium as it partially infiltrates into the breakwater, which is shown as in Figure 3 by comparing the horizontal velocity distributions between the experimental and numerical results at times of 1.5 s and 2 s. The numerical predictions of the flow velocities have slight discrepancies with the experimental measurements, which are attributed to the material assumptions made in the numerical model for the glass beads in the experimental setup.

Figure 1. The diagrammatic sketch of the numerical setup (non–scaled) (a) and the temporal evolution comparison of water surface between experimental and numerical results (b,c).

Figure 2. Water surface comparison between experimental and numerical results.

Figure 3. Comparison of horizontal velocity distribution between experimental and numerical results.

3.2. Forces on Pipeline

Another experimental test of a solitary wave impacting a pipeline was performed by Reference [34], which was also reproduced in this study for validation purposes. The calculation zone had dimensions of 40 m in length and 0.6 m in height. The solitary wave had a height of 0.0555 m with the initial water depth of 0.192 m. The pipe had a diameter of 0.048 m, which had a distance of 0.136 m over the bottom boundary of the model. A dense mesh consisting of 413,411 cells was employed with a mesh size of 0.1 mm around the pipe, which proved to be sufficiently fine through convergence studies. History of the horizontal and vertical forces, normalized by ρgL(πD²/4) with L as the unit length of 1 m, is compared between the experimental and numerical results (Figure 4). Both the peak values and the transient variations of the forces predicted by the numerical analysis converge to the measured values in the experimental test. The slight discrepancy between the numerical and experimental results at 2.5 s and 3.1 s, which may be induced by the error of the numerical model simulating the complicated turbulence behaviour, is acceptable in relation to the requirements of this study as our concern is mainly the peak values of forces.

Figure 4. Force comparison between the experimental and numerical results.

Therefore, the adopted numerical tool is sufficiently reliable to simulate the interactions between solitary waves and the permeable structure through the above validation cases.

4. Results and Discussion

Influence of the solitary waves on the performance of wrapped pipelines was investigated by considering different wave heights (H) and thicknesses (T) and wrapper porosities (n). The still water depth (h) was taken to be 6 m (Figure 5). The diameter of the porous medium was assumed to be 0.05 m. The pipeline diameter D was set at 1 m. In Figure 5 the variable G represents the gap between the permeable wrapper and the seabed. The scouring process had been completed before the simulation; therefore, the seabed boundary was taken as a rigid wall. The tandem pipelines had a distance of S between each other. The whole model had dimensions of 400 m in length and 12 m in height. The finest mesh around the pipeline was configured as 0.0025 m, which was verified to be sufficiently fine through trial calculations with finer meshes.

Figure 5. Layout for solitary wave impinging on the submarine pipeline encased in porous media.

4.1. Effect of Porous Wrapper

4.1.1. Wrapper Porosity

The pipeline was put on the seabed. The gap (G) between the wrapper and the seabed was considered to be zero. The height (H) of the solitary wave was considered to be 2 m. The porosity (n) was taken to be 0.0, 0.4, 0.6, and 1.0. Note that n = 0.0 indicates the impervious condition, while n = 1.0 corresponds to the non-wrapping condition. The thickness of the permeable wrapper remained at 0.5 m. In calculation, the wave approaches the pipe at around 6.3 s and departs from it at 10.2 s. When the wave approaches, the kinematic performance over the pipe is enhanced (Figure 6). Due to the wave disturbance, a number of small vortices are generated around the pipe (Figure 7). At the departure of the wave, the disturbance to the flow field seems to be more intense than that at its arrival, which further generates vortices around the pipeline. Without a wrapper, the pipe is fully exposed to the disturbance of the incoming wave, which maximises the velocity and vorticity values around the pipe. When the pipeline is wrapped by a porous medium, some water seeps into the wrapper, and the velocity in the wrapper is reduced to a very low value, which implies that the porous medium is capable of absorbing the dynamic energy of the flowing fluid. With an external coverage (n < 1.0), the disturbance is generated mainly at the outer surface of the wrapper. As the wrapper porosity increases, the domain of the low-speed flow underneath the pipeline expands.

Figure 6. The velocity contours of the flow fields under different porosities; (a) n = 0.0; (b) n = 0.4; (c) n = 0.6; (d) n = 1.0; left to right: arrival, departure.

Figure 7. The vorticity contours of the flow fields under different porosities; (a) n = 0.0; (b) n = 0.4; (c) n = 0.6; (d) n = 1.0; left to right: arrival, departure.

The peak velocity around the pipeline without a wrapper (1.9 m/s) is larger than that with a wrapper (1.6 m/s) (Figure 8). For pipes with wrappers, the peak velocities around them are similar to one another. In contrast, the velocity profiles at x = 23 m are quite different. When the pipeline has no wrapper (i.e., n = 1.0), the change in velocity is fairly moderate. When the pipeline has a wrapper, the porous wrapper causes a secondary fluctuation in the rear water body after the primary fluctuation due to the peak of the wave passing through the pipeline. This generates a series of velocity peaks. The secondary velocity peaks for a porosity coefficient of 0.4 are higher than those for a porosity coefficient of 0.6. Accordingly, the turbulent kinetic energy (TKE) also changes with the porosities, as shown in Figure 9. The TKE is expressed as

Figure 8. Comparison of horizontal and vertical velocities at front and rear of pipeline under different porosities.

Figure 9. Comparison of turbulent kinetic energy at front and rear of wrapper under different porosities.

With the propagation of the wave, the TKE increases gradually in front of the pipeline. The TKE value under the pipeline without a wrapper (n = 1.0) (0.0008 kJ) is nearly half of that with a wrapper (0.0015 kJ). In comparison, the TKE values for the wrapped pipelines (n < 1.0) are very close to each other. After the wave leaves the pipeline, the TKE in front of the pipeline decreases for around 50%. Then, the TKE in the rear of the pipeline with a porous wrapper increases intensively because the porous media perturb the flow field. Compared with the pipeline without the wrapper, the interaction between the wrapped pipeline with the flow field is more severe. Furthermore, the solid wrapper can cause a strong disturbance to the flow, but the interference of the solid wrapper (n = 0.0) in the rear flow is still weaker than the wrapper with the porosity of 0.4.
The hydrodynamic forces (F), including the pressure and shear stress, are normalized by ρgL(πD2/4) (Figure 10). With a fully solid (i.e., n = 0.0) wrapper, the pipeline tends to be unaffected by the external flow. Hence, the hydrodynamic forces are zero while the forces on the wrapper reach their maximum. With porous wrappers, water seeps into the wrapper, buffering the impact of the incoming waves on the pipe. As the porosity coefficient increases, the induced forces on the pipeline increase while those on the wrapper decrease. When the porosity coefficient is 0.4, the forces on the external wrapper become higher than that on the internal pipeline. Therefore, the porous wrapper is capable of protecting the pipeline. The smaller the porosity coefficient the better protection the wrapper provides to the pipeline. The pressure gradient and shear stress forces are also shown in Figure 11.

Figure 10. Comparisons of the maximum hydrodynamic forces on the pipeline and wrapper.

Figure 11. Decomposed pressure gradient force (a) and shear stress (b) force on the pipeline.

4.1.2. Thickness of Wrapper

Seven wrapper thicknesses are considered: T = 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, and 0.5 m. The porosity coefficient is taken to be 0.6. At the moment that the wave goes through the pipe, the transient evolution of the vorticity contours around the pipeline with a wrapper thickness of 0.25 m is depicted in Figure 12. A couple of vortices emerge on the upper and lower vertices of the pipeline as the wave approaches the pipeline. As the wave propagates, many vortices flow along the wrapper and then shed off. Compared with the top vortices, the bottom vortices are shed off faster for two reasons. Firstly, as the friction at the seabed is small, the bottom flow velocity is higher than that on the top. Secondly, when the wave peak departs from the pipeline, a strong disturbance by the water body occurs behind the pipeline, followed by the irregular swing and fall off of the vortices. After the wave travels far away, the water flow near the pipeline becomes weak, and the vortices are scattered around the pipeline.

Figure 12. Temporal evolutions of vorticity contours around pipeline with wrapper thickness of 0.25 m (a) 6.0 s (b) 6.6 s (c) 7.2 s (d) 7.8 s (e) 8.1 s (f) 8.7 s (g) 9.0 s (h) 10.2 s (i) 12.6 s.

Figure 13 shows a comparison of flow field stream traces and the velocity contours. When the fluid penetrates the wrapper, the streamline starts to diverge, which indicates that the free flow is hindered. Therefore, the flow becomes slower and the flow direction becomes non-uniform. For the fluid flows out of the wrapper, the stream traces are quite complex and chaotic. The reason is that the seeping fluid mixes with the bypass flows and causes strong interference in the water body behind the pipeline. The streamlines passing through the wrapper indicates frequent water exchange at the wrapper surface. Along with the small-attached vortices on the wrapper surface, more fluid passes over the wrapper and causes a large vortex behind the wrapper.

Figure 13. Comparisons of flow field streamtraces and velocity contours under different wrapper thicknesses; (a) T = 0.2 m; (b) T = 0.3 m; (c) T = 0.4 m; (d) T = 0.5 m.

The highest free surface elevations and velocities at the front and at the rear of the pipeline with different wrapper thicknesses are depicted in Figure 14. As the wrapper thickness increases, the highest elevation at the front of the pipeline seems to be quite stable, although the peak velocity increases by around 6%. At the moment that the wave bypasses the pipeline, the maximum elevation reduces with an increase in the wrapper thickness. This is because the pipeline blocks the wave propagation. However, due to the strong mixing effect of the seepage and bypass water, the maximum velocity rises to be higher than that in front of the pipe. The maximum forces on the wrapper and the pipeline for different wrapper thicknesses are shown in Figure 15. With an increase in the wrapper thickness from 0.2 to 0.5 m, the normalized forces on the wrapper are doubled as a larger interaction area is involved. In contrast, the vertical forces on the pipeline decrease by 12.5%. Therefore, the larger the thickness of the wrapper the safer the pipeline.

Figure 14. Comparisons of the maximum elevations and velocities in front and rear of the pipelines with different wrapper thicknesses; (a) free surface elevation (note: original water depth is 6 m); (b) velocity.

Figure 15. Hydrodynamic forces on the pipeline and wrapper.

4.2. Effect of Pipeline Structure

The in-situ pipelines may be under various suspended conditions since the seabed topography is often uneven. Some pipelines are also laid in tandem for the sake of the transportation efficiency. In order to examine the effects of porous wrappers on pipelines under different conditions, a study was carried out considering two scenarios, namely, suspended pipelines and pipelines in tandem. In the numerical models, the porosity coefficient (n) remained at 0.6, the thickness (T) of the wrapper was kept at 0.5 m, and the wave height (H) was assumed to be 2.0 m.

4.2.1. Suspended Pipelines

Six gaps (G) between the wrapper and the seabed (0.0, 0.2, 0.4, 0.6, 0.8, and 1.0 m) were considered [35,36,37]. The representative flow field at three points in time (6.3, 7.2, and 10.2 s) are shown in Figure 16. At the arrival of the wave at the pipeline (at 6.3 s), the flow is accelerated and the velocities over and beneath the pipe reach the maximum values due to the bypass effect of the fluid. At the moment that the wave peak is above the pipe (at 7.2 s), all the velocities around the pipe reach their highest values. After the wave passes over the pipe (at 10.2 s), the velocity decreases and several vortices are formed behind the pipeline. With a tiny wrapper-seabed gap, the velocity within the gap is very high while the flux is relatively small. An increase in the gap will result in an increase in the flux and a decrease in the velocity. A symmetric velocity distribution similar to a fisheye is observed behind the pipeline, which becomes more obvious when the gap increases (Figure 16c).

Figure 16. The velocity contours of the flow fields under different gaps; (a) G = 0.2 m; (b) G = 0.6 m; (c) G = 1.0 m. Left to right: 6.3 s, 7.2 s, and 10.2 s. Left to right: arrival, stay, departure.

With the bypass of the wave, the vortices generated around the pipeline become larger. The vorticity contours and the streamlines of the flow field are shown in Figure 17. As the solitary wave approaches, a pair of whirlpools shed off from the wrapper with a gap of 0.2 m. With an increase in the gap, the two whirlpools gradually disappear and are replaced with two smaller vortices. Due to the internal pores within the wrapper, the streamlines in the wrapper are dispersed, and it is hard for a vortex to be generated. With an increase in the gap, two anti-symmetric vortices shed off from the wrapper. Besides, some tiny vortices remain adhered to the wrapper due to the interaction by the seepage and the external flow. When the gap is very small, a few small vortices are generated between the wrapper and the seabed. In contrast to the interface of vortex from the flow around a solid cylinder, the vortex interface at the wrapper is not fully smooth. Because of the strong interactions of fluid over the wrapper surface, several small vortices mingle with the large shedding vortices. The flow direction also varies greatly according to the streamline mobilisation.

Figure 17. The vorticity contours of the flow fields under different gaps; (a) G = 0.2 m; (b) G = 0.6 m; (c) G = 1.0 m. Left to right: 6.3 s, 7.2 s and 10.2 s.

The gap is normalized by the pipeline diameter as β = G/D. With a small gap (β < 0.2), the horizontal forces on both the wrapper and the pipeline are slightly smaller than those on the wrapper and pipeline without a gap (Figure 18). With a further rise of the gap width, the horizontal forces are accordingly enlarged due to higher velocity around the pipeline as shown in Figure 16.

Figure 18. Comparisons of the maximum horizontal and vertical hydrodynamic forces on the pipeline and wrapper under different gaps.

In contrast, an increase in the gap width may inversely cause the reduction of vertical forces on both the wrapper and the pipeline. The vertical forces can be considered to consist of two parts. One is caused by the weight of the water body at the bypass of the wave from the pipeline, while the other can be caused by the velocity difference between the flow above and below the pipeline after the flow passes over. In summary, as the gap increases, the flow velocity within the gap initially increases when β < 0.2 and then decreases when β > 0.2. In contrast, the vertical forces caused by the wave’s weight always decrease with an increase in the gap.

4.2.2. Pipelines in Tandem

The hydrodynamic forces on pipelines in tandem are investigated considering five different distances (S) between the two pipeline centres (2.5, 3.0, 3.5, 4.0, and 4.5 m). The velocity and vorticity fields at 6.3, 7.2, and 10.2 s around the tandem pipelines with distances of 2.5, 3.5, and 4.5 m are depicted (Figure 19 and Figure 20). As the wave approaches the pipeline, the velocity within the pipeline gap is very small due to the blockage effect of the pipeline in front. As the distance increases, the velocity field within gap space is enhanced as more water flow is allowed. The velocity above the pipeline has its maximum value, and part of the high-speed fluid flows into the gap through the space underneath the pipeline. With a small distance, the vortices shedding off from the front pipeline impinge directly on the rear pipeline without any stretching. When the distance is increased, noticeable vortex shedding emerges in the middle space (Figure 20c). Similar vortex shedding behind the rear pipeline is observed for different distances. After the wave bypasses the pipeline, the increase in the distance between the pipelines will cause an increase in the velocity magnitudes in the space among the pipelines. As the distance increases, the flow becomes more chaotic due to the seepage from the wrapper and the limited flow space. In summary, influence of the distance between the pipelines over the whole kinematic field is not significant, although the local flow field around the pipelines is severely affected. When the wave bypasses the tandem pipelines, the largest forces on structures (i.e., the pipelines and wrappers) are shown in Figure 21, in which the distance ratio (θ) is calculated as θ = S/D.

Figure 19. The velocity contours of the flow fields under different spacings; (a) S = 2.5 m; (b) S = 3.5 m; (c) S = 4.5 m. Left to right: 6.3 s, 7.2 s and 10.2 s.

Figure 20. The vorticity contours of the flow fields under different porosities; (a) S = 2.5 m; (b) S = 3.5 m; (c) S = 4.5 m. Left to right: 6.3 s, 7.2 s and 10.2 s.

Figure 21. The maximum forces on the pipeline and wrapper under different distances.

As for the pipeline in front, as the distance ratio increases, the horizontal loads on the wrapper and pipeline increase slightly, while the vertical forces are almost doubled. As for the rear pipeline, as the distance reduces, the velocity in the gap becomes smaller and the forces on the pipelines and wrappers are also reduced, which is mainly attributed to the shield effect from the front pipeline. With an increase in the distance, the forces increase due to the increase in the turbulence energy in the gap.
Different ratios of the forces on the front and rear pipelines are depicted in Figure 22. The difference ratio is defined as ΔFn = (ff,max−fr,max)/ff,max, where ff,max and fr,max are the maximum forces on the pipeline or wrapper. It is found that the horizontal loads on the rear pipe and wrapper tend to be always higher than their counterparts on the front pipe. This means that a turbulent flow in the horizontal direction on rear pipe is more intense than that on the front pipe. For different distances, deviations for the forces on the pipelines and wrappers are also different. The deviation is found to be maximized at a distance of 1 m and this indicates that the pipeline is not well protected and needs to be avoided in engineering practice.

Figure 22. The deviation of the forces on the front and rear pipelines and wrappers under different distances.

4.3. Effect of Wave Height

Six groups of wave heights (H), i.e., 1.6, 1.8, 2.0, 2.2, 2.4, and 2.6 m, are selected to consider different marine environment. After bypassing the pipeline, the height of the wave decreases because of the blockage effect of the pipeline and the dissipation of the flow energy (Figure 23a). The deviation ratio of the wave heights before and after the wave passes over the pipeline is shown in Figure 23b and is defined as δ = (Hf,max − Hr,max)/Hf,max. The wave height attenuation becomes more significant as the wave height increases. This means that waves with larger heights are more easily affected by the pipelines.

Figure 23. Waves with different wave heights; (a) temporal evolutions; (b) attenuation deviation.

At the bypass of the wave through the pipe, the loads are increased until they reach the maximum values at the moment that the wave peak appears above pipeline (Figure 24). The forces gradually decrease as the wave propagates. Because of some reflux after the wave bypasses the pipeline, the flow is in the opposite direction to that of the wave propagation, resulting in a negative force. The vibration of the water body by the wave propagation induces oscillations of the forces on the pipeline and wrapper. When the wave height is larger, the force oscillation becomes fiercer and the maximum loads on the pipeline and the wrapper increase (Figure 25). The vertical forces on the pipeline are the largest compared with other forces under the same conditions. Besides, as the wave height increases, the increased amplitude of vertical forces on the pipeline is the most significant change since the weight of the water above the pipeline increases. Therefore, given that the wave height is very high, the protective function of wrapper on the pipeline tends to be weakened compared with that of the wrapper for a low wave height.

Figure 24. The temporal evolutions of forces on the pipeline and wrapper; (a) Horizontal maximum force on pipeline; (b) Vertical maximum force on pipeline; (c) Horizontal maximum force on wrapper; (d) Vertical maximum force on wrapper.

Figure 25. Variation of hydrodynamic forces on the pipeline and wrapper under different distances.

5. Conclusions

The effect of porous media on the dynamic performance of submarine pipelines under solitary waves was investigated. The porosity of the wrapper, the seabed topography, the structure of the pipeline, and the marine environment were considered. The study had a limitation of the model sizes due to the limited computational resource and the simplification of the solitary wave due to its mathematical complication, which will be tackled in future works. The following main conclusions have been made.

(1) When a pipe is wrapped by a porous medium, the velocity in the wrapper is relatively small because the porous medium can consume the water energy and weaken the flow. With an increase in the porosity, the range of the low-speed flow at the bottom of the pipeline expands. This indicates that the porous wrapper can slow down the flow and affect a wider region of the surrounding water. After the bypass of the wave through the pipe, the number and volume of the vortices behind the porous wrapper are larger than those for a pipeline with a solid wrapper or without a wrapper. As the porosity coefficient increases, the impact forces on the pipe increase, while those on the wrapper decrease. This implies that the porous wrapper is capable of protecting the pipeline.

With an increase in the wrapper’s thickness, the hydrodynamic forces on the wrapper tend to increase. In particular, the horizontal forces on the pipeline decrease with an increase in the thickness due to the protection of the wrapper, while the vertical forces are increased because of variations in the fluid’s stagnation point.

(2) For a wave bypassing a pipe with different heights, a symmetric speed change similar to a fisheye appears behind the pipeline, along with two antisymmetric vortices shedding off from the wrapper.

As the internal seepage interacts with the external fluid flow, several small vortices are still attached to the wrapper. The hydrodynamic vertical forces on both the wrapper and the pipeline decrease with the pipeline distance. With an increase in the suspension of the pipe, the velocity and TKE within the gap space increase and both the vortex intensity and the number of vortices increase. Therefore, the flow pattern appears to be chaotic. As for the front pipeline, an increase in the gap leads to a slight increase in the horizontal forces on both the wrapper and the pipeline, but a significant increase in the vertical forces. As for the rear pipeline, because of the shield function of the front pipeline, the velocity within the gap space and the forces on the pipes decrease with a decrease in the gap size.

(3)When the waves with different heights pass over the pipeline, the height of the wave is reduced because of the blockage function from the pipeline and the dissipation characteristic of the flow energy. When the wave height is increased, the velocity around the pipeline increases, inducing an increase in the TKE. As the wave height increases, all the maximum forces on the pipeline and wrapper also increase. Note that an increase in the vertical forces on the pipeline is the most significant change because the weight of the water above the pipeline increases, which implies that the protection function of the wrapper is enhanced by the reduction in the wave height.

From the above investigation, the mechanism of load transfer from the pipeline to the external wrapper has been presented. This encourages industrial experts and academic scholars to arrange more investigations of the functions and cost-efficiency of porous wrappers, which could form a new branch of the pipeline design practice.

Author Contributions

Contributor Roles Taxonomy: E.Z.: Conceptualization, Methodology, Validation, Investigation and Writing—Original Draft; Y.D.: Data Curation, Formal analysis; Y.D.: Visualization, Project administration; Y.D., L.C., Y.W. and Y.L.: Writing—Review & Editing. All authors have read and agreed to the published version of the manuscript.

Funding

The paper was supported by the National Natural Science Foundations of China (Grants No. 52001286 and No. 42272328), GuangDong Basic and Applied Basic Research Foundation (Grant No. 2022A1515240002) and Comprehensive Survey of Natural Resources in Huizhou-Shanwei Coastal Zone (Grant No. DD20230415).

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflict of interest.

References

1. Jeanjean, P.; Liedtke, E.; Clukey, E.C.; Hampson, K.; Evans, T. An operator’s perspective on offshore risk assessment and geotechnical design in geohazard-prone areas. In Proceedings of the International Symposium on Frontiers in Offshore Geotechnics, Perth, Australia, 19–21 September 2005; CRC Press: London, UK, 2005; pp. 115–143.
2. Fuhrman, D.R.; Baykal, C.; Sumer, B.M.; Jacobsen, N.; Fredsøe, J. Numerical simulation of wave-induced scour and backfilling processes beneath submarine pipelines. Coast. Eng. 2014, 94, 10–22.
3. Zhao, M.; Cheng, L. Numerical modeling of local scour below a piggyback pipeline in currents. J. Hydraul. Eng. 2008, 134, 1452–1463.
4. Eklund, T.; Høgmoen, K.; Paulsen, G. Ormen Lange Pipelines Installation and Seabed Preparation. In Proceedings of the Offshore Technology Conference, Houston, TX, USA, 30 April–3 May 2007. OTC-18967-MS.
5. Dong, Y.; Wang, D.; Randolph, M. Investigation of impact forces on pipeline by submarine landslide with material point method. Ocean Eng. 2017, 146, 21–28.
6. Dong, Y.; Wang, D.; Randolph, M. Quantification of impact forces on fixed mudmats from submarine landslides using the material point method. Appl. Ocean Res. 2020, 102, 102227.
7. Dong, Y.; Liao, Z.; Liu, Q.; Cui, L. Potential failure patterns of a large landslide complex in the Three Gorges Reservoir area. Bull. Eng. Geol. Environ. 2023, 82, 41–52.
8. Zhao, E.J.; Dong, Y.; Tang, Y.Z.; Sun, J.K. Numerical investigation of hydrodynamic characteristics and local scour mechanism around submarine pipelines under joint effect of solitary waves and currents. Ocean Eng. 2021, 222, 108553.
9. Sun, Q.L.; Wang, Q.; Shi, F.Y.; Alves, T.; Gao, S.; Xie, X.N.; Wu, S.G.; Li, J.B. Runup of landslide-generated tsunamis controlled by paleogeography and sea-level change. Commun. Earth Environ. 2022, 3, 244.
10. Rabinovich, A.B.; Titov, V.V.; Moore, C.W.; Eblé, M.C. The 2004 Sumatra tsunami in the southeastern Pacific Ocean: New global insight from observations and modeling. J. Geophys. Res. Ocean. 2017, 122, 7992–8019.
11. Madsen, P.A.; Fuhrman, D.R. Run-up of tsunamis and long waves in terms of surf-similarity. Coast. Eng. 2008, 55, 209–223.
12. Madsen, P.A.; Fuhrman, D.R.; Schäffer, H.A. On the solitary wave paradigm for tsunamis. J. Geophys. Res. Ocean. 2008, 113, C12012.
13. Fan, N.; Jiang, J.; Dong, Y.; Guo, L.; Song, L. Approach for evaluating instantaneous impact forces during submarine slide-pipeline interaction considering the inertial action. Ocean Eng. 2022, 245, 110466.
14. Dong, Y.; Cui, L.; Zhang, X. Multiple-GPU for three dimensional MPM based on single-root complex. Int. J. Numer. Methods Eng. 2022, 123, 1481–1504.
15. Xie, P.; Chu, V.H. The forces of tsunami waves on a vertical wall and on a structure of finite width. Coast. Eng. 2019, 149, 65–80.
16. Smith, L.; Kolaas, J.; Jensen, A.; Sveen, K. X-ray measurements of plunging breaking solitary waves. Eur. J. Mech. B Fluids 2019, 73, 112–121.
17. Lin, M.Y.; Liao, G.Z. Vortex shedding around a near-wall circular cylinder induced by a solitary wave. J. Fluids Struct. 2015, 58, 127–151.
18. Aristodemo, F.; Tripepi, G.; Meringolo, D.D.; Veltri, P. Solitary wave-induced forces on horizontal circular cylinders: Laboratory experiments and SPH simulations. Coast. Eng. 2017, 129, 17–35.
19. Van Gent, M.R.A. Wave Interaction with permeable coastal structures. Int. J. Rock Mech. Min. Sci. Geomech. 1996, 33, 227A.
20. Wu, Y.T.; Hsiao, S.C. Propagation of solitary waves over a submerged permeable breakwater. Coast. Eng. 2013, 81, 1–18.
21. Qu, K.; Sun, W.Y.; Deng, B.; Kraatz, S.; Jiang, C.B.; Chen, J.; Wu, Z.Y. Numerical investigation of breaking solitary wave runup on permeable sloped beach using a nonhydrostatic model. Ocean Eng. 2019, 194, 106625.
22. Jiménez, J.; Uhlmann, M.; Pinelli, A.; Kawahara, G. Turbulent shear flow over active and passive porous surfaces. J. Fluid Mech. 2001, 442, 89–117.
23. Liu, W.; Li, X.; Hu, J. Research on flow assurance of deepwater submarine natural gas pipelines: Hydrate prediction and prevention. J. Loss Prev. Process Ind. 2019, 61, 130–146.
24. Bruneau, C.; Mortazavi, I. Passive control of the flow around a square cylinder using porous media. Int. J. Numer. Methods Fluids 2004, 46, 415–433.
25. Akbari, H.; Pooyarad, A. Wave force on protected submarine pipelines over porous and impermeable beds using SPH numerical model. Appl. Ocean Res. 2020, 98, 102118.
26. Vestrum, O.; Kristoffersen, M.; Polanco-Loria, M.A.; Ilstad, H.; Langseth, M.; Børvik, T. Quasi-static and dynamic indentation of offshore pipelines with and without multi-layer polymeric coating. Mar. Struct. 2018, 62, 60–76.
27. Vestrum, O.; Langseth, M.; Børvik, T. Finite element analysis of porous polymer coated pipelines subjected to impact. Int. J. Impact Eng. 2021, 152, 103825.
28. Klausmann, K.; Ruck, B. Drag reduction of circular cylinders by porous coating on the leeward side. J. Fluid Mech. 2017, 813, 382–411.
29. Hirt, C.W.; Sicilian, J.M. A porosity technique for the definition of obstacles in rectangular cell meshes. In Proceedings of the 4th International Conference on Numerical Ship Hydrodynamics, Washington, DC, USA, 24–27 September 1985.
30. Coulombel, J.F.; Lagoutière, F. The Neumann numerical boundary condition for transport equations. arXiv 2018, arXiv:1811.02229.
31. Ding, D.; Ouahsine, A.; Xiao, W.; Du, P. CFD/DEM coupled approach for the stability of caisson-type breakwater subjected to violent wave impact. Ocean Eng. 2021, 223, 108651.
32. Wu, T. A Numerical Study of Three Dimensional Breaking Waves and Turbulence Effects. Ph.D. Thesis, Cornell University, Ithaca, NY, USA, 2004.
33. ANSYS. ANSYS FLUENT 14.0 Theory Guide. v.14.0.1; ANSYS, Inc.: Canonsburg, PA, USA, 2011.
34. Sibley, P.O. The Solitary Wave and the Forces It Imposes on a Submerged Horizontal Circular Cylinder: An Analytical and Experimental Study. Ph.D. Thesis, City University London, London, UK, 1991.
35. Zhao, E.J.; Shi, B.; Qu, K.; Dong, W.B.; Zhang, J. Experimental and numerical investigation of local scour around submarine piggyback pipeline under steady current. J. Ocean Univ. China 2018, 17, 244–256.
36. Zhao, E.J.; Qu, K.; Mu, L. Numerical study of morphological response of the sandy bed after tsunami-like wave overtopping an impermeable seawall. Ocean Eng. 2019, 186, 106076.
37. Zhao, E.J.; Sun, J.K.; Tang, Y.Z.; Mu, L.; Jiang, H.Y. Numerical investigation of tsunami wave impacts on different coastal bridge decks using immersed boundary method. Ocean Eng. 2020, 201, 107132.