# Improving Subseasonal Wind Speed Forecasts in Europe with a Nonlinear Model

Ganglin Tian, Camille Le Coz, Alexis Tantet, Anastase Charantonis, Naveen

Goutham, Riwal Plougonven

► To cite this version:

Ganglin Tian, Camille Le Coz, Alexis Tantet, Anastase Charantonis, Naveen Goutham, et al. Improving Subseasonal Wind Speed Forecasts in Europe with a Nonlinear Model. Monthly Weather Review, 2025, 153 (9), pp.1761-1779. 10.1175/MWR-D-24-0248.1. hal-05471373

HAL Id: hal-05471373

https://hal.science/hal-05471373v1

Submitted on 22 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Improving sub-seasonal wind-speed forecasts in Europe with a non-linear model

Ganglin Tian,ª Camille Le Coz,ª Anastase Alexandre Charantonis,ª,ᵇ Alexis Tantet,ª Naveen Goutham,ª,ᶜ and Riwal Plougonvenª

ª LMD/IPSL, École Polytechnique, Institut Polytechnique de Paris, ENS, Université PSL, Sorbonne Université, CNRS, Palaiseau, 91120, France

ᵇ INRIA, Paris, France

ᶜ EDF R&D, Palaiseau, France

<figure><img src="image_2.png" /></figure>

Corresponding author: ganglin.tian@lmd.ipsl.fr

1

**Early Online Release:** This preliminary version has been accepted for publication in *Monthly Weather Review*, may be fully cited, and has been assigned DOI 10.1175/MWR-D-24-0248.1. The final typeset copyedited article will replace the EOR at the above DOI when it is published.

© 2025 American Meteorological Society. This is an Author Accepted Manuscript distributed under the terms of the default AMS reuse license. For information regarding reuse and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).
ABSTRACT: Sub-seasonal wind speed forecasts provide valuable guidance for wind power system planning and operations, yet the forecast skills of surface winds decrease sharply after two weeks. However, large-scale variables exhibit greater predictability on this time scale. This study explores the potential of leveraging non-linear relationships between 500 hPa geopotential height (Z500) and surface wind speed to improve sub-seasonal wind speed forecast skills in Europe. Our proposed framework uses a Multiple Linear Regression (MLR) or a Convolutional Neural Network (CNN) to regress surface wind speed from Z500. Evaluations on ERA5 reanalysis indicate that the CNN performs better due to its non-linearity. Applying these models to sub-seasonal forecasts from the European Centre for Medium-Range Weather Forecasts, various verification metrics demonstrate the advantages of non-linearity. Yet, this is partly explained by the fact that these statistical models are under-dispersive since they explain only a fraction of the target variable variance. Introducing stochastic perturbations to represent the stochasticity of the unexplained part from the signal helps compensate for this issue. Results show that the perturbed CNN performs better than the perturbed MLR only in the first weeks, while the perturbed MLR's performance converges towards that of the perturbed CNN after two weeks. The study finds that introducing stochastic perturbations can address the issue of insufficient spread in these statistical models, with improvements from the non-linearity varying with the lead time of the forecasts.
# 1. Introduction

As renewable energy capacity, particularly wind energy, continues to grow, skillful sub-seasonal wind speed forecasts become important for pricing, production, transmission, and utilization of renewable energy resources. Energy producers depend on skillful sub-seasonal wind speed forecasts to plan and adjust operations of power plants, ensuring that wind turbines are active during optimal wind conditions and scheduling maintenance during low wind periods to prevent outages when demand is high, thereby enhancing efficiency in power production (Tawn et al. 2022). Grid operators require skillful sub-seasonal wind speed forecasts to maintain a stable power supply, especially in regions heavily dependent on wind power, by optimizing the integration of various energy sources, such as other renewable and fossil fuels, to meet consumer power demand (Cassola and Burlando 2012; Chang et al. 2014; White et al. 2017). Energy traders rely on skillful wind forecasts to estimate wind power availability and power demand, so that information on likely variations in wind speed and temperature contributes to anticipating likely changes in energy prices.

The predictability of surface variables like near-surface wind speed is generally lower than that of large-scale atmospheric circulation patterns on the sub-seasonal timescale. This is because large-scale circulation modes are closely linked with ocean-atmosphere interactions, planetary waves, and large-scale energy transfers (Wallace and Gutzler 1981; Zhang 2005), which evolve more slowly and are less sensitive to small-scale disturbances than surface variables. Existing global Numerical Weather Prediction (NWP) models (Hersbach et al. 2020; Haiden et al. 2023; WMO 2012) have limited ability to resolve small-scale processes that affect surface variables. Additionally, despite significant advancements in the parameterization of small-scale phenomena, uncertainties persist within these parameterization schemes (Hersbach et al. 2020; Haiden et al. 2023). The initialization of variables also presents a challenge: while large-scale atmospheric fields can be initialized accurately using direct satellite and radiosonde observations, the initialization of surface variables is often limited by the inconsistent quality and incomplete global coverage of in situ surface data (Hersbach et al. 2020; Bauer et al. 2015; Haiden et al. 2023).

Leveraging the predictable information in large-scale variables can significantly improve the skill of sub-seasonal forecasts for surface variables (Mariotti et al. 2020; Vigaud et al. 2017; Büeler et al. 2020). In-depth analysis of the relationship between large-scale phenomena and surface
wind speed has revealed both linear and non-linear statistics. Alonzo et al. (2017) and Goutham et al. (2023) have developed and refined methods such as Principal Component Analysis (PCA) and Redundancy Analysis (RDA) to capture these complex interactions between the geopotential height at 500 hPa (Z500) and surface wind speed in a linear way, to improve the forecast skill of wind speed. However, the relationships between large-scale variables and surface variables are often complex, containing both linear and non-linear components (Salameh et al. 2009; Wilby and Wigley 1997; Chen et al. 2011; Maraun et al. 2010; Sachindra et al. 2018). Neural networks effectively handle non-linear relationships in high-dimensional meteorological data, capturing complex interdependencies between large-scale and surface variables without requiring specific data distribution assumptions (Maraun et al. 2010; Sachindra et al. 2018; Pan et al. 2019; Rodrigues et al. 2018). This advantage is particularly valuable for meteorological applications, where data often follows complex, non-normal distributions. However, the primary focus has been on surface temperature and precipitation (Wilby and Wigley 1997; Sachindra et al. 2018; Pan et al. 2019; Rodrigues et al. 2018), with less attention given to surface wind speed. Capturing these complex relationships accurately is a potential approach to improving wind speed forecast skill. Yet, even with accurate representation of these relationships, the intrinsic uncertainty in atmospheric conditions on the sub-seasonal timescale poses ongoing challenges that need to be addressed. To address forecast uncertainty, meteorologists employ ensemble forecasting techniques. Sub-seasonal forecasting systems often suffer from under-dispersion, where ensemble spread inadequately captures the actual forecast uncertainty (Robertson et al. 2015; Bi et al. 2022; Kurth et al. 2023; Chen et al. 2023; Orth and Seneviratne 2014). Zhu et al. (2018) demonstrated that improving uncertainty quantification and representation can significantly enhance forecast skill.

While substantial progress has been made in medium-range forecast up to two weeks (Price et al. 2023; Rasp et al. 2023; Lam et al. 2022) and seasonal forecast up to 13 months (Alonzo et al. 2017), yet the sub-seasonal timescale, spanning three to six weeks ahead, presents distinct challenges as both initial atmospheric conditions and slowly evolving boundary forcings need to be considered. Advancing sub-seasonal wind speed forecast skill requires both a deep understanding of the complex relationships between large-scale variables and surface wind speed, and the accurate quantification of ensemble uncertainty. To improve sub-seasonal wind speed forecasts, this study aims to explore both the linear and non-linear relationships between large-scale atmospheric variables and surface
wind speed, while also addressing the under-dispersion of sub-seasonal forecasting ensembles.
Our research is guided by two key questions:

* Can we leverage non-linear regression relationships between large-scale variables and surface wind speed to improve sub-seasonal forecast skill over Europe, surpassing the performance of linear models?

* How can we effectively represent ensemble dispersion when applying the regression relationships to dynamical forecasts on sub-seasonal scales?

This study is structured as follows: section 2 introduces our methodology, starting with our framework, followed by a description of our linear and non-linear regression models and their training strategies. In section 3, we delve into the specifics of the cases studied, including the rationales behind the choice of input and target variables and data sets. Section 4 presents the performance of our regression models on a historical dataset. Subsequently, in section 5, we conduct probabilistic ensembles, applying the trained regression models to each ensemble member to assess whether the chosen non-linear model can outperform the linear model in improving forecast skill for wind speed on the sub-seasonal timescale. Finally, section 6 offers a summary.

## 2. Methodology

### a. Framework

We introduce a two-stage framework (as in figure 1) to improve the sub-seasonal forecast skill of a surface variable, incorporating a) training and validation of a regression model on historical data (reanalysis, see section 3b) and b) applying the regression model to ensembles. In the training and validation phase, as in figure 1a, we aim to establish a mapping between a large-scale variable $X$ and a surface variable $Y$ using the statistical information from historical deterministic data. On the training dataset, the inputs $X_{train}$ and the targets $Y_{train}$ are pre-processed using identical data pre-processing procedures (see appendix, section a). Subsequently, optimal parameters for the regression model are obtained by minimizing a loss between the estimated outputs $\hat{Y}_{train}$ and the targets $Y_{train}$ (see section 2b). Additionally, the differences between $\hat{Y}_{train}$ and $Y_{train}$ are used to derive the distribution parameters of their residuals (see section 2c).
<figure><img src="image_3.png" /><figcaption>FIG. 1. Our framework for improving surface-scale variable forecast skill comprises two main stages: (a) Training stage: Both deterministic inputs X and deterministic targets Y undergo identical data pre-processing (appendix, section a) before being used to train a regression model. The outputs of the training stage comprises the optimal model coefficients of the statistical model (Sec b) and the parameters of the residual distribution (Sec c). (b) Ensemble forecasting stage: The m-th member $X^m$ of an input ensemble, after pre-processing, is fed into the trained model to regress the corresponding regressed ensemble member $\hat{Y}^m$. Then we randomly sample P times from the residual distribution from the training stage to perturb the regressed member $\hat{Y}^m$, in order to obtain P perturbed members $\tilde{Y}^{m \times P}$ for this single regressed member.</figcaption></figure>

In the ensemble forecasting stage in figure 1b, in order to produce M-member forecasting ensembles $Y^M$ of the target variable, the trained regression model is applied to M-member forecasting ensembles $X^M$ of the input variable. The m-th member $X^m$ of the input variable independently feeds into the regression model to yield the corresponding regressed m-th member $\hat{Y}^m$. The outputs, M-member ensembles produced by the regression model, are termed the regressed ensembles $\hat{Y}^M$.

