*3.1. Hydrodynamics Calibration Test*

The hydrodynamic calibration process for the numerical model REEF3D is as follows. Figure 3 shows the results obtained by numerical modeling and those reported by Umeyama [3], based on experimental data, for a domain forced by waves and a mixed domain of wave and currents. The comparison variable corresponds to the instant surface elevation (η) nondimensionalized with the wave height (*H*) as a nondimensional time function ( *tT* ), with *T* as the period. The physical sense of the variable η*H* corresponds to the fraction of the increase or decrease in the water surface and, evidently, the maximum and minimum fractions are the descriptors of wave asymmetry. The variable *t T*corresponds to an indicator of the wave time fraction being simulated.

A comparison between the experimental data for cases where waves are acting alone (W1 to W3, illustrated with circles) and those for waves and currents (WC1 to WC3, illustrated with triangles) are shown in Figure 3, boxes 1.1 to 1.3. Here, it can be deduced that a codirectional current modifies the wavelength and amplitude of the wave. The numerical results are found in Figure 3, boxes 2.1 to 2.3, both for a flow with waves (blue line) and for waves and currents (red line). Thus, it can be verified that the wavelength and amplitude are modified between both flows, as shown in the experimental data (Figure 3, box 1.1 to 1.3). The results of the comparison between the numerical modeling REEF3D and the experimental data provided by Umeyama [3] are included in boxes 3.1 to 3.3 for waves acting alone, while waves and currents acting together are shown in boxes 4.1 to 4.2.

When analyzing numerical and experimental data associated to waves acting alone (W1 to W3), all model cases were able to adequately reproduce the maximum and minimum wave amplitude, as well as its temporal evolution within a wave time coinciding with the necessary time to reach the peak, the zero crossing time and the time to reach the minimum. Taking into account the difference in the estimation of the crest and trough, for the waves and currents (WC1 to WC3), it can be observed that the numerical model slightly underestimates the experimental data, but at a magnitude of less than one millimeter (around 2% error).

