Time-Dependent Upper Limits to the Performance of Large Wind Farms Due to Mesoscale Atmospheric Response

A prototype of a new physics-based wind resource assessment method is presented, which allows the prediction of upper limits to the performance of large wind farms (including the power loss due to wind farm blockage) in a site-specific and time-dependent manner. The new method combines the two-scale momentum theory with a numerical weather prediction (NWP) model to assess the “extractability” of wind, i.e., how high the wind speed at a given site can be maintained as we increase the number of turbines installed. The new method is applied to an offshore wind farm site in the North Sea to demonstrate that: (1) Only a pair of NWP simulations (one without wind farm and the other with wind farm with an arbitrary level of flow resistance) are required to predict the extractability. (2) The extractability varies significantly from time to time, which may cause more than 30% of change in the upper limit of the performance of medium-to-high-density offshore wind farms. These results suggest the importance of considering not only the natural wind speed but also its extractability in the prediction of (both long- and short-term) power production of large wind farms.


Introduction
The global wind power capacity has increased significantly over the last decade. In the UK, for example, the total wind capacity has grown from 5.4 GW in 2010 to 24 GW in 2019 [1]. Whilst this trend is very likely to continue in the foreseeable future, the upper limit to the growth of wind capacity is still unknown. This is partly because this limit will depend on the growth of our future economy and society in general, which is hard to predict, but also because existing theories and models of wind power generation are still inaccurate, especially for large wind farms.
One of the main difficulties in predicting the performance of large wind farms is to account for the power loss due to the so-called "wind farm blockage" [2]. For example, in 2019, a Danish wind farm developer Ørsted downgraded their prediction of power generation partly because this negative blockage effect had been previously underestimated [3]. Traditionally, wind farm models employed in the wind energy industry assumed that the wind speed upstream of an entire wind farm was unaffected by the presence of the wind farm itself. Now this assumption is known to be incorrect, and several correction methods to account for the wind farm blockage have been proposed in the last few years [4,5]. However, none of the existing correction methods are universal or robust enough to predict the wind farm blockage effect accurately under different conditions. This is essentially because, to correctly predict the wind farm blockage effect, we need to consider combined effects of many influential factors across a very wide range of scales, including not only "turbine-to-array-scale" factors, such as the layout, design and operating conditions of turbines, but also "farm-to-atmospheric-scale" factors, such as the size and location of the farm, nearby farms, and perhaps most importantly, atmospheric/weather conditions [6].
With wind energy at the forefront of the proposed "green economy", being able to predict the performance of large wind farms is essential. For the assessment of wind farm performance, power density is often a more useful metric than total power as it allows for comparison between different sized wind farms. Moreover, for better planning of future energy policies, it is crucial to predict upper limits to the power density of large wind farms. Miller et al. [7] employed a numerical weather prediction (NWP) model with a wind farm parameterization of Fitch et al. [8] over a very large onshore area of 10 5 km 2 in the central United States, and estimated the maximum power density to be about 1 W/m 2 . Although this estimation agrees fairly well with their simple theoretical prediction and some other studies [9,10], typical power densities of existing wind farms are higher than this estimation, mainly because their footprint areas are much smaller. Borrmann et al. [11] conducted a survey of existing wind farm data to empirically predict a realistic value of the output power density (or the "rated" power density multiplied by the capacity factor) of future offshore wind farms in the North Sea and Baltic Sea regions. Their survey suggests that a typical value would be around 3 W/m 2 until 2050, but the exact value would be site specific. More recently, Enevoldsen and Jacobson [12] also conducted a data investigation and reported that the output power densities of existing offshore wind farms in Europe were in the range of 1.15 to 6.32 W/m 2 (with a mean value of 2.94 W/m 2 ).
In this paper, we will present a prototype of a new physics-based approach to predict upper limits to the performance of large wind farms in a site-specific and time-dependent manner, including the wind farm blockage effect. This approach is based on the two-scale wind farm momentum theory proposed recently by Nishino and Dunstan [6] with the use of an NWP model. This theory allows us to predict, for example, the maximum power coefficient of ideal wind turbines in a large wind farm as a function of the effective array density (λ/C f 0 ), bottom friction exponent (γ) and the momentum response parameter for the atmosphere (ζ), as shown in Figure 1. As can be seen from the figure, the maximum power coefficient of ideal wind turbines decreases from the well-known 'Betz limit' (16/27 ≈ 0.593) to a lower value as the array density increases from zero to a larger value. However, the rate of decrease depends significantly on the value of ζ, which represents how strongly the atmosphere can sustain the wind against the resistance caused by the wind farm (as will be explained later). In this paper, we will demonstrate how this timedependent parameter ζ can be estimated using an NWP model. This will then allow us to predict the maximum power density of a given wind farm (with an arbitrary array density) in a time-dependent manner. Maximum power coefficient of ideal turbines (actuator discs) in a large wind farm predicted by the two-scale momentum theory [6] with C * T = 4α(1 − α) and M = 1 + ζ(1 − β). Solid and dashed lines are for γ = 2 and 1.5, respectively.
The significance of this new approach is that it may change the practice of wind resource assessment in the future. Unlike many traditional methods, this approach will allow us to assess the true wind resource that is determined not only by the natural wind speed but also by the quality or "extractability" of wind (at a given place and time).

