In the tidal current simulation, the TELEMAC-3D module (V8p4r0 version) is governed by the non-hydrostatic Navier–Stokes equations and the continuity equation. In the vertical direction, a σ-coordinate transformation is adopted, in which σ varies from 0 at the seabed to 1 at the free water surface. Assuming the fluid is incompressible [26], the effects of temperature and salinity are neglected in the simulation. The governing equations in the Cartesian coordinates are as follows:(1)∂U∂x+∂V∂y+∂V∂z=0(2)∂U∂t+U∂U∂x+V∂U∂y+W∂U∂z=−1ρwater∂Ppre∂x+v∆U+Fx(3)∂U∂t+U∂U∂x+V∂U∂y+W∂U∂z=−1ρwater∂Ppre∂y+v∆V+Fy(4)∂U∂t+U∂U∂x+V∂U∂y+W∂U∂z=−1ρwater∂Ppre∂z+v∆W+Fz where U, V, and W are the velocity components in the x, y, and z directions, respectively; F_x, F_y, and F_z represent the source terms of the fluid in the x, y, and z directions, respectively; Δ denotes the three-dimensional Laplacian operator; ρ_water is water density; t is the time; υ is the kinematic viscosity (υ = 1 × 10⁻⁶ m²/s); P_pre is the total pressure, which is the sum of the atmospheric pressure on the water surface (P_atm) and the hydrostatic pressure (P_d) and can be calculated as follows:(5)Ppre=Patm+ρwatergzs−z+Pd where z_s represents the free surface elevation.

In this study, the source terms include a structural drag source term used to represent the momentum-blocking effects of monopiles, which is described in detail in Section 2.2. The governing equations are solved using a semi-implicit numerical scheme.

The initial velocity field is initialized using the mapped velocity field derived from the TOPEX/Poseidon global tidal model (TPXO, TPXO10-atlas-v2 version), based on the Laplace tidal equations and satellite altimetry observations processed using the OSU Tidal Prediction Software [27]. The simulation focuses on barotropic tidal currents, and density stratification effects are neglected, while the Coriolis force is included in the momentum equations. To minimize the influence of initial transients, a spin-up period was applied [28]. Given that the TPXO forcing includes multiple tidal constituents, the shortest periods are semidiurnal. Accordingly, the 3–4-day spin-up corresponds to approximately 6–8 semidiurnal tidal cycles, which ensures that the flow field reaches a quasi-periodic state. Therefore, a 4-day spin-up was used for validation tests to provide a conservative margin, while 3 days were deemed sufficient for scenario simulations. The initial free-surface elevation is set to zero, and the model is allowed to adjust dynamically during the spin-up period.

The global tidal model TPXO provides the tidal forcing at the open-sea boundary and is among the most accurate available [29]. It includes 13 tidal constituents: M2, S2, N2, K2, K1, O1, P1, Q1, M4, Ms4, Mn4, Mm, and Mf, on a 1/6° spatial resolution grid. The land boundary is treated as a closed boundary, where both water and velocity fluxes are set to zero. During the simulation, the water level and horizontal velocity components (U and V) at the open boundary are expressed as the sum of harmonic constituents that vary over time and space. For each constituent, the water level or velocity component can be calculated as follows:(6)fx,y,t=∑ifi(x,y,t)(7)fix,y,t=fitAfix,ycos2πtTi−βfix,y+βi0+gi(t) where f denotes either the water depth or one component of the horizontal velocity (U or V); i is the index of the tidal constituent; T_i is the period of harmonic components; A_fi and β_fi are the amplitude and phase of variable f for the i-th constituent, respectively, while β_0i denotes the initial phase. The functions f_i(t) and g_i(t) represent the nodal factors that account for the lunar modulation of tidal effects. Both A_fi and β_fi depend solely on spatial location and remain constant in time.

In addition to the tidal boundary forcing described above, bottom friction is parameterized using a quadratic friction law to represent near-bed momentum dissipation. The structural drag is vertically distributed to represent the depth-integrated momentum loss induced by monopiles.

In offshore wind farms, strong interactions occur between wind turbine monopiles and the surrounding flow, resulting in local momentum dissipation and flow field redistribution. For large-scale numerical models, directly resolving the detailed flow around each pile is computationally prohibitive. Therefore, a parameterized drag source term is introduced to represent the drag effect of monopiles as a continuous body force distributed over the control volume, allowing the averaged resistance effect to be captured at the macro scale. This source term follows the classical drag force formulation for bluff bodies [30,31], and the structural drag force term can be expressed as(8)FD=−12ρwaterCDArefUinUin where C_D is the structural drag coefficient, A_ref is a reference area representing the projected area of the pile at a characteristic scale, and the source term acts opposite to the local flow velocity U_in, representing the momentum loss induced by the pile. |U_in| denotes the magnitude of the local flow velocity vector U_in.