However, the regression model captures only a portion of the variance of the target, as indicated by Mean Square Error (MSE) (detailed further in section 4). Consequently, when the model is applied to an ensemble dataset, the outputs exhibit reduced variance relative to the expected variance of the target. This leads to a systematic under-estimation of dispersion and variability of the target members. To maintain the total variance of the target and to attribute skill improvements more accurately to the representation of the predictable components of the signal, rather than to the vari-
ance of the unpredictable components, we employ a perturbed version of the model. For the *m*-th regressed output member $\hat{Y}^m$, *P* perturbations are randomly drawn from the residual distribution to perturb the $\hat{Y}^m$. This process produces *P* perturbed members for *m*-th member, denoted as $\tilde{Y}^{m\times P}$ (refer to section 2c).

In this study, the results obtained through the regression models are termed “statistical,” while those based solely on numerical weather prediction models are termed “dynamical.” To differentiate between model outputs when using reanalysis versus ensemble members as inputs, the term “prediction” is avoided. The outputs based on reanalysis are referred to as “regression outputs,” and the outputs deriving from ensemble members are called “regressed ensembles”. We use superscripts to differentiate between deterministic and probabilistic data. *X* and *Y* denote deterministic input and target dataset, respectively, while $X^M$ and $Y^M$ represent the input and target datasets comprised of *M* members, with the *m*-th member expressed as $X^m$ and $Y^m$. The notation “^” is employed to indicate values that are regressed. The ensembles and the models in which perturbations have been included are indicated with a tilde “~”. Lastly, all physical variables (such as wind speed and geopotential height) are represented as fields on latitude-longitude grids.

## b. Model architectures

As previously mentioned, Alonzo et al. (2017) and Goutham et al. (2023) used the linear relationships based on PCA and RDA to regress surface wind speed from Z500. This study explores whether the non-linearity between large-scale variables and surface wind speed contributes to further improving the ensemble forecast skill of surface wind speed, compared to linear models, rather than determining the most skillful model structure. Therefore, we opted for a simple Multiple Linear Regression model (MLR) to capture the linear relationship between inputs and a target. In our case, we developed the MLR using multiple grid points of Z500 over an input domain (more details can be found in section 3a) to regress U100 at a specific grid point.

Acknowledging the inherent spatial patterns in Z500 (Goutham et al. 2023), the sea-land distribution differences in U100, and the non-linearity from Z500 to U100 (Salameh et al. 2009), we adopted a Convolutional Neural Network (CNN) model to capture this non-linearity. CNNs, through multiple layers of convolution operations and non-linear activation functions, effectively learn spatio-temporal patterns and non-linear relationships within data. Building on this, we
adopted SmaAt-UNet (Trebing et al. 2021), a variant of UNet (Ronneberger et al. 2015). Compared to the standard UNet, which has been widely applied in the field of meteorology (Nguyen et al. 2022; Bouget et al. 2021), SmaAt-UNet integrates attention modules and depthwise separable convolutions, significantly reducing the number of parameters while maintaining model performance. The breadth of this model family is related to CNN architectures, and it is uncertain whether more intricate non-linear relationships could yield additional explanatory power. However, our experiments indicate that within the constraints of the available training set, augmenting complexity through the integration of additional layers and channels does not enhance the MSE of the regressed deterministic U100. The architectures of MLR and CNN, and their training configurations are detailed in appendix, section d.

### c. Stochastic forecasts

The reduction in variance associated with the deterministic regression models significantly affects the assessment of probabilistic ensembles. The skill of these ensembles, as quantified by Continuous Ranked Probability Skill Score (CRPS; see section 3d), is intrinsically linked to their reliability and includes their capacity to capture the variance (as evaluated by Spread Skill Ratio (SSR; see section 3d)). This lost variance from the deterministic regression does not reflect a decrease in forecast uncertainty; rather, it is of a purely statistical nature during the training stage in figure 1a. Consequently, one approach to preserve the variance that represents dynamical uncertainty is the introduction of stochastic perturbations. This stochastic perturbation approach should account for the portion of the variance that the deterministic regression model fails to explain. In other machine-learning-based ensemble forecasting studies (Bi et al. 2022; Kurth et al. 2023; Chen et al. 2023; Orth and Seneviratne 2014), their statistical models did not account for these systematic errors, potentially leading to the under-dispersion of ensembles. We represent the unexpressed fluctuations as residuals, which are the differences between $Y_{train}$ and $\hat{Y}_{train}$. These residuals, representing the uncertainty in the model's fit to a deterministic training dataset as systematic errors, are modeled with a Gaussian distribution to estimate their mean and variance as stochastic perturbations in figure 1a. As introduced at the beginning of this section, to express this uncertainty during the ensemble forecasting phase in figure 1b, for each regressed member $\hat{Y}^m$, P perturbations are randomly sampled from this distribution and then added to $\hat{Y}^m$ to obtain
$P$ perturbed members $\tilde{Y}^{m \times P} = \{\hat{Y}^{(m,1)}, \dots, \hat{Y}^{(m,P)}\}$. Our experiments (not shown) indicate that the ensemble skill converges when $P = 20$.

### d. Scores and significance test

We evaluate the regressed outputs $\hat{Y}$ and the ensembles $\hat{Y}^M$ and $\tilde{Y}^M$ against their corresponding reference $Y$. All scores are evaluated grid-point by grid-point separately and later spatially averaged with cosine-latitude weights over a domain of interest. The same approach is followed with respect to the lead times. During the regression validation phase in figure 1a, we evaluate the ability to capture the relationship between Z500 and U100 on reanalysis, quantified by the MSE between the ground truth $Y_{\text{validation}}$ and the regressed outputs $\hat{Y}_{\text{validation}}$ on the validation set of a historical reanalysis. Subsequently, in the ensemble forecasting phase in figure 1b, we measure the forecast skill of the regressed ensembles $\hat{Y}_{\text{test}}^M$ and the perturbed ensembles $\tilde{Y}_{\text{test}}^{M \times P}$ at each independent lead time by MSE of ensemble mean and CRPS (Wilks 2019; Matheson and Winkler 1976; Zamo and Naveau 2018). Additionally, ensemble reliability is quantified by SSR (Rasp et al. 2023). MSE quantifies the errors in deterministic forecasts, while CRPS assesses the accuracy of probabilistic forecasts by evaluating the alignment between ensemble distributions and reanalysis, while SSR measures the relative dispersion of ensembles. The formulas for these metrics are detailed in appendix, section b.

To assess the significance of the relative improvements of a statistical model $model_s$ exhibited over a benchmark $model_b$, for a given score, we employ a bootstrap technique similar to that described by Goddard et al. (2013). This involves randomly selecting identical samples of ensembles from $model_s$ and $model_b$, with replacement, to generate a new subset of samples, with as many samples as in the original samples. These samples are then used to calculate the aforementioned verification scores. We compute the relative difference in scores $\Delta_r Score = \frac{Score(model_s) - Score(model_b)}{Score(model_b)} \times 100(\%)$ between the two models and repeat this process 1000 times. For negatively oriented scores, such as MSE and CRPS, the proportion of $\Delta_r Score < 0$ serves as the $p$-value. Conversely, for positively oriented scores, the proportion of $\Delta_r Score > 0$ is used as the $p$-value. If the $p$-value is less than the predetermined significance level $\alpha$, the $model_s$ is considered significantly better than $model_b$ ($p < \alpha$).
## 3. Case and Data

### a. Case study

We apply the above framework to a specific case study. Here, the input variable is the geopotential height at 500 hPa, Z500, over the Europe-Atlantic domain (20°–80°N, 120°W–40°E). The output variable is the wind speed at 100-meter, U100, over the European region (34°–74°N, 13°W–40°E). These regions of interest are shown in figure 2.

#### 1) TARGET VARIABLE:

To study wind speed at hub height on a sub-seasonal timescale across the Europe domain in figure 2, we focus on U100, as in Goutham et al. (2023), obtained from the zonal (u) and meridional (v) components of 100-meter wind using the formula: $U100 = \sqrt{(u^2 + v^2)}$.

#### 2) INPUT VARIABLE:

The ensemble forecast skill of wind speed declines rapidly with increasing lead time; the U100 ensembles are skillful for less than about ten days (Buizza and Leutbecher 2015; Goutham et al. 2022). When selecting input variables, we aim to choose those that remain skillful within a two to six-week time scale. Z500, which reflects the large-scale atmospheric circulations in the mid-troposphere and is commonly associated with atmospheric fluctuations and 500 hPa geostrophic winds, is extensively used in the studies examining the regression relationship between large-scale variables and wind speed (Alonzo et al. 2017; Goutham et al. 2023; Liu et al. 2023). Moreover, the Z500 has higher predictability on the sub-seasonal timescale than the 10-meter wind (Toth and Buizza 2019). Buizza and Leutbecher (2015) demonstrated that Z500 remains more skillful than climatology up to approximately lead times of 22 days. In this study, Z500 is cropped into the Europe-Atlantic domain as depicted in figure 2. This specific domain is selected because the target domain for U100 is Europe, where downstream dependency of forecast errors typically occurs at mid-latitudes (Simmons and Hoskins 1979). Consequently, the domain extends further westward relative to the target domain than eastward to capture its large-scale circulations.
<figure><img src="image_4.png" /><figcaption>FIG. 2. The two domains of interest: the Europe-Atlantic domain (20°–80°N, 120°W–40°E) for the input variable Z500 and the Europe domain (34°–74°N, 13°W–40°E) for the target variable U100.</figcaption></figure>

### 3) THE SEASON OF INTEREST.

We validate our methodology during winter. This decision is based on the fact that Z500 exhibits the most significant anomalies during winter (Buizza and Leutbecher 2015). Additionally, the correlations between atmospheric circulations and surface variables vary seasonally and are strongest in this season (Laurila et al. 2021). Moreover, the higher frequency of low-pressure systems during winter leads to large wind speed variability, posing challenges for the energy sector and necessitating skillful estimations of U100.

## b. Data

### 1) HISTORICAL DETERMINISTIC DATASET.

We employ the ECMWF Reanalysis v5 (ERA5) high-resolution reanalysis dataset as our reference. Reanalysis combines historical meteorological observations with NWP models to generate physically consistent descriptions of past atmospheric states, widely serving as a reference in ensemble forecasting (Hersbach et al. 2018a,b). We source Z500, the u and v components of 100-meter wind reanalysis data spanning from December 1979 to March 2022, with a spatial resolution of 2.7 degree and temporal resolution of 6 hours from Climate Data Store (CDS). We then downsampled the Z500 and U100 reanalysis data to a weekly temporal resolution during the pre-processing step (discussed further in subsequent paragraphs).

