A correction to the enhanced bottom drag parameterisation of tidal turbines

Hydrodynamic modelling is an important tool for the development of tidal stream energy projects. Many hydrodynamic models incorporate the effect of tidal turbines through an enhanced bottom drag. In this paper we show that although for coarse grid resolutions (kilometre scale) the resulting force exerted on the flow agrees well with the theoretical value, the force starts decreasing with decreasing grid sizes when these become smaller than the length scale of the wake recovery. This is because the assumption that the upstream velocity can be approximated by the local model velocity, is no longer valid. Using linear momentum actuator disc theory however, we derive a relationship between these two velocities and formulate a correction to the enhanced bottom drag formulation that consistently applies a force that remains closed to the theoretical value, for all grid sizes down to the turbine scale. In addition, a better understanding of the relation between the model, upstream, and actual turbine velocity, as predicted by actuator disc theory, leads to an improved estimate of the usefully extractable energy. We show how the corrections can be applied (demonstrated here for the models MIKE 21 and Fluidity) by a simple modification of the drag coefficient.

tidal turbines through an enhanced bottom drag. In this paper we show that although for coarse grid resolutions (kilometre scale) the resulting force exerted on the flow agrees well with the theoretical value, the force starts decreasing with decreasing grid sizes when these become smaller than the length scale of the wake recovery. This is because the assumption that the upstream velocity can be approximated by the local model velocity, is no longer valid. Using linear momentum actuator disc theory however, we derive a relationship between these two velocities and formulate a correction to the enhanced bottom drag formulation that consistently applies a force that remains closed to the theoretical value, for all grid sizes down to the turbine scale. In addition, a better understanding of the relation between the model, upstream, and actual turbine velocity, as predicted by actuator disc theory, leads to an improved estimate of the usefully extractable energy. We show how the corrections can be applied (demonstrated here for the models MIKE 21 and Fluidity) by a simple modification of the drag coefficient.