Theory
The two-scale momentum theory [6] uses the assumption of scale separation between large-scale fluctuations of the atmospheric flow (with periods of more than about an hour) and small-scale fluctuations due to turbulence (with periods of less than a few minutes). This assumption allows us to split the problem of wind-farm aerodynamics into "internal" (turbine/array-scale) and "external" (farm/atmospheric-scale) sub-problems; the former may capture the small-scale fluctuations and the latter the large-scale fluctuations. Even with the scale separation assumption, the two sub-problems are still loosely coupled to each other, and the relationship between the two can be derived from the law of momentum conservation as: where C * T is the (average) local thrust coefficient of all turbines in the farm (defined using the "farm-average" wind speed, U F , instead of an upstream wind speed), λ is the array density, which is the ratio of the total rotor swept area to the farm area, C f 0 is the natural bottom friction coefficient (defined using the natural or "undisturbed" farm-average wind speed that would be observed without the farm, U F0 ), β = U F /U F0 is the farm wind-speed reduction factor, γ is the bottom friction exponent, and M is the momentum availability factor defined as:

M =
Momentum supplied by the atmosphere to the farm site with turbines Momentum supplied by the atmosphere to the farm site without turbines (2) In general, C * T and γ can be obtained from the internal sub-problem (using CFD or reduced-order wake models, for example) whereas M can be obtained from the external sub-problem (using an NWP model), but the two sub-problems need to be solved simultaneously in such a way that the values of β observed in the two sub-problems agree and satisfy Equation (1). However, for the assessment of wind resource (or the maximum power of an "unknown" wind farm) it is convenient to use the following two approximations: where α = U T /U F is the local wind-speed reduction factor, U T is the average wind speed over the rotor swept area, and ζ is the momentum response parameter. The first approximation was originally proposed in [13] using an analogy to the classical actuator disc theory, and this seems to provide a good approximation to the upper limit of C * T for wide ranges of α and λ [6,14]. The second approximation was proposed in [6] with limited information on its validity, and this is properly investigated for the first time in this paper. The point of this approximation is that M is a function of β but ζ is not. This means that now the external sub-problem can be solved to obtain ζ (instead of M) without solving the internal sub-problem simultaneously, to assess the wind resource at a given site. By substituting Equations (3) and (4) into Equation (1), we obtain: which can be solved to obtain β for a given set of parameters α, λ/C f 0 , γ and ζ. If we define the power coefficient of ideal wind turbines (or actuator discs) in the farm as: then, together with Equation (5), we can obtain C P as a function of α, λ/C f 0 , γ and ζ. Hence, by calculating C P for a range of α between 0 and 1, we can find the maximum power coefficient, C Pmax , for a given set of λ/C f 0 , γ and ζ as shown earlier in Figure 1. Note that the exact value of γ is unknown without solving the internal sub-problem for a given wind farm; however, γ = 2 is deemed a good approximation for the prediction of upper limits of power production [13,15]. It is also worth noting that the aforementioned simple theoretical prediction of the upper limit by Miller et al. [7] can be reproduced by the present model with λ/C f 0 → ∞, γ = 2 and ζ = 0 (corresponding to an infinitely large and dense wind farm; see [16] for further details). We therefore assume γ = 2 in this study for simplicity. Once the value of C P(max) has been obtained from the above, we can also calculate the (maximum) power density of a wind farm immediately as (Maximum) power density = (Maximum) total power Wind farm area Hence, if we obtain U F0 , C f 0 and ζ from an NWP model in a time-dependent manner, the above theoretical model, namely Equations from (5) to (7), will provide a timedependent (and array-density-dependent) upper limit to the power density at a given site. As this new approach accounts for the power loss due to wind farm blockage, the maximum power density predicted by this approach is a better indicator of the true wind resource than traditional predictions using U F0 only.
All non-dimensional parameters used in the original and approximated farm momentum equations are summarized in Table 1 and Table 2, respectively. Table 1. Summary of the non-dimensional parameters used in the original farm momentum Equation (1); see [6] for further details.