### 2) ENSEMBLES DATASET.

We use the ECMWF extended-range hindcasts (Vitart et al. 2019) as the dataset for ensemble forecasting. Hindcasts, generated through NWP models for retrospective periods, help calibrate
and assess forecasting models' performance. The ECMWF hindcasts are initialized from ERA5 reanalysis data and generated with a control member and ten perturbed members under slightly different initial atmospheric and oceanic conditions with stochastic parameterizations. In this study, we use only these ten perturbed members, which are expected to sample part of the errors that grow from uncertain initial conditions and part of the uncertainty associated with model errors.

Besides hindcasts, ECMWF also provides forecast ensembles, which are estimates of the atmospheric state over specific future periods based on current observations and NWP models. Unlike hindcasts, forecasts are initialized from analysis and consist of 50 perturbed members and a control member. ECMWF extended-range forecasts are available from 2015 and thus cover a shorter period than the hindcasts. This study primarily uses hindcasts for verification owing to their extended temporal coverage, encompassing a range of meteorological phenomena, including extreme weather events and seasonal variations, thereby facilitating a comprehensive evaluation of model performance. Unless stated, the dynamical ensembles mentioned hereafter refer to hindcasts, while the results for forecasts are also provided in appendix, section g. The forecasts and hindcasts of Z500 and U100 were obtained from ECMWF via Meteorological Archival and Retrieval System (MARS). We downloaded the hindcasts initialized in December, January, and February (DJF) from 1995 to 2021 and the forecasts in DJFs from 2015 to 2021 from MARS, covering 128 initializations for forecast ensembles and 2560 initializations for hindcast ensembles. The Integrated Forecasting System (IFS) of ECMWF has undergone several updates during these periods; however, the statistical differences between these various versions are expected to be marginal (Goutham et al. 2022). In this study, we employ hindcasts only up to March 2022, due to data availability.

### 3) SPATIO-TEMPORAL RESOLUTION:

Predictability depends on spatial and temporal scales (Vitart 2014). Small-scale features fluctuate more frequently than large-scale features, making them less predictable over sub-seasonal timescales (Vitart and Robertson 2019). A simple method to reduce unpredictable noise is averaging data over time and space (Buizza and Leutbecher 2015). Here, bilinear interpolation is applied to achieve a desired resolution of 2.7×2.7 degree (approximately 300 km) for both the U100 and Z500 ensembles, leading to (22×59) grid points for Z500 over Europe-Atlantic domain and (14×19) grid points for U100 over Europe domain. These ensembles are averaged weekly to
<figure><img src="image_5.png" /><figcaption>FIG. 3. Nested cross-validation process consisting in two layers: the outer layer enclosed in solid red lines and the inner layer in dashed red lines. In the outer layer, each rectangle represents a period of 9 years. The gray rectangles indicate the training reanalysis folds, the blue rectangles signify the test reanalysis folds, and the rectangles with blue hatching depict the test ensembles' folds. In the inner layer, each rectangle spans 3 years, with gray rectangles representing the training reanalysis folds and light yellow rectangles indicating the validation reanalysis folds.</figcaption></figure>

focus on the potentially predictable parts of the signal. The reanalysis are also aggregated at the same spatio-temporal resolution to maintain coherence between different datasets.

4) NESTED CROSS-VALIDATION:

Nested cross-validation is extensively employed for cross-validation on small datasets to prevent overfitting. We reserve 15 years of reanalysis as climatology for data pre-processing (in appendix, section a), thus limiting the size of data available for neural network training from December 1995 to March 2022. This relatively short period may lead to overfit the neural networks; therefore, we use nested cross-validation for such small datasets to optimize model selection and hyperparameter configurations while ensuring the generalization capability. Our nested cross-validation comprises outer and inner cross-validation layers, as illustrated in figure 3. For the outer layers, we divide the 27 years of reanalysis data, from December 1995 to March 2022, into three 9-year folds. For each outer split, one fold is reserved as a test set (illustrated in blue in figure 3) to verify the performance of deterministic outputs $\hat{Y}_{test}$ obtained by regressing using the Z500 reanalysis as input. The remaining two folds are used as the training set (shown in gray in figure 3), further split by the inner split for inner training and validation. Therefore, within each outer fold's 18-year training dataset, we further divide it into six 3-year inner folds for subsequent cross-validation.
<figure><img src="image_6.png" /><figcaption>FIG. 4. The spatially averaged MSSS, CRPSS and the SSIM are displayed as a function of lead week averaged across the Europe-Atlantic domain for Z500 and the Europe domain for U100.</figcaption></figure>

Here, one fold is kept as a validation set (depicted in light yellow in figure 3), and the other five serve as the training set (again shown in gray in figure 3). During the inner cross-validation, various hyperparameter combinations are trained and validated to identify the optimal configuration. The selected hyperparameter combination is then evaluated on the corresponding outer fold’s reanalysis test set and the paired hindcasts test fold. To prevent data leakage, we partition the 27 years of hindcast data $\hat{Y}_{test}^M$ based on their initialization dates to align with the temporal segments of each outer test fold (represented by the hatched blue rectangles in figure 3). We then average the skills from these three test folds to represent the final skill of our methodology.

### c. Skill of the dynamical ensembles for Z500 and U100

Toth and Buizza (2019) have noted that large-scale/low-frequency variables possess a longer forecasting horizon compared to small-scale/high-frequency variables, consistently with expectations from geophysical fluid dynamics. Here, we revisit this by quantifying how much more skillful is Z500 relative to U100 in the case considered. Since Z500 and U100 represent different physical variables, the predictability of these variables is compared using relative skill scores (Wilks 2019), such as Mean Squared Skill Score (MSSS) in figure 4a and Continuous Ranked Probability Skill Score (CRPSS) in figure 4b, for Z500 and U100 with respect to their respective 15-year rolling climatology reference. The formulas for MSSS and CRPSS are detailed in appendix, section 3. MSSS assesses the deterministic skill of the ensemble mean for Z500 and U100 hindcasts, which decreases over time. In figure 4a, U100 converges by lead week 4, while Z500 continues to
exhibit a declining trend up to lead week 6. During this period, Z500 demonstrates significantly higher MSSS than U100, indicating that, the ensemble mean of Z500 is more predictable than that of U100. CRPSS evaluates the probabilistic skill of the ensemble members for Z500 and U100 dynamical hindcasts, which also decrease over time. In figure 4b, up to lead week 6, Z500 consistently shows higher CRPSS than U100, suggesting that the distribution of the Z500 hindcasts aligns more closely with the corresponding reanalysis. These findings from the MSSS and the CRPSS reflect the inherent challenges in forecasting U100, which is less predictable than Z500. However, these scores do not account for spatial correlations for each variable.

We employ Structural Similarity Index (SSIM) as a relative metric to assess the spatial similarity between the dynamical ensemble members of Z500 and U100 and their corresponding reanalysis. SSIM, which accounts for variations in luminance, contrast, and local structure of two images, measures the spatial mean similarity, spatial variance similarity, and normalized spatial covariance, as shown in equation (A11) in appendix, section b. Predominantly utilized in computer vision (Wang et al. 2004), SSIM values range from 0 to 1, where 1 denotes perfect identity between images. SSIM quantifies the spatial covariance between the hindcasts and their respective reanalysis, facilitating the understanding of the spatial variations in the atmospheric variables across various lead times. Figure 4c presents SSIM for Z500 and U100 across lead times. The SSIM values for both Z500 and U100 exhibit a decreasing trend, with a sharp decline after lead week 2, and a convergence after lead week 4. This trend implies a progressive decrease in the dynamical forecast skill of the ECMWF hindcasts over time. Up to 6 weeks, Z500 always exhibits a higher SSIM compared to U100, indicating greater predictability for Z500.

To gain a deeper understanding of the impacts of the spatial mean, the spatial variance and the spatial covariance on the Z500 and U100 hindcasts, we delineate the three components of SSIM in appendix, section e. As evident from figure A2a in appendix, the structure component exhibits a similar declining trend to the SSIM and serves as the primary factor on the SSIM decreasing trend. Whereas the variations in the luminance component and the contrast component remain relatively minor (as in figures A2b and A2c in appendix), indicating a modest change in the dynamical spatial-averaged bias and the dynamical spatial variability, respectively. In other words, these dynamical hindcasts do not become more blurred, and the primary loss in predictability stems from changes in the spatial structure of dynamical hindcasts, particularly the misrepresentation of
features, such as local pressure lows and associated fronts. When comparing the relative skills
of Z500 to U100, it is evident that Z500 consistently exhibits higher SSIM and structure. Due
to dynamical factors, it is anticipated that Z500 is a larger-scale and more slowly evolving field
compared to wind. Furthermore, a mid-tropospheric field, characterized by its large scale, tends to
be more predictable than a surface field. That Z500 exhibits the higher MSSS, CRPSS and SSIM
values compared to U100 is encouraging to improve U100 forecasts by downscaling information
from Z500. However, improvements are only possible if the sufficient information about U100 can
be regressed from Z500.

4. Regression from reanalysis input

Whether more information from U100 can be regressed from Z500 using a non-linear model
compared to a linear model is tested here by comparing U100 targets $Y_{test}$ to U100 regressed
outputs $\hat{Y}_{test}$ from Z500 inputs $\hat{X}_{test}$ on ERA5 reanalysis. Here, we employ rolling climatology
as our benchmark. As mentioned before, this study defines rolling climatology as an interannual
mean of the same calendar dates over past 15 years. We use the deterministic metric MSE to
evaluate performance on the reanalysis test sets (figure 3), with figure 5a presenting the MSE maps
for the climatology, the MLR, and the CNN across Europe. A lower MSE indicates a more accurate
reconstruction of U100 reanalysis. As indicated at the top right corner of each subplot, the spatially
averaged MSE across Europe for the climatology, the MLR, and the CNN are 4.21, 2.39, and 2.04
(m/s)², respectively. A conspicuous feature is the land-sea contrast, with higher MSE values at
oceanic grid points and lower at terrestrial ones, due to higher wind speed and greater variability
at sea.

According to the definitions of the relative improvements in section 3d, we dis-
play the relative improvements of MSE between models (Δ<sub>r</sub>MSE(MLR, Climatology),
Δ<sub>r</sub>MSE(CNN, Climatology) and Δ<sub>r</sub>MSE(CNN, MLR)) in figure 5b. It is evident that both
the MLR and the CNN significantly outperform the climatology across all the grid points over Eu-
rope, with the spatially averaged Δ<sub>r</sub>MSE(MLR, Climatology) and Δ<sub>r</sub>MSE(CNN, Climatology)
of 40.39% and 47.00%, respectively. Especially in Western Europe, the MLR and the CNN ex-
hibit approximately 50% improvements. To explore the non-linearity in modeling the regression
relationship between Z500 and U100 reanalysis, we also compared the Δ<sub>r</sub>MSE(CNN, MLR) in
figure 5b, averaging around 11.10% spatially, with the significant improvements observed across most grid points.

