The inventory of observations (Table 1) covers the German Bight and nearby land area (Figure 2). The sites are categorised into atmospheric (“atmos”) sites (measurements of wind speed, U or wind direction, D), “waves” sites (measurements of significant wave height, H _s , peak direction, θ _p or mean wave direction, θ _m ) and ‘ocean’ sites (measurements of sea surface temperature, SST, or water temperature, t _w ). The data sets differ in the periods covered, their temporal resolutions and averaging times. To align the different measurements time-wise, first the buoy measurements of EMODnet (European Marine Observation and Data network, Emodnet, 2020a) and FINO1 and FINO3 (Forschungsplattformen In Nord- and Ostsee 1 and 3, BSH, 2020) are re-indexed to the nearest full 10 min, secondly the measurements are filtered for suspicious values and thirdly all measurements are hourly averaged to align the temporal frequency of the different data sets. The filtering is done following recommended filter criteria (Supplementary Table S1) for the different data sets (Table 1) and afterwards outliers are detected using the IOOS (Integrated Ocean Observing System) QARTOD (Quality Assurance/Quality Control of Real-Time Oceanographic Data) method (IOOS, 2021). After the automatic checks, a manual visual quality control is performed.

We created a base sample using the quality controlled data sets, which is used to ensure that sufficient observations are available for the entire period to validate the simulations. Thus, the base sample should have an availability above 98% during the 30 year period. This criterion was only fulfilled by the 10-m wind measurements at the 12 stations around the German Bight (circles with black edge color in Figure 2), which are therefore used as the basis for base sample. The available dates were filtered so that only dates with at least 22 h on each of two consecutive days at all stations are available. A simulation period of two consecutive days is used to balance computational overhead and sufficient meteorological variability between the days. In total 77% of the 30-year period was included in the base sample. The monthly frequency indicates that the availability after the filtering varies slightly throughout the year with lower availability in summer and higher in winter (Supplementary Figure S1). Nevertheless, the minimum summer availability is still about 70% of the 30-year period and the filtering does not introduce a systematic seasonal bias.

The base sample (Section 2.1.1) is used to select representative days. For that 10, 20, … , 190 random pairs of two consecutive days N(t) = 20, 40, … 380 are randomly selected from the base sample while ensuring that no date t is selected twice. This random sampling procedure is repeated R = 1000 times to capture the scatter of individual samples.

To derive the representativity of the group of selected days (N(t)), we apply the skill score by Perkins (PSS, Perkins et al., 2007) following the approach in Boettcher et al. (2015). The PSS evaluates how well two discrete PDFs, Z ^c (1) and Z N ( t ) s ( i ) , overlap for each bin i: P S S l , t = ∑ i = 1 n min Z c i , l , Z N t s i , t , l . (1)

Here, Z ^c (i, l) is the discrete daily climatic PDF at a certain location l for all hourly measurement times T, i.e. Z c ( i , l ) ≡ Z T s ( i , l ) ≡ Z s ( i , T , l ) . Z N ( t ) s ( i , t , l ) is the daily sample PDF for a certain number of sample days N made up of t days within T. The value of PSS is close to zero for almost no overlap of the two PDFs and approaches one for a perfect match. Since the PSS operates on discrete PDFs, the variables are binned according to Table 2. An example for a climatic PDF Z ^c and sample PDFs Z N ( t ) s for a collection of dates t is shown in Supplementary Figure S2.

Figure 3 shows the combined PSS for all base stations (black circles in Figure 2) and all R = 1000 resamples of different collections of dates t for the different sample sizes N in form of boxplots for (A) U and (B) D. Following Perkins et al. (2007) and Boettcher et al. (2015) we classify an agreement between two PDFs as good for PSS > 0.8 (blue line) and as near-perfect for PSS > 0.9 (green line). The red lines indicate the mean and maximum data availability over all stations. Since Figure 3 only shows the results for the base sample, which was selected based on its high data availability (Section 2.1.1, Table 2), the availability is close to 1.