Parameter Definition
Array density λ = nA/S F Average local thrust coefficient Table 2. Summary of the non-dimensional parameters in the approximated farm momentum Equation (5).

Parameter Definition Note
Array density

Methodology
In this section, the methodology of the study is presented, including the setup of an NWP model and how U F0 , C f 0 and ζ are computed. In this study, we consider an offshore wind farm site in the North Sea as an example to demonstrate how the new approach works.

Setup of NWP Simulations
The NWP model used in this study is the UK Met Office's Unified Model (UM). The simulations were carried out using its nesting suite, which allows a high-resolution local area model to be nested in a lower resolution global model. The global model configuration used is GA6.1 [17] at n768 resolution (∼15 km). The local area model is a 200 km × 200 km domain on a rotated-pole lat-lon grid with uniform horizontal (x and y) grid spacing at ∼1 km resolution. In the vertical (z) direction, 70 grid points are used non-uniformly with the first point at z = 5 m and the 70th at z = 40 km. The center of the local area model is at 56.44035 N, 2.241083 W, corresponding to a region approximately 15 km off the east coast of Scotland. The science configuration for the inner domain is RA1M [18]. The sea surface roughness length z 0(sea) has a squared dependence on the surface friction velocity u * , following Charnock's formula [19]. Following initialization the local model is free-running, meaning it is only constrained by the lateral boundary conditions taken from the global model.
The NWP model is used to run "twin" simulations or a pair of simulations that have identical initial and boundary conditions, with the only difference being that one accounts for the effect of a wind farm and the other does not, as shown in Figure 2. In this study we model the effect of an offshore wind farm by increasing the bottom (sea surface) roughness length for the entire farm area. We employ this simple wind farm parameterization since our main aim is to demonstrate how ζ can be computed from twin NWP simulations to predict the maximum power density in conjunction with the two-scale momentum theory. A more sophisticated farm parameterization, such as [8,20], may yield better predictions of ζ in future studies. We consider a hypothetical offshore wind farm located at the center of the 200 km × 200 km domain. The farm has a circular footprint area (which is expected to simplify the effect of wind direction on ζ) with a diameter of 20 km, covering a large offshore area of 314 km 2 . To put this into context, two existing large wind farms in the North Sea, the London array (630 MW) and Gemini Wind Farm (600 MW) cover areas of approximately 100 km 2 and 68 km 2 , respectively. Comparing the surface areas, the wind farm considered in this study is substantially larger. However, this is still much smaller than the scale of a typical atmospheric circulation, which usually has a horizontal length scale of much larger than 20 km. This farm size is also smaller than the scale above which the maximum power density is expected to start decreasing due to the farm-size effect that is governed by a timescale inversely proportional to the Coriolis parameter [21]. We therefore speculate that the effect of farm size on ζ is relatively minor in this range, although this needs to be investigated further in future studies.