The improvements of the CNN over the MLR (negative Δᵣ*MSE(CNN, MLR)*) demonstrate the presence of non-linear components in the Z500-U100 statistical relationship on the ERA5 reanalysis, which the MLR does not capture. These improvements manifest distinct spatial patterns; for instance, the degradation is seen over complex terrain like the Pyrenees, the Alps and the Scandinavian Mountains, whereas more substantial improvements occur over regions with greater wind variability, such as the North Sea and the Bay of Biscay.

Overall, the CNN is more effective at extracting the information from Z500 reanalysis, thereby better representing the spatial variability of U100 over Europe. Together with the results of section 3c, these improvements motivate to use the non-linearity between Z500 and U100 reanalysis to regress U100 ensembles from Z500 ensembles.

**5. Improving ensemble forecasts using non-linear regression**

In the previous section, we demonstrated that non-linearity can improve the reconstruction of U100 from Z500 reanalysis. In this section, we explore whether this non-linearity could be exploited to improve the sub-seasonal skill of U100 hindcasts compared to the linear model in figure 1b. First, we apply the deterministic MLR and CNN as trained in the previous section to the 10-member dynamical Z500 hindcasts $X_{test}^M$, where $M = 10$, from ECMWF to regress statistical ensembles following the workflow in figure 1b. We evaluate the forecast skill of these regressed ensembles $\hat{Y}_{test}^M$ from the MLR and the CNN on the sub-seasonal time scale, using the ECMWF hindcasts $Y_{test}^M$ as a benchmark. Subsequently, we assess the effectiveness of the stochastic perturbations $\tilde{Y}_{test}^{M \times P}$, where $P = 20$, following the same approach, when we need to attribute in the skill improvements to a better representation of the predictable part of the signal or to a better representation of the uncertainty.

a. Forecasts from deterministic regression

As outlined in section 3, the ensemble uncertainty of Z500 hindcasts increases with lead time due to the chaotic nature of the atmosphere. In this section, we explore whether the non-linear CNN, compared to the MLR, can further improve U100 skill in the presence of hindcast uncertainty. In
<figure><img src="image_7.png" /><figcaption>FIG. 5. Maps of MSE (top) and Δ<sub>r</sub>MSE (bottom) for the climatology, the MLR and the CNN estimated from validation reanalysis data. Δ<sub>r</sub>MSE(model, bench) denotes the difference between the MSE of model and the MSE of benchmark over the MSE of the latter. Thus, negative (positive) values correspond to improved (worsened) skill. To assess the robustness of Δ<sub>r</sub>MSE, gray scatter points marked with "+" on the grid points denote statistically significant (0.01 significant level) improvements or degradations at those grid points. The spatial mean μ of the score and the Δ score are provided at the top of each subfigure.</figcaption></figure>

this subsection, ensembles from the ECMWF, the MLR, and the CNN each consist of 10 members, allowing us to conduct verifications with equal member size to avoid unfair comparison due to differing member size (Zamo and Naveau 2018).

Figure 6 presents the spatial mean of MSE, CRPS and SSR of the ECMWF, the MLR, the CNN and the perturbed version of the MLR (MLR) and the CNN (CNN), as functions of lead weeks. We first focus on the skill of the ECMWF, the MLR and the CNN. Only the medians derived from the bootstrap method (section 3d) are provided to enhance the readability of these figures, while significant improvements in MSE and CRPS across models for different lead weeks are displayed in table 1. The MSE and the CRPS of these models increase with lead time and stabilize from the
<figure><img src="image_8.png" /><figcaption>FIG. 6. The spatial mean of MSE, CRPS and SSR as a function of lead week across the Europe domain for U100 hindcasts from the ECMWF, the MLR, the CNN, the *MLR* and the *CNN*.</figcaption></figure>

third week most likely due to increasing uncertainty. After lead week 4, there is a slight decrease in both MSE and CRPS for these models, likely because the hindcasts initialized in February include end dates in spring for lead week 5 and 6, which tend to be more predictable due to the lower U100 variability compared to winter dates (Cortesi et al. 2019). At lead weeks 1 and 2, the ECMWF exhibits skillful MSE and CRPS. At these lead times a significant fraction of the information from the initial conditions remains and the evolution of the state of the atmosphere is relatively well captured by the ECMWF IFS, making it difficult for the MLR and the CNN models, which only use Z500 hindcasts as inputs, to surpass. The CNN outperforms the MLR in terms of MSE up to lead week 2, demonstrating that non-linearity can improve the deterministic skill in the first lead weeks. However, this advantage diminishes as we progress into longer lead times. Starting from lead week 3, the MSE curves for the MLR and the CNN converge towards the ECMWF benchmark. This convergence reflects a fundamental aspect of sub-seasonal forecasting, where the input-output relationships become increasingly dominated by unpredictable components at longer lead times. In such conditions, the complex non-linear mappings learned by CNN become less effective, and the simpler linear approximation provided by MLR appears to be equally capable of capturing the remaining predictable signals. Table 1 further reveals that, at lead week 3, the $\Delta_r \text{MSE}(\text{MLR}, \text{ECMWF})$ and $\Delta_r \text{MSE}(\text{CNN}, \text{ECMWF})$ are marginal, remaining within a 1%, but statistically significant. This observation suggests that increased model complexity alone cannot overcome the fundamental predictability barriers inherent in sub-seasonal forecasting.
<table><caption>TABLE 1. Spatially averaged $\Delta_r$ MSE (top) and $\Delta_r$ CRPS (bottom) comparing the improvements in the MLR over the ECMWF, the CNN over the ECMWF, the CNN over the MLR, and their perturbed version over the Europe domain for weeks 3 to 6. Negative values denote improved skill, while positive values signify a deterioration. Superscripts $^a$, $^b$, and $^c$ indicate statistical significance: $^a$ for $0 < p < 0.01$, $^b$ for $0.01 < p < 0.05$, $^c$ for $0.05 < p < 0.1$, where $p$ is the p-value estimated from the bootstrap method mentioned in section 3d.</caption><thead><tr><th>Δ<sub>r</sub> MSE (%)</th><th>Week 3</th><th>Week 4</th><th>Week 5</th><th>Week 6</th></tr></thead><tbody><tr><td>(MLR,ECMWF)</td><td>-0.84 <sup>a</sup></td><td>-0.42 <sup>a</sup></td><td>-0.6 <sup>a</sup></td><td>+0.7 <sup>a</sup></td></tr><tr><td>(CNN,ECMWF)</td><td>-0.05 <sup>a</sup></td><td>+0.71</td><td>+1.01 <sup>b</sup></td><td>+1.65 <sup>a</sup></td></tr><tr><td>(CNN,MLR)</td><td>+0.85 <sup>b</sup></td><td>+1.17 <sup>a</sup></td><td>+1.63 <sup>a</sup></td><td>+0.92 <sup>b</sup></td></tr><tr><td>(<span style="text-decoration: overline;">MLR</span>,ECMWF)</td><td>-0.56 <sup>a</sup></td><td>-0.14 <sup>a</sup></td><td>-0.33 <sup>a</sup></td><td>+0.97 <sup>a</sup></td></tr><tr><td>(<span style="text-decoration: overline;">CNN</span>,ECMWF)</td><td>-0.55 <sup>a</sup></td><td>+0.2 <sup>a</sup></td><td>+0.62 <sup>a</sup></td><td>+1.52 <sup>a</sup></td></tr><tr><td>(<span style="text-decoration: overline;">CNN</span>,<span style="text-decoration: overline;">MLR</span>)</td><td>+0.05</td><td>+0.39 <sup>b</sup></td><td>+1.0 <sup>a</sup></td><td>+0.59 <sup>c</sup></td></tr></tbody><thead><tr><th>Δ<sub>r</sub> CRPS (%)</th><th>Week 3</th><th>Week 4</th><th>Week 5</th><th>Week 6</th></tr></thead><tbody><tr><td>(MLR,ECMWF)</td><td>+5.11 <sup>a</sup></td><td>+4.81 <sup>a</sup></td><td>+4.63 <sup>a</sup></td><td>+4.97 <sup>a</sup></td></tr><tr><td>(CNN,ECMWF)</td><td>+2.52 <sup>a</sup></td><td>+2.77 <sup>a</sup></td><td>+2.89 <sup>a</sup></td><td>+3.07 <sup>a</sup></td></tr><tr><td>(CNN,MLR)</td><td>-2.43 <sup>a</sup></td><td>-1.92 <sup>a</sup></td><td>-1.64 <sup>a</sup></td><td>-1.8 <sup>a</sup></td></tr><tr><td>(<span style="text-decoration: overline;">MLR</span>,ECMWF)</td><td>-2.87 <sup>a</sup></td><td>-2.72 <sup>a</sup></td><td>-2.75 <sup>a</sup></td><td>-2.29 <sup>a</sup></td></tr><tr><td>(<span style="text-decoration: overline;">CNN</span>,ECMWF)</td><td>-2.57 <sup>a</sup></td><td>-2.19 <sup>a</sup></td><td>-1.94 <sup>a</sup></td><td>-1.68 <sup>a</sup></td></tr><tr><td>(<span style="text-decoration: overline;">CNN</span>,<span style="text-decoration: overline;">MLR</span>)</td><td>+0.32 <sup>b</sup></td><td>+0.56 <sup>a</sup></td><td>+0.85 <sup>a</sup></td><td>+0.64 <sup>a</sup></td></tr></tbody></table>

When analyzing the CRPS in figure 6b, the CNN exhibits a consistently lower CRPS throughout the sub-seasonal time scale compared to the MLR. However, both the MLR and the CNN still underperform relative to the ECMWF. Table 1 demonstrates a significant improvement in CRPS for the CNN over the MLR by 2.43%, suggesting that the non-linearity enhances the CRPS performance, yet the CNN remains 2.52% worse than that of the ECMWF.