If the body force is applied only at a single grid point, it may cause numerical instability and unrealistic local momentum losses. To avoid this, a spatial smoothing kernel is introduced to distribute the drag effect continuously around the pile, forming an effective drag cloud in the computational domain. The kernel radius typically ranges from several times the pile diameter to several grid spacings, depending on the local resolution and numerical stability. The kernel function generally takes the following form:(9)Wr=exp−rRkernel2 where r denotes the horizontal distance from the pile center, and R_kernel defines the spatial decay range of the drag field.

To ensure conservation of total momentum dissipation in the numerical sense, the kernel function is normalized as shown in Equation (10):(10)∫VdWrdVol=1 where V_ol denotes the control volume over which the kernel function is integrated, V_d is the drag distribution volume.

For the overall effect of a monopile array, structural drag is represented as a volumetrically distributed drag force source term constructed from kernel functions, thereby capturing the cumulative effect of the array at the model scale.

A self-implemented structural drag source-term module, following the approach of [18], is implemented within the TELEMAC-3D framework as a dynamically evaluated drag force source term. In the numerical implementation, the structural drag is introduced as an additional source term in the momentum equations. It is treated explicitly within the semi-implicit hydrodynamic solver of TELEMAC-3D, with the source term evaluated dynamically at each time step.

For each horizontal grid node, the horizontal distance to the center of each monopile is computed. A Gaussian-based kernel is adopted to enhance numerical robustness in the near-pile region. To account for non-uniform mesh geometry, the kernel weights are area-weighted and subsequently normalized, ensuring conservation of the total applied drag force over the horizontal mesh. For each monopile, the drag force is evaluated based on a representative inflow velocity, computed as a horizontally weighted spatial average of the surrounding flow within a prescribed neighborhood around the pile center. This representative velocity is then distributed across the vertical layers in proportion to the local layer thickness, ensuring consistency with the depth-integrated representation of the pile-induced drag. The resulting drag force associated with each pile is distributed to the surrounding mesh nodes according to the corresponding kernel weights. The kernel radius R_kernel is prescribed as a model parameter and determines the spatial extent over which the pile-induced drag is applied.

In the present simulations, the kernel radius is set equal to the monopile diameter, which controls the spatial extent of the applied drag force and ensures numerical robustness of the explicit source-term treatment. Additionally, kernel support is truncated beyond a prescribed cutoff distance for computational efficiency. At each horizontal grid point, the total structural drag is obtained by accumulating the contributions from all piles whose kernel support overlaps the local mesh, thereby representing array-scale momentum dissipation in a spatially distributed and dynamically coupled manner.

Within the three-dimensional hydrodynamic framework, the structural drag source term is coupled directly to the momentum equations of TELEMAC-3D. After being constructed on the horizontal plane, the drag source term is vertically distributed uniformly throughout the water depth. This treatment ensures that the depth-integrated drag force remains consistent with the pile-based parameterization while allowing the pile-induced resistance to be projected into the full three-dimensional domain.

In the discrete form, the horizontal components of the structural drag source term can be written as Equation (11):(11)Sx,y=−12ρwaterCDDpileAkernUinUinWr where D_pile is the pile diameter, and A_kern represents the integrated area of the truncated Gaussian kernel within its influence radius.

The drag source term is applied consistently with the model time-stepping procedure without altering the original structure of the hydrodynamic solver. By directly incorporating directionally resolved drag forces into the three-dimensional momentum equations, this coupling strategy enables accurate reproduction of the modulation effects of the monopile array on tidal momentum and flow field distribution, without explicitly resolving local pile-scale flow structures.

The nearshore region off the northern Shandong Peninsula, China, was selected as the numerical simulation domain (Figure 1). This region lies at the junction of the Yellow Sea and the Bohai Sea, with geographical coordinates of approximately 120° E–123° E and 36.0° N–38.5° N. Influenced by the peninsula’s complex topography, the coastline is highly indented, forming diverse local hydrodynamic patterns. The nearshore water depth is generally less than 10 m and gradually increases to about 70 m offshore, with isobaths extending seaward in a fan–shaped pattern. The region is characterized by a regular semidiurnal tide, driven by the combined influence of the East Asian monsoon and the tidal regime. As a representative shallow coastal area in the northern Yellow Sea, this region is characterized by strong tidal currents and increasing offshore wind energy development, making it particularly suitable for investigating monopile array-induced hydrodynamic modifications at both local and regional scales.