The representativity of the group of sampled days for the 12 base sample stations combined (Figure 3) increases with increasing number of sample days. However, depending on the randomly selected sample days, the representativity varies considerably, especially for smaller sample sizes. The spread decreases with increasing sample size and for sample sizes with more than 140 days, all samples provide at least a good representativity at all 12 stations.

The selected sample days should represent the wave climate and hub-height wind climate in accordance with the aim of this study. Thus, a PSS is also derived for all other variables and stations using the bins as in Table 2. Only those stations from Figure 2 are taken into account that cover a minimum period within the 30 years (Table 2). The availability thresholds were chosen to maximize spatial and temporal coverage of the stations and variables, while disregarding stations that are not representative for a long-term climate. The PSS is calculated for those stations using a similar equation as Eq. (1) but now only those randomly selected days that lie within the measurement climate period (m) at that station can be used for comparison. Thus, instead of evaluating the agreement between Z N ( t ) s and Z ^c , the agreement is evaluated between Z N ( t ) s and Z m c = Z T m s , i.e. including all hourly measurement times T that fall in m. Thus PSS ^mc is P S S m c t , l = ∑ i = 1 n min Z m c i , l , Z N t s i , t ∩ m , l (2)

As a consequence of using mc, the number of days included in the calculation of PSS ^mc varies 1) for each station, as the measurement climate period varies, and 2) for each resample R, as a different number of sample days N(t) can lie within the measurement climate period. Hence, when interpreting the results of PSS ^mc one has to keep in mind that even if 180 days are sampled, the number of days included to derive PSS ^mc might be smaller. In addition, Z ^mc might not fully represent Z ^c , since the measurement period covers less than 30 years.

The PSS ^mc is shown in Figure 4 for the different variables described in Table 2. The mean and maximum data availability shows the spread of data availability per station and variable. The H _s at 3 hour resolution is used as a separate variable, since those measurements are available for a longer period, as indicated by the higher mean and max availability in Figure 4F compared to Figure 4E.

Similar to the results for the base sample (Figure 3), PSS ^mc increases with increasing number of sample days. However, the number of days required to represent the long-term conditions in a good-to-near-perfect way, differs for each variable. The long term wind climate is better represented close to the surface compared to hub height, especially for D _hub . This is due, among other factors, to the relatively low availability of measurements (50%, red line with dots), which means that in practice only about half of the sample days are used to calculate PSS ^mc . However, even for D _hub , the median is close to near perfect for 180 or more samples. This indicates that the surface and hub height wind climate can be well represented using 180 days.

The wave climate is not as well represented as the wind climate, as indicated by the lower median PSS ^mc . However, especially the 3-h H _s shows near perfect agreements with increasing number of samples. The low mean availability of the data indicates that on average less than half of the number of days could actually be evaluated. However, even under this constrain the wave climate should be reasonably represented, if actually at least 180 days are used.

In conclusion, using a sample of 180 days can well represent the 30 years 10-m wind climate around the German Bight and the long-term climate of other variables around the German Bight.

While most combinations of 180 sample days provide good or near perfect representativity depending on the different variables, the R = 1000 resamples indicate that there is still a lot of spread. Thus, to find the best possible resample for both wind and wave climate, the maximum sum of the PSS ^mc for all variables in Figure 4 for 180 sample days, is used as a criterion. This best overall sample is shown in Figures 3, 4 as a purple star. Supplementary Figure S1 shows that the sample days spread relatively equally over the different months and thus no seasonal bias is evident.

To confirm the representativity of the selected samples, two plausibility checks are performed: Firstly, we use European Centre for Medium-Range Weather Forecasts Reanalysis fifth Generation ERA5 reanalysis data (Hersbach et al., 2018), which can provide a 30 years time series even at stations where the measurement climate is short and provides an opportunity to confirm the spatial pattern of the sampled period. Secondly, we investigate the two-dimensional representativity, i.e. the representativity with respect to variable combinations.