To understand why the CRPS of these statistical models is inferior to that of the ECMWF and to determine the contribution of non-linearity in improving CRPS, we employ SSR to assess ensemble reliability. Recall from section 3d that an SSR value of 1 indicates a reliable ensemble, while values less than 1 suggest that the ensemble is under-dispersive, while values greater than 1 indicate an over-dispersive ensemble. In figure 6c, across the time scales of interest, the ensembles of the ECMWF, the MLR, and the CNN (shown in solid lines) are consistently under-dispersive and their dispersion increases but insufficiently after the first two weeks. The ECMWF approaches an SSR of 0.92, while the MLR and the CNN converge to 0.55 and 0.68, respectively, indicating these models lack spread, and so underestimate the forecast uncertainty. The low SSR of the MLR
means its regressed members are close to each other, reflecting its ensemble mean. Notably, in the first two weeks, the dispersion of the MLR and the CNN rapidly increases, implying a swift growth in the uncertainty of the Z500 ensemble members used as inputs, with the CNN displaying greater dispersion than the MLR. Thus, even though non-linearity does not improve the MSE of ensemble mean, it contributes to improving the probabilistic metrics (SSR and CRPS) by better representing the regression model uncertainty.

Aggregating these skills spatially obscures the spatial patterns, which are valuable for wind energy applications. In order to analyze the dependence of the skills on different grids, figure 7 displays $\Delta_r$MSE and $\Delta_r$CRPS for lead week 3. Given the consistency in spatial patterns of the $\Delta_r$MSE and $\Delta_r$CRPS across various lead weeks, only the improvements for lead week 3 are presented. The maps detailing the $\Delta_r$MSE and $\Delta_r$CRPS from weeks 4 to 6 can be found in appendix (figures A3 and A4). For lead week 3, while the spatially averaged $\Delta_r$MSE in table 1 shows marginal improvement for the MLR and the CNN relative to the ECMWF, distinct spatial patterns are observable in figure 7. The MLR and the CNN exhibit the more positive $\Delta_r$MSE in some grids of Northern Europe than the ECMWF, indicating poorer performance in these regions. Conversely, in Western and Southeastern Europe, both the $\Delta_r$MSE(MLR, ECMWF) and $\Delta_r$MSE(CNN, ECMWF) are negative, with about a 5% improvement in MSE noted in these grids relative to the ECMWF. Similar to the findings in section 4, when applying the CNN to hindcasts, the CNN still struggles to effectively reconstruct wind speed in the grids with complex topography, resulting in poorer MSE and CRPS compared to the ECMWF and the MLR in these grids. The improvements of the CNN over the MLR also display a spatial distribution; the CNN performs better than the MLR in the North Sea, likely due to the greater wind speed variability in this region during winter and the CNN's non-linearity that more effectively reconstruct wind speed here (see section 4). However, the CNN's MSE is approximately 4% worse than the MLR's in Northern Europe.

For the maps of $\Delta_r$CRPS in figure 7b, a spatial pattern similar to $\Delta_r$MSE is observed, where regions with negative $\Delta_r$MSE correspond to lower $\Delta_r$CRPS, even though most grid points across Europe underperform relative to the ECMWF. When considering $\Delta_r$CRPS(CNN, MLR), the marginal but significant improvements (2.43%) are noted in most parts of Europe, except for Northern Europe and the grids with complex topography. This improvement is attributed to the
<figure><img src="image_9.png" /><figcaption>FIG. 7. Δ<sub>r</sub>MSE (top) and Δ<sub>r</sub>CRPS (bottom) at lead week 3 comparing the improvements in the MLR over the ECMWF, the CNN over the ECMWF, and the CNN over the MLR. The significance of the results is represented in the same way as in figure 5.</figcaption></figure>

CNN's ability to better fit the signal than the MLR, resulting in a higher variance of the regressed signal during model training. Therefore, when we use ensembles as inputs, the CNN can express more variance, thereby increasing ensemble dispersion.

Overall, the non-linearity indeed contributes to the improvement of MSE and CRPS relative to the linear model during the first two lead weeks. However, the advantage of non-linearity in the deterministic skill measured by MSE diminishes after lead week 3. Despite this, in certain regions, such as Northwestern Europe, the non-linear model still outperforms the linear model. The CNN exhibits superior performance over the linear model in terms of CRPS by better representing the dispersion of the ensembles within the sub-seasonal timescale. Nevertheless, both the MLR and the CNN exhibit under-dispersion.
<figure><img src="image_10.png" /><figcaption>(a) $\Delta_r$ MSE map at lead week 3 for the perturbed models</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) $\Delta_r$ CRPS map at lead week 3 for the perturbed models</figcaption></figure>

FIG. 8. Same as figure 7. But comparing the improvements in the $\overline{MLR}$ over the ECMWF, the $\overline{CNN}$ over the ECMWF, and the $\overline{CNN}$ over the $\overline{MLR}$.

b. Forecasts from stochastic perturbations

To address the issue of under-dispersion in statistical ensembles, we employ the method discussed in section 3c to perturb the statistical ensembles of the MLR and the CNN. This subsection compares the perturbed versions of the MLR and the CNN, denoted as $\overline{MLR}$ and $\overline{CNN}$, respectively, which integrate model uncertainty during training, examining both the temporal evolution and spatial improvements of ensemble skills.

Dashed lines in figure 6 display the temporal evaluation of forecast skill of the perturbed ensembles $\tilde{Y}^{M \times P}$ from the $\overline{MLR}$ and $\overline{CNN}$. In addition to the results for the non-perturbed models presented in section 5a, table 1 also details the improvements at lead week 3 for these perturbed models. As expected, the perturbed ensembles from the $\overline{MLR}$ and the $\overline{CNN}$ do not exhibit great improvements in MSE compared to the MLR and the CNN, where improvements by these perturbed models over their direct models are less than 1% and not statistically significant (not shown
here), which is not surprising considering that MSE is not a probabilistic score and does not derive
benefits from the stochastic perturbations. However, the CRPS shows significant improvement
from lead week 1 to 6. By lead week 3, the $\overline{MLR}$ and the $\overline{CNN}$ show over a 4% improvement over
their non-perturbed versions, represented by $\Delta_r \text{CRPS}(\overline{MLR}, \text{MLR})$ and $\Delta_r \text{CRPS}(\overline{CNN}, \text{CNN})$,
which are not provided in table 1; however, they can be calculated using the table information such
as $\Delta_r \text{CRPS}(\overline{MLR}, \text{MLR}) = \Delta_r \text{CRPS}(\overline{MLR}, \text{ECMWF}) - \Delta_r \text{CRPS}(\text{MLR}, \text{ECMWF})$ (same for
the $\overline{CNN}$). In the first two weeks, the $\overline{CNN}$ reports a lower CRPS than the $\overline{MLR}$, though
the gap between the $\overline{MLR}$ and the $\overline{CNN}$ narrows after week 2, converging towards parity (the
$\Delta_r \text{CRPS}(\overline{CNN}, \overline{MLR})$ is less than 1% in table 1). From week 3 onwards, the CRPS of the $\overline{MLR}$
and the $\overline{CNN}$ surpass that of the ECMWF, with an improvement of over 2% compared to the
ECMWF (in table 1). The SSR (figure 6c) for the $\overline{MLR}$ and the $\overline{CNN}$ approaches 0.84,
indicating that much of the under-dispersion has been corrected. Notably, the $\overline{MLR}$ and the $\overline{CNN}$,
which incorporate the stochastic perturbations, exhibits greater reliability than the MLR and the
CNN. This suggests that representing uncertainty from the training model is key to enhance the
probabilistic skill of statistical ensembles.

The grid-level comparisons of $\Delta_r \text{MSE}$ and $\Delta_r \text{CRPS}$ for the $\overline{MLR}$ and the $\overline{CNN}$ are illustrated
in figure 8, with a specific emphasis on week 3. Similarly, the maps for $\Delta_r \text{MSE}$ and $\Delta_r \text{CRPS}$ from
week 4 to week 6 can be found in figures A3b and A4b. The improvements in MSE for perturbed
ensembles are temporally consistent, corroborating the findings presented in figure 8a. Across
the European domain, except for the $\overline{MLR}$ in the North Sea, there are improvements in CRPS for
both models. The introduction of the perturbations reduces the $\Delta_r \text{CRPS}(\overline{CNN}, \overline{MLR})$ of 2.43%
in figure 7b to the $\Delta_r \text{CRPS}(\overline{CNN}, \overline{MLR})$ of 0.32% in figure 8b, rendering them statistically
insignificant over most grids. However, the slight improvement in CRPS for the $\overline{CNN}$ over the
$\overline{MLR}$, although less than 1%, is significant and predominantly observed in the North Sea and
Southern Europe.

Overall, the introduction of uncertainty quantification in these regression models has significantly
improved the skills of the MLR and the CNN. Furthermore, for the perturbed ensembles, the
spatially averaged improvements from the non-linearity are marginal.
# 6. Conclusion

This study investigates the potential of non-linearity in improving the skill of U100 ensemble forecasts on the sub-seasonal timescale during European winter, while also exploring methods to better represent ensemble spread when applying regression models to sub-seasonal dynamical forecasts. Building on the works of Alonzo et al. (2017) and Goutham et al. (2023), we employed both linear (MLR) and non-linear (CNN) models to compare their performance in modeling the regression relationship from Z500 reanalysis and hindcasts to surface variable U100.

Our results demonstrate that the CNN exhibits superior capability in capturing non-linear features compared to the MLR on the reanalysis dataset, thus enabling more accurate reconstruction of U100 over Europe, from knowledge of Z500 over the Europe-Atlantic domain, as indicated by the spatially averaged MSE. Specifically, the MSE for the MLR is $2.39 (m/s)^2$, while the CNN achieves an MSE of $2.04 (m/s)^2$. This represents approximately 11.10% improvement in performance attributed to the non-linear characteristics on the reanalysis.

When applied to the ECMWF weekly Z500 hindcasts, the improvements of forecast skill attributed to the non-linearity diminishes over time in terms of spatial averages. However, within specific regions (such as North Sea), the non-linearity still contributes to both MSE of ensemble mean and CRPS.

We also observed that the U100 ensembles generated by both the MLR and the CNN exhibited under-dispersion, representing an underestimation of forecast uncertainty, leading to poorer skills compared to the ECMWF dynamical ensembles. In other data driven weather forecast models(Robertson et al. 2015; Bi et al. 2022; Kurth et al. 2023; Chen et al. 2023; Orth and Seneviratne 2014), it has also been found that the forecasts focused on the predictable part of the signal, and hence lacked dispersion. To address this issue, we introduced stochastic perturbations, which substantially improved the reliability and spread of the ensembles. As quantified by SSR, the MLR and the CNN exhibit the SSR values of 0.55 and 0.68, respectively. When the stochastic perturbations are introduced, the SSR increases to 0.84 for both models.