To achieve a reasonable balance between computational accuracy and efficiency, a grid sensitivity analysis was conducted to determine the appropriate horizontal resolution. The computational domain was discretized using an unstructured triangular mesh to accommodate the complex shoreline geometry and variations in bathymetric gradients. To better resolve the hydrodynamic interactions around the array, a local mesh refinement was applied in the offshore area. Except for the horizontal grid resolution, all other model settings were kept identical throughout the sensitivity tests, including the time step, vertical layering, and boundary conditions, with only the horizontal mesh scale being varied. A vertical resolution was selected based on previous regional tidal modelling studies in similar shallow-water environments, providing a reasonable compromise between resolving vertical shear and computational efficiency.

As shown in Table 1, four candidate grids (G1–G4) were constructed with progressively refined resolutions to systematically evaluate the influence of horizontal resolution on the simulation results. The sensitivity analysis indicates that the resolutions of G1 and G2 were relatively coarse, leading to larger deviations between the simulated and observed tidal elevations and current velocities. With increasing resolution, the deviations gradually decreased and became nearly convergent between G3 and G4. The differences in the main evaluation metrics between these two grids were small, and the spatial distribution patterns of tidal currents were also largely consistent. However, the computational cost of G4 increased significantly, with the number of mesh elements being approximately 3.17 times that of G3, and the total computation time about 4.7 times longer, while the improvement in accuracy was relatively limited.

Considering both computational accuracy and efficiency, grid G3 was finally selected as the unified mesh for all subsequent simulations. This grid provides a good representation of shorelines, bays, and bathymetric gradients, achieving a balance between model accuracy and computational cost, and ensuring the comparability of results among different physical scenarios.

The open boundary was located in the central Yellow Sea, where tidal forcing was driven by the main harmonic constituents provided by the TPXO global tidal database. The initial water level was obtained from the superposition of harmonic components, while the initial velocity field was set to zero. The land boundaries were specified as impermeable (no-flow) conditions. Bottom friction was parameterized based on the median sediment size, with the Nikuradse roughness height defined as Equation (12):(12)ks=2.5d50 where d₅₀ is the median sediment grain size.

The model used a σ-coordinate system with six vertical layers, and the time step of 10 s satisfied the CFL stability criterion. All subsequent simulations were conducted on the G3 mesh to ensure consistency and comparability among different scenarios.

To ensure the accuracy of the tidal model, data from five tidal elevation stations, three depth-averaged velocity stations, and two vertically stratified velocity stations are used for model validation. Figure 1 and Table 2 show the details of these measurements.

Figure 2 compares the observed tidal elevations with the simulated results at five tidal stations (E1–E5). Stations E1–E4 include both summer (orange squares) and winter (blue circles) measurements, allowing assessment of the model’s seasonal performance under varying hydrodynamic conditions. Station E5, although not separated by season, spans a relatively long observation period (green triangles), providing additional validation of the model’s long-term stability. Overall, the simulated tidal elevations are in reasonable agreement with the observations. The model captures both the phase and amplitude of the dominant tidal constituents well, with only minor discrepancies at certain peaks and troughs, indicating satisfactory overall accuracy of the tidal-elevation simulation.

Figure 3 compares the simulated and measured tidal velocity and direction at the three observation stations (V1–V3). The model reproduces the temporal variations of tidal currents in good agreement with the field measurements. At all stations, the simulated velocity aligns with the observed data in both magnitude and phase, and the ebb-flood cycles are accurately represented under spring and neap tidal conditions. The simulated current directions also capture the observed rotational and reversal patterns throughout each tidal cycle, indicating that the model can effectively describe the periodic variations in flow orientation. Minor discrepancies occur near several turning points or extreme values, likely due to local bathymetric influences or measurement uncertainties. Overall, the results confirm that the model provides a reliable representation of the tidal current dynamics in the study area.

Furthermore, Figure 4 and Figure 5 compare the simulated and measured tidal velocity and direction at three vertical layers in the water column at Stations V4 and V5. The model captures the main vertical and temporal variations of tidal currents, with simulated amplitudes and phases closely matching the field observations. The periodic reversals of current direction are also captured, indicating that the model can reasonably describe the three-dimensional tidal dynamics. Minor discrepancies remain, mainly manifested as slight underestimation near the bottom at V4 and weaker surface currents with small phase lags at V5, which are likely related to near-bed friction and surface turbulence. Nevertheless, these deviations are minor and do not affect the overall agreement, confirming that the model can reasonably simulate the three-dimensional tidal dynamics in the study area.

In summary, the validation results confirm that the model can adequately capture the regional tidal regime and resolve the vertical structure of tidal currents, ensuring a robust hydrodynamic framework for subsequent analyses.

To validate the applicability of the structural drag source-term approach at the monopile scale, an idealized monopile test was conducted following the experimental conditions reported by Roulund et al. [6]. The computational domain was configured as a rectangular flume with a length of 35 m, a width of 3 m, and a uniform water depth of 0.54 m. A constant inflow velocity of U₀ = 0.326 m/s was prescribed, and the monopile diameter was set to D = 0.536 m. The monopile was represented by an equivalent body force applied at the pile center, enabling an assessment of the model’s capability to reproduce wake characteristics and associated momentum dissipation.