Introduction
One of the key advantages of tidal energy as a renewable energy source, is the predictable nature of the resource. Methods for the detailed prediction of tidal dynamics using hydrodynamic numerical models have been developed over many years and have been applied for many different purposes. Less well understood is how the placement of tidal energy converters in the flow will modify the existing tidal currents at both local and regional scales [1]. The challenge here is that the detailed flow around a turbine is a three-dimensional phenomenon comprising far smaller length scales than those of the underlying tidal resource. A typical approach therefore is to model the turbine scale flow in a three-dimensional CFD simulation based on a actuator disc, blade element, or actuator-line model (see e.g. Sun et al. [22], Harrison et al. [10], Batten et al. [2], Malki et al. [15], Churchfield et al. [3]). The effects of the turbine in a large scale hydrodynamic model are then parameterised, based on properties extracted from the CFD model.
The main property of the turbine that needs to be parameterised is the amount of thrust force exerted by the turbine on the flow (and vice-versa) as a function of the flow speed. This also determines the amount of energy taken out of the flow. Thrust curves typically take the form of a quadratic function of current speed with a non-dimensional thrust coefficient, and can be derived as described above in a high-resolution CFD model, or in lab experiments.
Turbine specific properties such as cut-in and rated speeds however, mean that the curve does not necessarily follow a quadratic and therefore the coefficient is not constant but itself varies as a function of flow speed.
It is important to note, that the turbine properties derived in e.g. a CFD model, or from lab experiments, typically consider the placing of a single turbine in uniform background flow. Speed dependent properties are then expressed in terms of the background velocity, which, because the velocity is slowed down in the presence of a turbine, is available as the undisturbed upstream velocity. In a finite width channel, blockage effects may also affect the resulting thrust curve but can be corrected for (see e.g. Garrett and Cummins [9], Whelan et al. [26]) to derive the thrust curve for an idealised free-standing turbine. In addition, the results may be dependent on the turbulent properties of the flow.
An approach followed in many models is to implement the thrust in the form of an equivalent drag force term. For depth-averaged models this effectively comes down to an increased bottom drag [13,23,19,8,16] Three-dimensional models may implement the drag as a force over the entire water column [5], or if the vertical resolution allows it the drag can be applied over a vertical cross section (e.g. Roc et al. [20]), i.e. an idealised actuator disc.
Since the thrust force is given as a function of the upstream velocity, it is important to consider what velocity to use for the equivalent drag force in the model. One option is to probe the numerical velocity solution somewhat upstream of the turbine location. This however brings with it various difficulties such as the question of how far upstream is appropriate, or the fact that the flow upstream might not actually return to the uniform background flow condition that was considered in the CFD model, due to bathymetric changes or the presence of other turbines. Additionally, the use of a non-local velocity is not desirable for numerical and computational purposes: it makes it hard to treat the term implicitly (in the time-integration sense), potentially leading to time step restrictions for stability, and memory access outside of a fixed numerical stencil, or across sub-domains in domain-decomposed parallel model, is computationally inefficient.
When enough mesh resolution is available, both in the horizontal and vertical dimensions, to resolve the flow through the turbine the relationship between the upstream velocity and the turbine velocity can be predicted using Linear Momentum Actuator Disc Theory (LMADT). Using this relationship the quadratic drag law can be reformulated into a function of the local velocity, thus overcoming the difficulties and ambiguities mentioned above. This is the approach followed in Roc et al. [20]. The typical width of a tidal turbine, order 20m, can however be orders of magnitude smaller than the spatial scales of the tidal flow so that resolving an individual turbine may become prohibitively expensive even in large-scale unstructured mesh models that allow for the efficient focusing of mesh resolution. Also this approach requires the alignment of the mesh with the position and direction of the turbine, thus limiting the flexibility to quickly evaluate different turbine positionings and angles.
If the mesh resolution available is such that computational cells are much larger than the turbine scale, the drag force is necessarily applied over a larger area. In a typical implementation a constant drag is applied over a single cell (the cell that contains the turbine). If the cell size is in fact large enough it may be expected (this will be further investigated in this paper), that the local velocity is not actually affected greatly by the presence of the drag term since the drag force is "smeared" out over a large area and the local cell velocity represents an average of the velocity in a large area around the turbine. In that case the difference between the undisturbed background flow and the local cell velocity may be neglected and the turbine can be implemented using a simple quadratic drag law, function of the local velocity.
As will be shown in this paper however, when the mesh resolution is refined such that mesh distances are closer to the turbine scale, this approximation is no longer tenable as the difference between upstream and local velocity becomes too large. As long as individual turbines are not resolved however, the approach in Roc et al. [20] is also not valid as the local velocity is still larger than the theoretical turbine velocity predicted by linear momentum actuator disc theory.
In particular, for depth-averaged models the local velocity will remain higher than the actual turbine velocity even when the horizontal scales are sufficiently resolved. This is due to the fact that the drag acts on the entire water column and thus the depth-averaged model velocity will represent an average of the actual turbine velocity and a higher by-pass velocity above and below the turbine.
Even in three-dimensional models the drag force is often applied over the entire water column [5], or limited to one or only a few layers [27,11], and does not necessarily give an accurate representation of the actual turbine cross-section and thus the model velocity where the drag is applied is not necessarily equal to the real turbine velocity.
Here we demonstrate how the actuator disc computation may be modified to include the fact that the drag force numerically is applied over a different cross section than the actual turbine. Thus again an analytical relationship can be derived between the undisturbed upstream flow and the local cell velocity, and similarly the drag force can be reformulated as a drag law dependent on the local cell velocity. Like the approach in Roc et al. [20], this leads to a correction to the drag law, which in this case depends on the local cell width, but that nonetheless can easily be implemented in existing models, as will be demonstrated here for the Fluidity and MIKE 21 models.
An alternative method for the parameterisation of turbines in large-scale hydrodynamic models that also makes extensive use of actuator disc theory, is the line momentum sink method [7,21]. Actuator disc theory is used to express the effect of turbines, and more specifically an entire fence of turbines, as a relative head loss across the whole near-field flow pattern starting from the assumed uniform upstream flow at one end , down to the end of individual turbine wakes at the point where uniform flow is again achieved (within the far-field wake of the fence). This head loss is then applied as a jump condition across an edge, or multiple aligned edges within the computational grid using a Discontinuous Galerkin discretisation of the depth-averaged shallow water equations. The advantage of this method that it incorporates a detailed LMADT treatment of the near-field effects, including blockage effects for multiple turbines in a fence.
It does require however that these effects are treated at the sub-grid level, and is therefore only appropriate for hydrodynamic models with grid sizes larger than the length scale of the near-field/turbine wake (typically 10-20 turbine diameters) [7].
For any numerical modelling study it is important to look at the effect of changing the grid resolution on the results of interest. In the modelling guidelines for tidal resource assessments in [14], a range of grid resolutions is rec-ommended depending on the stage of the resource assessment, ranging from kilometre scale for regional studies, down to a range of 500 m to 50 m for specific site feasibility studies. Since the wake of a turbine is a three-dimensional phenomenon, it is not expected that an accurate description of the near-field flow can be obtained with a depth-averaged model. Nevertheless, such models should be capable of studying far-field effects. This relies on the correct forces and their effect on the large-scale flow being modelled correctly. As this paper shows however, the results of the standard enhanced bottom drag parameterisation of the turbine thrust force will deteriorate as the mesh resolution falls below that of the near-field/wake length scale (≈ 200 − 300 m for a typical turbine). The correction proposed in this paper ensures that consistent results can be obtained with grid resolutions smaller than the length scale of the turbine wake, all the way down to the turbine scale.