Although the perturbed models ($\overline{MLR}$ and $\overline{CNN}$) still demonstrated some under-dispersion, their reliability improved relative to the original MLR and CNN. Notably, the dispersed ensemble had significantly enhanced probabilistic skill, with a corresponding improvement of approximately 4% in CRPS for both the $\overline{MLR}$ and the $\overline{CNN}$ over their non-perturbed versions. Thus, better
quantification of uncertainty on the sub-seasonal timescale is the key to improve sub-seasonal
forecast skill.

This study represents the first application of non-linear reconstruction of U100 from Z500 to
improve sub-seasonal skill in U100 forecasts, as well as the introduction of the unexplained part
of the signal by the regression model to improve ensemble reliability. Using historical data and
only forecasting information of the large-scale part of the flow, as described by Z500, our approach
combining statistical models with stochastic perturbations demonstrates improved probabilistic
forecast skill over dynamical forecasts after lead week 3, particularly in CRPS, while maintaining
comparable MSE. While our research focuses on Z500-U100 relationships, this methodology has
potential applications to other variables with regression relationships (Orth and Seneviratne 2014;
Goutham et al. 2023).

In terms of methodological choices, we utilized established MLR and CNN models. While
more complex CNN architectures may yield further improvements, we found that designing a more
realistic representation of regression uncertainty is crucial for significantly enhancing ensemble
skill. Our approach did not account for spatial correlation among the stochastic perturbations,
which presents an avenue for future research. Additionally, our calibration method (Mean Variance
Adjustment, in appendix, section c) faces consistency challenges across temporal and spatial
dimensions, as it calibrates dynamical forecasts independently at each lead time and grid point,
potentially disrupting the inherent spatial and temporal dependencies.

Future research directions include exploring whether the correlations in the regressed spatial
patterns can be improved through a more sophisticated representation of stochastic perturbations.
Additionally, it is worth investigating whether representing forecast uncertainty with ensemble
spread could be improved by inputting the full ensemble distribution into the machine learning
model, rather than applying the model separately to each ensemble member. Another direction for
future research could examine whether learning directly from hindcasts/forecasts, rather than from
reanalysis, and including the dynamical forecasts of target variables as inputs can help improve
ensemble forecast skill. Furthermore, explainable AI methods could help elucidate which input
variables and which hidden units are actually responsible for skill improvements on average and
geographically.
*Acknowledgments.* This work has been carried out at the Energy4Climate Interdisciplinary Center (E4C) of IP Paris and Ecole des Ponts ParisTech, which is in part supported by 3rd Programme d'Investissements d'Avenir [ANR-18-EUR-0006-02], and by the Foundation of Ecole polytechnique (Chaire "Défis Technologiques pour une Énergie Responsable" financed by TotalEnergies).