Figure 6 shows that the implemented drag source-term method captures the primary hydrodynamic features induced by a monopile at mid- to upper-water depths. Upstream of the monopile, the model shows the continuous deceleration of the incoming flow resulting from monopile blockage, with streamwise velocity variations in good agreement with experimental measurements. However, when the measurement location approaches the bed (z < 0.02 m), the agreement between simulated and measured velocities upstream of the monopile deteriorates. In the downstream region, discrepancies are also evident in the near-wake zone. While the experiments exhibit pronounced flow recirculation and strong vortex structures, the present simulations yield a comparatively smoother representation of these local flow features. The agreement between the model and reference data is further quantified using the mean absolute error (MAE), as summarized in Table 3. It should be noted that the MAE is computed using the experimental data of Roulund et al. [6] as the reference and the numerical results as the predictions.

To further assess the validity of the proposed approach within a numerical modelling framework, the simulated wake velocity distributions and their spatial evolution were compared with results obtained using a geometric masking (geo-mask) method, which was also implemented within the TELEMAC framework. A reasonable agreement between the numerical approaches indicates that, in large-scale hydrodynamic modelling, the structural drag source-term method can reasonably capture the main features of the wake structure and associated momentum dissipation. The MAE values range from 0.009 to 0.0166 m/s (Table 4), providing an objective measure of the consistency between the results of this study and the geo-mask data (numerical simulation result). It should be noted that the discrepancies observed in the downstream near-wake region and in the near-bed layer are not attributable to numerical instability or implementation errors. Instead, they arise from the parameterized nature of the structural drag source-term approach, which represents the monopile through an equivalent momentum source without explicitly resolving fine-scale processes such as flow separation, vortex shedding, and shear-layer development. Consequently, a close quantitative match with high-resolution CFD simulations or laboratory measurements in terms of near-field flow details is not expected. These discrepancies, therefore, reflect the inherent limitations of large-scale parameterized models when applied to strongly three-dimensional and boundary-layer-dominated flows. Accordingly, the drag-related parameters used in the structural drag source-term approach, including the drag coefficient and the kernel width, were selected based on the idealized monopile validation against laboratory measurements, to ensure that the overall momentum loss induced by the monopile is reasonably represented in an integrated sense rather than to match detailed local velocity profiles.

Overall, by preserving the monopile-induced momentum loss in an integrated sense, the proposed method provides a physically consistent representation of the hydrodynamic response at its intended modelling scale. Although the local-scale validation is conducted using a single monopile configuration, it primarily serves to verify the physical consistency and numerical stability of the structural drag source-term formulation. At the array and regional scales, where the flow response is governed by the spatial superposition of identical source terms and large-scale momentum balance. The validated single-structure behavior provides a necessary physical basis for its extension to the wind farm array. At these larger spatial scales, where the hydrodynamic evolution is primarily governed by large-scale momentum balance and dissipation rather than instantaneous vortex dynamics, the structural drag source-term approach offers a theoretically coherent and computationally efficient modelling strategy. The extension from a single monopile to an array configuration is physically justified because the present model represents pure drag-induced momentum loss from individual monopiles, without involving any energy conversion process. The drag force acting on each monopile is governed by local flow conditions and remains fundamentally consistent, even in the presence of array-induced flow interactions. Consequently, the drag coefficient calibrated from the single monopile validation can be directly applied in the array-scale simulations. In this study, C_D is set to 0.3, which falls within the typical range reported for cylindrical monopiles under high Reynolds number flow conditions. This value accounts for the spatial spreading of the drag force through the kernel function, avoiding an overestimation of momentum loss that would occur if using the full-scale drag coefficient of an isolated monopile. It should be noted that the present approach is not intended to replace high-resolution CFD simulations for detailed near-field flow analysis, but rather to efficiently represent structural effects within large-scale hydrodynamic models.

Based on the above validation results, the structural drag source-term approach shows good applicability and physical consistency across scales. In the idealized monopile test, the simulated flow field agrees well with experimental measurements, confirming its reliability at the local scale. In the field validation, the model simulates tidal elevations and current variations with good accuracy, indicating good performance under complex boundary conditions and time-varying tidal forcing. Overall, the method exhibits consistent hydrodynamic behavior, supporting its application in the wind turbine monopile array simulations and regional-scale impact assessments.

