2.1. Dataset

We use the dataset introduced by Grothe et al. (Reference Grothe, Kächele and Watermeyer2022), comprising 40 years of hourly meteorological and synthetic power data for 29 European offshore wind farms. The meteorological data is derived from ERA5, a global reanalysis dataset produced by the European Centre for Medium-Range Weather Forecasts (ECMWF) that provides consistent hourly estimates of certain atmospheric variables. By combining the ERA5 reanalysis data with turbine-specific power curves, the authors generate realistic hourly power production time series, designed for research purposes. While this dataset relies on generated power curves and grid-averaged reanalysis data rather than high-resolution operational SCADA measurements, it provides a controlled and consistent baseline across all target farms. Avoiding the missing values, sensor errors, and turbine degradation typical of real-world SCADA data allows us to isolate the theoretical efficacy of our proposed transfer learning framework.

2.2. Selecting source wind farms

To capture a comprehensive representation of meteorological patterns across Northern Europe, source wind farms are selected to maximize meteorological and operational diversity. Exposure to such diversity allows the models to learn from a broader range of atmospheric patterns, enhancing their ability to generalize to new sites that have distinct meteorological characteristics.

For each of the 29 available wind farms, a 38-dimensional vector is extracted from data recorded between January 1, 2018, and December 31, 2019. The feature set included wind speed statistics (mean, variance, and percentiles), sea surface roughness distributions, wind directional consistency, power output metrics (including capacity factor and efficiency), and the frequency of extreme weather events. Dimensionality reduction is performed using Uniform Manifold Approximation and Projection (McInnes et al., Reference McInnes, Healy, Saul and Großberger2018), configured to consider 15 nearest neighbors for each data point and with a minimum distance of 0.1. This is followed by agglomerative hierarchical clustering (Ward, Reference Ward1963), which partitions the farms into six distinct clusters based on both their meteorological and operational profiles. From each cluster, one farm is selected to represent the cluster as a source wind farm (Figure 2).

Figure 2.

Location and cluster assignment for available wind farms with the selected source wind farms.

Figure 2. Long description

A two-panel geographical map set against a grid of latitude (46 degrees North to 58 degrees North) and longitude (2.5 degrees West to 12.5 degrees East).

Panel a, titled Distribution of wind farm clusters, shows colored dots representing different groupings.

- In the North, near Scotland, are two dark green dots.

- Along the East coast of England and the central North Sea are clusters of red and yellow dots.

- The English Channel and the West coast of France feature a string of blue dots.

- The Dutch and German coasts show a cluster of dark green dots.

- Further East, near Denmark and the Baltic Sea, are clusters of orange and purple dots.

Panel b, titled Wind farms selected for source models, identifies six specific locations marked with red dots and text labels.

- Beatrice is located in the North off the coast of Scotland.

- Walney is located to the West in the Irish Sea.

- Hornsea is positioned in the central North Sea off the English coast.

- Gode Wind is located off the German coast.

- Baltic Eagle is situated further East in the Baltic Sea.

- Îles d’Yeu et de Noirmoutier is located to the South off the West coast of France.

2.3. Identifying and clustering weather patterns

To identify distinct weather patterns across the six source wind farms, hourly meteorological data from all farms is segmented into continuous, nonoverlapping time periods of length $ p $ hours. Each period $ i $ is represented as a multivariate time series $ {\mathbf{X}}_i\in {\mathrm{\mathbb{R}}}^{p\times d} $ , with $ d=6 $ meteorological features: wind speed, sea surface roughness, wind direction components ( $ \sin \theta $ and $ \cos \theta $ for a wind direction $ \theta $ ), and horizontal and vertical wind velocity components at 100 m above sea level ( $ {u}_{100} $ , $ {v}_{100} $ ). To ensure the forecasting models are trained on the most distinct weather patterns, the time-period length $ p $ and cluster count $ K $ are optimized via grid search, maximizing the separation between identified weather patterns (see Section 2.4).