Figure 5A,C show, respectively, the 30-years average wind speed at 100 m and significant wave height based on ERA5 reanalysis. Coloured diamonds show the average U ₁₀₀ and H _s at the different measurement stations based on the best 180 sample days. The figures show that the average climate of both hub-height wind speed and significant wave height can be well represented with the sample days and that the spatial pattern around the German Bight is met. Figure 5B,D shows as boxplot for each station the measurement climate (blue), the measurement distribution from the sample days (orange) as well as the station climate based on ERA 5 (green) and based on ERA5 distribution from the sample days (red). The plot confirms that not only the mean climate is well met by the selected days, but also the distributions as a whole, especially within the inner quantiles. Extreme values cannot be captured with this statistical downscaling method, and those are outside the goal of this study. To capture extremes a different downscaling method was developed by Larsén et al. (2019). The spatial distributions are also well matched for other variables as shown in Supplementary Figures S3–S6.

So far, the selection of representative days was based on individual variables. However, these variables depend on each other. For instance higher waves in the German Bight are usually associated with a wave direction from the open sea and not from the land area. Therefore, the two-dimensional representativity of the selected 180 days should be verified.

Due to a phenomenon, Richard Bellman introduced as “the curse of dimensionality” (Banks and Fienberg, 2003), the number of samples required to replicate a n-dimensional histogram increases exponentially. Also the limiting values of a ‘good’ and ‘near-perfect’ result do not apply in higher dimensions. For a valid assessment of the two-dimensional distribution, the two-dimensional PDF of two variables v ₁ and v ₂, Z v 1 , v 2 , is reduced to one dimension by using the projection of the distribution of one variable v ₁ onto the axis of v ₂: P S S ∑ v 2 m c t , l = ∑ i 1 = 1 n 1 min ∑ i 2 = 1 n 2 Z v 1 , v 2 m c i 1 , i 2 , l , ∑ i 2 = 1 n 2 Z N t , v 1 , v 2 s i 1 , i 2 , t ∩ m , l (3)

Both P S S ∑ v 1 m c and P S S ∑ v 2 m c can be evaluated using the limiting values of a ‘good’ and ‘near-perfect’.

This study focuses on the wind and wave climate, thus, we assess joint distributions for 10-m wind speed, U ₁₀, and direction, D ₁₀, near-surface wind speed, U _NS , and direction, D _NS , (lower mast heights in Table 1), significant wave height, H _s and peak wave direction, θ _p . The same availability filters as in the previous analysis are applied (Table 2) and only stations where both variables are measured are taken into account. This corresponds to 24 stations for U ₁₀ versus D ₁₀ and two stations, FINO1 and FINO3, for the other two comparisons.

Figure 6 shows the PSS ^mc for different aforementioned joint distributions (rows) for both axes (columns). The two-dimensional distributions of U ₁₀ and D ₁₀ (Figure 6A,B), are well met for both axis projections. The joint distribution of U _NS and H _s (Figure 6C,D) is well matched when projected onto the H _s -axis but less so, when projected onto the U _NS -axis. This is due to the relatively coarse binning (Table 2) of the wave height compared to the wind speed. The projection onto the wave axis, i.e. the sum over U _NS , is much narrower and therefore easier to match with the sample days compared to the wider projection onto the wind speed axis (sum over H _s ). The joint distribution for wind and wave direction (Figures 6E,F) are similarly met for both projections. Overall relatively high skill scores are also reached for the two-dimensional distributions based on the available days out of the 180 days and in particular for the chosen sample days (purple star).

The evaluation of different matrices shows that about 180 days should be considered to match the probability density functions of the long-term wind and wave climate in the German Bight. To reduce the computational costs to simulate those days, they were sampled as a 48-h consecutive period. Unfortunately, four 48-h periods out of the 90 periods could not be performed with the set-up described in Section 2.2. Thus, those 8 days could not be taken into account for the analysis in Section 3. However, the blue diamonds in Figures 3, 4, 6 for 172 days shows that the performance in terms of PSS is still comparable to 180 days and thus the simulated sample is still climatically representative.