Figure 7 presents the spatial distributions of tidal elevation and current velocity fields, with flow directions superimposed, at representative instants during spring and neap tides. Overall, the tidal current intensity during spring tides is significantly higher than that during neap tides. Specifically, as shown in Figure 7a,b, the spring flood and ebb stages are characterized by relatively strong currents, with instantaneous domain-averaged velocities of 0.478 m/s and 0.470 m/s, respectively. During neap tides (Figure 7c,d), the tidal current intensity is substantially reduced, with the instantaneous domain-averaged velocities decreasing to 0.308 m/s and 0.281 m/s, respectively. Similarly, tidal elevation shows a larger amplitude during spring tides than during neap tides, whereas its spatial distribution remains broadly similar across tidal stages. This indicates that tidal elevation primarily acts as a regional background water-level variation, without introducing pronounced spatial heterogeneity across the study area.

Figure 8 shows the tidal current velocity fields at both the array scale and representative monopile locations. At the array scale (Figure 8a), a clear spatial gradient of background tidal current velocity is observed, increasing gradually from the nearshore region toward offshore areas, with localized low-velocity wake regions forming downstream of individual monopiles along the dominant flow direction. The black square marks the representative area selected for detailed local inspection (zoomed in). Figure 8b–e presents the local velocity fields at a representative monopile location at characteristic instants selected near the peak current velocities during the spring flood, spring ebb, neap flood, and neap ebb stages, respectively. The results show that, at these representative near-peak instants, monopile wakes appear as coherent low-velocity regions extending downstream under different tidal stages.

During spring tides, although the absolute velocities within the wake regions remain relatively high, the velocity deficits relative to the incoming flow are more pronounced, making wake structures more clearly identifiable and spatially extensive. In contrast, under neap tide conditions, the weaker background currents result in reduced relative velocity deficits and weaker wake-induced disturbances, causing wake features to appear less distinct and more spatially confined. In addition, during flood and ebb phases, the downstream extension of the wake aligns with the dominant tidal flow direction. Overall, these observations indicate that the prominence of monopile wake structures is influenced not only by the absolute velocity levels within the wake regions, but also by the relative contrast with the background inflow.

The array-scale wake structures discussed in Section 4.1.1 primarily illustrate the spatial patterns of tidal current velocity under different tidal conditions, but do not provide a direct quantitative measure of the magnitude of array-induced flow modification relative to the free-field condition. To facilitate a consistent and quantitative comparison among different tidal regimes, a normalized velocity metric is introduced. This metric is defined as the relative deviation of the flow velocity under array conditions from the corresponding free-field velocity and is expressed in percentage form. Positive and negative values indicate relative velocity enhancement and reduction, respectively. This formulation provides a unified, normalized measure of array-induced flow modulation. The normalized velocity change is defined as follows:(13)∆Unorm%=Uarray−UfreeUfree×100% where U_array and U_free denote the flow velocity under wind turbine monopile array and free-field conditions at the same spatial location, respectively. In this study, the velocity magnitude used in the calculation is time-averaged over the spring tide period and the neap tide period, respectively.

As shown in Figure 9, both spring and neap tidal conditions exhibit clear and consistent response of the flow field to the presence of the array. Specifically, under both tidal regimes, the normalized velocity change reaches comparable percentage levels, with localized extrema of approximately −10% occurring in the wake regions downstream of individual monopiles. Despite the similarity in extreme values, the normalized velocity response under both spring and neap tide conditions remains primarily concentrated in the near-wake regions downstream of individual monopiles. The quantitative results corroborate the qualitative observations presented in Section 4.1.1, indicating that the two tidal regimes differ primarily in the wake contrast with respect to the background flow rather than in their spatial extent.

In general, positive changes on the order of 0% to +3% are widely distributed within the inter-row and peripheral areas of the array under spring tide conditions, resulting in a heterogeneous spatial pattern of coexisting positive and negative responses. Under neap tide conditions, the normalized velocity change within the array domain exhibits a more uniform and unidirectional pattern, with large portions of the array area being continuously dominated by negative changes, typically ranging from −3% to −8%. In contrast, regions exhibiting positive changes are significantly reduced in both spatial extent and frequency. Overall, spring and neap tides differ not in the peak normalized velocity change, but in the spatial organization and coherence of the response, as well as the relative dominance of positive and negative changes within the array domain.

The geometric layout of the array, together with the numbering scheme and sampling definitions, is illustrated in Figure 10. The left and right sides correspond to fixed geographical directions in the map view and represent purely geometrical positions in the layout, rather than being defined with respect to the dominant flow or downstream direction. This array consists of 8 transverse rows and 10 columns, and the rows are labeled from Line A to Line H from top to bottom. Monopiles are numbered starting from the upper-right corner of the array in the map view as monopile No. 1 and continue sequentially along the column direction.