How to Compute U F0 , C f 0 and ζ
To compute U F0 , C f 0 and ζ, we need to decide on the nominal farm-layer height, H F , across which the farm-average wind speed is computed. Choosing the value of H F is similar to choosing a reference wind speed for non-dimensionalization; i.e., we can choose the value of H F arbitrarily as long as the same H F value is used in both internal and external sub-problems to calculate C * T , γ, C f 0 and ζ [6]. In the present study we do not solve the internal sub-problem but employ Equation (3) instead, which means that the value of H F needs to be decided in such a way that Equation (3) is approximately valid. This height is usually about 2-3 times the turbine hub-height (depending on the turbine design and the wind profile) and we use 2.5 times the turbine hub-height as recommended in [6]. We also assume that the turbine hub-height is 100 m, which is similar to many of the offshore wind turbines currently used in the North Sea. Hence, H F = 250 m in this study.
It should be noted that the exact value of H F with which Equation (3) is satisfied depends (only moderately) on the stability condition of the atmosphere [22]. This means that using a fixed H F value would result in an over-or under-prediction of the wind resource at a given time (depending on the stability condition at that time). However, this is not a major concern as the results are not very sensitive to the exact value of H F . For example, we tend to observe only less than 2% of changes in the farm-average wind speeds as we change the assumed value of H F by 10%.
To compute the farm-average wind speeds (U F and U F0 for the cases with and without wind farm, respectively) we first calculate the farm-average wind direction or the "streamwise" direction using a horizontal average over the entire farm area at the turbine hub-height, z = 100 m. Note that this "streamwise" direction may change in time and may also differ slightly between the two cases; hence, this direction needs to be calculated at a given time instant in each simulation. Then U F and U F0 can be computed from a volume average (i.e., average horizontally over the entire farm area and vertically from z = 0 to 250 m) of the "streamwise" velocity in the corresponding simulation. The natural bottom friction coefficient C f 0 can also be computed immediately as where τ w0 is the bottom shear stress averaged over the entire farm area (for the case without wind farm). The momentum response parameter ζ can be obtained from Equation (4) with M and β computed from the twin simulations (at each time instant, as both M and β may change in time). Since the effect of the wind farm is modelled in this study using an enhanced bottom roughness length, we can compute M simply as where τ + w is the farm-average bottom shear stress computed from the simulation with the wind farm model (i.e., this τ + w represents the sum of the actual bottom shear stress on the sea surface and the resistance due to wind turbines). It should be noted that both M and β depend on the roughness length (z 0 ) given to represent the wind farm, and the exact relationship between z 0 and the wind farm conditions (such as the array density and turbine operating conditions) is not known a priori. However, as will be shown in the next section, our simulation results for a wide range of z 0 values suggest that the linear approximation of M given by Equation (4) is satisfactory for most of the cases tested. This means, importantly, that only a single set of twin simulations (with only a single z 0 value to represent an "unknown" wind farm) needs to be conducted to compute ζ approximately; then, together with computed U F0 and C f 0 , we can theoretically predict the power density of any wind farms (i.e., with an arbitrary array density λ and an arbitrary local wind-speed reduction factor α) by solving Equations from (5) to (7).

Results
In this section, we first present NWP simulation results for a selected day with five different wind-farm z 0 values, to investigate the validity of the linear approximation of M and to assess the amount of error in ζ due to this approximation. We then present results for 10 different days with a single z 0 value of 0.1 m, to demonstrate how this approach may be used in future studies to develop a dataset (histogram) of ζ for the assessment of true wind resource at a given offshore site. Finally, we demonstrate, using computed data of U F0 , C f 0 and ζ for two different days, how this approach can be used to predict time-dependent upper limits to the power density of three hypothetical offshore wind farms with low, medium and high array densities.