Enhanced bottom drag formulation
In this section we will describe the enhanced bottom drag parameterisation of turbines used in many models [19,6,16,27] and demonstrate some issues with mesh dependency. We will do this within the framework of MIKE 21 [25], a depth-averaged hydrodynamics model widely used in the marine renewable industry, and an equivalent drag-based implementation in Fluidity, an open source, finite element modelling package [18,12]. By comparing results between the two models we verify that the implementation in the closed source model MIKE 21 is indeed based on the same theory that underlies our implementation in Fluidity, and that the same issues are observed.
The aim of the turbine parameterisation is to represent the drag force of the turbine on the flow, which is typically given as: here u is the flow velocity, ρ the density of sea water, C t the dimensionless drag or thrust coefficient, and A t the effective cross-sectional area of the turbine in the flow. The drag coefficient C t may itself be a function of speed due to turbine properties such as rating, pitch control and the use of a cutin speed. As discussed in the introduction the drag law, often derived from a small-scale three-dimensional CFD model, is typically expressed as a function of the undisturbed background flow velocity, which in the case of an idealised domain corresponds to the uniform velocity upstream of the turbine.
The depth-integrated shallow water equations (in conservation form) are given by where H is the total water depth between bottom and free surface, elevated at a level z = η, u is the depth-averaged velocity, g the gravitational acceleration and c b is the bottom friction coefficient.
A local momentum balance in a fixed local horizontal area A is derived by integrating (2) over this area, multiplied by ρ: The second term represents momentum flux through the boundary ∂A. The third term can be rewritten as an integral of hydrostatic pressure around the three-dimensional water column below A. The last term represents a momentum sink term due to bottom friction.
To implement the turbine thrust force through an enhanced bottom friction, we need the additional momentum sink to be equal to the force, (1). To address the question of which velocity u is used to compute F ( u), in a first attempt we simply employ the local, depth-averaged velocity and average the force over the area A. Thus, we require that Combined with (1) , it readily follows that the enhanced bottom drag coefficient c t in this case should be set to: Since we consider the parameterisation of turbines in hydrodynamic models where mesh distances are larger than the size of an individual turbine, the force is applied over the smallest area possible, typically the area of a single mesh cell.
Thus the area A in (6) corresponds to the cell area over which the enhanced drag coefficient is applied. In models where the cell area is much larger than the turbine cross section A t , the additional drag is small and therefore the presence of the turbine will not have a large effect on the numerical solution for u in that cell. As an example, for typical values of C t = 0.6, a mesh distance ∆x = 200 m and turbine diameter D = 18 m, if the drag is applied over a single square computational cell of ∆x × ∆x, we get which is only half of a typical value of c b = 0.0025 for the background bottom friction coefficient.
Since the effect of the additional drag is relatively small it is to be expected that the assumption that the local velocity within the cell is close to the undisturbed background flow is valid for relatively coarse resolution models, and can therefore be used in the averaged force in the right-hand side of (5). As the resolution is increased however and the mesh distances become closer to the turbine scale, the drag is applied over a smaller area and the reduction in local flow speed may become much larger. Because of the quadratic dependency of the drag force on the flow speed, this may have a significant impact on the force that is applied in the model.

Local velocity drop in idealised channel
We investigate the mesh-dependent reduction in local flow speed in more de-     figure 1 it can be seen that the obtained velocity is indeed highly mesh-dependent, and drops with increased mesh resolution. Since the square of this velocity is used to implement the drag term, a 10% drop in the local velocity leads to an approximate 20% drop in the drag force.
In a model of a fully resolved turbine the local velocity is expected to drop.
After all, the velocity through the turbine is known to be smaller due to momentum exchange with the turbine, whereas the bypass flow around the turbine is expected to accelerate. The deceleration of the flow through the turbine can be estimated using linear momentum actuator disc theory (LMADT, see [9] for an application of this theory to tidal turbines). The theory assumes inviscid flow and a uniform upstream velocity u 0 . Furthermore, it defines a velocity u 1 through the turbine, and velocities u 3 and u 4 in respectively the wake and bypass flow (see figure 2). It also defines pressures: p 0 for the upstream pressure, p 1 and p 2 directly on either side of the turbine, and a uniform pressure p 4 downstream where the velocities u 3 and u 4 are defined. At the same downstream location, the cross-sectional area of the wake flow is defined as A 3 . In addition, it defines the known cross sections A c for the total channel cross section and A t for the turbine cross section.
Through selective application of the continuity equation, momentum conservation and Bernoulli's principle, seven equations can be derived for the un- For our idealised channel case considered above, we may compute u 1 = 2.49 ms −1 . As we will see however in the next section, the difference between upstream and turbine velocity is smaller in the numerical results because the drag force is spread out over an area with a larger width. As already discussed, this means that at very coarse mesh resolutions the velocity in the drag cell is hardly different from the upstream velocity. As the mesh resolution is refined however, and the force can thus be applied over a smaller width, this velocity will drop. This decrease in local velocity will continue with increasing mesh resolution, and only when the resolution is sufficient that the drag force can be applied numerically over exactly the same cross section as that of the turbine, e.g. in a three-dimensional model, should we expect this velocity to have reached the value of u 1 computed from standard actuator disc theory.