The longitudinal midpoint line of the outermost monopile column is taken as a reference baseline. The sampling distance R_kernel is defined as the lateral distance from the reference baseline to a given sampling location. It represents the characteristic spatial scale of the kernel function used to distribute the drag force induced by individual structures into the surrounding flow field. This parameter controls the lateral extent over which the structural resistance is applied, thereby ensuring a continuous representation of the array-induced drag effect. In the present study, R_kernel is set to 10 m, corresponding to the characteristic diameter of the monopile (e.g., 15 MW, D = 10 m) [33]. This choice ensures that the spatial extent of the distributed drag force remains physically consistent with the structural scale, thereby localizing the momentum loss around each monopile. For brevity, R is used to denote R_kernel in the following contents. The downstream distance p is measured from the same baseline along the dominant flow direction.

The dimensionless distance p/R is then used to describe the relative spatial positions of different downstream sampling locations, with p/R = 1, 5, 10, and 20 corresponding to near-field to far-field regions, respectively. These distances approximately represent the near-wake region (p/R ≈ 1), transition zone to the far-wake (p/R ≈ 5), the far-wake recovery region (p/R ≈ 10) and the far field region (p/R ≈ 20). This selection ensures that the full spatial evolution of the flow disturbance induced by the array is captured. Moreover, the chosen range is broadly in line with wake recovery behavior reported in studies of flow past bluff bodies and offshore wind turbine arrays [34,35,36,37]. It should be noted that under realistic tidal conditions, the flow direction reverses between the flood and ebb phases. Consequently, the upstream and downstream directions are not fixed in space, and a given geometrical side of the array may correspond to an upstream position during one phase and a downstream position during another. Accordingly, both geometrical sides of the array are considered to ensure a comprehensive characterization of the array-induced velocity response under reversing tidal flows. Point O denotes the geometric center of the array, while points M and N are located at representative positions along the outermost monopile on either side of the array. These points are selected as characteristic sampling locations for subsequent velocity time-series analysis.

Figure 11 and Figure 12 illustrate the variations in time-averaged velocity magnitude along the array transverse rows (Lines A–H) at different downstream distances (p/R = 1, 5, 10, and 20) under spring and neap tidal conditions, respectively. For a given downstream distance, the velocity variation patterns along the array rows on the left and right sides of the array are generally similar along lines A–H, indicating a similar overall trend in the array-induced response. At the row-wise local scale, pronounced velocity attenuation is observed at the rows corresponding to monopile locations. As the incoming flow passes through a monopile row, the velocity decreases rapidly due to wake effects, while smoother and gradually recovering velocity profiles are observed in the inter-row regions. Quantitatively, under spring tide conditions, the velocity deficit (defined relative to the no-turbine case at the same location) at p/R = 1, 5, 10, and 20 on the left side is 37.9%, 17.2%, 10.3%, and 6.9%, respectively; on the right side, the deficits are 22.2%, 11.1%, 5.5%, and 1.4%, respectively. Under neap tide conditions, the corresponding deficits on the left side are 25%, 15.6%, 9.4%, and 3.1%, and on the right side are 22.2%, 8.3%, 4.4%, and 1.1%. This behavior reflects momentum dissipation by the monopiles and the subsequent velocity recovery driven by turbulent mixing and lateral transport. Therefore, the spatial structure of the velocity modulation is primarily governed by the overall array configuration rather than by lateral asymmetry in the flow field.

A comparison between spring and neap tidal conditions indicates that, although the overall velocity patterns remain similar, spring tide exhibits larger modulation amplitude. For example, at the inter-row region at p/R = 1, the velocity is 0.28 m/s under spring tide versus 0.16 m/s under neap tide, a 75% increase. This suggests that tidal intensity primarily affects the strength of the array-induced modulation, without altering its fundamental spatial structure. To quantify downstream recovery, the recovery length scale defined as the distance where the velocity deficit falls below 10%. This threshold is adopted as a practical criterion to enable consistent comparison of recovery distances under different tidal conditions. Using the linear interpolation method, the recovery distances under spring tide are p/R ≈ 11 and p/R ≈ 6 on the left and right side, respectively. While under neap tide conditions it reduces to about p/R ≈ 9.5 (left) and p/R ≈ 4.5 (right). These results indicate a faster recovery of velocity on the right side of the array and under weaker tidal forcing. The shorter recovery length under neap tide conditions is primarily attributed to the smaller initial velocity deficit associated with weaker tidal forcing, which allows the velocity deficit to fall below the threshold over a reduced downstream distance.