Linearity of the Atmospheric Response
Five roughness lengths (z 0 = 0.05, 0.1, 0.35, 0.7 and 1.4 m) are tested to investigate the linearity of M against β. As will be shown below, the highest roughness length (z 0 = 1.4 m) tends to yield a low β value of about 0.75, i.e., about 25% reduction in the farm-average wind speed, which could be achieved only by a very dense wind farm. For example, the CFD simulation results presented in [2] for three existing wind farms suggest that the hub-height wind speed reduction in these wind farms may be up to about 16%. Hence, we consider that the range of z 0 tested here is wide enough.
We also test three different data-sampling methods here, referred to as "one-day instantaneous", "two-day instantaneous" and "two-day hourly average" methods. In the first method, only the selected day of interest (23 March 2015) is simulated by the UM local area model (with the initial flow conditions obtained from the UM global model) and instantaneous flow field data are used to compute M and β at a given time instant over the 24 h. This is the simplest method requiring the minimum computational resources, but two potential issues are: (i) there is no spin-up time for non-physical initial effects (such as those due to the wind farm roughness added abruptly at the beginning of the selected day of interest) to disappear; and (ii) instantaneous flow data may contain small-scale fluctuations that should be averaged out for the computations of M and β. Hence, in the "two-day instantaneous" method, we introduce a spin-up time of 24 h, i.e., we start the local area model simulation 24 h ahead of the selected day of interest, and in the "two-day hourly average" method we also perform time-averaging of results to use hourly averaged flow field data to compute M and β. The choice of this filter size (one hour) is based on the assumption that a typical time scale for the atmosphere to respond to the resistance of a large wind farm should be within a range from 10 min to an hour, i.e., smaller-scale fluctuations should be irrelevant to the response characteristics of the atmosphere that we aim to quantify here. Figure 3 shows the time variations of β values computed with the five different roughness lengths, using the "one-day instantaneous" and "two-day hourly average" methods. Overall, the two methods (and also the "two-day instantaneous" method not shown here for brevity) capture the general trend that β decreases as we increase z 0 . However, the "one-day instantaneous" method yields some noticeable abrupt changes in β from 12:00 to 16:00, whereas such abrupt changes are largely eliminated by the "two-day hourly average" method. Moreover, the "two-day hourly average" method yields similar patterns of the time variation of β for all five roughness lengths.
Next, we present in Figure 4 the relationship between the values of (M − 1) and (1 − β) computed using the "two-day hourly average" method. Note that there are five data points at each hour, taken from the five separate simulations with the five different roughness lengths. It can be seen that, at each time instant, the relationship between these two parameters is approximately linear for a wide range of (1 − β) from 0 up to about 0.25 at least. This suggests that the linear approximation model of M given by Equation (4) is useful for most practical cases. As discussed in [6] and confirmed in [23], the main reason why M becomes above 1 is because the pressure difference between the upstream and downstream of the entire wind farm tends to increase as the farm-average wind speed decreases. However, it can also be seen from the figure that the relationship between (M − 1) and (1 − β) is not perfectly linear, especially for a short period of time between 14:00 and 16:00. Hence, we now assess the amount of errors resulting from the linear approximation.   Here the results are presented for the three data-sampling methods, i.e., these errors represent the discrepancies of the values of ζ for z 0 = 0.05 m from the "correct" ζ values for a higher z 0 (computed using the same data-sampling method). For z 0 = 0.1 m, the errors are very small (less than a few percent) for most of the time, especially with the "two-day hourly average" method, which gives a large error only at 15:00 (under-prediction of about 15% as can be visually confirmed from the first two data points for 15:00 in Figure 4). This large error at 15:00 suggests that a certain type of short-time weather event responds non-linearly to different roughness lengths. Nevertheless, these results largely support the validity of the linear approximation of M across this range of z 0 .
For a much higher z 0 value of 1.4 m, the errors are larger even with the "two-day hourly average" method (around 10% for the majority of time, Figure 5). These results suggest that a quadratic approximation (with two model parameters, such as ζ 1 and ζ 2 ) may provide a better description of M with respect to β. However, we anticipate that the linear approximation would be more useful in practical applications due to its simplicity. It should also be noted that these percentage errors in ζ lead to smaller percentage errors in the maximum power coefficient predicted. For example, an error of 10% in ζ tends to result in only a few percent of error in C Pmax (as can be estimated from Figure 1). Figure 6 summarizes 24 h averages of the errors in ζ. It can be seen that the 24 h average error also tends to increase with the roughness length, but the largest error is less than 10%, and the "two-day hourly average" method always yields smaller average errors compared to the other two methods. Hence, we employ "two-day hourly average" method in the rest of the paper.