Predicting the reduced velocity in the enhanced drag cell
To simplify matters, we first consider the case where the turbine is represented by a square area of enhanced drag instead of a triangle. Additionally, When we neglect variations in the streamwise-direction (here denoted as the x-direction), the model results should correspond to those of an infinitely thin actuator disc as is considered in actuator disc theory. The actuator disc modelled by this shallow water model has a cross-sectional area of ∆yH. Here, and in the rest of the paper, ∆y is the width of the drag area, in the crossstream direction. In this section in particular ∆y = ∆x. Since we consider mesh resolutions where ∆y > D, and additionally H > D, this "numerical" cross section will be much larger than the actual cross section A t . Therefore, to predict the results of the shallow water model using actuator disc theory, we should be careful to use the cross section ∆yH applicable to this model. If we neglect any variations within the horizontal square, we may then hope to predict the velocity within the square as the disc velocity u 1 from this modified actuator disc theory calculation.
Following the assumption made above (5), the magnitude of the force applied in the enhanced bottom drag approximation is given by: Note that here we need to use the actual turbine cross section A t as that is the user input in this formulation to calculate the enhanced drag c t in (6). Further we assume that the velocity that is used to compute the force in this approximation, which is simply the local velocity in the drag cell, will be accurately predicted as the velocity u 1 in the modified actuator disc theory that follows below.
Following the steps in the derivation of (8), (A.22) in the appendix, but now applied to an actuator disc of cross sectionÂ t = ∆yH, we first define a modified thrust coefficient (cf. (A.16) in the appendix): Following the same derivation of (8), we then obtain a relationship between the local model velocity u 1 and the upstream velocity u 0 if in (8) we replace C t withĈ t . This gives an expression for the ratio u 1 /u 0 than can be substituted in (10), to give: After some algebraic manipulation 1 , this can be reworked tô 1 the authors made use of SymPy, a python library for symbolic mathematics: www.sympy.
org Finally, the relationship between the local velocity u 1 within the cell that the enhanced drag is applied in, and the upstream velocity u 0 is given by Figure 3 shows that the speed predicted by (13) closely follows that computed with Fluidity. Note that the results here differ from the Fluidity results in figure 1. This is because in figure 1, the drag is applied over an arbitrary triangle in an unstructured, triangular mesh generated by the mesh generator with a characteristic edge length set to the value of ∆x on the x-axis. The meshes used for the results here, figure 1, are also unstructured, triangular and use the same characteristic edge lengths, but incorporate a square, consisting of two triangles, with dimensions ∆x × ∆x over which the drag is applied. Comparing the two figures, it can be observed that the drag being applied over a square area, aligned with the flow direction, leads to a different relationship between the upstream velocity and the velocity within the drag area, than when the drag is applied over a triangle. For this reason, in the following we will derive two different corrections to the enhanced bottom drag formulation for these two cases.