As Figure 13 shows, from the perspective of downstream evolution, pronounced periodic fluctuations in velocity are observed in the near-field region (p/R = 1) along the transverse direction, with distinct velocity deficits about the troughs occurring at locations corresponding to the monopiles. The array-induced velocity deviations exhibit clear negative peaks with magnitudes on the order of 10⁻² m/s. This denotes that the local flow field is primarily influenced by the direct effect of wakes generated by individual monopile rows. As the downstream distance increases to p/R = 5 and 10, the amplitude of velocity modulation gradually decreases, and the overall distribution becomes smoother, with reduced transverse gradients and weakened periodic fluctuations. Meanwhile, the recovery process exhibits spatial non-uniformity, with faster velocity recovery occurring in the inter-row regions compared with the monopile-row locations, indicating a progressive weakening of wake effects during downstream propagation. When the downstream distance further increases to p/R = 20, the velocity variations along the transverse direction are substantially reduced. The residual velocity deviation is confined to a very small range within ± 0.005 m/s. This indicates that the influence of the array on the velocity structure is significantly weakened at this downstream distance. However, weak periodic features associated with the array geometry can still be identified, suggesting that residual array wakes retain spatial organization in the mid- and far-field.

The time-averaged results show that the velocity deviation distributions on the left and right sides of the array remain consistent at all downstream distances. This symmetry reflects that, under tidal current forcing, instantaneous upstream-downstream asymmetries induced by the array are effectively mitigated through flow reversal, resulting in an overall symmetric spatial pattern of velocity modulation. Nevertheless, small local differences between the two sides can still be observed. These residual asymmetries are likely attributable to inherent differences between flood and ebb tidal conditions, including variations in peak velocity, duration, and acceleration characteristics, as well as background flow non-uniformity. Such differences remain minor compared with the dominant symmetric response induced by flow reversal.

To further examine the temporal characteristics of the array-induced velocity response, velocity time series at several representative locations are analyzed. These locations correspond to the representative points previously defined in Figure 10, namely the array center (O) and two characteristic positions (M and N) along the boundary monopile columns.

The velocity time series at the representative locations exhibit consistent tidal periodicity and aligned phases across multiple consecutive tidal cycles (Figure 14), indicating that the deployment of the array does not alter the fundamental temporal structure of the tidal flow. Instead, the array influence is expressed as a stable modulation of velocity amplitude, representing a persistent background effect rather than a transient response. Although spring tides are associated with higher velocity levels and stronger modulation amplitudes than neap tides, the relative temporal behavior among the representative locations remains consistent. From a single-point time-series perspective, the array-induced modulation therefore appears more pronounced during spring tides, whereas under neap tidal conditions, it is smoother and partially masked by the weaker background flow.

Although some differences in velocity magnitude can be observed among the representative locations, the overall temporal variations exhibit similar tidal oscillations with aligned trends. However, the magnitude differences remain minor, as the selected locations are not situated within the immediate near-wake cores of individual monopiles. Moreover, under tidal current forcing, periodic flow reversal further smooths local asymmetries over time, reducing the apparent strength of array-induced effects in single-point velocity time series. As a result, the temporal response at individual locations reflects a comparatively weak but consistently present modulation induced by the array.

The flow kinetic energy density is calculated based on the time-averaged value over the selected period. It is defined as(14)KE=12ρwateru¯2+v¯2¯ where u¯ and v¯ denote the depth-averaged horizontal velocity components, the overbar indicates time averaging.

To quantitatively characterize the overall energy state of the array region, two metrics are considered. The area-averaged kinetic energy density and total kinetic energy of the array region, both derived from the time-averaged kinetic energy density field. The former represents the mean kinetic energy density within the array coverage, while the latter quantifies the total kinetic energy obtained by spatial integration over the same region. Their expressions are given as follows:(15)KEarray¯=1Aarray∫AarrayKEdA(16)KEarraytot=∫AarrayKEdA where KE denotes the time-averaged kinetic energy density, A_array represents the horizontal area covered by the array, and the integration is performed over the array region in the horizontal plane.

Based on the definitions given above, the area-averaged and the total kinetic energy of the array region are quantified for spring and neap tidal conditions. As Figure 15 shows, in tidal spring conditions, the time-averaged kinetic energy density within the array-covered region is approximately 141.057 J/m³, and the corresponding area-integrated kinetic energy is about 4.48 × 10¹² J/m, indicating that both the overall energy level and the energy scale within the array region remain relatively high. In contrast, under neap tide conditions, the mean kinetic energy density within the array-covered region decreases to approximately 27.396 J/m³, with a corresponding area-integrated kinetic energy of about 8.71 × 10¹¹ J/m, suggesting a pronounced reduction in both the energy level and the energy scale within the array region under neap tidal conditions compared with spring tides.

From a spatial distribution perspective, despite the pronounced differences in absolute kinetic energy density levels between spring and neap tide conditions, the array-induced patterns of kinetic energy density remain structurally consistent. Low-energy regions are predominantly associated with the monopile locations and their surrounding wake areas, exhibiting a clear spatial organization aligned with the monopile layout and forming localized deficits with elongated wake-like features along the dominant flow direction. Because the absolute magnitude of kinetic energy density varies considerably between tidal conditions, different color scales are used to highlight the internal spatial distribution in each case, and these scales are not intended for direct comparison of absolute values.