Histogram of the Response Parameter
Next, we demonstrate how this "twin simulation" approach can be used to obtain a histogram of ζ. Such a histogram, if used together with a histogram of the natural wind speed and the two-scale momentum theory, would allow a better estimate of the annual energy production (AEP) of a large wind farm at a given site. Here we present a mini case study with twin simulation results for 10 different days (i.e., 10 different two-day twin simulations, generating 240 hourly averaged ζ values). These 10 days were selected arbitrarily to cover different types of weather conditions from different seasons of the year, so that the results provide a realistic histogram of ζ. However, for a proper prediction of the AEP, it would be necessary to obtain a histogram of ζ from a much larger set of twin simulations.
In the rest of the paper, we employ z 0 = 0.1 m as a representative roughness length for an "unknown" wind farm and thus compute ζ from twin NWP simulations. It should be noted that the value of z 0 for a "given" wind farm may be estimated using, e.g., an analytical model of Calaf et al. [24], which shows that z 0 = 0.1 m corresponds to a relatively low array density of λ ≈ 0.006 (for a typical thrust coefficient of 0.75). Of importance here, however, is that the current approach using ζ does not require the specification of the wind farm (such as the array density) a priori. In other words, ζ is computed from twin NWP simulations using a certain z 0 value (0.1 m in this case), but the obtained ζ value can be used to estimate the performance of different wind farms (with different z 0 values) due to the linear relationship between M and β as discussed in the previous subsection. Figure 7a shows a histogram of all 240 hourly averaged ζ values obtained. It can be seen that most of the ζ values lie between 5 and 25, but there are some outliers shown by light-green bars (and there is an extreme outlier of ζ = −321, which is outside the range of the graph). All these outliers (10 data points out of 240) were found at rare occasions when β < 0.8 or β > 0.9, whereas a typical β value for the roughness length tested here (z 0 = 0.1 m) is about 0.85 (i.e., about 15% of reduction in the farm-average wind speed). If we remove these outliers using the arbitrary cut-off conditions of β < 0.8 and β > 0.9, the histogram will be as in Figure 7b. Another (slightly less arbitrary) method to remove these outliers is to remove the data points that correspond to low wind speeds, since in reality wind turbines can be operated to produce power only when the wind speed is above their cut-in speed. If we remove the 27 data points at which U F < 5 m/s, for example, the histogram will be as in Figure 7c. There is still one outlier (ζ = −39) but most outliers were found to be removed by this method. This outlier, at which the value of β is 1.07, is a consequence of the interaction of unstable air masses with subtle changes to the regional pressure field induced by the wind farm. This causes complex non-local effects of the wind farm on the regional flow field as shown in Figure 8b. For most cases, the effects of the wind farm are localized, as shown in Figure 8a. In future studies, it would be beneficial to develop a more systematic method to remove such rare anomalies from the data set.
A summary of the three data sets of ζ is presented in Table 3. It can be seen that whilst the medians of the three data sets are very similar (within 0.3 of each other), the means vary slightly more (within 1.6 of each other). The mean value of ζ obtained from the second and third data sets is about 14, which is a little higher than (but still fairly close to) the value of about 11 reported in [25] using Reynolds-averaged Navier-Stokes (RANS) simulations of a steady neutral ABL flow over 625 actuator discs for an idealized offshore wind farm. A more important finding here is that the values of ζ obtained from an NWP model are time-dependent and wide-spread (between about 6 and 24 in this case). Such a large change in ζ (which is much larger than the expected error in ζ due to the linear approximation discussed earlier) results in a substantial change in the maximum performance of wind turbines in a large wind farm (as can be expected from the large gap between the orange and purple curves in Figure 1). This suggests the importance of taking into account the variability of ζ in future wind resource assessments. (a) Commonly observed regional flow patterns (b) Very rare occasion where U F exceeds U F0 Figure 8. Comparison of hourly averaged regional flow fields: (left) For the case with wind farm (center) for the case without wind farm (right) wind speed difference between the two cases. The wind speed (V and V 0 for the cases with and without wind farm, respectively) has been averaged vertically from the lowest model level up to the 7th model level (233 m). Table 3. Comparison of the three data sets of ζ plotted in Figure 7.