Turbine correction for square cells
We have shown that actuator disc theory, using the width of a square enhanced drag cell and water depth, can accurately predict the relationship between upstream and local cell velocities. This can be used to reformulate the drag force applied in the cell to be a function not of the local cell velocity, but effectively of the upstream velocity. Instead of applying the force in (9), which is based on neglecting the difference between upstream and local velocity, we want to apply the force where u 0 is the upstream velocity which is not readily (and locally) available.
This expression is the same as in standard actuator disc theory, except now we  The decreasing cell speed can be accurately predicted using (13) derived from actuator disc theory.
need to take into account that this force is not applied over the cross-section A t but over a cross-sectionÂ t = ∆yH. Thus we obtain a modified thrust coefficient:Ĉ Note that this modified thrust coefficient differs from the one in the previous section, used to predict the results in the unmodified enhanced drag formulation, as we now assume that the correct force is applied.
Assuming u 1 is an adequate estimate for the local velocity in the cell with enhanced drag c t and cell area A, the force applied by the enhanced drag is given by (cf. the left-hand side of (5)): After updating (8) to use the modified thrust coefficientĈ t , we can substitute it here to make F a function of the upstream velocity u 0 : To obtain the appropriate value of c t we simply equate this expression with the desired force in (14). This leads to: In comparison with (6) from the standard enhanced bottom drag formulation, we have obtained an additional factor that corrects for the fact that we are using the local cell velocity instead of the upstream velocity. For coarse resolution runs, we have A t /Â t → 0, and thus we fall back, as expected, to the unmodified enhanced drag formulation, since the cell velocity is close to the upstream velocity. As we have seen for finer resolutions, still coarser than the turbine scale, the difference between cell and upstream velocities becomes significant.
The correction derived above can also be applied to three-dimensional simulations with a resolved turbine, where the drag force is applied in threedimensions over a vertical cross-sectional area (actuator disc) withÂ t = A t and thereforeĈ t = C t . The correction factor then simplifies to exactly that given in Roc et al. [20]. For the unresolved case however, both in two and three dimensions, the correction derived here not only corrects for the difference between upstream and turbine velocity, but also for the difference between the actual turbine cross-section and the cross-section over which the drag is applied numerically.
Returning to our idealised channel case, in figure 4 it is shown how the force in the standard enhanced bottom drag formulation applied to a square decreases with increasing mesh resolution. It is to be noted that the relative drop in drag force is larger than the relative drop in speed, due to the quadratic dependency of the force on the speed. Adjusting the drag formulation according to (18), the applied force is not only more accurate at coarse resolution, but also remains much closer to that computed from the upstream velocity directly as the mesh resolution approaches the turbine scale.   (6), the applied force is a quadratic function of the local velocity in the drag cell (here, a square area of ∆x × ∆x). As the mesh resolution increases, the local velocity drops, and therefore the force that is applied within the model decreases. Using the correction in (18) however, the same force can be maintained more or less independent of resolution.

Turbine correction for triangular cells
We now return to the case where the enhanced drag formulation is applied to a single triangular cell, not necessarily aligned in any way with the flow. Again we may approximate the applied drag by an actuator disc spanning the width of the triangle. In this case however, if we thus collapse the applied drag force to a single line, the amount of drag varies along the disc.
We assume here that the streamlines run parallel through the triangle and use a local coordinate system where x is in the streamwise direction and 0 ≤ y ≤ ∆y in the transversal direction, where ∆y is the largest width of the triangle. We may subdivide the triangle into a number of streamtubes of infinitesimal width dy, which can be considered as rectangles ∆x × dy, whose length ∆x = ∆x(y) is a function of y. When approximating this situation with actuator disc type theory, we make the following assumptions: 1. The drag in each streamtube, which in the numerical model is applied over a length ∆x(y), is collapsed in the x-direction and applied at a single point along the streamline, representing an infinitesimal actuator disc with cross section Hdy.

The results in each of the streamtubes are independent of one another.
This means that we take no blockage effects into account and assume laminar flow.
For simplicity we first consider a triangle that is oriented in such a way that it is at its widest at y = ∆y, in other words its top edge is aligned with the streamline at y = ∆y (see figure 5). Furthermore, we have ∆x(y = 0) = 0 in the bottom vertex, and ∆x(y) varies linearly for 0 ≤ y ≤ ∆y. Its area can be computed as A = 1 2 ∆x(∆y)∆y. The function ∆x(y) is therefore given by: The force applied in each streamtube is given by where u 1 (y) is the velocity through the streamtube. Similar to (10), we apply actuator disc theory where we assume that this force is applied over a cross section Hdy and obtain a modified thrust coefficient: Following the same steps as in equations (10)- (13) we may derive the following relation between u 1 (y) and the upstream velocity u 0 : The varying width ∆x(y) thus leads to a variation of the velocity u 1 (y) for 0 ≤ y ≤ ∆y. In the computer models the accuracy of this variation is limited by the numerical approximations employed.
In MIKE, the underlying discretisation is based on a piecewise-constant velocity in each cell. To estimate the cell average obtained in the model we therefore evaluate (22) in the centroid at y = 2 3 ∆y, which gives: For the case where the triangle does not have one of its edges aligned with a streamline, we may consider splitting the triangle into two triangles that share an edge that is aligned along the streamline (see figure 5). The length of this shared edge is the maximum width ∆x max of the triangular drag cell in the streamwise direction. The area of either of the two triangles that the cell is split into, can be computed as A 1,2 = 1 2 ∆x max ∆y 1,2 , where ∆y 1,2 is the height of either triangle. Therefore for each of the triangles we have A 1,2 /∆y 1,2 = actual model where the original, non-aligned triangular drag cell is not split, we can use the same equation (23) for the estimated average velocity of the entire cell as we did for the aligned case.
Using this estimated average, the force applied in the model is then: By equating this to the desired force (1), we may derive a quadratic expression In Fluidity, the P1 DG -discretisation prescribes a linear variation for velocity.
Thus we approximate (22) Equating with the desired force in (1) this time results in a cubic expression for c t : In case the triangular drag cell does not have an edge that is aligned with the streamlines, we may again consider splitting it into two triangles with a shared edge that is aligned with the flow. Here however, (27) does not predict the same linear function for u Fluidity 1 in both triangles, since although A/∆y = 1 2 ∆x max is the same, the value for ∆y in the denominator of y/∆y is different for both triangles, and due to the different orientation of the top triangle, the sign of the gradient of u Fluidity 1 with respect to y will be opposite. The combined piecewise solution is therefore not supported by the underlying discretisation. However, we did find that when using the value of c t found by solving (30), the discrete model gave results that varied only slightly for different orientations of the triangular cell.
The results in figure 6 indicate that again the force applied in the unmodified enhanced drag implementation, in Fluidity and MIKE reduces significantly with increasing mesh resolution. A modification to the enhanced bottom drag c t was derived in this section, solving for c t in (26) and (30) for MIKE and Fluidity respectively, that is shown here to lead to a force that remains close to the desired value. The correction in MIKE was implemented by first finding the value for c t from (26) and then working back from (6) to compute what value of C t should be entered in the GUI to achieve this value in MIKE.