The 2020 built-out of OWFs clearly has an impact on the available wind resources. Figure 9 shows the simulated wind climate at 100 m height for A) the uncoupled simulation using WRF-only and D) the coupled simulation using the WRF and SWAN models. The relative differences in wind speed (Figures 9B,C,E) due to the presence of wind farms as seen in the different scenarios by subtracting the no-wind-farm scenario from the respective wind farm scenario and normalizing the difference using the no-wind-farm scenario.

Wind speed reductions are on average about 10% or about 2 m s⁻¹ close to the wind farm and are reduced with increasing distance. This is because the wind speed is reduced at the grid cells where turbines are present, but the wind deficits recover with increasing distance though turbulent mixing. Some groups of wind farms act like a wind farm cluster, meaning that there is a continuous wind speed deficit > 1  % linking the individual farms. For FIT three such cluster areas are visible in the German Bight, while for EWP five clusters are simulated.

The magnitude of the wind speed deficit depends on the chosen WFP and the model complexity (coupled or uncoupled). It amounts to a maximum of −21.7% for WRF + FIT, −13.3% for WRF + EWP and −22.1% for WRF + SWAN + FIT. Thus, the maximum wind deficit for EWP is just 61% of the maximum wind deficit for FIT at 100 m height. The lower wind speed deficit for EWP is in line with previous studies (Pryor et al., 2020; Shepherd et al., 2020; Fischereit et al., 2021b). The maximum deficits for WRF + SWAN + FIT is 1% smaller than for WRF + FIT, which is discussed in more detail in Section 3.2.2.

The extent of the wake affected area is also important for wind energy planning, since it indicates the impact of one wind farm on the surrounding farms. Following Pryor et al. (2020), we define the wake affected area A _2% as the area where the velocity deficit is greater than 2%. The respective areas for the scenarios are A _{2%,(WRF+FIT)−WRF} = 9,052 km², A _{2%,(WRF+EWP)−WRF} = 6,488 km² and A _{2%,(WRF+SWAN+FIT)−(WRF+SWAN)} = 9,024 km². Thus, the A _{2%,(WRF+EWP)−WRF} is only 72% of A _{2%,(WRF+FIT)−WRF} . This ratio is slightly larger than the ratio between EWP and FIT that Pryor et al. (2020) estimated for the US mid-west. However, the same trend of larger wake-affected areas for FIT than for EWP is seen. The higher value in our study could be both due to the offshore location or because the simulations in Pryor et al. (2020) do not take into account the bug-fix by Archer et al. (2020).

The influence of the model complexity is smaller on average than the influence of the chosen WFP on the wake affected area. The wake affected areas of both the coupled and uncoupled simulations are comparable, although it is slightly smaller in the coupled simulations. The reason for smaller wake affected area is the higher mixing in the coupled simulations, which is reflected in up to 6% higher values of TKE in 100 m height in the coupled simulations (TKE ₁₀₀ , Figure 10E,F). Those higher TKE ₁₀₀ are a consequence of higher TKE at 10 m (Supplementary Figure S8), which is related to the higher surface roughness lengths calculated in the coupled simulations (Section 3.1).

Defining A _5% and A _10% accordingly to A _2%, results in A _{5%,(WRF+FIT)−WRF} = 3,028 km², A _{10%,(WRF+FIT)−WRF} = 1092 km² and A _{5%,(WRF+SWAN+FIT)−(WRF+SWAN)} = 3,036 km², A _{10%,(WRF+SWAN+FIT)−(WRF+SWAN)} = 1112 km², respectively. Thus, the wake affected area decreases with increasing wind speed deficit threshold more in the coupled simulations than in the uncoupled ones. This agrees with the previous finding that the maximum wind speed deficit is larger in the coupled simulations and will be discussed in more detail in Section 3.2.2.

The climatic mean impact of neighbouring wind farms on wind speed at 100 m (Figure 9) indicates the expected magnitude of losses from other farms. These losses vary with time as shown by the standard deviation of the wind speed difference over the 86 48-means (Supplementary Figure S9). Like the mean difference, the standard deviation of the differences also decreases with increasing distance to the farms.