All Data Points
Only for 0.8 < β < 0.9 Only for U F > 5 m/s where n is the number of turbines in the farm, and S F is the total area of the farm. These three λ values were selected as most offshore wind farms in the North Sea have an average array density within this range. Even with the highest λ value considered here (0.027), the optimal value of β (that is required to maximize C P and the power density) is found to be around 0.75, which suggests that the linear approximation of M discussed earlier is still valid in this range of λ. Figure 9 shows time-histories of ζ, C Pmax (obtained theoretically for a given set of λ/C f 0 and ζ at each hour), U F0 , and the maximum power density (obtained theoretically for a given set of λ, C Pmax and U F0 at each hour) for the two selected days. In Case A (windy day), ζ changes only moderately (between about 11.5 and 14) and thus the C Pmax of ideal turbines in a given wind farm also changes only moderately throughout the day, i.e., C Pmax depends mostly on the array density ( Figure 9c). The natural (undisturbed) farm-average wind speed U F0 is high throughout the day, varying between about 15 m/s and 19.5 m/s (Figure 9e), and thus the theoretical upper limit to the power density of a given wind farm (with ideal turbines, Equation (7)) also varies in a very similar manner throughout the day (Figure 9g).
Note, however, that these high power density values predicted in Figure 9g are unlikely to be achieved by real wind turbines, mainly because the rated wind speed of real wind turbines is usually lower than the wind speed available in this example case. For reference purposes, three dashed lines in Figure 9g show "rated" power density values for the three hypothetical wind farms operating ideal "Betz turbines" with an assumed rated wind speed of 12 m/s (i.e., power density calculated from Equation (7) with fixed C Pmax = 16/27 and U F0 = 12 m/s). It can be seen that, for the low-and medium-density farms, the rated power density values (about 2 and 8 W/m 2 , respectively) are always below the theoretical maximum power density predicted, implying that the output power would be capped at the rated power, i.e., the capacity factor would be 100% throughout the day. Interestingly, for the high-density farm, the theoretical maximum value becomes lower than the rated value (about 17 W/m 2 ) from 9:00 to 17:00 due to the high array density causing reductions in β and thus C Pmax , demonstrating that the capacity factor of such a high-density farm may not reach 100% even on a very windy day. In Case B (less windy day) the natural wind speed is around 12-13 m/s and does not change substantially (Figure 9f). However, ζ increases substantially throughout the day, from about 8 at 5:00 up to about 14 at 23:00 (Figure 9b). This leads to a noticeable increase in C Pmax towards the end of the day, especially for the high-and medium-density farms, for which the C Pmax values at 23:00 are about 35% and 20% higher than their daily average values shown by dotted lines, respectively. This change in C Pmax also affects the theoretical maximum power density plotted in Figure 9h, where the three dashed lines again show the rated power densities of the three wind farms (same as those plotted in Figure 9g). As can be seen from the figure, the theoretical maximum values for the highand medium-density farms are always much lower than their rated values in this case, and only the low-density farm may achieve a capacity factor of close to 100%. A more important observation here, however, is that the time dependency of the maximum power density is a little different from that of U F0 in this case. For example, the maximum power density for the high-density farm at 23:00 is about 25% higher than that at 17:00 (Figure 9h) whilst U F0 values at 17:00 and 23:00 are identical (Figure 9f), highlighting the importance of taking into account the time dependency of ζ (and thus C Pmax ) in this case.