Power production
The correction to the enhanced drag formulation, derived in this paper, is to ensure that the correct amount of momentum is extracted from a shallow water model. This means that the force F applied by the enhanced drag in the drag cell (or region) is an accurate approximation of the real thrust exerted by the turbine on the flow. The amount of energy taken out of the flow within the cell is given by: where u 1,model is the velocity in the enhanced drag cell. As we have seen however, in the case where turbines are not fully resolved this velocity will be larger than the real velocity u 1,turbine that goes through the turbine (as predicted by actuator disc theory). Therefore, the real power production P turbine = F u 1,turbine will be smaller than the amount of power P cell taken out of the model in the drag cell.
This discrepancy can be explained from the fact that part of the mixing losses are not modelled explicitly within the model, but occurs at the sub-grid scale. Following the analysis of Vogel et al. [24], the total amount of power taken out of the flow can be split as follows: where P mixing takes account of the mixing losses due to a.o. shear between the wake and bypass flows. The total power can be computed as [24]: Therefore, as long as the model applies an accurate representation of the thrust force F , using the correction presented in this paper, and an accurate value for the upstream velocity u 0 , the total power extracted from the flow in the model will be accurate as well. The fact that the power P cell extracted within the drag cell, according to (31), is larger than P turbine means that the mixing loss that occurs in the model (outside the drag cell) must be smaller than the real P mixing predicted by actuator disc theory. Therefore part of the mixing loss occurs within the drag cell itself. Thus P cell accounts for both the power P turbine taken out by the turbine itself and additional losses that happen at the sub-grid level.
Vogel et al. [24] considers the case where the drag of an entire farm is smeared out over an enhanced drag region, with the assumption that all mixing losses actually occur within this region. In that case it may be assumed that the total power extraction in the model is a good approximation of the total power extraction predicted by actuator disc theory, so that the available usefully extracted power can be computed as a fraction of that using the same theory.
For the case, considered in this paper, where individual turbines are modelled but are not necessarily fully resolved, part of the mixing losses are modelled explicitly. As argued above however, using the power extracted from the flow by the turbine parameterisation still leads to an overprediction of the usefully extractable energy. It is to be noted that in a shallow water model, even if an individual turbine is resolved in the horizontal mesh, with a minimum mesh distance smaller or equal than the turbine diameter D, the effective cross-section A t = ∆yH will still be larger than the actual turbine cross-section A t . This is because the actual cross-section does not span the entire depth of the water.
Thus, the velocity at the turbine in the model should be interpreted as a depthaveraged velocity that averages between the velocity through the turbine, and the bypass velocity above and below the turbine. This velocity is therefore expected to be higher than the real turbine velocity itself, and therefore the power extraction by the depth-averaged turbine-parameterisation will always be an overprediction of the actual power available to the turbine. The difference between these power values roughly corresponds to vertical mixing losses that are not explicitly modelled in the depth-averaged model. In the next section we will explain how the relationship between the upstream velocity and the local velocity in the model, derived in this paper, can also be used to predict the usefully extractable energy, excluding mixing losses, more accurately.

Implementation details
In this section we summarise, how the analysis derived in this paper can be practically applied in existing models, in order to ensure that the correct force is applied on the flow and an accurate estimate of the available turbine power can be made.

