The Influence of Intra-Array Wake Dynamics on Depth-Averaged Kinetic Tidal Turbine Energy Extraction Simulations

Assessing the tidal stream energy resource, its intermittency and likely environmental feedbacks due to energy extraction, relies on the ability to accurately represent kinetic losses in ocean models. Energy conversion has often been implemented in ocean models with enhanced turbine stress terms formulated using an array-averaging approach, rather than implementing extraction at device-scale. In depth-averaged models, an additional drag term in the momentum equations is usually applied. However, such array-averaging simulations neglect intra-array device wake interactions, providing unrealistic energy extraction dynamics. Any induced simulation error will increase with array size. For this study, an idealized channel is discretized at sub 10 m resolution, resolving individual device wake profiles of tidal turbines in the domain. Sensitivity analysis is conducted on the applied turbulence closure scheme, validating results against published data from empirical scaled turbine studies. We test the fine scale model performance of several mesh densities, which produce a centerline velocity wake deficit accuracy (R2) of 0.58–0.69 (RMSE = 7.16–8.28%) using a k-ε turbulence closure scheme. Various array configurations at device scale are simulated and compared with an equivalent array-averaging approach by analyzing channel flux differential. Parametrization of array-averaging energy extraction techniques can misrepresent simulated energy transfer and removal. The potential peak error in channel flux exceeds 0.5% when the number of turbines nTECs ≈ 25 devices. This error exceeds 2% when simulating commercial-scale turbine array farms (i.e., >100 devices).