*Data availability statement.* The ERA5 reanalysis data used in this study are available from the Climate Data Store (https://cds.climate.copernicus.eu). The ECMWF extended-range forecast data were obtained through Meteorological Archival and Retrieval System with institutional license. Researchers interested in sub-seasonal forecast and hindcast can access related datasets through the ECMWF S2S project (https://www.ecmwf.int/en/research/projects/s2s). All source codes have been made available on our GitHub repository (https://github.com/TIANGANGLIN/s2s-wind-Non-linearity). The repository also contains the trained parameters for both the MLR and CNN models, as well as their perturbed versions.

# APPENDIX

## a

Table A1 defines the symbols used in this article. For instance, $y_{g,n}$ represents the deterministic value of the input at the *t*-th of *T* initialization times at a specific grid point, indicated by $g \in G$. Similarly, $y_{g,n,l}^m$ symbolizes the *m*-th member of a *M*-member ensemble at the *t*-th initialization time (of *T* total times) for a given lead time *l* at the corresponding grid point. For ensembles, the notation $\overline{\hat{y}}^M = \frac{1}{M} \sum_m^M y^m$ is used to denote the ensemble mean. The ensembles and the models that are perturbed are denoted with a tilde "~".

### a. Pre-processing

Our data pre-processing module involves sequential steps of spatial and temporal averaging and normalization. Reanalysis and ensembles are initially downsampled to a spatial resolution of 2.7 degrees (approximately 300 km) and a temporal resolution of 7 days. Before normalization, we perform deseasonalization and detrending to isolate the underlying sub-seasonal signals. First, we remove the seasonal cycle by subtracting the mean climatology for each initialization date *t*. Unlike traditional grid wise deseasonalization, we calculate the spatial mean of climatology for each *t* and subtract this seasonal pattern from all grid points to preserve spatial coherence in the
<table><caption>TABLE A1. Symbols in this study</caption><thead><tr><th>Symbol</th><th>Range</th><th>Description</th></tr></thead><tbody><tr><td>x</td><td></td><td>Input</td></tr><tr><td>y</td><td></td><td>Target</td></tr><tr><td>ŷ</td><td></td><td>Output</td></tr><tr><td>o</td><td></td><td>Ground Truth</td></tr><tr><td>~</td><td></td><td>Perturbed version</td></tr><tr><td>t</td><td>1, ..., T</td><td>Initialization time index</td></tr><tr><td>l</td><td>1, ..., L</td><td>Lead time index</td></tr><tr><td>g</td><td>1, ..., G</td><td>Grid point index</td></tr><tr><td>m</td><td>1, ..., M</td><td>Ensemble-member index</td></tr></tbody></table>

anomalies. Second, we remove the long-term trend in the anomalies by fitting a linear model to the anomalies time series of spatially-averaged values, and then remove this estimated trend from all grid points, ensuring consistency in spatial patterns. Subsequently, the normalization of an input $x$ at the grid point $g$ is performed as $\text{Normalization}(x_{t,g}) = (x_{t,g} - \mu_g(X_{train})) / \sigma_g(X_{train})$, where $\mu_g(X_{train})$ and $\sigma_g(X_{train})$ respectively denote the mean and standard deviation along the initialized time dimension of the training dataset. These pre-processing operations are also applied independently to each forecast member $m$ and each lead week $l$ for forecasts and hindcasts.

## b. Scores

We employ the symbols defined in Table A1 to calculate various verification scores, using the target variable $y$ as an example. The Mean Squared Error (MSE), the Continuous Ranked Probability Score (CRPS), the Spread Skill Ratio (SSR), the Mean Squared Skill Score (MSSS) and the Continuous Ranked Probability Skill Score (CRPSS) are initially calculated at the grid level to generate spatial distribution maps, then averaged spatially with cosine-latitude weights to illustrate the evolution of forecast skill over lead weeks.

### 1) MEAN SQUARED ERROR (MSE)

MSE is evaluated independently at each grid point for both the deterministic and probabilistic datasets. Specifically, for the regressed deterministic 100uv, represented as $ŷ$, MSE is computed as $\text{MSE}(y_{t,g}, \hat{y}_{t,g})$. For the regressed probabilistic 100uv ensembles $\hat{y}^M$, MSE is calculated using the ensemble mean $\hat{y}^M$ at each grid point $g$ and for each lead week $l$, denoted as $\text{MSE}(y_{t,l,g}, \overline{y}^M_{t,l,g})$.
$$MSE(y_{t,g}, \hat{y}_{t,g}) = \frac{1}{T} \sum_{t}^{T} (y_{t,g} - \hat{y}_{t,g})^2 \quad (A1)$$

$$MSE(y_{t,l,g}, \overline{\hat{y}}_{t,l,g}^{M}) = \frac{1}{T} \sum_{t}^{T} (y_{t,l,g} - \overline{\hat{y}}_{t,l,g}^{M})^2 \quad (A2)$$

$$\overline{\hat{y}}^{M} = \frac{1}{M} \sum_{m}^{M} \hat{y}^{m} \quad (A3)$$

## 2) CONTINUOUS RANKED PROBABILITY SCORE (CRPS)

CRPS is used to assess the probabilistic skill of ensembles. Our benchmark dynamical hindcasts ensembles from ECMWF comprise merely 10 members, whereas our perturbed ensembles consist of 200 members. Zamo and Naveau (2018) and Goutham et al. (2022) indicate that an increase in ensemble member size improves the CRPS, making ensembles more skillful. To compare ensembles with varying member size, Zamo and Naveau (2018) suggests using a fair version of CRPS (Ferro 2014) that estimates the CRPS as the member size approaches infinity. However, unlike the work of Zamo and Naveau (2018), our perturbed ensembles do not satisfy the Fair CRPS assumption of exchangeability. This is because the perturbed ensembles are sampled based on each regressed member, which does not meet the requirements of being independently and identically distributed. Therefore, to ensure a fair comparison between the 10-member benchmark and the 200-member perturbed ensembles, following Zamo and Naveau (2018)'s approach, we downsample the 200-member perturbed ensembles to 10 members. To preserve the Cumulative Distribution Function (CDF) of the 200 members, we do not randomly select 10 members; instead, we use 10 quantiles of the 200 members as 10 ensemble members. Consequently, our comparison involves the 10-member ECMWF hindcasts and the 10 quantiles of the perturbed ensembles. Zamo and Naveau (2018) refers to the ensembles produced through this method as “optimal ensembles,” which achieve the minimal CRPS.
We employ a discrete version of CRPS (Zamo and Naveau 2018), calculating $CRPS_{l,g}$ for each lead time and each grid point.

$$CRPS_{l,g} = \frac{1}{t} \sum_{t}^{T} \left( \frac{1}{M} \sum_{m}^{M} |\hat{y}_{t,l,g}^{m} - y_{t,l,g}| - \frac{1}{2M^2} \sum_{m}^{M} \sum_{n}^{M} |\hat{y}_{t,l,g}^{m} - \hat{y}_{t,l,g}^{n}| \right) \quad (A4)$$

### 3) MEAN SQUARED SKILL SCORE (MSSS) AND CONTINUOUS RANKED PROBABILITY SKILL SCORE (CRPSS)

MSSS and CRPSS are used as relative scores. The subscript *ens* denotes the ensembles of interest, while *ref* refers to the reference with which comparisons are made. In this study, the 15-year rolling climatology is treated as reference.

$$MSSS = 1 - \frac{MSE_{ens}}{MSE_{ref}} \quad (A5)$$

$$CRPSS = 1 - \frac{CRPS_{ens}}{CRPS_{ref}} \quad (A6)$$

### 4) SPREAD SKILL RATIO (SSR)

SSR is the ratio of ensemble spread to the Root MSE of ensemble mean (Rasp et al. 2023), with ensemble spread represented by $Spread_{l,g} = \sqrt{\frac{1}{T} \sum_{t}^{T} var_{m} (y_{t,l,g}^{m})^2}$, where $var_{m}$ is the variance across ensemble members.

$$SSR_{l,g} = \frac{Spread_{l,g}}{\sqrt{MSE_{l,g}}} \quad (A7)$$
## 5) STRUCTURAL SIMILARITY INDEX MEASURE (SSIM)

SSIM comprises luminance L, contrast C and structure S (Wang et al. 2004), with $\mu_g$, $\sigma_g$, and $cov_g$ denoting the spatial mean, the standard deviation, and the covariance, respectively.

$$SSIM_l = L_l * C_l * S_l \qquad (A8)$$

$$L_l = \frac{1}{T * M} \sum_{t}^{T} \sum_{m}^{M} \frac{2\mu_g(\hat{y}_{t,m,l})\mu_g(y_{t,m,l})}{\mu_g^2(\hat{y}_{t,m,l}) + \mu_g^2(y_{t,m,l})} \qquad (A9)$$

$$C_l = \frac{1}{T * M} \sum_{t}^{T} \sum_{m}^{M} \frac{2\sigma_g(\hat{y}_{t,m,l})\sigma_g(y_{t,m,l})}{\sigma_g^2(\hat{y}_{t,m,l}) + \sigma_g^2(y_{t,m,l})} \qquad (A10)$$

$$S_l = \frac{1}{T * M} \sum_{t}^{T} \sum_{m}^{M} \frac{cov_g(\hat{y}_{t,m,l}, y_{t,m,l})}{\sigma_g(\hat{y}_{t,m,l})\sigma_g(y_{t,m,l})} \qquad (A11)$$

### c. Calibration

we apply a simple and effective bias adjustment method, known as the Mean-Variance Adjustment (MVA). MVA is commonly utilized for the calibration of seasonal and sub-seasonal ensembles (Goutham et al. 2022, 2023; Manzanas et al. 2019). The mean and variance of ensembles are adjusted towards the mean and variance of reference data. We calibrate the forecasts and hindcasts using the same MVA as in Goutham et al. (2022)'s work.

### d. Models

#### 1) MULTIPLE LINEAR REGRESSION (MLR)

This study employs a multiple linear regression (MLR) model to capture the linear relationship between Z500 and U100 in a spatially independent way. For each grid point g, a separate regression model is applied:

$$\hat{Y}_g = \beta_{0,g} + \sum_{i=1}^{G} \beta_{i,g} \cdot X_i + \epsilon_g$$

where G denotes the total number of Z500 input grid points across the Europe-Atlantic domain, $\hat{Y}_g$ is the regressed U100 at grid point g, $\beta_{0,g}$ is the intercept term specific to grid point g, $\beta_{i,g}$ represents the regression coefficient linking input point i to output point g, $X_i$ is the Z500 value at input grid point i, and $\epsilon_g$ is the error term.
<figure><img src="image_12.png" /><figcaption>FIG. A1. SmaAt-UNet architecture for regressing U100 from Z500. The diagram shows the network structure with Convolutional Block Attention Modules (CBAM), depthwise-separable convolutions (DSC), batch normalization (BN), and other key components that enable the model to capture the non-linear relationships between geopotential height and wind speed fields.</figcaption></figure>

## 2) CONVOLUTIONAL NEURAL NETWORK (CNN) - SMAAt-UNET

The SmaAt-UNet in figure A1 is employed to capture the non-linear relationships between Z500 and U100. Our experiments showed that increasing model complexity through additional layers or more advanced architectures did not improve MSE performance on hindcasts, suggesting we had reached the complexity limit supported by our available data. The input is the Z500 across the European-Atlantic domain (22×59 grid), while the output is the regressed U100 over the Europe domain (14×19 grid). The encoder of SmaAt-UNet includes four downsampling stages, each comprising: (1) two consecutive 3×3 depthwise-separable convolutions (DSC), followed by batch normalization (BN) and ReLU activation; (2) a 2×2 max-pooling layer to reduce spatial dimensions; (3) Convolutional Block Attention Modules (CBAM) are placed after the first double convolution and after each encoder to amplify important features. As the network deepens, the number of feature channels increases (from 16 to 128), while the spatial resolution decreases. The decoder also comprises four stages, reconstructing U100 through bilinear interpolation. Key characteristics at each layer include: (1) skip connections from corresponding encoder layers to preserve detailed information; (2) two consecutive depthwise-separable convolution operations identical to those in
<figure><img src="image_13.png" /><figcaption>FIG. A2. Same as figure 4, but for Structure, Luminance and Contrast.</figcaption></figure>

the encoder; (3) bilinear upsampling operations to double the spatial resolution of feature maps. Ultimately, a 1×1 convolution operation is applied to map the feature space and a crop operation is applied to truncate it to the target resolution of a 14×19 grid for U100.

Both MLR and CNN are implemented on Pytorch (Paszke et al. 2019) and trained using Adam optimizer to minimize MSE, with hyperparameters (learning rate and weight decay) optimized via Optuna (Akiba et al. 2019).

### e. Forecast skills of z500 and 100uv

The three components of SSIM, calculated using the formula presented in appendix, section b, are displayed in figure A2. Luminance and Contrast exhibit small variation across the lead weeks, indicating that the spatial mean and the variance of ensemble members remain relatively stable over time. As illustrated in figure 4c, the decrease in SSIM across the lead weeks primarily stems from a decrease in Structure. This suggests significant variations in the covariance between ensemble members and the corresponding reanalysis, inaccurately representing the spatial patterns of the atmospheric states.

### f. Skill improvements of the 100uv ensembles from the ECMWF, the MLR, the CNN, the MLR, the CNN from week 4 to week 6

In the main text, we presented the $\Delta_r$MSE and $\Delta_r$CRPS for the models during lead week 3. Here, we extend our analysis from lead week 4 to lead week 6, demonstrating $\Delta_r$MSE and $\Delta_r$CRPS for this period in figures A3 and A4. Notably, during lead weeks 4 and 5, from
the $\Delta_r \text{MSE}(\text{MLR}, \text{ECMWF})$ and the $\Delta_r \text{MSE}(\text{CNN}, \text{ECMWF})$ both the MLR and the CNN exhibit significant improvements over the ECMWF primarily in France, Spain and Eastern Europe; however, by lead week 6, a notable degradation is observed in the North Atlantic. The comparison between the $\Delta_r \text{MSE}(\text{MLR}, \text{ECMWF})$ and the $\Delta_r \text{MSE}(\text{CNN}, \text{ECMWF})$ highlights these trends. The improvements attributed to non-linearity are primarily concentrated in the North Sea, Germany, and Eastern Europe. However, the CNN exhibits worse MSE compared to the MLR in regions with complex terrain and Southern Europe. Similarly, the perturbed versions, $\Delta_r \text{MSE}(\widetilde{\text{CNN}}, \widetilde{\text{MLR}})$ in figure A3b demonstrates similar spatial distribution patterns as $\Delta_r \text{MSE}(\text{CNN}, \text{MLR})$.

Figure A4 illustrate the $\Delta_r \text{CRPS}$ from lead week 4 to lead week 6, where the MLR and the CNN demonstrate worse CRPS than the ECMWF. However, the benefits of non-linearity, indicated by more negative values in CRPS as analyzed in the main text, are evident. For the $\widetilde{\text{MLR}}$ and the $\widetilde{\text{CNN}}$, the introduction of the stochastic perturbations consistently aids CRPS across the sub-seasonal timescales, although the non-linear assistance diminishes with the increasing lead weeks.

g. Spatially averaged skills of the 100uv ensembles from the ECMWF, the MLR, the CNN, the $\widetilde{\text{MLR}}$, the $\widetilde{\text{CNN}}$ on forecasts

Similar to figure 6, we present the spatially averaged MSE, CRPS, and SSR for the ECMWF, the MLR, the CNN, the $\widetilde{\text{MLR}}$ and the $\widetilde{\text{CNN}}$ across the forecasts dataset in figure A5. It is noteworthy that the size of the forecasts dataset is considerably smaller than that of the hindcasts, leading to less representative results from the forecasts. Nonetheless, the conclusions drawn from the forecasts are consistent with those from the hindcasts dataset.

Compared to figure 6, the performance of the MLR, the CNN, the $\widetilde{\text{MLR}}$ and the $\widetilde{\text{CNN}}$ on forecasts aligns with their performance on hindcasts. For MSE, within the initial two weeks, the ECMWF benefits from its complex dynamics, exhibiting better MSE of the ensemble mean. From lead week 2, the MLR and the CNN demonstrate better MSE than the ECMWF, and the $\widetilde{\text{MLR}}$ and the $\widetilde{\text{CNN}}$ do not further improve the MSE of the MLR and the CNN. Regarding CRPS, the non-linear CNN maintains a lower CRPS in the first four weeks, underscoring the benefits of non-linearity in the enhancing probabilistic forecast skill. When the stochastic perturbations are introduced, the impact of non-linearity on CRPS diminishes with increasing lead weeks. As for SSR, the MLR
<figure><img src="image_14.png" /><figcaption>(a) $\Delta_r$MSE maps of MLR and ECMWF, CNN and ECMWF, CNN and MLR at lead week 4, 5 and 6</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) $\Delta_r$MSE maps of $\overline{MLR}$ and ECMWF, $\overline{CNN}$ and ECMWF, $\overline{CNN}$ and $\overline{MLR}$ at lead week 4, 5 and 6</figcaption></figure>

FIG. A3. Same as figures 7 and 8, but for the $\Delta_r$MSE from week 4 to week 6
<figure><img src="image_16.png" /><figcaption>(a) $\Delta_r$CRPS maps of MLR and ECMWF, CNN and ECMWF, CNN and MLR at lead week 4, 5 and 6</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) $\Delta_r$CRPS maps of $\tilde{MLR}$ and ECMWF, $\tilde{CNN}$ and ECMWF, $\tilde{CNN}$ and $\tilde{MLR}$ at lead week 4, 5 and 6</figcaption></figure>

FIG. A4. Same as figure 7 and 8, but for the $\Delta_r$CRPS from week 4 to week 6
<figure><img src="image_18.png" /><figcaption>FIG. A5. The same as figure 6, but on forecasts</figcaption></figure>

and the CNN's SSR are under-dispersive. The addition of the stochastic perturbations significantly
improve the reliability of forecasts, but insufficiently.

References

Akiba, T., S. Sano, T. Yanase, T. Ohta, and M. Koyama, 2019: Optuna: A next-generation hyperparameter optimization framework. *Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining*, 2623–2631.

Alonzo, B., H.-K. Ringkjob, B. Jourdier, P. Drobinski, R. Plougonven, and P. Tankov, 2017: Modelling the variability of the wind energy resource on monthly and seasonal timescales. *Renewable energy*, **113**, 1434–1446.

Bauer, P., A. Thorpe, and G. Brunet, 2015: The quiet revolution of numerical weather prediction. *Nature*, **525 (7567)**, 47–55.