Turbine drag applied over a rectangular area
For models where the turbine parameterisation consists of an enhanced bottom drag applied over a fixed, rectangular area A (e.g. [23]), we may use the analysis presented in section 5. Where existing models typically make no distinction between upstream and local turbine velocity, they calculate the enhanced drag coefficient as c t = C t A t /2A. Such implementations can be improved using the correction given by (18). The extra factor at the end of (18) can easily be included by the user in either C t or A t , without the need for code modification, if these are the input parameters to the model.
An additional complexity arises if C t itself is not a constant. This occurs for example if a cut-in speed and/or rating are applied to the turbine. In this case, C t is typically given as a function (thrust curve) of the upstream velocity u 0 .
In the model however only the local velocity u 1 is available. Using the formula however, it is straight-forward to transform a lookup table that gives the thrust coefficient for different values of u 0 , into a lookup table that is a function of u 1 , by computing u 1 for the given values of u 0 as a pre-processing step.
For the computation of the power available to the turbine, we may use (A.23). Here, again we use (34) to derive the upstream velocity u 0 from the local cell velocity u 1 . Combining these two equations, we derive: Again, in the case that C t is not a constant, a lookup table may be used to obtain the correct value of P turbine for each value of u 1 .

Turbine parameterisation in an arbitrary triangular mesh
For models such as MIKE 21 and Fluidity that employ triangular meshes and which implement turbines through an increased drag applied within a single triangle, the theory presented in section 6 can be applied. In triangular mesh models where the drag force is based on a cell-averaged velocity, the value for the enhanced drag coefficient can be found by solving (26) for c t . Models that use a linear interpolation of velocities stored in the vertices, such as Fluidity should use the value of c t found by solving (30). The same approach could also be followed to implement a turbine in a single drag cell in Telemac 2D, where its Finite Element modus is expected to behave in a similar manner as Fluidity, using a linear representation of the velocity within a cell.
In models, like MIKE, where the applied drag force and the associated coefficient c t are not explicitly prescribed, the same effect can be achieved by modifying the value of C t . This is done by assuming the implementation is equivalent to the standard enhanced bottom drag formulation according to equation (6).
Indeed the results in figure 1 where the standard drag implementation of Fluidity is compared with results in MIKE show that this is true to at least a good approximation. By providing MIKE with a modified value of C t we can therefore create the effect of applying a value of c t obtained from (26) without modifying the code. Note, that in equation (26) we use the original value of C t for the real turbine.
For non-constant C t that is given as a thrust curve, MIKE (and similar models) use the local cell velocity u 1 instead of the upstream velocity to look up the value of C t . This can be corrected by converting the upstream values u 0 in a u 0 → C t look-up table into cell velocities u 1 using equation (23).
To compute the power that can be usefully extracted by the turbine we again use (A.23) this time combined with (23), giving: For finite element models, such as Fluidity, that consider a linear variation of the velocity within the cell we can use (27) which predicts the relationship between the upstream velocity and the velocity in the cell as a function of y.

Support structure
The drag exerted on the flow by the support structure, e.g. pylons or tripods, can typically also be parameterised as a force that depends quadratically on the upstream velocity u 0 : where A s and C s are the cross-sectional area and the drag coefficient of the support structure. To include this drag in the form of an enhanced bottom drag coefficient, we have to deal with the same issue of expressing this force in terms of a local velocity u 1 . Although the theory derived so far can be straightforwardly applied to any force in this quadratic form, we cannot simply derive the enhanced bottom drag coefficient that represents the support drag, denoted by c s , independently of c t and add them up. This is because the local velocity u 1 is the depth-averaged velocity that will be slowed down by the two sources of drag simultaneously.
For the drag parameterisation applied over a square area, the correct value for u 1 is obtained from (34) by usinĝ We can then derive a combined enhanced drag coefficient that represents both turbine and support drag (cf. equation (18)). The effectively extracted power is still given by (35) using the combined value ofĈ t from (40), but only using the values for the turbine itself for C t and A t .
For models where the drag is applied over a triangular cell, the relation between u 0 and u 1 is expressed in terms of the actual enhanced bottom drag coefficient c t . Thus if we include both turbine and support drag in c t , we can maintain equations (23) and (27) for models with cell-wise constant, and piecewise linear velocities respectively. The actual combined value of c t can then be found from the quadratic and cubic equations (26) and (30) respectively, by Finally, the extracted power is still given by (37), using only the turbine values for C t and A t , but using the combined value of c t .