Introduction
Internationally coordinated efforts to set renewable energy targets have enhanced interest in the marine renewable energy (MRE) sector [1,2].Tidal energy converter (TEC) technologies that transfer kinetic energy of tidal stream currents into electrical power are moving from the testing and feasibility stage [2] to commercial scale array deployments.Notably, NOVA Innovation have three 100 kW devices deployed in Shetland, and Meygen have four 1.5 MW devices in the Pentland Firth, UK.Prior to large-scale development, it is crucial that we understand the physical feedbacks of energy extraction with the resource itself and the surrounding environment (e.g., [3,4]).Ocean models can aid this process, although uncertainties in the way that energy extraction is represented in models should be robustly quantified and minimized.A TEC array contains multiple turbines that operate independently based on the local flow passing each turbine rotor.A technique to reduce feedbacks between devices (i.e., increase performance) is to adopt a staggered array design.Common array configurations (e.g., [5]) lack any technical basis and fail to account for device and site-specific features, such as intra-array wake effects and misalignment to the flow [6,7].Improved understanding of such feedbacks and impacts through accurate model simulations could strengthen investor confidence in the sector [1,8].
Tidal stream energy extraction is often implemented as an enhanced momentum sink or decelerating seabed stress term in coastal modelling studies (e.g., [4,[9][10][11][12][13][14][15][16][17]), and Neill and Hashemi [18] provided a detailed review of techniques and applications.TECs are commonly parameterized in depth-averaged models as an enhanced seabed frictional force applied over a defined area (e.g., [13,19,20]).As this stress term denotes a decelerating momentum sink, valid results are produced when the drag parameterization is directly related to the total system energy lost, similar in effect to applying spatially varying seabed roughness to account for variations in seabed composition.As drag is implemented as a quadratic function of flow velocity, the area over which the force is applied is crucial, and will influence subsequent flow calculations.
Present guidance regarding resource assessment suggests incorporating energy extraction following resource characterization, where potential installed capacity may exceed 10 MW or extraction is greater than 2% of the available theoretical resource [21].Many shelf-scale studies typically use a mesh size in the order of hundreds of meters or kilometers (e.g., [14]), which would fail to resolve the dynamics of typical individual turbines (<20 m).Industry standard resource assessment (stage 2a) suggests 200-500 m at pre-feasibility stage [5].For full-feasibility (stage 2b) and design development (stage 3) assessments, which include device feedbacks, models become more computationally expensive and care must be taken when discretizing to prevent instability [22].The inclusion of turbine dynamics enables stage 3 assessment and the prediction of potential impacts based on development size.However, the derivation of enhanced terms for modelling energy extraction must be sufficiently justified [21].
In situ observations or computational fluid dynamics (CFD) modelling (e.g., [23]) should inform intra-array feedbacks.Assessments of mid to far-field impacts can only be addressed realistically using a numerical coastal modelling framework.Thus, a staged approach to specification of technique is suggested by the International Electrotechnical Commission (IEC [21]), involving:

•
Actual energy captured by TECs requiring consideration of device 'water to wire' efficiency factor that accounts for hydrodynamic, mechanical and electrical losses.

•
Intra-array spacing of devices relative to the plan area of model cells.

•
Support structure form drag associated with each TEC.

•
Any associated energy losses in downstream wake mixing and free stream velocity interaction between devices.
The combination of these factors ensures that macro-scale system feedbacks are captured.This enables confident determination of hydrodynamic interaction curves that link inflow velocities to levels of energy extraction and array capacity potential within a system.Simulations of energy extraction that average the enhanced stress term across the seabed area of an entire array fail to capture the intra-array wake effects and misrepresent the power density delivered to each turbine rotor within the array by as much as 21% [24].
An issue noted by Shives et al. [25] is that when considering an array-averaging approach versus a device-scale approach, there is difficulty defining the free stream velocity profile applicable to each device.The velocity experienced at the face of the array will differ from that in the lee of the array, where flow retardation over a large area is significant due to the decelerating force.Therefore, wake 'recovery', or mixing length, is misrepresented using an array-averaging approach, because energy extraction by a single device will induce a flow wake effect that disturbs the currents available to downstream devices, thus influencing the performance of surrounding devices.This phenomenon is ignored in array-averaging studies and ultimately results in uncertainties in simulations [7,26].Further complications arise due to a lack of availability in validation data for mid to far field wakes from full-scale device arrays; therefore, CFD models may help improve coastal modelling simulations [27].The ability to upscale this information to simulate regional scale environmental feedbacks provides essential detailed calibration and validation of methods.
Individual TECs influence the local flow field and adjacent devices as flow velocity reduces immediately in the lee of a turbine, is enhanced at the rotor periphery, and creates a backpressure or thrust effect immediately upstream [28].TEC wake velocity reduction is a function of rotor operating state values, such as thrust; therefore, the rate of downstream velocity recovery will depend upon mixing between flow recirculation and the device wake.In unbounded flow, the velocity remains close to that of the incident regime; however, in the presence of bounding surfaces formed by the seabed, free surface and wakes from other devices, a pressure drop forms across each turbine.The fluid will accelerate in the regions of restriction, causing localized velocity increases in accordance with Bernoulli principles [28].
In this study, we investigate whether coastal models discretized at individual device-scale, sufficiently capture more hydrodynamic feedbacks compared with array-averaged models.To achieve this, the mesh density should, as a minimum, sufficiently resolve each 'installed' device area, and resolve the flow in-between devices.This approach requires sensitivity analysis of mesh size and careful consideration of the applied model turbulence closure scheme (TCS).Furthermore, appropriate modification of applied formulae is needed to account for the proportional influence of the device on depth varying flow [21,29].
We implement a self-modified energy extraction method, using a depth-averaged coastal model (Telemac-2D) that solves the Shallow Water Equations (SWE).The use of 2D modelling strategies is justified by the IEC [21].However, in regions of rapid flow divergence and convergence, where vertical velocities are greatly enhanced, the hydrostatic assumption is no longer valid.The use of 3D models may provide better approximations in these cases, particularly those with higher order turbulence resolving schemes.In either case, authors should report and justify in sufficient detail the derivation of the additional extraction term applied [21].The aim here is to ascertain a validated method of energy extraction that improves upon an array-averaged approach in depth-averaged simulations.The usefulness of this 2D model application is tested in anticipation of future work involving coupled wave and sediment dynamics across large spatial and temporal scales that would prove difficult using 3D modelling of arrays at device-scale.
Our objectives are: a.
To incorporate tidal energy extraction at device scale.b.
To test the performance of the modelled hydrodynamics.c.
To compare array-averaging with device-scale results.d.
To quantify associated uncertainty.
Objectives are achieved by utilizing published empirical scaled turbine tank test and CFD studies to apply a validation of the methodology.We then utilize the validated device-scale approach to implement multiple devices, for various small-scale array configurations and compare this approach with an equivalent array-averaging technique.A determinant error threshold is suggested, beyond which it becomes important to account for intra-array turbine wake dynamics in coastal environmental impact studies caused by device-wake interactions.In doing this, we present a validated approach to shelf-scale turbine array implementation at device-scale.The inclusion of intra-array dynamics improves confidence in hydrodynamic resource assessment.This becomes significant when project feasibility studies scale up from single devices (<3 MW) and small (3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20) or medium  arrays, to full-scale (>50 MW) commercial array farms.

Materials and Methods
Telemac is an open-source numerical model that solves the Reynolds Averaged Navier-Stokes (RANS) equations non-hydrostatic mode using Boussinesq approximations on an unstructured triangular computational mesh.The RANS equations include the eddy viscosity concept into the Navier-Stokes equations and the Boussinesq hypothesis parameterizes the value of eddy viscosity by linearizing the inverse of density, or buoyancy effects.In depth-averaged (2D) calculations, the SWE or Saint-Venant equations resolve the vertically averaged RANS equations using a finite element or finite volume algorithm [30].Here we utilize Telemac-2D v7p2 in finite element mode, run in parallel on Supercomputing Wales [31], a high-performance computing system.

Seabed Stress Term
Bottom friction or bed shear stress, τ 0 (kg m −1 s −2 or N m −2 ), is a decelerating force exerted per unit area of seabed due to the flow.The implementation of turbines enhances or retards the total seabed shear stress within the local region.Numerical 2D models represent this process via an additional stress term in the momentum equations.A simple representation of seabed shear stress, assuming a gently sloping seabed is represented by the quadratic friction law [32]: where the stress term is related to the depth-averaged flow speed, U (m s −1 ), by the quadratic friction law based on the drag coefficient; C D .ρ is water density (kg m −3 ); h is water depth (m); g is acceleration due to gravity (9.81 m s −2 ); and n (s m −1/3 ) is a Manning coefficient value that can be related to bed roughness and grain size.Bed shear stress can be used to determine sediment transport rates, based on predetermined flow thresholds for sediment resuspension and settlement.For some model applications, it may be important to accurately parameterize both spatial and temporal variations in seabed roughness and, hence, bed shear stress.TEC array developments, for example, both affect and are influenced by differences in substrate type and water column turbidity levels [33].

Blockage Ratio
Any validation method must be independent of blockage effects.This occurs when the blockage ratio, β of the validation data differs significantly from that of the simulated device.
In accordance with flume tank experiment values used for validation as published by Gaurier et al. ([34] p. 89), where R is turbine rotor radius, b is channel width and h is water depth, we calculate a nominal empirical experimental value of β = 0.048.Based on this value and a full-scale turbine rotor diameter of 20 m (at the extreme end of the range for present-day deployments), a uniform channel of width 130 m and depth 50 m is discretized in order to validate the methodology.

Mesh Resolution Dependency
For this study, an idealized irregular unstructured domain containing a regular sub-mesh in the central channel region (4 km long) is discretized.Various mesh densities of sub 10 m horizontal scale (10 m, 5 m and 1 m; Figure 1) were analyzed.A regular sub-mesh is used because irregular node spacing causes inconsistency of applied force magnitude per node within a specified area [29].We vary the mesh resolution to investigate the influence of mesh dependency on simulated velocity deficit.A model specific algorithm calculates seabed area based on the number of nodes (and their associated element size) that fall within a user-defined polygon, related to rotor diameter.As mesh resolution increases, the associated error in model calculated area increases, affecting subsequent enhanced bed stress calculations and hence the momentum sink.Computational inaccuracy related to mesh resolution is therefore a potential source of error in simulations and is expressed as a percentage (Error = 100[|approx − exact|/exact]) in Table 1.Plan view of the discretized tidal channel.Zoomed section illustrates various mesh density resolutions applied in order to implement a single TEC (solid black polygon) in the center of the domain.The associated system calculated, seabed area (shaded dashed polygon) over which the enhanced stress term is applied is thus element mesh size dependent (Table 1).

Turbulence Closure Schemes
The TCS of a numerical coastal model provides a means to close a system of time averaged flow equations [35].TCS available in RANS simulations have limited capability and will not resolve microscale turbulence [36]; however, for commercial array feasibility studies, macrophysics in the mid-to far-field are resolved sufficiently to meet industry assessment needs [37].Classical closure schemes based on time averaged RANS equations include: The number of equations denotes the additional partial differential equations to be solved.Most first-order models (zero, one and two equation schemes) resolve micro-scale eddy simulation as space-filtered, time-dependent equations to calculate larger circulations.Second-order models (Algebraic Stress Model, Reynolds Stress Model) resolve smaller-scale effects on the flow pattern at sub-grid scales [38].For environmental numerical modelling applications, grid resolution and time step are usually too coarse to resolve sub-grid turbulent processes and scales of motion.Only when the variables are not to be temporally and spatially averaged is there a need for filtering, and a more appropriate closure scheme assumption is required.
Telemac-2D offers four schemes of varying complexity: Plan view of the discretized tidal channel.Zoomed section illustrates various mesh density resolutions applied in order to implement a single TEC (solid black polygon) in the center of the domain.The associated system calculated, seabed area (shaded dashed polygon) over which the enhanced stress term is applied is thus element mesh size dependent (Table 1).

Turbulence Closure Schemes
The TCS of a numerical coastal model provides a means to close a system of time averaged flow equations [35].TCS available in RANS simulations have limited capability and will not resolve micro-scale turbulence [36]; however, for commercial array feasibility studies, macrophysics in the mid-to far-field are resolved sufficiently to meet industry assessment needs [37].Classical closure schemes based on time averaged RANS equations include:

•
One equation models (Eddy viscosity concept, Bradshaw et al. etc.) The number of equations denotes the additional partial differential equations to be solved.Most first-order models (zero, one and two equation schemes) resolve micro-scale eddy simulation as space-filtered, time-dependent equations to calculate larger circulations.Second-order models (Algebraic Stress Model, Reynolds Stress Model) resolve smaller-scale effects on the flow pattern at sub-grid scales [38].For environmental numerical modelling applications, grid resolution and time step are usually too coarse to resolve sub-grid turbulent processes and scales of motion.Only when the variables are not to be temporally and spatially averaged is there a need for filtering, and a more appropriate closure scheme assumption is required.Constant viscosity (CV, the default in Telemac) 2.
Smagorinski model These turbulence closure models operate on the assumption that turbulent viscosity is either constant or is directly dependent on known calculated parameters.In 2D models, the horizontal components of eddy viscosity and diffusivity should also contain a contribution due to the shear or vertical variation of horizontal flow.Where the horizontal velocities remain considerably greater than vertical velocities, a depth-averaged simulation approach can provide a reasonable solution.In this study, we apply two of the four schemes to simulations: the default (1), and that suggested by Telemac experts as the most applicable for energy extraction methods (3).Hervouet [35] discusses in more detail all four TCS available.

Constant Viscosity Scheme
This scheme equates a constant viscosity coefficient that represents all forms of mixing including molecular viscosity, turbulent viscosity and dispersion.The Telemac-2D velocity diffusivity coefficient has a model default value of 10 −6 , corresponding to the molecular viscosity of water.When this approach is applied to simulations, the optimal velocity diffusivity coefficient value must be determined.Values for horizontal velocity diffusivity (m 2 s −1 ) applicable to coastal applications are difficult to source, although values of between 0.1 and 50 are suggested for fluvial environments [39].Haverson et al. [40,41] use the CV TCS and apply coefficient values of 10 and 10 −6 , respectively, whereas Fallon et al. [42] use a value of 1.
Low coefficient values tend to dissipate smaller eddies whereas higher values will dissipate larger circulation systems.The selected value will be dependent upon recirculation dissipation and the mean angular velocity of recirculating flow.Please note that a coefficient value that results in dissipation of eddies smaller than two mesh sizes has negligible impact on the computation.This approach proves to be sufficient when the flow is governed primarily by pressure gradients and advection, such as in tidal flow regimes and particularly for modelling oceanographic scale circulation.However, its suitability within the context of fine scale energy extraction may be limited; therefore, we test its applicability by varying the velocity diffusivity coefficient in the range 10 −6 to 1.This sensitivity test is calibrated against the empirical data (Appendix B).

Depth-Averaged K-Epsilon Scheme
An extended version of the classic k-ε scheme adapted for the Saint-Venant equations is suggested by Rastogi and Rodi [43] based on the vertically integrated RANS equations that calculates dispersion terms from non-uniformity in vertical velocity profiles.The vertically averaged k and ε values are defined as: Hervouet [35] explains how these terms are formulated into production terms using horizontal velocity components and how shear forces in the vertical subsequently affect them.The shear velocity being related to several coefficients based on one-directional flow and a selected friction law, such as Manning's.Telemac experts suggest that the k-ε TCS provides the best means to resolve TECs in model simulations.It should be noted that the default coefficient value (10 −6 ) should always be applied when using this scheme.

Energy Extraction Term
A 2D model allows parameterization of the available energy in a tidal channel when considering seabed mounted turbines deployed under typical constraints.For first generation deployments, this applies to typical devices operating at water depths <50 m, operational velocities <4 m/s and hub heights within the lower 1/5th to 3/5th of the water column to negate changes in vertical shear (Figure 2).The equivalent stress term for each turbine is the sum total of the rotor thrust and structural drag forces applied to the flow (Figure 3).The drag due to the structure of a turbine may be significant (approximately 20% of the total rotor thrust force), depending upon inflow turbulent intensity levels [7,44].Importantly, structural drag is not capped at high velocities, unlike rotor thrust.For simplicity, we assume the support structure of a device to be a single cylindrical vertical monopole for this study, having a structural drag coefficient value (CD) of 0.6 [13,20].Published coastal modelling studies have used an increased (i.e., rougher) value of CD, for example Martin-Short et al. [19] use 0.7, whereas other studies [13,20,41] use 0.9.Previous research has also assumed that biofouling on structures might increase drag by up to 50% [45].However, further research may be required in order to justify such an assumption.
Implementation of simulated energy extraction is based on the generic TEC parameters presented in Table 2, along with the determination of thrust characteristic coefficient (CT) values that vary with change in velocity (Figure 3).Below device cut-in and above cut-out velocity, thrust (or flow backpressure when the turbine rotates and generates electrical power) is zero (Equation ( 5)).Between cut-in and rated velocity, CT is continuous at 0.85 (Equation ( 6)), when transitioning from rated to cut-out velocity, CT follows a polynomial curve fit based on the values presented by Baston et al. [46] (Equation ( 7)).This simulates feathering of the turbine blades to maintain optimal generation and rated power as flow velocity increases.The equivalent stress term for each turbine is the sum total of the rotor thrust and structural drag forces applied to the flow (Figure 3).The drag due to the structure of a turbine may be significant (approximately 20% of the total rotor thrust force), depending upon inflow turbulent intensity levels [7,44].Importantly, structural drag is not capped at high velocities, unlike rotor thrust.For simplicity, we assume the support structure of a device to be a single cylindrical vertical monopole for this study, having a structural drag coefficient value (C D ) of 0.6 [13,20].Published coastal modelling studies have used an increased (i.e., rougher) value of C D , for example Martin-Short et al. [19] use 0.7, whereas other studies [13,20,41] use 0.9.Previous research has also assumed that biofouling on structures might increase drag by up to 50% [45].However, further research may be required in order to justify such an assumption.2, along with the determination of thrust characteristic coefficient (C T ) values that vary with change in velocity (Figure 3).Below device cut-in and above cut-out velocity, thrust (or flow backpressure when the turbine rotates and generates electrical power) is zero (Equation ( 5)).Between cut-in and rated velocity, C T is continuous at 0.85 (Equation ( 6)), when transitioning from rated to cut-out velocity, C T follows a polynomial curve fit based on the values presented by Baston et al. [46] (Equation ( 7)).This simulates feathering of the turbine blades to maintain optimal generation and rated power as flow velocity increases.∅ R The influence of a TEC is approximated in Telemac-2D using a modified version of the DRAGFO subroutine for enhanced bed shear stress calculations altered by Joly et al. [47] to simulate tidal energy extraction.The thrust and structural drag forces, F T and F S (kg m s −2 or N) are: and The total drag force exerted by the TEC on the components of flow velocity in the x and y direction becomes: where U r (m s −1 ) is a reference upstream velocity and θ ( • ) is the angle of orientation of the central axis of the TEC to the flow, having a positive clockwise value when rotated from the x-axis as defined by Joly et al. [47].A T is the turbine rotor swept area and A S is the exposed area of the TEC support structure to flow (Equation ( 10)).When implementing the decelerating force terms (Equations ( 11) and ( 12)) into the model, we calculate the defined area of seabed, A, containing each turbine (i) based on the number of mesh nodes (j) and their associated element size.For this reason, regular grid elements are utilized that facilitate equidistant spacing in order to minimize errors associated with Energies 2018, 11, 2852 9 of 21 mesh density [48].The calculated total stress is thus distributed evenly across the total number of mesh nodes within the defined device area (i.e., 20 m 2 , see Table 1).We thus define the stress term (N m −2 ) to be: To solve these equations, adapted stress term calculations are necessary that include a proportional weighting of device drag relative to changes in h.We therefore adapt Equation ( 13) and ( 14) so that calculations conform to the desired unit base within the subroutine (s −1 ): Both water density and x, y components of velocity (u, v) are applied in the higher-order Navier-Stokes calculations of Telemac.Modification of the associated externally formatted data input file for TEC characteristics must also be performed (Appendix A), where Joly et al. [47] nominally include: The modified input file should include variables containing device dependent cut-in, rated and cut-out velocity, efficiency at cut-in and rated flow and a monopile structure hub height and diameter.Removal of the original simple, fixed power coefficient (C P ) value occurs as this is calculated within the subroutine in a similar manner to C T , i.e., to enable more realistic estimation of time varying outputs with change in velocity.The standard turbine power, P (W) equation is: When turbine efficiency is assumed to increase linearly from cut-in to rated velocity, we define a C P value for each device as: Energies 2018, 11, 2852 10 of 21 where P RAT is the maximum power generated at rated velocity.C P is assumed to be the water to wire coefficient value that includes all losses due to energy transfer (i.e., mechanical, electrical, hydrokinetic etc.).It should be noted that for this study, we utilize the same generic turbine coefficient values for all TECs simulated.In addition, whereas array-averaged simulations exhibit an inherent inability to account for device power alteration in relation to rotor orientation to flow direction, using a device-scale approach does not.Published studies [6,16] have considered the effect that device rotor face offset, or yaw misalignment to flow has in combination with other effects, such as flow asymmetry and ambient turbulent intensity.The influence of device orientation offset will affect available power generated at each rotor due to changes in channel flow power density caused by intra-array turbulence.
A convenient measure of change in available array power is capacity factor (CF), which is the average power generated over a given period, divided by the rated peak power.It is, therefore, not possible to precisely determine CF using an array-averaged approach.

Method Validation
A simple idealized channel case study is analyzed, whereby tidal energy extraction is implemented based on the method presented above.Sensitivity analysis is performed to understand the effects of varying mesh density and TCS on results.Calibration and validation are achieved using published results from scaled physical turbine rotor tank tests [7,[49][50][51][52][53] and CFD modelling studies [54,55].These empirical experiments, conducted at various device scales and using different turbine designs, therefore provide determinant mean centerline velocity deficit and transverse downstream wake profiles.Thus, we quantify the ability to simulate altered hydrodynamics when both mesh density and TCS are altered.
Minimal bottom roughness is applied (Manning coefficient of 0.01) and sidewall friction ignored (free slip) to emulate Perspex flume tank test conditions.The channel is forced at the two boundaries by surface elevation change with a free surface at the input boundary and a fixed elevation at the output (i.e., a constant pressure gradient that drives a constant along-channel flow).Simulated maximum velocity magnitude reaches steady state above cut-in and close to rated velocity.Simulations are performed using 16 processors on a parallel supercomputing platform, with a nominal water density of 1025 kg m −3 .A time step of 0.1 s was chosen based on the Courant stability criterion of the finest mesh (although clearly a larger time step could be used for the coarser mesh simulations).

Axial Wake Centreline Velocity Deficit Profile
Downstream flow in lee of the axial turbine wake centerline is analyzed to produce a longitudinal or along-channel velocity deficit profile (U def ): where U 0 is a measured undisturbed velocity taken at a distance equal to five rotor diameters upstream of the rotor position.U x (x = 1−20) are a series of measurements taken equidistantly one rotor diameter apart in the downstream axial centerline wake (Figure 4).A polynomial fit and nearest neighbour extrapolation is applied to approximate the mean profile from published data (Figure 5).Error estimates (±1 s.d) about the mean provide a realistic range of values across downstream positions.This type of validation accounts for the effects that axial turbine thrust [49], turbine proximity to boundaries and turbulent intensity levels [49,53,56] impart on measurements for various rotor designs (e.g., two or three blades).
rotor diameter apart in the downstream axial centerline wake (Figure 4).A polynomial fit and nearest neighbour extrapolation is applied to approximate the mean profile from published data (Figure 5).
Error estimates (±1 s.d) about the mean provide a realistic range of values across downstream positions.This type of validation accounts for the effects that axial turbine thrust [49], turbine proximity to boundaries and turbulent intensity levels [49,53,56] impart on measurements for various rotor designs (e.g., two or three blades).

Transverse Wake Profile
Transverse profiles (across the wake y plane) of along wake velocity are almost symmetric in nature for this type of assessment and should therefore follow a Gaussian profile.Immediately behind the device, an enhanced bypass flow at rotor peripheries will occur, in line with Bernoulli principles.The velocity decreases from ambient in the lee of the rotor.The velocity deficit is, by nature, at its maximum in the near wake and decreases in the mid to far wake as the flow converges.Velocity profiles from a free shear layer can be normalized, i.e., expressed in self-similar form that follow a Gaussian profile given by [54,55]: ∆ =  −  ( … ) (23) where ∆ represents the transverse velocity deficit profile at any defined downstream position, and is the nominal maximum profile value.The half-width,  / equates to the distance from axial centreline at which the velocity deficit has decayed to half its nominal maximum value.

Transverse Wake Profile
Transverse profiles (across the wake y plane) of along wake velocity are almost symmetric in nature for this type of assessment and should therefore follow a Gaussian profile.Immediately behind the device, an enhanced bypass flow at rotor peripheries will occur, in line with Bernoulli principles.The velocity decreases from ambient in the lee of the rotor.The velocity deficit is, by nature, at its maximum in the near wake and decreases in the mid to far wake as the flow converges.Velocity profiles from a free shear layer can be normalized, i.e., expressed in self-similar form that follow a Gaussian profile given by [54,55]: Energies 2018, 11, 2852 where ∆U represents the transverse velocity deficit profile at any defined downstream position, and is the nominal maximum profile value.The half-width, y 1/2 equates to the distance from axial centreline at which the velocity deficit has decayed to half its nominal maximum value.

Array-Scale Evaluation
Differences in flow flux (m 2 s −1 ) were calculated using a 5 km wide discretized channel domain that incorporates applied energy extraction for various array layouts (Figure 6).The differential channel width integrated flux between an upstream and downstream section spaced a kilometer either side of the array central position is analyzed.This highlights changes caused by intra-array effects that alter the flow through the channel, where:  In each case, device spacing is 2.5 ∅ R laterally and 10 ∅ R longitudinally, as suggested by [5], with the rotor face of each simulated device aligned perpendicular to flow.When implementing an array-averaged approach, simulations are conducted such that Equations ( 15) and ( 16) are adapted; whereby A becomes the equivalent area occupied by all simulated devices based on the suggested device spacing (i.e., dashed polygons in Figure 6).Thus, the resultant enhanced bed stress term is determined from the product of the adapted original term, given the total number of devices simulated, n TECs .It is hypothesized that existential difference in channel flow will be exhibited due to the nuances of the applied methodologies (Figure 7).

Single Device Validation
Mesh resolution dependency is assessed using the k-Ɛ scheme.Velocity deficit profiles downstream of the simulated TEC are compared to the mean wake profile from physical tank test studies (Figure 8).Regression analysis provides coefficient of determination (R 2 ) and root mean squared error (RMSE) values (Table 3).Determinant accuracy of fit ranges by 11% (R 2 ) and RMSE by 1.12%, when altering mesh size.The finest grid discretization provides the most accurate energy extraction determined from axial centerline deficit profiles (R 2 = 0.69, RMSE = 7.16%).Coarser grids

Single Device Validation
Mesh resolution dependency is assessed using the k-ε scheme.Velocity deficit profiles downstream of the simulated TEC are compared to the mean wake profile from physical tank test studies (Figure 8).Regression analysis provides coefficient of determination (R 2 ) and root mean squared error (RMSE) values (Table 3).Determinant accuracy of fit ranges by 11% (R 2 ) and RMSE by 1.12%, when altering mesh size.The finest grid discretization provides the most accurate energy extraction determined from axial centerline deficit profiles (R 2 = 0.69, RMSE = 7.16%).Coarser grids provide a less accurate fit (R 2 = 0.58 and 0.61, RMSE = 8.28% and 8.02%, respectively).Maximum velocity deficit deviation from ambient flow conditions occurs immediately behind a turbine.The mean axial centerline velocity deficit (Figure 5) is approximately 52% one rotor diameter downstream, and 17% at 10 rotor diameters.In comparison, maximum simulated velocity deficit (30-38%) occurs 2-4 rotor diameters downstream (Figure 8).provide a less accurate fit (R 2 = 0.58 and 0.61, RMSE = 8.28% and 8.02%, respectively).Maximum velocity deficit deviation from ambient flow conditions occurs immediately behind a turbine.The mean axial centerline velocity deficit (Figure 5) is approximately 52% one rotor diameter downstream, and 17% at 10 rotor diameters.In comparison, maximum simulated velocity deficit (30-38%) occurs 2-4 rotor diameters downstream (Figure 8).Beyond 10 rotor diameters downstream, available empirical data is scarce; therefore, simulated transverse profiles are not considered in this region.Half-width is found to be ±0.54Y/∅ R .The transverse curve profile exhibits a Gaussian fit, transitioning from decelerated flow at the axial centerline to increased peripheral flow outside ±0.6 Y/∅ R .In general, simulated energy extraction is conservative in nature compared to both the empirical data (Figure 8) and the ideal normalized transverse curve from CFD modelling (Figure 9).Beyond 10 rotor diameters downstream, available empirical data is scarce; therefore, simulated transverse profiles are not considered in this region.Half-width is found to be ±0.54Y/∅ .The transverse curve profile exhibits a Gaussian fit, transitioning from decelerated flow at the axial centerline to increased peripheral flow outside ±0.6 Y/∅ .In general, simulated energy extraction is conservative in nature compared to both the empirical data (Figure 8) and the ideal normalized transverse curve from CFD modelling (Figure 9).

Array-Scale Analysis
Comparison of array-averaged versus device-scale methods reveals the potential to misrepresent flow conditions in simulations when the former approach is used (Figure 10).Effects increase with device thrust as energy extraction increases, during the power generation phase.Peak difference occurs as flow velocity approaches the rated value.For both inline and staggered TEC array configurations analyzed here (Figure 6), the resultant estimated induced simulated peak channel flux error is less than 0.2%.A determinant nominal number of turbines that can be effectively simulated using the array-averaged approach before peak flux difference exceeds 0.5% is thus estimated using linear interpolation and extrapolation.The nominal number of simulated devices, n TECs = 24 for inline turbine rows and n TECs ≈ 26 devices for staggered turbine rows.
When a single TEC is considered in the domain and rotor orientation to flow is altered from zero offset to 10 • and 20 • of misalignment, a reduction in simulated power is induced (Figure 11a).This reduction increases as misalignment to flow increases (0.04% to 1%).When scaled up to any larger array configuration (inline or staggered layout, see Video S1), there is an inherent change in array CF that is dependent upon both device misalignment and array configuration (Figure 11b, Table 4).
difference occurs as flow velocity approaches the rated value.For both inline and staggered TEC array configurations analyzed here (Figure 6), the resultant estimated induced simulated peak channel flux error is less than 0.2%.A determinant nominal number of turbines that can be effectively simulated using the array-averaged approach before peak flux difference exceeds 0.5% is thus estimated using linear interpolation and extrapolation.The nominal number of simulated devices, nTECs = 24 for inline turbine rows and nTECs ≈ 26 devices for staggered turbine rows.When a single TEC is considered in the domain and rotor orientation to flow is altered from zero offset to 10° and 20° of misalignment, a reduction in simulated power is induced (Figure 11a).This reduction increases as misalignment to flow increases (0.04% to 1%).When scaled up to any larger array configuration (inline or staggered layout, see Video S1), there is an inherent change in array CF that is dependent upon both device misalignment and array configuration (Figure 11b, Table 4).When a single TEC is considered in the domain and rotor orientation to flow is altered from zero offset to 10° and 20° of misalignment, a reduction in simulated power is induced (Figure 11a).This reduction increases as misalignment to flow increases (0.04% to 1%).When scaled up to any larger array configuration (inline or staggered layout, see Video S1), there is an inherent change in array CF that is dependent upon both device misalignment and array configuration (Figure 11b, Table 4).

Discussion
Simulated tidal energy extraction at individual device-scale (i.e., representing a typical horizontal axis turbine) is analyzed in a 2D model using an idealized channel domain.Discretization is designed to replicate physical tank test conditions and thus emulate published empirical experimental scaled turbine data for validation.Analysis of centerline velocity deficit and transverse wake profiles are considered.Other modelling studies use similar approaches, for example, Roc et al. [57] compared simulated results to porous disc analysis.How analogous this method is in replicating actual turbine dynamics is unclear.Our results show that the Telemac k-ε TCS provides simulated axial velocity deficit profiles with an accuracy ranging from 0.58 to 0.69 (R 2 ).The subsequent RMSE is 7.16%, 8.28% and 8.02% dependent upon mesh resolution (1, 5 and 10 m respectively).Simulations are based on a 20 m rotor diameter, ensuring that a minimum of two computational nodes fall within the defined device seabed area for the coarse mesh and twenty at finest resolution.
The depth-averaged simulation of near wake dynamics, i.e., within a few meters in lee of a device, exhibits greatest deviation from both empirical values (14-22%) and CFD simulations (Figures 8 and 9).Where the magnitude of vertical velocity increases and becomes similar to horizontal velocity, the hydrostatic approximations are no longer valid.Furthermore, low-order TCS in 2D models fail to effectively resolve complex flow divergence and convergence.This is an inherent restriction for coastal modelling simulations utilizing this technique; Masters et al. [23] report similar findings.This is highlighted in Figure 8, where simulated velocity deficit within two rotor diameters downstream is almost half that expected, based on mean empirical values.However, in the mid to far wake profile the simulated momentum sink improves.At quasi steady state values approximately 20 rotor diameters downstream the velocity deficit is within ±1 s.d for all mesh dependent simulations using the k-ε TCS (Table 3) and is achieved for the CV TCS on a 1 m grid when a velocity diffusivity coefficient value of 0.5 is applied (Appendix B).By utilizing the generic device subroutine formulation and validation methodology presented here, the inherent accuracy of simulations can therefore be implied when upscaling to commercial array scales.Therefore, misrepresentation of phenomena such as wake turbulence in simulations that add to inaccuracies is both minimized and quantified, as both mesh dependency and TCS coefficient can be calibrated and optimized where necessary.
Analysis suggests that intra-array effects remain negligible when seabed shear stress and device structural drag are the contributing factors to flow deceleration.The influence of turbine feedbacks on dynamics is significantly enhanced (approximately four times greater) following device cut-in as power generation begins, and thus channel power density take-off increases.The impact peaks at rated speed (Figures 3c and 10).This implies that turbine thrust (C T ), (i.e., the backpressure exerted on flow) and the associated intra-array device wake interactions that result, greatly influence turbulent mixing and energy dissipation and transfer within the hydrokinetic system.Comparison of the difference between undisturbed inflow conditions to the downstream, disturbed flux for device-scale and array-averaged simulations, illustrates that wake interactions enhance impacts as the number of devices in the array increases.
To consider difference in hydrokinetic system losses, the error induced into simulations when applying array-averaged energy extraction versus the device-scale method presented here was analyzed based on differential along channel flux (Figure 10).In each configuration, simulated devices were aligned perpendicular to the incoming flow direction (i.e., simulating a non-yawing device that has no misalignment to flow).When linear interpolation and extrapolation is applied, the difference in channel flux (and hence flow velocity and subsequent power density) increases to 0.5% when n TECs ≈ 25 devices for inline and staggered rows of turbines.
Power extracted from flow reduces when device rotor face is offset to flow direction [6,16].The ability to orientate the angle of individual devices to ambient flow direction when device-scale approaches are adopted allows models to capture this effect.Results based on CF suggest that an initial 0.04% and 1.00% reduction in power captured by a single device occurs, due to a 10 • and 20 • offset to flow, respectively.This CF change alters when device numbers increase due to intra-array flow changes.The results are configuration-dependent, with an enhanced reduction in CF of 0.18% and 1.71% when five inline devices are simulated and 0.42% and 1.81% for fifteen (Table 4, Video S1).Staggering rows of TECs takes advantage of enhanced flow at the periphery of upstream devices to counter act power availability losses.Configuring turbines in this manner effectively negates (<0.02% for five devices) or improves (−0.11% and −0.45% for fifteen devices) any change in CF.However, ultimately, power take off is related to velocity change via C P .Flow velocity is affected by system feedbacks from energy extraction terms based on drag (C D ) and thrust (C T ) parameterization.Therefore, individual device power capture is independent of the overall associated energy loss.
Given the limitations of the TCS available in Telemac, the approach presented here offers a robust improvement to array-averaged simulations by including intra-array wake effects.Improvements to this study might include extending the simulated transverse and longitudinal velocity deficit profile comparisons to full-scale deployed coastal observations where available.The use of effective validation against empirical observations and CFD simulations is key to improving confidence in results.Implementing this approach when up-scaling to large arrays therefore ensures induced model errors are minimized.The method presented provides a generalized approach that accounts for axial turbine thrust [49], turbine proximity to boundaries and turbulent intensity levels [49,53,56] for various rotor design types.Further improvements necessitate device field measurements, combined with CFD analysis to alter turbulence terms and improve 2D wake dynamics in coastal models.Reliable assessment of environmental resource impacts (e.g., combining sediment dynamics that include morphological changes [12] and wave-tide interactions [58]) can be made using 2D ocean-scale models that include this intra-array wake technique.

Conclusions
The aim of this study was to establish a device-scale approach to energy extraction in depth-averaged simulations that improves upon present techniques.The goal was to create an automated, validated and calibrated methodology, which improves upon array-averaging simulations when implementing energy extraction as an enhanced seabed friction term.This modelling study has implemented energy extraction at individual device-scale into a depth-averaged coastal numerical model (Telemac-2D).It was hypothesized that array-averaged energy extraction simulations do not account for intra-array wake effects and hence misrepresent device-resource interaction and feedbacks.This was investigated by adapting an irregular unstructured model mesh of an idealized tidal channel, capable of resolving individual turbines using a fine resolution regular sub-mesh (10, 5 and 1 m).In doing so, we determine that results are mesh size-dependent; however, the adapted formulae, produces repeatable results.Using the k-ε TCS, a centerline axial wake velocity deficit profile accuracy of 0.58-0.69(R 2 ) and 7.16-8.28%(RMSE) was achieved when compared with empirical, scaled tank tests.If using a constant viscosity TCS, a velocity diffusivity coefficient of 0.5 must be assigned to achieve similar accuracy.
The estimated TEC power captured from the flow alters with device offset to flow orientation and array layout configuration.An increase in CF difference from 0.04% to 0.18% and 0.42% occurs when single rows of five and fifteen turbines are offset by 10 • to the flow.This increases from 1.00% to 1.71% and 1.81% when the offset is 20 • .However, power captured by an array is independent of energy extraction governed by drag and thrust terms, but dependent upon intra-array dynamics and associated velocity manipulation from device-wake feedbacks.Difference in CF is negated or improved (0.02%, 0.01% and −0.11%, −0.45% for five and fifteen devices, respectively) when devices are configured in staggered rows.
The calibrated and validated device-scale energy extraction methodology is compared to an array-averaged approach based on channel flux differential.Various configurations are analyzed, based on generic turbine characteristics and array configurations, with turbines simulated perpendicular to flow (i.e., no misalignment offset).Linear interpolation and extrapolation of these results identifies that the flux differential increases to >0.5% when n TECs ≈ 24-26 devices.This suggests that quantification of the uncertainty associated with utilizing an array-averaging approach, particularly when upscaling to commercial tidal array configurations, may be pertinent.Our results indicate that the uncertainty associated with upscaling must necessarily be established for environmental assessments that negate intra-array effects, as this will exceed 2% error when n TECs ≥100 devices.

Energies 2018 , 21 Figure 1 .
Figure 1.Plan view of the discretized tidal channel.Zoomed section illustrates various mesh density resolutions applied in order to implement a single TEC (solid black polygon) in the center of the domain.The associated system calculated, seabed area (shaded dashed polygon) over which the enhanced stress term is applied is thus element mesh size dependent (Table1).

Figure 2 .
Figure 2. Depth-averaged velocity (dashed vertical line) values typically represent those occurring at 2/5th of the vertical water column, regardless of the applied power law.Even distribution of the profile occurs approximately 1/5th to 3/5th of the vertical water column.

Figure 2 .
Figure 2. Depth-averaged velocity (dashed vertical line) values typically represent those occurring at 2/5th of the vertical water column, regardless of the applied power law.Even distribution of the profile occurs approximately 1/5th to 3/5th of the vertical water column.

Figure 3 .
Figure 3. (a) TEC thrust coefficient curve; (b) Resultant thrust and structural drag force; (c) Formulated normalized stress for the parameterized generic device.

•
number of TECs • center position coordinates of each TEC • the length and width of the grid area occupied by each TEC • orientation of center axis of TEC to x axis • TEC rotor radius • TEC upstream reference velocity distance • TEC structural drag coefficient

Figure 4 .
Figure 4. Plan view of simulated steady state flow and subsequent wake velocity profiles for (a) 1 m (b) 5 m and (c) 10 m mesh density discretization.Direction of flow is shown and the undisturbed velocity, U0 is taken at five rotor diameters upstream.Downstream along channel measurement points, Ux (red dots) are equidistantly spaced at one rotor diameter apart and transverse wake measurements (black dashed lines) are analyzed at 1, 2, 4, 6, 8 and 10 rotor diameters downstream.

Figure 4 . 21 Figure 5 .
Figure 4. Plan view of simulated steady state flow and subsequent wake velocity profiles for (a) 1 m (b) 5 m and (c) 10 m mesh density discretization.Direction of flow is shown and the undisturbed velocity, U 0 is taken at five rotor diameters upstream.Downstream along channel measurement points, U x (red dots) are equidistantly spaced at one rotor diameter apart and transverse wake measurements (black dashed lines) are analyzed at 1, 2, 4, 6, 8 and 10 rotor diameters downstream.Energies 2018, 11, x FOR PEER REVIEW 11 of 21

Figure 5 .
Figure 5. Mean downstream axial centerline velocity deficit profile (solid line) ascertained from published empirical scaled physical tank test modelling studies.Associated standard deviation error bars provide a realistic range of values.

Figure 6 .
Figure 6.Configuration plan providing an indication of seabed area occupied and used for calculations when individual turbines (solid polygons) and array-averaged (dashed polygons) approaches are implemented.

Figure 6 .
Figure 6.Configuration plan providing an indication of seabed area occupied and used for calculations when individual turbines (solid polygons) and array-averaged (dashed polygons) approaches are implemented.

Figure 6 .
Figure 6.Configuration plan providing an indication of seabed area occupied and used for calculations when individual turbines (solid polygons) and array-averaged (dashed polygons) approaches are implemented.

Figure 7 .
Figure 7.An example of simulated channel flow for array-averaged and individual device resolving energy extraction methods using configuration layout E.

Figure 7 .
Figure 7.An example of simulated channel flow for array-averaged and individual device resolving energy extraction methods using configuration layout E.

Figure 8 .
Figure 8. Simulated Telemac-2D along channel axial centerline velocity deficit profiles using the k-ε turbulence closure scheme for various mesh densities.Solid line represents mean empirical test results, grey shading highlights ±1 s.d.(see Figure 5).

Figure 8 .
Figure8.Simulated Telemac-2D along channel axial centerline velocity deficit profiles using the k-Ɛ turbulence closure scheme for various mesh densities.Solid line represents mean empirical test results, grey shading highlights ±1 s.d.(see Figure5).

Table 3 .
Evaluation of k-Ɛ turbulence closure scheme fit to mean observed physical modelling data.

Figure 10 .
Figure 10.(a) Comparison of channel-width integrated flux differential when energy extraction between array-averaged and device-scale configurations (Figure 6) are analyzed; (b) Number of TECs simulated, nTECs, versus peak percentage error in simulations for both single and multiple turbine rows.

Figure 11 .
Figure 11.(a) Power curve of a single generic turbine simulated with zero, 10° and 20° rotor misalignment to flow; (b) CF for inline (circles) and staggered (crosses) array configurations both with and without TEC offset to flow.

Figure 10 .
Figure 10.(a) Comparison of channel-width integrated flux differential when energy extraction between array-averaged and device-scale configurations (Figure 6) are analyzed; (b) Number of TECs simulated, n TECs , versus peak percentage error in simulations for both single and multiple turbine rows.

Figure 10 .
Figure 10.(a) Comparison of channel-width integrated flux differential when energy extraction between array-averaged and device-scale configurations (Figure 6) are analyzed; (b) Number of TECs simulated, nTECs, versus peak percentage error in simulations for both single and multiple turbine rows.

Figure 11 .
Figure 11.(a) Power curve of a single generic turbine simulated with zero, 10° and 20° rotor misalignment to flow; (b) CF for inline (circles) and staggered (crosses) array configurations both with and without TEC offset to flow.

Figure 11 .
Figure 11.(a) Power curve of a single generic turbine simulated with zero, 10 • and 20 • rotor misalignment to flow; (b) CF for inline (circles) and staggered (crosses) array configurations both with and without TEC offset to flow.

Table 1 .
Mesh resolution dependent Telemac geometrical approximation of seabed area compared with user defined area for a single TEC.

Table 1 .
Mesh resolution dependent Telemac geometrical approximation of seabed area compared with user defined area for a single TEC.

Table 2 .
Generic model parameterization used for simulated tidal energy extraction.

Table 3 .
Evaluation of k-ε turbulence closure scheme fit to mean observed physical modelling data.

Table 4 .
Change in capacity factor due to farm array layout and device orientation to flow.

Table A2 .
Evaluation of CV turbulence closure scheme fit to mean observed physical modelling data for various velocity diffusivity coefficient values.