The vertical profiles for the instant velocity associated with each of the simulated cases in the present investigation (which were experimentally registered by Umeyama [3]) are compared in Figure 4 for waves acting alone and in Figure 5 for waves and currents combined. Both figures illustrate four instants of time for each simulated case, ordered from left to right and corresponding to *t*/*T* = 0.00, 0.25, 0.50, 0.75, and 1.00.

**Figure 3.** Phase-average surface displacements comparison between experimental data (Umeyama [3]) and numerical results.

**Figure 4.** Instantaneous horizontal velocity profile comparison between the experimental data (Umeyama [3]) and numerical results for different time steps, for cases with waves alone.

**Figure 5.** Instantaneous horizontal velocity profile comparison between the experimental data (Umeyama [3]) and numerical results for different time steps, both for waves alone and for current cases.

The results of the comparison of numerical and experimental data considered in the presence of waves alone (Figure 4) reflect a high consistency between vertical and temporal behavior, as the model is capable of adequately representing the flood direction (negative *u*/*C*) and the ebb direction (positive *u*/*C*). The experimental and numerical data show equivalent vertical structures for the velocity profile, with a slight increase toward the surface, which is evidenced in a greater proportion when analyzing the case with the highest wave (W3).

Near the bed, the current magnitudes determined based on the numerical model coincide with those determined by experimental means, showing slight differences between the simulated and instrumental data.

The detected differences for *t*/*T* = 0.25 and 0.75 show a low magnitude. However, the flood and ebb conditions reached different magnitudes described, as follows. For the W1 case, the maximum difference in the nondimensional velocity (*u*/*C*) obtained by the ebb direction was 0.004 (*t*/*T* = 1.00), while for the flood direction it was −0.001. For the W2 case, the differences fluctuated between −0.005 and −0.011.

Greater differences between the numerical model and experimental data were found for the W3 case for *t*/*T* = 0.00, which is produced at the water surface. Meanwhile, near the bed, the greatest difference in *u*/*C* was 0.015 for *t*/*T* = 1.00.

In Figure 4, the maximum and minimum velocities are not at 0.25 and 0.75 *t*/*T* respectively, because the wave phase effect on the initial condition was adjusted to represent the same oscillatory flow that Umeyama [3] reported in his research.

When incorporating a codirectional current to waves, Figure 5 shows that the vertical velocity distribution along the channel only shows the flood direction because the currents control the hydrodynamics. This can be confirmed by calculating the flow relative velocity proposed by Sumer and Fredsøe [26] (*Ucw*), which offers results equal to 0.84, 0.70, and 0.60, for WC1, WC2, and WC3, respectively.

In general terms, the numerical model adequately captured the behavior of the velocity profile, showing a greater similarity between the currents near the bed and those obtained toward the free surface. Compared to the hydrodynamic scenarios, where only waves were present, the differences found for the dimensionless velocity (*u*/*C*) have a greater magnitude when the current is incorporated in the centre, reaching 0.010 for WC1, −0.017 for WC2, and −0.037 for WC3.

The previous results correspond to the scenarios of interactions between waves and currents without the incorporation of a circular pile blocking the flow. However, they allow us to confirm that the numerical model is capable of representing the complex hydrodynamics resulting from the combined action of both components (waves and currents). To strengthen this analysis, the results of the model for a circular pile in the flow are shown next.

The results of the comparison of the numerical model with the experimental data obtained from Qi and Gao [26] are presented in Figure 6, which considers the total flow velocity (*Uc* + *Uw*) non-dimensionalized with the characteristic velocity (*C*). Based on this comparison, it was observed that the numerical model was able to represent the dynamic behavior of the combined flow velocity, for the different characteristics of the simulated waves and currents. In test case C01, it can be observed that the numerical model predicted slightly higher velocities in the trough located between 1.00 < *t*/*T* < 1.50, (*Uc* + *Uw*)/*C* (equal to 0.02). This result, however, was not observed in cases C02 and C03. Hence this case does not correspond to the numerical model configuration.

It is important to emphasize that the total flow velocity measurements experimentally obtained by Qi and Gao [26] were registered upstream from the pile at a distance of 20*D*. Ergo, the effects of the opaque structure on the velocities field would not be shown in their behavior. Therefore, this comparison (Figure 6) complements that previously shown in Figure 3 for the instant surface elevation. Thus, the numerical model is capable of representing the wave and current interactions in a freestream.

**Figure 6.** Total flow velocity comparison between the experimental data (Qi and Gao [26]) and numerical results for waves alone and for current cases in the presence of a cylindrical pile.

The effects on the velocity fields caused by the cylindrical pile in a flow field with waves perpendicular to the currents obtained from the numerical model were compared with data provided by Miles et al. [51]. The results for which waves and currents come from perpendicular directions are shown in Figure 7 (associates with case C04). When analyzing monitoring station P1, a high consistency is observed between the velocity profiles modeled (red line) and experimentally obtained data (circles), highlighting that the model is capable of representing the mean velocity near the bed and in the vertical direction as well. Equivalent performance was also verified for P2 and P8, which corresponds to the monitoring stations exposed to the current direction.

**Figure 7.** Mean velocity profiles comparison between the experimental data (Miles et al. [51]) and numerical results associated with case C04, for waves perpendicular to the current cases in presence of a cylindrical pile.

In station P4, it is noted that the numerical model may have slight differences in its vertical velocity distributions for elevation range, here *z*/*h* is equal to −0.98 to −0.96, which is not clearly presented in the other velocity monitoring stations. This performance detected in station P4 may be caused by the combined wake effect due to currents and waves, since geometrically this location is completely downstream of both forcings (associates with case C04).

#### *3.2. Hydrodynamics Behaviour of the Flow Around a Cylindrycal Pile*

The results obtained for the spatio–temporal evolution of the velocity and vorticity fields for case E07 are illustrated in Figure 8. These results consider four instants in time that allow us to visualize the interaction between the flow and the pile. All charted variables have been nondimensionalized by test characteristic scales. For example, the longitudinal and the vertical axes have been made dimensionless with the pile diameter, and, the beginning of the coordinated system has been placed at the center of the pile. The velocities (*u* and *w*) have been made dimensionless with the characteristic velocity (*C*). Vorticity (<sup>ω</sup>*x*, <sup>ω</sup>*y*, and <sup>ω</sup>*z*), on the other hand, was nondimensionalized by the integral temporal scale of the experiment and corresponds to the proportion *D*/*C*.

In Figure 8, in the first instant of time (Time 1 in Figure 8), it can be observed that the free surface is on a trough nearby the pile where the longitudinal velocity is mainly an ebb type, which usually form a horseshoe vortex when interacting with the incoming boundary layer. Descending velocities are developed in the vertical axis and near the pile, which may be associated with the down-flow.

As time progresses and the crest approaches the pile (Time 2 in Figure 8), the longitudinal velocity starts to diminish its magnitude (nearing zero), thus reflecting the oscillating effect, since this hydrodynamic behavior is an indicator that both ebb and flow can be present based on the phase of the wave. Contrary to the previous instant, the vertical velocity near the pile shows a positive magnitude, which would not produce down-flow and would consequently modify the horseshoe vortex behavior, as follows.

At time instant 3 (time 3 in Figure 8), the wave crest interacts directly with the pile, developing longitudinal velocities mainly oriented downstream the pile, while the vertical component of the velocity is once again oriented toward the bed, thus producing a down-flow. This is maintained for as long as the wave interacts with the pile during the crest phase (time 4 in Figure 8), while the longitudinal velocity changes to an ebb direction.

The association of the main component of the vorticity with the cross vorticity (<sup>ω</sup>*yD*/*C*) is obtained from the vorticity development near the pile and the adjacent bed for the four illustrated instants of time in Figure 8. Thus, that the minimum intensities occur under the trough, while in the crest there is a general tendency to maximize these intensities. When the crest is approaching the pile, it is observed that in the bed (upstream), a cross vorticity structure with a positive magnitude approaches the pile, and, conversely, the cross vorticity considerably diminishes its magnitude when the trough moves through the pile. The geometrical characteristics of this cross vorticity structure mainly present a greater longitudinal than vertical development.

The vorticity behavior upstream from the pile (described above) corresponds to the horseshoe vortex, and its intermittence would be conditioned by direction changes of the vertical and longitudinal velocity.

The downstream sector under the pile showed a vorticity behavior (<sup>ω</sup>*xD*/*C* and <sup>ω</sup>*zD*/*C*) with alternate structures (positives and negatives) and a greater development in the vertical axis than in the horizontal axis, compared to the cross vorticity. This spatio–temporal development could be associated with vortex shedding, resulting from structure-fluid interactions.

Based on a general analysis of the information obtained from the numerical modeling (a visual inspection of the results), there are no significant differences in the mean spatio–temporal behavior of the streamlines, vorticities, and velocity field; hence, these element can be broadly described by the specification of one of the eight case simulations, with scenario E01 selected for this effect.

**Figure 8.** Temporal behavior for velocities and vorticities near to the pile, waves and current are coming from the left to the right.

The characteristic results of scenario E01 are presented in Figure 9. Unlike Figure 8, which illustrates the behavior from the free surface to the bed, in this section, the analysis focused on the bottom to obtain the characteristics of the hydrodynamics that might be responsible for sediments transport and, consequently the scour around the pile.

**Figure 9.** General description of the mean velocities and vorticities near the pile for the longitudinal and cross profiles (associated to the scenario E01). Waves and currents are move from the left to the right.

Figure 9 is divided into two main boxes. The upper box illustrates the flow characteristics in the longitudinal profile of the channel, while the lower box illustrates the transverse section. Both have included the streamlines and total vorticity (vector addition, <sup>ω</sup>*TD*/*C*), the mean velocity along the channel (*u*/*C*), the vertical mean velocity ( *w*/*C*), the mean vorticity along the channel (<sup>ω</sup>*xD*/*C*), the transversal mean vorticity ( <sup>ω</sup>*yD*/*C*), and the vertical mean vorticity ( <sup>ω</sup>*zD*/*C*) in a nondimensionalized manner.

The total vorticity in the longitudinal profile and flow lines reflect the presence of a horseshow vortex as the main structure upstream the pile (*x*/*D* < 0), this structure approximately centered at *x*/*D* = −2.6 and *z*/*D* = 0.3. This vorticity system also reflects the average behavior of the flow longitudinal velocity, as in the zone equivalent to the horseshoe vortex location, *u*/*C* has negative values that are driven by the vortex counter clockwise rotation. Additionally, *w*/*C* also reflects the horseshoe vortex effects, since negative magnitudes of the vertical velocity (known as the down-flow) can be revealed near the pile. Meanwhile, at a distance lower than *x*/*D* = −2.6, the down-flow becomes positive (mainly driven by the counter clock turn of the vortex).

The pile downstream area (*x*/*D* > 0) in the mean streamline condition showed a rotational centered structure at around *x*/*D* = 2.4 and *z*/*D* = 1.6, which is related to the vortex shedding and could produce alterations in the mean field of the vertical and longitudinal velocities. Longitudinally, near the pile, a flow could be produced toward it, mainly caused by the momentum balance that would be developed at an approximate distance of *x*/*D* = 1.2. Then, a positive direction of the flow velocity up to a distance of *x*/*D* = 6 would be present, to subsequently re-adopt a negative velocity; since this behavior is associated with vortex shedding, as stated previously.

Vertical velocities in the downstream area show that, nearby the pile, the flow moves upwards. Meanwhile, an *x*/*D* = 2.7 distance would make the vertical velocity negative. This and the flow line behavior could be caused by clockwise rotation and related to the formation of vortex shedding.

Analyzing the transverse section of the flow in Figure 9, clear that in the mean streamlines, a horseshoe vortex system develops around the pile, with a longitudinal component of mainly positive velocities near the pile. In the case of a vertical component, these velocities would be negative. This means that: a downstream flow component is present once the horseshoe vortex surrounds the pile, and it remains near the bed pushed by the flow vertical velocities.

The presence of structures concordant with the horseshoe vortex around the pile is shown by the transverse vorticity ( <sup>ω</sup>*yD*/*C*). Thus, the vertical development of these structures is limited to an approximately height of *z*/*D* = 0.5 from the bed. The mean vorticity fields obtained for vertical and longitudinal components are congruen<sup>t</sup> with the classical description of the flow around the cylindrical piles for a permanent flow, which has been widely addressed in the literature [1–4].

The mean velocity vertical profiles for monitoring stations P1 to P8 are presented in Figure 10 for each modeled scenario. From this, it can be observed that the mean characteristics of the codirectional currents and wave velocity profiles are not significantly di fferent under opposite flow conditions in either of the proposed scenarios (those controlled by currents or those controlled by waves).

Clear evidence of this phenomenon is shown in P5 (Figure 10), which corresponds to a downward station for codirectional currents and waves and an upstream station for waves in opposite flow scenarios. In each case, the station presented the smallest flow magnitudes, indicating that for all the simulated scenarios of station P5, currents dominate the flow, according to the flow relative velocity (*UCW*) proposed by Sumer and Fredsøe [25], this result indicates that a combined regimen is present. This aspect will be further addressed in the analysis of the results.

It is important to note that the currents are symmetrical based on the pile geometrical characteristics (cylinder) and the studied flow. Figure 10, shows that the velocity profiles of station pairs P2 with P8; P3 and P7; and P4 and P6 are concordant not only in their magnitudes but also in their vertical axis. Thus, describing only one of the stations associated in a pair is su fficient to understand the flow around the pile.

*Water* **2019**, *11*, 2256

The velocity profiles for both codirectional and opposite cases reached their maximum value in station P3 (P7), mainly due to the contraction of flow lines producing accelerations, thereby increasing the velocities from station P1 towards P3 (P7), followed by a gradual decrease from station P3 (P7) towards P5.

**Figure 10.** Mean streamwise velocity vertical profile for all the simulated cases.

In light of the results described above, and for the simulated scenarios, the hydrodynamics around a pile for combined flow of waves and currents (codirectional and opposite) should behave like a normal flow around a cylinder due to its steady state flow, which is widely described in the literature.

The amplification of the mean shear stresses (ατ) made dimensionless with an undisturbed bed shear stress around a cylindrical pile is shown in Figure 11 for each of the simulated cases, in which the left column shows codirectional currents and waves cases and the right column shows the opposite current and wave cases.

The maximum amplifications of the main bed shear stress were produced in the pile lateral edge and reached magnitudes of 2.5 for codirectional and opposite currents and waves cases. Nevertheless, when waves act on the current in an opposite direction, the amplification zone coverage is reduced (smaller area) compared to the codirectional current and wave cases. As a general trend, the amplification obtained shows that shear stress gradually decreases from hydrodynamic scenarios dominated by currents (E01 and E02) to those dominated by waves (E07 and E08).

Analyzing the results of E01 and E02, the differences found between the spatial distributions of the bed shear stress are not significant when the waves act codirectionally with currents or when they are opposite. No significant movements were noticed in the maximum amplification localization. This means no direct influence of the waves flow or ebb was found for the development of the shear stress in the bottom.

**Figure 11.** Bed shear stress made dimensionless with undisturbed bed shear stress amplification for each simulated case.

A behavior equivalent to that described for E01 and E02 was identified in E03, E04, E05, and E06, that is, no influence of the waves in the mean bed shear stress distribution was evidenced, even though according to *Ucw* waves domain. In general hydrodynamics, should become more significant when transiting from scenario E01 to E07 and E08. In these last scenarios (E07 and E08), the lowest amplifications of the mean bed shear stress were obtained, which were 1.5.

## *3.3. Scour around a Cylindrycal Pile*

This section presents the analysis results of the scour around the pile. Figure 12 presents the dimensionless scour time series (*S*/*D*) obtained experimentally by Qi and Gao [26] and the series resulting from the numerical model, for codirectional (E01) and opposite (E04) waves and currents.

**Figure 12.** Maximum scour development comparison between the experimental data (Qi and Gao [26]) and numerical simulation results for cases E01 and E04.

The information for the codirectional flow (E01) illustrated in Figure 12 shows a strong agreemen<sup>t</sup> between the experimental and simulated data, in its temporal evolution (curve form) and in the magnitude reached by the dimensionless scour. For example, after 10 minutes of simulation the model reached a magnitude of *S*/*D* = 0.173, and the experimental data reached a magnitude of *S*/*D* = 0.169 (2.4% of the relative error).

The comparison of numerical data versus experimental data for the dimensionless scour with and opposite flow (E04) illustrated in Figure 12 (as the previously described for E01) showed a strong correspondence in both its the temporal evolution and the magnitude reached. For example, within 20 minutes, the experimental data shows the dimensionless scour would be 0.146, while the numerical model showed that the scour would be 0.139 (4.8% of the relative error).

The generality of the results obtained and presented in Figure 12 indicates that the scour would be greater when the flow acts codirectional to the waves and currents than in the opposite case. The latter is supported by the results summarized in Table 7 column *S*1/*D*, which is the dimensionless scour obtained in the last time step of the numerical model. Additionally, Table 7 lists results for the adjustment coefficients of the equilibrium scour equation from Sheppard et al. [54] (*<sup>a</sup>*1 to *a*4), the equilibrium dimensionless scour (*St*/*D*), and the relative scour factor (*St*/*S*1). A similar method was used by Qi and Gao [26] but with experimental data.


**Table 7.** Adjustment parameter for equilibrium scour estimate and the results obtained from the numerical simulation.

The relative scour factors indicated in Table 7 were greater in cases where the flow velocity defined at a distance of *z* = *D*/2 (*Uc*) was greater. This mean that the scour obtained from the numerical model of a 25 minutes sediment transport and resulting morphodynamics evolution of the bed, is very different from the equilibrium in cases whit greater flow magnitude.

The latter may be associated with major currents acting on the center, producing major scours and subsequently requiring greater action times in order to reach scour equilibrium [5]. Nonetheless, when magnitudes of the dimensionless equilibrium scour estimated by the equation developed by Sheppard et al. [54] and based on the data obtained here for the 25 minutes simulation are compared with the experimental data obtained by third parties, equivalent results and estimates within an acceptable range of experimental variability are obtained, as shown in Figure 13.

**Figure 13.** The equilibrium scour estimated from numerical data and its comparison with the equilibrium scour obtained by other authors.

Figure 13 presents a comparison between the numerical data obtained in this article (blue circles for codirectional scenarios and red circles for opposite flows), and experimental ones found in the literature [21,25,26,55,56]. This graph was formed similarly to that presented by Sumer and Fredsøe [25], i.e., the dimensionless equilibrium scour as a function of the Keulegan–Carpenter (*KC*) number for di fferent relative velocity ranges between waves and currents ( *UCW*).

From the general analysis of Figure 13, it is observed that, for the range 0.10 < *UCW* < 0.40, the data obtained from the numerical model implemented in this article would be close to the data obtained by Raaijmakers and Rudolph [21] and Sumer et al. [55] for similar Keulegan–Carpenter numbers (around five) for both the codirectional flow condition and the opposite. This situation repeats in the comparison between ranges 0.40 < *UCW* < 0.50 and 0.50 < *UCW* < 0.80, meaning that the data obtained from the results projection of the numerical model toward the equilibrium scour based on the equation of Sheppard et al. [54] allows us to gather magnitudes that can be compared with the experimental records obtained by other researchers, for similar hydrodynamic characteristics.