Conclusions
In order to accurately estimate the resource available to tidal turbines and to assess their impact on the hydrodynamics, it is important to accurately represent the drag force exerted by the turbines on the flow. In depth-averaged, and more generally under-resolved hydrodynamic models, one should keep in mind that the local model velocity at the turbine is different from both the upstream and the actual velocity passing through the turbine. The relationship between them is dependent on the mesh resolution, and in the case of depth-averaging, the ratio between the actual turbine cross section and the flow cross section spanning the entire depth. Therefore, although the use of the local velocity for the implementation of the drag force is computationally attractive, it is required to take these relationships into account to avoid spurious and mesh-dependent results. In addition, a better understanding of the relation between local and upstream velocity is necessary for an accurate estimate of the power available to the turbine.
Here we have presented the theory for a single, isolated turbine, and demonstrated that a correction based on linear momentum actuator disc theory taking into account the actual numerical cross section that the force is applied over in the model, can be used to obtain results that are consistent over a range of grid scales. It was shown that the standard enhanced bottom drag formulation results in a drag force that decreases with decreasing grid lengths, in particular when the grid size falls below the length scale of the turbine wake (roughly 10-20 turbine diameters). With the correction the applied force can be kept constant to a large degree, thus ensuring that the effect of the turbine on the large scale flow is correctly modelled.
The analysis for single, isolated turbines may be sufficient for sparsely populated turbine sites which see little interaction between turbines. It is generally recognised however, that in order to achieve the maximum available energy from certain sites, one needs to consider turbine configurations that benefit from local and global blockage effects [17], e.g. fence structures. The analysis in this paper could be extended to include blockage effects. Here again one should make a distinction between the influence of blockage on the relation between upstream and turbine velocities in reality, and the influence of blockage on the relation between upstream and local velocities in the model, in particular taking into account the difference in effective cross sections between reality and the model.
With more closely packed turbines the representation of turbine wake structures and wake recovery also becomes much more important. In addition, the turbulence characteristics may have a great impact on the performance of the turbines. As mentioned in the introduction, depth-averaged models will not be sufficient to accurately model these three-dimensional near-field effects. In further work we would like to explore however, how well these effects can still be approximated in depth-averaged models, possibly through parameterisation and tuning of horizontal turbulence models. Nonetheless, we recognise that in general it may no longer be possible to simply extrapolate from the results of a single isolated turbine, and it may be required to study the effects of combining multiple turbines in detailed three-dimensional CFD calculations and lab experiments.

Acknowledgements
The authors would like to kindly acknowledge the UK's Engineering and

Appendix A. Linear Momentum Actuator Disc Theory
In this appendix we briefly review the main steps in the derivation of the actuator disc theory used in tidal turbine calculations. This is so we can refer to the relevant equations when the modifications, that take into account the numerical implementation details of the enhanced bottom drag formulation, are derived in the main text. These results can be found in e.g. Garrett and Cummins [9], or Whelan et al. [26].
We consider a channel of cross-sectional area A c in which a turbine is located with cross section A t . We assume a uniform flow across the channel upstream of the turbine with velocity u 0 , the flow through the turbine is u 1 . Further downstream we define u 3 to be the velocity in the wake, and u 4 the bypass velocity. Furthermore we assume that at the point down-stream where u 3 and u 4 are defined we have a uniform water level η 4 . The water level upstream is denoted by η 0 , and the water levels just upstream and downstream of the turbine, associated with the pressure drop across the turbine are denoted by η 1 and η 2 .
First we formulate the conservation of mass for the flow through the turbine and in the bypass flow where A 3 is the cross-sectional area of the wake at the location where u 3 is defined. Here we neglect the influence of the water level on the cross sections, so that the cross-sectional area of the bypass flow is given by A c − A 3 . Inclusion of the dependency of cross section on the water level is only significant for high Froude numbers, with details given in [26].
The force F exerted by the turbine on the flow (and vice-versa), can be related to a conservation of momentum principle in the entire channel, or to the pressure drop across the turbine: where g is the gravitational acceleration. Finally, applying Bernoulli's principle along streamlines: 1) from upstream, where u 0 is considered uniform, to just before the turbine, where water level η 1 is defined; 2) from just after the turbine, where water level η 2 is defined, to downstream where a uniform water level η 4 is defined; and 3) in the bypass flow from upstream to downstream. This yields three more equations: Assuming boundary conditions for u 0 and η 4 , and an expression for F as a function of u 0 , we have seven equations for seven unknowns: u 1 , u 3 , u 4 , η 0 , η 1 , η 2 , and A 3 .

General solutions
The Bernoulli equations (A.5) to (A.7) can be rewritten as expressions for water level differences: We can rearrange (A.12) and (A.13) in the following manner, respectively: We introduce the additional definitions, Note that we do not have to assume that F is actually quadratic in u 0 , so that C t is not necessarily a constant; it may still be dependent on u 0 . With these we can derive the following quartic polynomial in u 4 from (A.14) and (A.15): which in combination with (A.12), gives: (A.20)

Zero blockage limit
From the above, it follows that in the limit → 0: u 4 → u 0 and thus η 4 → η 0 .
In this limit, (A.18) becomes The energy yield then becomes: The maximum yield as a function of C t is obtained by: . (A.25) Thus the maximum power (assuming no blockage) is P max = 16 27 · 1 2 A t ρu 3 0 ≈ 0.59 · 1 2 A t ρu 3 0 . (Betz limit) (A. 26)