A t-test is performed to derive where the difference is statistically significant. Supplementary Figure S10B,E,H,K shows exemplary for four sites (Supplementary Figure S10C,F,I,L) the difference distributions between WRF + FIT and WRF along with the derived p-values for the t-test described in Section 2.3. The selected sites illustrate the narrowing of the difference distribution with increasing distance from the farms and the corresponding increase in p-values. Supplementary Figure S11 shows the area of significant differences based on p < 0.01 as dotted hatched area. Regardless of the chosen WFP, the area of significant differences connects all offshore farms in the German Bight. This indicates that already with the 2020 built-out, the wind farms in the German Bight influence each other’s wind resources.

Hub-height wind speed is reduced by on average 1% (or about 0.15 m s⁻¹) above the sea (Figure 11B), with a more realistic representation of wind-wave interaction in the coupled simulations. This is due to the higher roughness length z ₀ in the wind speed range 5 m s⁻¹ to 20 m s⁻¹ calculated through the WBLM in the SWAN model compared to the Charnock-based relationship used in the WRF model (Section 3.1). Thus, compared to the impact of neighbouring wind farms (Section 3.2.1), the mean impact of waves is smaller. Also the standard deviation of the 48-h mean wind speed difference between the uncoupled and coupled simulation (Supplementary Figure S12) is smaller than the difference between simulation with and without wind farms (Supplementary Figure S9). However, the impact of a more realistic representation of wind-wave interaction in the coupled simulations is not concentrated in a particular area, but spread out over the entire model area. Thus, average differences are also significant over larger areas, in particular offshore (Supplementary Figure S13).

Despite the average reduction in Figure 11 being small, there is a huge spread in wind speed difference due to coupling for each 10 min-period. Figure 12 shows the histogram of grid point differences between WRF + SWAN and WRF over the analysed area in Figure 9 for every 10 min-period of the 172 days. Simulated wind speeds between the coupled and uncoupled simulations can differ by more than ±20% or ±2 m/s on a 10 min time scale, respectively. The mean difference between WRF + SWAN and WRF over the entire analysed area is almost zero, but slightly skewed towards negative values due to the increased drag by the waves (Section 3.1). Differences are positive more often onshore than offshore, which means that for the coupled simulations wind resources at hub height are more often increased onshore than offshore. However, those more positive 10-min-fluctuations do not manifest in statistically significant differences of the climatic means onshore (Supplementary Figure S11) and thus should be treated with caution. The relative differences are strongly influenced by outliers due to small wind speeds in WRF, which are used for normalization. Thus, the positive means for the relative differences for “onshore” and “all” are artifacts of those outliers. Medians, being less affected by outliers, are negative for all areas (dotted lines).

Wind-wave-wake interactions lead to an additional reduction in hub-height wind speed on average in the coupled simulations compared to the uncoupled ones in most wind farms (Figure 11D). Although some areas outside of the wind farms also show reductions larger than 1.75%, the additional reduction is most systematic within the wind farms areas. The additional relative reduction is visible throughout the lower atmosphere, as shown by vertical profiles in and outside wind farms for (WRF + SWAN + FIT)-(WRF + FIT) Supplementary Figure S14C. The absolute difference for (WRF + SWAN + FIT)-(WRF + FIT) is larger outside the farms than inside up to 50 m height, when the difference outside the farms becomes larger (Supplementary Figure S14B). The additional reduction in the wind farms stems mostly from the non-linear dependence of the thrust coefficient (c _T ) on wind speed: The slope of the average thrust coefficient over all farms correlates well with the difference in wind speed at 100 m height for (WRF + SWAN + FIT)-(WRF + FIT) between the mean over the center point of all wind farms and the mean over three sites outside wind farms (Supplementary Figure S15). Thus, the slightly smaller wind speeds in the coupled simulations due to the higher surface roughness (Section 3.1) cause higher c _T -values for the coupled simulations compared to the uncoupled simulations.