Discussion and Conclusions
In this paper, we have presented a new type of approach to wind resource assessment using the two-scale momentum theory [6] together with an NWP model. As noted briefly in the introduction, the significance of this new approach is that it allows wind resources to be assessed not only by the natural wind speed observed at a given site, but also by the "extractability" of wind, i.e., how high the wind speed at the farm site tends to be maintained as we increase the number of turbines installed. As described by the two-scale momentum theory, this "extractability" can be predicted by modelling the momentum availability factor M as a function of the farm wind-speed reduction factor β, and here we have shown using NWP simulations that the simple linear approximation model given by Equation (4) seems satisfactory for this purpose.
We have also shown that this "extractability" (represented by the momentum response parameter ζ) seems to vary significantly from time to time, with a typical value of ζ lying between 5 and 25 for an offshore wind farm site in the North Sea. As can be seen from Figure 1, the impact of this variability is significant. For example, for a medium-density offshore wind farm with a typical λ/C f 0 value of 5, the value of C Pmax increases from about 0.28 to 0.37 (i.e., 32% increase) as ζ increases from 10 to 20. It has also been shown that such a large change in ζ may occur even within a day, causing an hourly average C Pmax value to be more than 20% higher than the daily average C Pmax for medium-to high-density wind farms. As there is no strong correlation between ζ and the natural wind speed, the maximum power density may also vary in a different way from the natural wind speed; i.e., it is possible for the maximum power density to increase in time while the natural wind speed stays constant or even decreases. These results highlight the importance of considering the "extractability" of wind in any type of wind resource assessment (including both long-and short-term power predictions) for large wind farms in the future.
In future studies, it would be useful to develop a more reliable histogram of ζ from a much larger set of NWP simulations with a higher-fidelity wind farm parameterization. Although our preliminary study using RANS simulations of a steady neutral ABL flow [25] has shown that the value of ζ computed using a simple roughness model (as in the present study) is similar to that computed using a large array of actuator discs (less than 2% difference in ζ), this should be confirmed with an NWP model for further validation. It would also be useful to investigate how the histogram of ζ changes with the location of the wind farm (e.g., how it changes with the distance from the nearest coastline). Another challenge would be to develop an empirical model of ζ, which would allow prediction of ζ from existing metocean data (such as atmospheric stability parameters and ocean wave heights) without running NWP simulations.
Looking further ahead, the proposed approach based on the two-scale momentum theory may be extended in future studies such that the approach predicts not only the theoretical maximum but also the actual power output of real large wind farms. To do this, it will be necessary to develop a model of the local thrust coefficient C * T for real wind turbine arrays. In the present study we employed a simple analytical model of C * T given by Equation (3) as an approximation to the upper limit [6], but in reality C * T should be a function of the layout, design and operating conditions of turbines as well as the wind conditions, including the wind direction and atmospheric stability conditions. If such a general C * T model is developed to replace Equation (3), the proposed approach would become a powerful tool not only for wind resource assessments as demonstrated in this paper but also for the optimization of design and operating conditions of large wind farms in the future.

Data Availability Statement:
The data presented in this study are available on reasonable request from the corresponding author.