For a fixed time-period length $ p $ , we train a Variational Autoencoder (VAE) (Kingma and Welling, Reference Kingma and Welling2014) to learn a compressed representation of each weather pattern, which we model as a Gaussian distribution $ {q}_{\phi}\left(\cdot |\mathbf{X}\right) $ with mean $ {\mu}_{\phi}\left(\mathbf{X}\right) $ and variance $ {\sigma}_{\phi}\left(\mathbf{X}\right) $ . It simultaneously learns a decoder, a Gaussian distribution $ {p}_{\theta}\left(\cdot |\mathbf{z}\right) $ with mean $ {\mu}_{\theta}\left(\mathbf{z}\right)=\hat{\mathbf{X}} $ and fixed variance to perform the inverse reconstruction. The encoder $ {q}_{\phi } $ is parametrized with an architecture containing two 1D convolutional layers with kernel sizes of 5 and 3, designed to extract short-term local features. These feature maps are processed by a Bidirectional Long Short-Term Memory (BiLSTM) network, which captures global temporal dependencies by traversing the sequence in both forward and backward directions (Schuster and Paliwal, Reference Schuster and Paliwal1997). The final concatenated hidden states are projected via fully connected layers to parametrize the variational posterior, $ {\mu}_{\phi}\left(\mathbf{X}\right)\in {\mathrm{\mathbb{R}}}^8 $ and $ \log {\sigma}_{\phi}^2\left(\mathbf{X}\right)\in {\mathrm{\mathbb{R}}}^8 $ . The input multivariate time series are thus mapped to latent variables that are an eight-dimensional latent vector $ \mathbf{z}\in {\mathrm{\mathbb{R}}}^8 $ . The latent variables are subsequently sampled via the reparametrisation trick $ \mathbf{z}={\mu}_{\phi}\left(\mathbf{X}\right)+{\sigma}_{\phi}\left(\mathbf{X}\right)\odot \boldsymbol{\varepsilon} $ , where $ \boldsymbol{\varepsilon} \sim \mathcal{N}\left(\mathbf{0},\mathbf{I}\right) $ and $ \odot $ denote element-wise multiplication. The decoder, $ {p}_{\theta}\left(\mathbf{X}|\mathbf{z}\right) $ , performs the inverse reconstruction by employing a unidirectional LSTM for sequence generation. The latent vector $ \mathbf{z} $ is first expanded via a fully connected layer and repeated across the temporal dimension $ p $ to initialize the sequence. This sequence is processed by the unidirectional LSTM to reconstruct the temporal dynamics, followed by transposed convolutional layers to recover the original multivariate dimensions and local structure of the time series $ \hat{\mathbf{X}} $ . The architecture of the model is inspired by Mehmood et al. (Reference Mehmood, Yuan, Gong, Li, Hu, Li and Abbas2026).

The model parameters $ \theta $ and $ \phi $ are optimized by minimizing the $ \beta $ -VAE loss function (Higgins et al., Reference Higgins, Matthey, Pal, Burgess, Glorot, Botvinick, Mohamed and Lerchner2017) with mean squared error reconstruction loss and KL divergence regularization:

$$ \mathrm{\mathcal{L}}\left(\theta, \phi; \mathbf{X}\right)={\unicode{x1D53C}}_{q_{\phi}\left(\mathbf{z}|\mathbf{X}\right)}\left(\parallel \mathbf{X}-\hat{\mathbf{X}}{\parallel}_2^2\right)+\beta \hskip0.1em {D}_{\mathrm{KL}}\left({q}_{\phi}\left(\mathbf{z}|\mathbf{X}\right)\parallel \mathcal{N}\Big(\mathbf{0},\mathbf{I}\Big)\right). $$

Training was run for up to 200 epochs using AdamW optimization (Loshchilov and Hutter, Reference Loshchilov and Hutter2019). To mitigate posterior collapse, we employ a monotonic annealing schedule for $ \beta $ , linearly increasing from $ 0.01 $ to $ 1.0 $ . For inference, each period was deterministically encoded as the variational posterior mean, $ {\mathbf{z}}_i={\mu}_{\phi}\left({\mathbf{X}}_i\right) $ providing a concise representation of the meteorological temporal dynamics. The latent vectors are L2-normalized to prioritize directional similarity over magnitude, then partitioned into $ K $ distinct weather clusters using Hierarchical Agglomerative clustering with Ward’s linkage (Ward, Reference Ward1963) and Euclidean distance.

2.4. Optimal configuration selection

The optimal configuration ( $ p $ , $ K $ ) is selected via grid search by maximizing a composite quality score $ Q $ . Drawing upon the multi-objective framework of Rana et al. (Reference Rana, Koprinska and Agelidis2016), the composite score aggregates six metrics, designed to balance statistical rigor in the latent space with physical interpretability and temporal coherence:

$$ Q=0.2S+0.2\tilde{D}+0.2\tilde{C}+0.2\tilde{M}+0.1T+0.1H. $$

Structural quality is evaluated using the Silhouette score $ S $ (Rousseeuw, Reference Rousseeuw1987), which measures cluster cohesion and separation; the normalized Davies–Bouldin index $ \tilde{D} $ (Davies and Bouldin, Reference Davies and Bouldin1979), which evaluates intra-cluster similarity relative to inter-cluster separation; and the normalized Calinski–Harabasz index $ \tilde{C} $ (Caliński and Harabasz, Reference Caliński and Harabasz1974), which quantifies the ratio of between-cluster to within-cluster dispersion. Physical interpretability is validated using meteorological separability $ \tilde{M} $ , calculated as the log-normalized mean ANOVA F-statistic across wind speed, wind variability, and power output, while robustness is assessed through temporal coherence $ T $ , defined as the average within-cluster autocorrelation (Warren Liao, Reference Warren Liao2005), and distributional consistency ( $ H=1- JSD $ ), which uses the Jennen–Shannon Divergence (Lin, Reference Lin1991) to penalize shifts in cluster proportions between training and testing sets. As detailed in Supplementary Table 4, a sensitivity analysis reveals a trade-off between these objectives: while structural metrics ( $ S $ , $ \tilde{D} $ , $ \tilde{C} $ ) prefer compact representations with fewer clusters, physical and temporal metrics ( $ \tilde{M} $ , $ T $ ) improve with more clusters. To prioritize identifying distinct weather clusters, the composite weights were heuristically calibrated into two tiers. Primary structural and meteorological metrics ( $ S $ , $ \tilde{D} $ , $ \tilde{C} $ , $ \tilde{M} $ ) drive the core clustering behavior and were assigned dominant weights (0.2 each). Meanwhile, temporal coherence and distributional consistency ( $ T $ , $ H $ ) act as secondary regularization terms to prevent fragmented or unstable clusters, receiving lower weights (0.1 each). These generally pointed to similar results when clustering without a composite score, which is why they were combined to agree on final cluster parameters.

2.5. Gaussian process models

Following the assignment of meteorological periods to the $ K $ identified weather clusters, we model the relationship between atmospheric conditions and power generation using cluster-specific Gaussian Process regression (Williams and Rasmussen, Reference Williams and Rasmussen2005). GPs are particularly advantageous for tackling our “cold-start” problem because they operate nonparametrically and are highly data efficient, allowing them to adapt dynamically to limited target data. Crucially, unlike standard deterministic models, GPs inherently generate well-calibrated predictive distributions. For offshore WPF, this uncertainty quantification is an operational necessity, providing grid operators with the reliable risk margins required to securely manage reserve requirements and maintain system stability (Pinson, Reference Pinson2013). For each weather cluster $ k $ , a separate GP $ {f}_k\left(\mathbf{x}\right)\sim \mathcal{GP}\left({m}_k\left(\mathbf{x}\right),{k}_k\left(\mathbf{x},{\mathbf{x}}^{\prime}\right)\right) $ is trained to learn the cluster-dependent power curve, enabling each model to specialize in the distinct dynamics characterizing that specific weather type.

To capture both the overall period trajectory and immediate local conditions, we construct a 20-dimensional input vector $ \mathbf{x}\in {\mathrm{\mathbb{R}}}^{20} $ for each hourly forecast. This vector concatenates the 8-dimensional latent representation $ \mathbf{z} $ with 12 physical predictors: current wind speed and sea surface roughness, autoregressive lags $ \left(t-1,t-2\right) $ for both power and wind speed, and cyclical encodings ( $ \sin, \cos $ ) of the hour, month, and wind direction.

To map this input to the target power output, the covariance function (kernel) must be tailored to the multi-scale physical reality of wind power generation. Standard applications of GPs in WPF frequently utilize the infinitely differentiable Squared Exponential (RBF) kernel, which effectively captures the smooth, theoretical aerodynamic relationship between wind speed and power (Rogers et al., Reference Rogers, Gardner, Dervilis, Worden, Maguire, Papatheou and Cross2020). However, real-world wind power curves contain high-frequency stochastic variations induced by local turbulence and wake effects, for which recent literature highlights the one-differentiable Matérn- $ \nu =3/2 $ is better suited (Prakash et al., Reference Prakash, Tuo and Ding2023). To explicitly capture both phenomena, we construct an additive composite kernel. The RBF component isolates the smooth, macroscopic baseline of the ideal power curve, while the Matérn- $ \nu =3/2 $ component captures the rougher, unpredictable atmospheric fluctuations:

$$ k_k(\mathbf{x}, \mathbf{x}') = \sigma_{f,k}^2 \left( k_{\mathrm{RBF}}(\mathbf{x}, \mathbf{x}'; \boldsymbol{\ell}_k^{\mathrm{RBF}}) + k_{\mathrm{Mat\acute{e}rn}}(\mathbf{x}, \mathbf{x}'; \boldsymbol{\ell}_k^{\mathrm{Mat}}, \nu = 1.5) \right), $$

where $ {\sigma}_{f,k}^2 $ controls the output variance. Both components utilize Automatic Relevance Determination (ARD), where $ {\mathbf{\ell}}_k\in {\mathrm{\mathbb{R}}}^{20} $ represents a vector of length scales. This allows the GP to automatically determine the relevance of each input feature for each specific cluster, for example, weighting lags more heavily in stable conditions versus turbulence features during storms.

Observation noise is modeled via a Gaussian likelihood with variance $ {\sigma}_{n,k}^2 $ , assuming $ y={f}_k\left(\mathbf{x}\right)+\varepsilon $ , where $ \varepsilon \sim \mathcal{N}\left(0,{\sigma}_{n,k}^2\right) $ . Consequently, the covariance of the noisy observations is given by $ {\mathbf{K}}_{y,k}={k}_k\left(\mathbf{x},{\mathbf{x}}^{\prime}\right)+{\sigma}_{n,k}^2I $ . The complete hyperparameter set $ {\theta}_k=\left\{{\sigma}_{f,k}^2,{\sigma}_{n,k}^2,{\mathbf{\ell}}_k^{\mathrm{RBF}},{\mathbf{\ell}}_k^{\mathrm{Mat}},{c}_k\right\} $ is optimized by maximizing the exact marginal log-likelihood (MLL) of the observations. Using the noisy covariance $ {\mathbf{K}}_{y,k} $ , the MLL is defined as:

$$ \log p\left(\mathbf{y}|X,{\theta}_k\right)=-\frac{1}{2}\left({\mathbf{y}}^{\top }{\mathbf{K}}_{y,k}^{-1}\mathbf{y}\right)-\frac{1}{2}\log \mid {\mathbf{K}}_{y,k}\mid -\frac{N}{2}\log \left(2\pi \right). $$

2.6. Transfer learning

To evaluate the cross-site generalization capability of the proposed framework, we apply a transfer learning strategy where the library of GP models is adapted to eight unseen target farms. This adaptation utilizes the fixed VAE encoder $ {\phi}_S $ and the set of cluster-specific GP priors, proceeding in three stages: latent projection, cluster alignment, and GP fine-tuning. First, the target meteorological time series $ {\mathbf{X}}_t $ are projected into the source latent space using the fixed source encoder $ {\mathbf{z}}_T={\mu}_{\phi_S}\left({\mathbf{X}}_T\right) $ . This enforces that target weather patterns are interpreted through the learned representation of the source domain. Each target period is then assigned to one of the $ K $ source weather clusters via a nearest-neighborhood classification in the latent space, utilizing the centroid from the source training set.

For each identified weather pattern $ c\in \left\{1,\dots, k\right\} $ , we initialize a target GP $ {f}_T^{(c)} $ using the hyperparameters of the corresponding source GP, $ {f}_S^{(c)} $ . This initialization provides a strong prior on the temporal covariance structure of the power output. The kernel hyperparameters $ {\theta}_k $ are then optimized by maximizing the exact marginal log-likelihood on the available target training data $ {\mathcal{D}}_{\mathrm{train}} $ :

$$ {\hat{\theta}}_T=\underset{\theta }{\mathrm{argmax}}\log p\left({\mathbf{y}}_{\mathrm{train}}|{\mathbf{X}}_{\mathrm{train}},\theta, {\sigma}_{n,\mathrm{fixed}}^2\right). $$

To assess data efficiency, we conducted experiments varying the data availability of the 2-year target datasets. Each target dataset is split into training fractions $ \gamma \in \left\{\mathrm{0.1,0.2,0.3,0.4,0.5}\right\} $ , simulating a “cold start” scenario where a newly commissioned wind farm has been operational for only a short period (ranging from ~2.4 to 12 months). The remaining $ 1-\gamma $ is reserved for testing.