Bi, K., L. Xie, H. Zhang, X. Chen, X. Gu, and Q. Tian, 2022: Pangu-weather: A 3d high-resolution model for fast and accurate global weather forecast. *arXiv preprint arXiv:2211.02556*.

Bouget, V., D. Béréziat, J. Brajard, A. Charantonis, and A. Filoche, 2021: Fusion of rain radar images and wind forecasts in a deep learning model applied to rain nowcasting. *Remote Sensing*, **13 (2)**, 246.

Büeler, D., R. Beerli, H. Wernli, and C. M. Grams, 2020: Stratospheric influence on ecmwf sub-seasonal forecast skill for energy-industry-relevant surface weather in european countries. *Quarterly Journal of the Royal Meteorological Society*, **146 (733)**, 3675–3694.
Buizza, R., and M. Leutbecher, 2015: The forecast skill horizon. *Quarterly Journal of the Royal Meteorological Society*, **141** (693), 3366–3382.

Cassola, F., and M. Burlando, 2012: Wind speed and wind energy forecast through kalman filtering of numerical weather prediction model output. *Applied energy*, **99**, 154–166.

Chang, W.-Y., and Coauthors, 2014: A literature review of wind forecasting methods. *Journal of Power and Energy Engineering*, **2** (04), 161.

Chen, J., F. P. Brissette, and R. Leconte, 2011: Uncertainty of downscaling method in quantifying the impact of climate change on hydrology. *Journal of hydrology*, **401** (3-4), 190–202.

Chen, L., X. Zhong, F. Zhang, Y. Cheng, Y. Xu, Y. Qi, and H. Li, 2023: Fuxi: A cascade machine learning forecasting system for 15-day global weather forecast. *arXiv preprint arXiv:2306.12873*.

Cortesi, N., V. Torralba, N. González-Reviriego, A. Soret, and F. J. Doblas-Reyes, 2019: Characterization of european wind speed variability using weather regimes. *Climate Dynamics*, **53** (7), 4961–4976.

Ferro, C., 2014: Fair scores for ensemble forecasts. *Quarterly Journal of the Royal Meteorological Society*, **140** (683), 1917–1923.

Goddard, L., and Coauthors, 2013: A verification framework for interannual-to-decadal predictions experiments. *Climate Dynamics*, **40**, 245–272.

Goutham, N., R. Plougonven, H. Omrani, S. Parey, P. Tankov, A. Tantet, P. Hitchcock, and P. Drobinski, 2022: How skillful are the european subseasonal predictions of wind speed and surface temperature? *Monthly Weather Review*, **150** (7), 1621–1637.

Goutham, N., R. Plougonven, H. Omrani, A. Tantet, S. Parey, P. Tankov, P. Hitchcock, and P. Drobinski, 2023: Statistical downscaling to improve the subseasonal predictions of energy-relevant surface variables. *Monthly Weather Review*, **151** (1), 275–296.

Haiden, T., M. Janousek, F. Vitart, Z. Ben-Bouallegue, and F. Prates, 2023: 911. Evaluation of ecmwf forecasts, including the 2023 upgrade. ECMWF, URL &nbsp;; https://doi.org/10.21957/d47ba5263c.
Hersbach, H., and Coauthors, 2018a: Era5 hourly data on pressure levels from 1979 to present. *Copernicus climate change service (c3s) climate data store (cds)*, **10** (10.24381).

Hersbach, H., and Coauthors, 2018b: Era5 hourly data on single levels from 1979 to present. *Copernicus climate change service (c3s) climate data store (cds)*, **10** (10.24381).

Hersbach, H., and Coauthors, 2020: The era5 global reanalysis. *Quarterly Journal of the Royal Meteorological Society*, **146** (730), 1999–2049.

Kurth, T., and Coauthors, 2023: Fourcastnet: Accelerating global high-resolution weather forecasting using adaptive fourier neural operators. *Proceedings of the Platform for Advanced Scientific Computing Conference*, 1–11.

Lam, R., and Coauthors, 2022: Graphcast: Learning skillful medium-range global weather forecasting. *arXiv preprint arXiv:2212.12794*.

Laurila, T. K., V. A. Sinclair, and H. Gregow, 2021: Climatology, variability, and trends in near-surface wind speeds over the north atlantic and europe during 1979–2018 based on era5. *International Journal of Climatology*, **41** (4), 2253–2278.

Liu, Y., S. Feng, Y. Qian, H. Huang, and L. K. Berg, 2023: How do north american weather regimes drive wind energy at the sub-seasonal to seasonal timescales? *npj Climate and Atmospheric Science*, **6** (1), 100.

Manzanas, R., J. M. Gutiérrez, J. Bhend, S. Hemri, F. J. Doblas-Reyes, V. Torralba, E. Penabad, and A. Brookshaw, 2019: Bias adjustment and ensemble recalibration methods for seasonal forecasting: A comprehensive intercomparison using the c3s dataset. *Climate Dynamics*, **53**, 1287–1305.

Maraun, D., and Coauthors, 2010: Precipitation downscaling under climate change: Recent developments to bridge the gap between dynamical models and the end user. *Reviews of geophysics*, **48** (3).

Mariotti, A., and Coauthors, 2020: Windows of opportunity for skillful forecasts subseasonal to seasonal and beyond. *Bulletin of the American Meteorological Society*, **101** (5), E608–E625.
Matheson, J. E., and R. L. Winkler, 1976: Scoring rules for continuous probability distributions. *Management science*, **22** (10), 1087–1096.

Nguyen, B. M., G. Tian, M.-T. Vo, A. Michel, T. Corpetti, and C. Granero-Belinchon, 2022: Convolutional neural network modelling for modis land surface temperature super-resolution. *2022 30th European Signal Processing Conference (EUSIPCO)*, IEEE, 1806–1810.

Orth, R., and S. I. Seneviratne, 2014: Using soil moisture forecasts for sub-seasonal summer temperature predictions in europe. *Climate dynamics*, **43**, 3403–3418.

Pan, B., K. Hsu, A. AghaKouchak, and S. Sorooshian, 2019: Improving precipitation estimation using convolutional neural network. *Water Resources Research*, **55** (3), 2301–2321.

Paszke, A., and Coauthors, 2019: Pytorch: An imperative style, high-performance deep learning library. *Advances in neural information processing systems*, **32**.

Price, I., and Coauthors, 2023: Gencast: Diffusion-based ensemble forecasting for medium-range weather. *arXiv preprint arXiv:2312.15796*.

Rasp, S., and Coauthors, 2023: Weatherbench 2: A benchmark for the next generation of data-driven global weather models. *arXiv preprint arXiv:2308.15560*.

Robertson, A. W., A. Kumar, M. Peña, and F. Vitart, 2015: Improving and promoting subseasonal to seasonal prediction. *Bulletin of the American Meteorological Society*, **96** (3), ES49–ES53.

Rodrigues, E. R., I. Oliveira, R. Cunha, and M. Netto, 2018: Deepdownscale: A deep learning strategy for high-resolution weather forecast. *2018 IEEE 14th International Conference on e-Science (e-Science)*, IEEE, 415–422.

Ronneberger, O., P. Fischer, and T. Brox, 2015: U-net: Convolutional networks for biomedical image segmentation. *Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18*, Springer, 234–241.

Sachindra, D., K. Ahmed, M. M. Rashid, S. Shahid, and B. Perera, 2018: Statistical downscaling of precipitation using machine learning techniques. *Atmospheric research*, **212**, 240–258.
Salameh, T., P. Drobinski, M. Vrac, and P. Naveau, 2009: Statistical downscaling of near-surface wind over complex terrain in southern France. *Meteorology and Atmospheric Physics*, **103 (1)**, 253–265.

Simmons, A. J., and B. J. Hoskins, 1979: The downstream and upstream development of unstable baroclinic waves. *Journal of the Atmospheric Sciences*, **36 (7)**, 1239–1254.

Tawn, R., J. Browell, and D. McMillan, 2022: Subseasonal-to-seasonal forecasting for wind turbine maintenance scheduling. *Wind, 2 (2)*, 260–287.

Toth, Z., and R. Buizza, 2019: Weather forecasting: What sets the forecast skill horizon? *Sub-Seasonal to Seasonal Prediction*, Elsevier, 17–45.

Trebing, K., T. Stanczyk, and S. Mehrkanoon, 2021: Smaat-unet: Precipitation nowcasting using a small attention-unet architecture. *Pattern Recognition Letters*, **145**, 178–186.

Vigaud, N., A. W. Robertson, and M. K. Tippett, 2017: Multimodel ensembling of subseasonal precipitation forecasts over north america. *Monthly Weather Review*, **145 (10)**, 3913–3928.

Vitart, F., 2014: Evolution of ecmwf sub-seasonal forecast skill scores. *Quarterly Journal of the Royal Meteorological Society*, **140 (683)**, 1889–1899.

Vitart, F., and A. W. Robertson, 2019: Introduction: Why sub-seasonal to seasonal prediction (s2s)? *Sub-Seasonal to Seasonal Prediction*, Elsevier, 3–15.

Vitart, F., and Coauthors, 2019: 854. Extended-range prediction. ECMWF, URL https://www.ecmwf.int/node/19286, https://doi.org/10.21957/pdivp3t9m.

Wallace, J. M., and D. S. Gutzler, 1981: Teleconnections in the geopotential height field during the northern hemisphere winter. *Monthly weather review*, **109 (4)**, 784–812.

Wang, Z., A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, 2004: Image quality assessment: from error visibility to structural similarity. *IEEE transactions on image processing*, **13 (4)**, 600–612.

White, C. J., and Coauthors, 2017: Potential applications of subseasonal-to-seasonal (s2s) predictions. *Meteorological applications*, **24 (3)**, 315–325.

Wilby, R. L., and T. M. Wigley, 1997: Downscaling general circulation model output: a review of methods and limitations. *Progress in physical geography*, **21 (4)**, 530–548.
Wilks, D. S., 2019: *Statistical Methods in the Atmospheric Sciences*. 4th ed., Elsevier, Amsterdam, https://doi.org/10.1016/C2017-0-03921-6.

WMO, A., 2012: Guidelines on ensemble prediction systems and forecasting. *World Meteorological Organization Weather Climate and Water*, **1091**.

Zamo, M., and P. Naveau, 2018: Estimation of the continuous ranked probability score with limited information and applications to ensemble weather forecasts. *Mathematical Geosciences*, **50 (2)**, 209–234.

Zhang, C., 2005: Madden-julian oscillation. *Reviews of Geophysics*, **43 (2)**.

Zhu, Y., and Coauthors, 2018: Toward the improvement of subseasonal prediction in the national centers for environmental prediction global ensemble forecast system. *Journal of Geophysical Research: Atmospheres*, **123 (13)**, 6732–6745.