To quantitatively characterize the influence of the array on flow kinetic energy, a radial analysis framework centered on the array centroid is established in the horizontal plane, based on which a characteristic array scale is introduced.

As shown in Figure 16, the array characteristic radius R_array is defined as half of the diagonal length of the geometric envelope enclosing all monopiles in the array. The selected bin width represents a reasonable balance between spatial resolution and statistical robustness, enabling the capture of local flow features while maintaining sufficient samples within each annular bin. Specifically, half of the inter-turbine spacing is adopted to ensure that the annular averaging scale remains smaller than the turbine spacing, thereby resolving the gradual radial variation of array-induced effects without excessive spatial smoothing. Based on this characteristic scale, a core region is defined within the array to represent a representative source region of momentum dissipation. This core region is located inside the geometric boundary of the array and is selected to characterize the superposition effects of multiple monopile foundations. By focusing on the interior of the array, where edge-related bypass flow and ambient-fluid compensation are minimized, allowing for a more representative and approximately additive assessment of array-induced momentum dissipation. Accordingly, a radial extent of 0.3 R_array is adopted to define the core area.

As illustrated in Figure 16, a series of concentric annular bins centered on the array centroid is constructed in the horizontal plane to evaluate array-induced effects as a function of radial distance statistically. This annular averaging framework enables a direction-independent quantification of the spatial modulation imposed by the array as an integrated structure, without introducing assumptions regarding a preferred flow direction. It is therefore particularly suitable for tidal environments with temporally varying flow directions. However, it should be noted that tidal currents exhibit pronounced directional characteristics associated with alternating flood and ebb phases, and wake structures may differ under different flow directions. In this context, the radial averaging method should be interpreted as a statistical spatial averaging approach rather than a representation of directional flow variations. Although this method inevitably weakens directional wake features, it remains effective for characterizing array-scale averaged effects. Accordingly, within each annular bin, the array-induced energetic response is quantified using the radially averaged relative kinetic energy variation, defined as Equation (17):(17)KE¯cri=1Ari∫AriKE¯cx,ydA where A_ri denotes the area of the annular region corresponding to the i-th radial interval centered at r_i, c denotes the test condition, corresponding to the array and free-field cases.

In this study, the radial interval is set to 850 m based on the characteristic spacing between adjacent monopiles in the array, with the typical transverse inter-monopile distance being approximately 1700 m. It should be noted that the selected annular-bin width is introduced solely for statistical analysis and does not represent a physical decay length of the array-induced effects.

Based on radially averaged kinetic energy defined (Equation (17)), the relative change in kinetic energy density with respect to the free-field condition, ΔKE, is defined in Equation (18). This metric is employed to quantify the relative energetic modulation imposed by the array on the flow system.(18)∆KE=KE¯arrayri−KE¯freeriKE¯freeri where KE¯(array) and KE¯(free) corresponds to the array and free-field cases, respectively.

Figure 17 presents the variation of ΔKE as a function of the normalized radial distance from the array center, x/R_array. ΔKE increases with radial distance and reaches a maximum at approximately x/R_array ≈ 2.2 during spring tides and ≈2.7 during neap tides. This behavior suggests that the relative reduction in flow kinetic energy induced by the array is not confined to the source region associated with array-induced momentum dissipation, but progressively manifests during the radial expansion of array-induced disturbances. This radial range reflects the location where array-induced disturbances accumulate during radial expansion, resulting in the maximum deviation from the free-field condition. With further increases in radial distance, ΔKE exhibits an overall decreasing trend. This indicates that the relative energy perturbations induced by the array are progressively compensated through lateral momentum exchange and wake diffusion, and that the flow field gradually recovers toward the free-field state. It should be noted that the radial location of the maximum ΔKE does not necessarily coincide with the core region, as the core region is defined for statistical representativeness rather than to fully capture all monopile-induced momentum dispersion. This definition avoids the influence of bypass and compensating flows near the array boundary, enabling a clearer assessment of the array-scale energetic modulation.

A comparison between different tidal conditions shows that the radial variation of ΔKE exhibits a generally consistent pattern, while differences arise in the decay rate and the spatial extent of the array influence. Under spring-tide conditions, ΔKE decreases more rapidly with radial distance following its peak, even after normalization by the free-field kinetic energy. This behavior is consistent with a more efficient attenuation of array-induced disturbances under stronger background flows, where advection and lateral momentum redistribution by the ambient flow are more pronounced. In contrast, under neap-tide conditions, the peak value of ΔKE is relatively higher and remains significant over a larger radial range, suggesting that the relative modulation imposed by the array on the flow system is more pronounced under weaker flow conditions, leading to an expanded spatial influence.