Bivariate Gaussian models for wind vectors in a distributional regression framework

TL;DR

This paper introduces a novel probabilistic post-processing method for wind vectors using bivariate Gaussian models within a distributional regression framework, improving forecast accuracy by modeling all distribution parameters simultaneously.

Contribution

The study presents a new method that models all wind vector distribution parameters jointly, including correlation, with flexible regression splines, and incorporates ensemble forecast data for enhanced accuracy.

Findings

Improved predictive skill for wind vectors across various terrains.

Rotation-allowing model outperforms baseline in all tested sites.

Small benefits observed from modeling correlation in alpine valleys.

Abstract

A new probabilistic post-processing method for wind vectors is presented in a distributional regression framework employing the bivariate Gaussian distribution. In contrast to previous studies all parameters of the distribution are simultaneously modeled, namely the means and variances for both wind components and also the correlation coefficient between them employing flexible regression splines. To capture a possible mismatch between the predicted and observed wind direction, ensemble forecasts of both wind components are included using flexible two-dimensional smooth functions. This encompasses a smooth rotation of the wind direction conditional on the season and the forecasted ensemble wind direction. The performance of the new method is tested for stations located in plains, mountain foreland, and within an alpine valley employing ECMWF ensemble forecasts as explanatory variables for all distribution parameters. The rotation-allowing model shows distinct improvements in terms of predictive skill for all sites compared to a baseline model that post-processes each wind component separately. Moreover, different correlation specifications are tested and small improvements compared to the model setup with no estimated correlation could be found for stations located in alpine valleys.

Tables

Table 1. Table 1 : Overview of bivariate Gaussian model specifications. For the ‘baseline model’ ( BLM-0 , see Sect.  2.2 ) and the ‘rotation-allowing model’ ( RAM-0 , see Sect.  2.3 ) no correlation is employed, i.e., fixed at zero. For all tested correlation specifications, the RAM-0 setup is employed for the location and scale part (see Sect.  2.4 ). In all setups for each distribution parameter, a seasonally varying intercept effect is estimated. For the BLM-0 , a seasonally varying slope coefficient is fitted for the two wind components in the location and scale parts. For the setup RAM-0 , the slope coefficients are additionally dependent on the wind direction. In the correlation model RAM-ADV , the wind speed is modeled conditional on the wind direction. The ditto symbol ‘–  "  –’ implies that the same configuration as in the line above is employed.

Name	Location part	Scale part	Correlation part
BLM-0	${\mu }_{\star }\sim {\text{vec}}_{\star ,mean}$	${\sigma }_{\star }\sim {\text{vec}}_{\star ,log.sd}$	$\rho =0$
RAM-0	${\mu }_{\star }\sim {\text{vec}}_{1,mean},{\text{vec}}_{2,mean}$	${\sigma }_{\star }\sim {\text{vec}}_{1,log.sd},{\text{vec}}_{2,log.sd}$	$\rho =0$
RAM-EMP	– " –	– " –	$\rho ={\text{vec}}_{\star ,corr}$
RAM-IC	– " –	– " –	$\rho \sim 1$
RAM-DIR	– " –	– " –	$\rho \sim {\text{dir}}_{mean}$
RAM-ADV	– " –	– " –	$\rho \sim {\text{dir}}_{mean},{\text{spd}}_{mean}$

Equations

$$L(\mathbf{\mu },\Sigma |\mathbf{y})=\frac{1}{\sqrt{(2\pi {)}^2|\mathbf{\Sigma }|}}\exp \left(-\frac{1}{2}(\mathbf{y}-\mathbf{\mu }{)}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{y}-\mathbf{\mu })\right),$$

$$\mathbf{\Sigma }=\left(\begin{array}{cc}{\sigma }_1^2& \rho {\sigma }_1{\sigma }_2\\ \rho {\sigma }_1{\sigma }_2& {\sigma }_2^2\end{array}\right),$$

$$\begin{array}{rl}{\mu }_{\star }& =\underset{\text{intercept}}{\underbrace{{\alpha }_0+f_0(\text{doy})}}+\underset{\text{slope coefficient}}{\underbrace{({\alpha }_1+f_1(\text{doy}))}}\cdot {\text{vec}}_{\star ,mean},\\ \text{log}({\sigma }_{\star })& =\underset{\text{intercept}}{\underbrace{{\beta }_0+g_0(\text{doy})}}+\underset{\text{slope coefficient}}{\underbrace{({\beta }_1+g_1(\text{doy}))}}\cdot {\text{vec}}_{\star ,log.sd},\end{array}$$

$$\begin{array}{rl}{\mu }_{\star }={\alpha }_0+f_0(\text{doy})& +({\alpha }_1+f_1(\text{doy})\cdot f_2({\text{dir}}_{mean}))\cdot {\text{vec}}_{1,mean}\\ & +({\alpha }_2+f_3(\text{doy})\cdot f_4({\text{dir}}_{mean}))\cdot {\text{vec}}_{2,mean},\\ \text{log}({\sigma }_{\star })={\beta }_0+g_0(\text{doy})& +({\beta }_1+g_1(\text{doy})\cdot g_2({\text{dir}}_{mean}))\cdot {\text{vec}}_{1,log.sd}\\ & +({\beta }_2+g_3(\text{doy})\cdot g_4({\text{dir}}_{mean}))\cdot {\text{vec}}_{2,log.sd},\end{array}$$

$$\text{rhogit}(\rho )$$

$$\eta ={\alpha }_0+\underset{\text{linear effect for x1}}{\underbrace{{\alpha }_1\cdot {\text{x}}_1}}+\underset{\text{pot. non-linear effect for x2}}{\underbrace{f_1({\text{x}}_2)}},$$

$$\eta$$

$$\text{skill score}=1-\frac{{\text{score}}_{\text{forecast}}}{{\text{score}}_{\text{reference}}}.$$

$$\text{ES}(F,\mathbf{y})=\frac{1}{m}\sum _{i=1}^m\Vert {\mathbf{X}}_i-\mathbf{y}\Vert -\frac{1}{2m^2}\sum _{i=1}^m\sum _{j=1}^m\Vert {\mathbf{X}}_i-{\mathbf{X}}_j\Vert ,$$

$$\text{LS}(F,\mathbf{y})=\log (f(\mathbf{y})),$$

Full text

Bivariate Gaussian models for wind vectors in a distributional regression framework

Moritz N. Lang

Universität Innsbruck &Georg J. Mayr

Universität Innsbruck &Reto Stauffer

Universität Innsbruck &Achim Zeileis

Universität Innsbruck

moritz.lang@uibk.ac.at

\Plaintitle

Bivariate Gaussian models for wind vectors

\PlainauthorMoritz N. Lang, Georg J. Mayr, Reto Stauffer, Achim Zeileis \Abstract A new probabilistic post-processing method for wind vectors is presented in a distributional regression framework employing the bivariate Gaussian distribution. In contrast to previous studies all parameters of the distribution are simultaneously modeled, namely the means and variances for both wind components and also the correlation coefficient between them employing flexible regression splines. To capture a possible mismatch between the predicted and observed wind direction, ensemble forecasts of both wind components are included using flexible two-dimensional smooth functions. This encompasses a smooth rotation of the wind direction conditional on the season and the forecasted ensemble wind direction. The performance of the new method is tested for stations located in plains, mountain foreland, and within an alpine valley employing ECMWF ensemble forecasts as explanatory variables for all distribution parameters. The rotation-allowing model shows distinct improvements in terms of predictive skill for all sites compared to a baseline model that post-processes each wind component separately. Moreover, different correlation specifications are tested and small improvements compared to the model setup with no estimated correlation could be found for stations located in alpine valleys. \Keywordsbivariate Gaussian distribution, distributional regression, regression splines, wind vectors, ensemble forecasts \Plainkeywordsbivariate Gaussian distribution, distributional regression, regression splines, wind vectors, ensemble forecasts \Address Moritz N. Lang

Department of Statistics

Department of Atmospheric and Cryospheric Sciences

Universität Innsbruck

6020 Innsbruck, Austria

E-mail:

1 Introduction

Accurate forecasts of wind speed and direction are of great importance for decision making processes and risk management in today’s society and will likely become more important in the future. Not only because of the rapid change in climate and the resulting increase of severe storms (e.g., Kunkel et al., 2012; Vose et al., 2013), but also due to the change of the society itself and its technical revolution. As an example, the European Union is aiming to increase the amount of wind energy by 2030 to 35%, which would be more than double the capacity installed at the end of 2016 (WindEurope, 2017). In the field of aviation and air traffic control for instance, more flexible landing procedures with a so-called time-based separation are currently tested at Heathrow Airport and are planned to go operational in the near future (EuropeanCommission, 2018). In both cases, wind (power) forecasts are of fundamental importance; probabilistic wind forecasts are in particular advisable as they permit optimal risk assessment and decision making (Gneiting, 2008).

Probabilistic weather forecasts are usually issued in form of ensemble predictions. To account for the underlying uncertainty in the atmosphere, numerical ensemble prediction systems (EPSs) provide a set of weather forecasts using slightly perturbed initial conditions and different model parameterizations (Palmer, 2002). Despite recent advances in the development of EPSs, the resulting forecasts still often show displacement errors and usually capture only part of the forecast uncertainty, especially when comparing EPS forecasts and point measurements (Buizza et al., 2005; Gneiting and Katzfuss, 2014). This often results from structural model deficiencies, and insufficient resolution or unresolved topographical features. To remove systematic biases and to provide corrected variance information, statistical post-processing methods are often employed. For wind, various ensemble post-processing methods have been proposed over the last decade, mainly focusing on wind speed. For a single location, parametric examples are non-homogeneous regression (Thorarinsdottir and Gneiting, 2010; Lerch and Thorarinsdottir, 2013; Baran and Lerch, 2015, 2016), kernel dressing methods with similarities to Bayesian model averaging (Sloughter et al., 2010; Courtney et al., 2013; Baran, 2014), and extended logistic regression (Messner et al., 2014a, b). A non-parametric approach based on quantile regression forests was applied by Taillardat et al. (2016). On a regular grid, ensemble post-processing based on non-homogeneous regression was performed by Scheuerer and Möller (2015).

To account for the circular characteristics of wind or utilizing information of wind speed and direction, an intuitive post-processing approach is to model a bivariate process for the zonal and meridional wind components. Gneiting et al. (2008) suggested to use a bivariate Gaussian response distribution for the wind components, which Pinson (2012) implemented. He estimates a dilation and translation factor for the individual ensemble members utilizing the empirical correlation structure of the EPS. This procedure can be seen as a variant of the ensemble copula coupling (ECC) method introduced by Schefzik et al. (2013). With the ECC, both wind components are calibrated with univariate approaches and a discrete sample drawn from each univariate predictive distribution is rearranged in the rank order structure of the raw ensemble. The method introduced by Schuhen et al. (2012) also fits a bivariate Gaussian distribution for the wind components; however, in their approach the post-processed probabilistic forecast consists of a fully specified predictive distribution instead of a discrete ensemble. As their analyses show that the correlation mainly depends on wind direction, they model the correlation as a trigonometric function of the ensemble mean wind direction, separately for different wind sectors. In addition, an extra group is formed for cases with low wind speeds unconditionally on their wind direction. The estimation of the correlation parameter is done offline in a pre-processing step for a separate year either for all stations combined or for each station individually. However according to Schuhen et al. (2012), the fitting can be critical for individual stations since wind sectors may only contain few data points.

For stations in complex terrain, a possible drawback of the bivariate post-processing approaches introduced by Pinson (2012) and Schuhen et al. (2012) is that the models are not able to correct for a systematic distortion in the wind directions due to discrepancies between the model and real topography. Especially, when the respective valley orientations differ a meridional wind component might be partially rotated into a zonal wind component and vice versa. In the field of post-processing deterministic weather forecasts, Glahn and Lowry (1972) already suggested to employ both forecasted wind components for the calibration of the zonal or meridional wind component, which is partly able to correct for systematic distortions in wind directions in a linear manner. However, to our knowledge, this approach has not yet been further pursued in the field of bivariate probabilistic post-processing of wind vectors.

Alternatively to bivariate calibration methods, wind direction can also be employed in univariate settings. In a post-processing approach for wind speed, Eide et al. (2017) suggest to utilize the potentially non-linear information of the wind direction by a generalized additive model (GAM; Hastie and Tibshirani, 1986). GAMs were first applied in the meteorological context by Vislocky and Fritsch (1995), and provide a powerful statistical model framework which can capture potential non-linear relationships between the covariates and the response by smooth functions or splines. Eide et al. (2017) employ wind direction as an additional covariate for the estimation of wind speed, by accounting for its cyclic and potential non-linear characteristics utilizing thin-plate regression splines.

In this study, we directly model the zonal and meridional wind components employing the bivariate Gaussian distribution as suggested by Gneiting et al. (2008), and performed by Pinson (2012) and Schuhen et al. (2012). However, we capture all distribution parameters, namely the location and scale parameters for both wind components, and also the correlation coefficient between them in a single flexible model. In the estimation of the two-dimensional location and scale parameters the information value of both ensemble wind components is utilized to allow for a smooth rotation of the forecasted wind direction accounting for unresolved topographical features. To consider the correlation characteristics detected in Schuhen et al. (2012) and to allow for possible non-linear effects, as e.g., suggested by Eide et al. (2017), we model the correlation as a function of wind speed and direction utilizing cyclic regression splines. To account for potential time-varying effects, all linear predictors use a time-adaptive intercept and time-adaptive slope coefficients based on cyclic smooth splines.

The paper is structured as follows: Section 2 introduces the employed statistical models. The underlying data of this study are shortly described in Sect. 3. Within Sect. 4, first a model comparison and validation is presented for two weather stations with different site characteristics, followed by aggregated scores for station sites located in plains, mountain foreland, and within an alpine valley. The article ends with a brief discussion and a conclusion given in Sect. 5.

2 Methods

In Sect. 2.1, the bivariate Gaussian distribution is reviewed and briefly presented in a distributional regression framework. Subsequently, three broad model classes are introduced all of which are based on a time-adaptive training scheme but employ different specifications for the location, scale, and correlation parameters of the bivariate distribution. First, the baseline model is presented in Sect. 2.2 that serves as a benchmark and simply combines two univariate heteroscedastic regression models that post-process each wind component separately. Second, the baseline model is extended in Sect. 2.3 by always adding both EPS wind components as regressors using smooth splines and thus allowing for potential misspecifications in the EPS wind direction. Finally, Sect. 2.4 also considers models with estimated correlation coefficients based on various regression specifications. Table 1 provides a synoptic summary of all bivariate Gaussian model specifications tested within this study.

2.1 Distributional regression for a bivariate Gaussian response

The zonal and meridional components of the horizontal wind vector are represented by a bivariate Gaussian distribution. Its likelihood function $L$ is given by

[TABLE]

where $\boldsymbol{y}=(y_{1},y_{2})^{\top}$ are bivariate observations and $\boldsymbol{\mu}=(\mu_{1},\mu_{2})^{\top}$ the distributional location parameters with $\mu_{\star}\in\mathbb{R}$; the subscript asterisk acts as a placeholder for the zonal and meridional wind component from here on. The covariance matrix is defined as

[TABLE]

with correlation parameter $\rho\in[-1,1]$ and scale parameters $\sigma_{\star}>0$. In the framework of distributional regression, the location parameters $\mu_{\star}$, the scale parameters $\sigma_{\star}$, and the correlation parameter $\rho$ are linked to additive predictors by an identity, logarithm and rhogit link, respectively (Klein et al., 2014).

To be able to utilize the information of cyclic covariates as e.g., wind direction in addition to linear covariates, we follow Eide et al. (2017) and fit a GAM but utilize cubic smooth functions with cyclic constraints for all cyclic covariates (Wood 2017; see Appendix A.1). In the context of distributional regression, additive models with smooth effects are typical referred to as ‘generalized additive models for location, scale and shape’ (GAMLSS, Rigby and Stasinopoulos, 2005). In this study, we utilize GAMLSS in a Bayesian framework, which allows us to examine the estimated effects based on Markov chain Monte Carlo (MCMC) simulations. A comprehensive summary of the method is given in Umlauf et al. (2018).

2.2 Baseline model (BLM-0)

The baseline model (BLM-0) combines two univariate heteroscedastic regression models that post-process each wind component separately with correlation fixed at zero. Hence, for the location and scale part, it uses its direct counterparts of the EPS as covariates, namely EPS-forecasted zonal wind information ($\text{vec}_{1}$) to model the zonal component of the bivariate response, and EPS-forecasted meridional wind information ($\text{vec}_{2}$) to model the meridional component:

[TABLE]

where, $\alpha_{\bullet}$ and $\beta_{\bullet}$ are regression coefficients, and $f_{\bullet}(\text{doy})$ and $g_{\bullet}(\text{doy})$ employ cyclic regression splines conditional on the day of the year (doy). The subscripts $mean$ and $log.sd$ refer to mean and log standard deviation of the ensemble wind components, respectively. We follow Gebetsberger et al. (2017) and use the logarithm transformation for the standard deviation of the ensemble members to assure positivity which is preferable for the estimation process.

Equation (3) specifies a time-adaptive training scheme (with further details in Appendix A.2), where the linear predictors consist of a global intercept and slope coefficient plus a seasonally varying deviation. Thus, the intercept and slope coefficients can smoothly evolve over the year in case that the bias or the covariate’s skill varies seasonally. If there is no seasonal variation, the non-linear effects become zero and Eq. (3) simplifies to a regression model with a constant intercept and slope coefficient ($\mu_{\star}=\alpha_{0}+\alpha_{1}\cdot\text{vec}_{\star,mean};~{}\text{log}(\sigma_{\star})=\beta_{0}+\beta_{1}\cdot\text{vec}_{\star,log.sd}$).

2.3 Rotation-allowing model without correlation (RAM-0)

In the second model, labeled as rotation-allowing model (RAM-0), we extend the setup BLM-0 by employing the zonal and meridional wind information of the ensemble for the linear predictors of all location and scale parameters. That means we use the ensemble information of both the zonal and meridional wind components for the two components of the response (cf. Glahn and Lowry, 1972). In case of a perfect EPS the zonal wind predictions are non-informative covariates for the meridional wind component and vice versa. However, if, e.g., the model topography is not sufficiently resolved or in case of local shadowing effects, both EPS wind components may contain valuable information for the zonal and meridional wind components of the response. Especially in a mountain valley, when the model and real valley orientation differs, both wind components of the ensemble can potentially contain information about both location and scale parameters, respectively. Thus, we propose to employ seasonally varying effects depending on the ensemble wind direction, which allows the model to rotate the forecasted wind direction if necessary. To do so, we obtain a two-dimensional smooth function with a respective cyclic constraint for the day of the year (doy) and for the mean ensemble wind direction ($\text{dir}_{mean}$):

[TABLE]

where, as before, $\alpha_{\bullet}$ and $\beta_{\bullet}$ are regression coefficients, and $f_{\bullet}$ and $g_{\bullet}$ employ cyclic regression splines. From a more physical perspective, the two-dimensional smooth effects rotate the ensemble wind components conditional on the day of the year and the ensemble wind direction.

2.4 Rotation-allowing models with correlation

By explicitly modeling the correlation, we further extend the RAM-0 setup within this section. For the estimation of the correlation structure different model specifications are tested. The most advanced specification, RAM-ADV, assumes that the correlation mainly depends on the mean ensemble wind direction ($\text{dir}_{mean}$) and speed ($\text{spd}_{mean}$) by modeling a linear interaction between these two covariates:

[TABLE]

with $\text{rhogit}(\rho)=\rho/\sqrt{(1-\rho^{2})}$; $\gamma_{0}$ is the global intercept and $h_{0}(\text{doy})$ the seasonally varying intercept. The effect $h_{1}(\text{dir}_{mean})$ estimates the dependence of the correlation given the wind direction and $(\gamma_{1}+h_{2}(\text{dir}_{mean}))\cdot\text{spd}_{mean}$ employs a varying effect of wind speed conditional on the wind direction. The estimation of the underlying correlation structure is in accordance to results of Schuhen et al. (2012), who employ wind direction and an offset of wind speed as informative covariates in the estimation of the correlation parameter.

Other implementations tested for the correlation parameter are an intercept-only model (RAM-IC), a model with a cyclic effect solely depending on wind direction (RAM-DIR), the univariate model RAM-0 of Sect. 2.3, and a univariate model using the empirical correlation ($corr$) of the raw ensemble (RAM-EMP). A synoptic table of all models tested in this study is given in Table 1.

3 Data

3.1 Observational data

The validation and comparison of the different model specifications is performed for $15$ measurement sites located across Austria, Germany, and Switzerland. The sites are chosen to investigate the influence of different underlying topographies or varying discrepancies between the real and the model topography on the post-processing. The stations are divided into three groups representing sites located in plains, mountain foreland, and within an alpine valley. An overview of the stations is given in Fig. 1. The results for the stations Hamburg and Innsbruck, which are labeled in Fig. 1, are discussed in more detail in Sect. 4. At all meteorological sites, wind speed and direction measurements are reported for the $10$ m height level. The data are $10$-minute averages for the period January 26, 2010 to March 7, 2016, yielding a total of 2233 days.

3.2 Ensemble prediction system

Covariates are derived from the global 50-member EPS of the European Centre for Medium-Range Weather Forecasts (ECMWF). These EPS forecasts have a horizontal resolution of approximately $30$ km (T639) for the time between January 2010 and March 2016 and are bilinearly interpolated to the measurement sites. Covariates employed in this study are the zonal and meridional wind components as well as the derived quantities wind speed and direction valid at $10$ m above ground. For all these variables, two statistics are computed over the 50 perturbed ensemble members, namely the $mean$ and the logarithm of the standard deviation ($log.sd$). Additionally, the empirical correlation ($corr$) is computed from the raw ensemble members to capture their rank dependence structure. Forecasts are taken from the EPS run initialized at $00$ UTC for forecast steps ranging from $+12$ h to $+72$ h ahead on a $12$ hourly temporal resolution. Figure 2 shows the empirical wind distributions of the observed and predicted winds for Innsbruck and Hamburg for forecast steps $+12$ h and $+24$ h corresponding to $12$ UTC and $00$ UTC.

4 Results

This section presents the results of the statistical post-processing models. The structure is as follows: First, the estimated effects of the baseline model BLM-0 (Sect. 4.1) and the rotation-allowing model RAM-0 (Sect. 4.2) are shown for two stations representative for one alpine valley site and for a measurement site in the plains. For both models a constant correlation of zero is employed, and their predictive performance is discussed in Sect. 4.3. Afterwards, a model comparison (Sect. 4.4) and validation (Sect. 4.5) of the different correlation specifications is given for the two representative stations. In Sect. 4.6, the overall performance of the model setups is evaluated for three groups of stations classified as topographically plain, mountain foreland, and alpine valley sites.

The model estimation is performed on data of the first $4$ years, leaving an out-of-sample validation data set ranging from February 24, 2014 to March 7, 2016.

4.1 Baseline model (BLM-0)

For BLM-0, the cyclic seasonal effects for the stations Hamburg and Innsbruck are shown in Fig. 3 as solid and dashed lines with the respective $95$% credible intervals. The estimated effects are on the scale of the additive predictor; i.e., on the linear scale for the location parameters $\mu_{\star}$ and on the log-scale for the scale parameters $\sigma_{\star}$. Each of the four distribution parameters is described by a (potentially) seasonally varying effect for the intercept (left column) and the slope coefficient (right column) as specified in Eq. (3).

For Hamburg, for both location parameters $\mu_{\star}$, the intercept effect is almost zero (Fig. 3a, c) and the effect for the slope coefficient is close to one (Fig. 3b, d) with very little seasonal variability. This means apparently no bias correction is necessary and the ensemble mean wind components are mapped almost one-to-one to the location parameters. Similarly, barely any seasonal variation exists for the scale parameters $\sigma_{\star}$ (Fig. 3e–h); however, here the intercept and slope coefficients actually post-process the EPS variances of the wind components (rather than a one-to-one mapping only), leading to an increase of the scale parameters compared to the under-dispersed ensemble. The $95$% credible intervals indicate a higher uncertainty of the estimated scale parameters compared to the location parameters. In summary, the EPS performance for Hamburg is almost constant over the year and no time-adaptive training scheme seems to be necessary.

On the contrary, for Innsbruck the estimated effects show a distinct annual cycle for the location parameters $\mu_{\star}$, which indicates a varying information content of the predictor variables $\text{vec}_{\star,mean}$ and the need of some adaptive training scheme. For the location parameter $\mu_{1}$, the intercept is rather large during winter (Fig. 3a) while, at the same time, the slope coefficient (Fig. 3b) is close to zero due to an apparently low skill of the EPS. For the location parameter $\mu_{2}$ (Fig. 3c, d), the higher slope coefficients during spring and autumn suggest a higher information content of the raw EPS in the transitional seasons than for the rest of the year. For the scale parameters (Fig. 3e–h), the estimated effects show high variability; this indicates a seasonally varying skill of the EPS variance information. In summary, for the weather station in Innsbruck, the information content of the ensemble wind predictions seem to be rather low. This is in accordance with the clearly different pattern of observations and EPS forecasts shown in Fig. 2.

4.2 Rotation-allowing model without correlation (RAM-0)

Figure 4 shows the estimated mean effects of the setup RAM-0 in comparison with BLM-0 on the wind direction at Innsbruck and Hamburg for the forecast steps $+12$ h and $+24$ h valid at $12$ UTC and $00$ UTC, respectively. The marginal effects are non-centered and shown for the mean covariates within 10 degree wind sectors conditional on the day of the year. The BLM-0 model for Innsbruck shows a distinct seasonal dependency of the post-processed wind direction for both times of the day (Fig. 4a, c). During winter at $12$ UTC and $00$ UTC, mainly down-valley winds (approximately $280$ degree) are predicted, whereas over the rest of the year, the EPS mainly forecasts up-valley wind directions. This pattern is more pronounced during night ($00$ UTC) and has less variability in summer. In general, the setup BLM-0 seems to mainly capture the climatological mean wind direction; this leads to little variations between the different wind directions issued by the EPS. In contrast, the rotation-allowing setup RAM-0 has the flexibility to post-process the wind directions conditional on the forecasted EPS wind direction, which is apparent for the station Innsbruck at both $12$ UTC and $00$ UTC: For $12$ UTC the seasonal dependency leads to either up-valley or down-valley wind conditional on the ensemble wind direction and on the day of the year (Fig. 4b); whereas at $00$ UTC almost no seasonal variation exists and the predicted wind direction solely depends on the issued ensemble wind forecasts (Fig. 4d). At Hamburg, a completely different picture can be seen: Almost no post-processing conditional on the ensemble wind direction or the day of the year is visible for both time steps and model setups (Fig. 4e–f). In other words, the predicted ensemble wind direction fits the observed wind direction quite well and only little statistical correction is needed.

4.3 Predictive performance – models without correlation

To investigate the predictive performance of the two competing setups, Fig. 5 shows the discretized logarithmic score (LS; see Appendix B) and the energy score (ES; see Appendix B) for the forecast steps from $+12$ h to $+72$ h on a $12$ hourly temporal resolution. In addition, skill scores are shown with the raw EPS as reference or comparing the different setups to each other. Both multivariate scores are proper scores (Gneiting and Raftery, 2007) and evaluate the full predictive distribution returned by the statistical models. The scores for the different forecast horizons show an overall better predictive performance at Hamburg than at Innsbruck. For both stations, the forecasts valid at $00$ UTC have more skill than those for $12$ UTC with higher diurnal variations at Innsbruck. In terms of the ES, the improvements of the BLM-0 model over the raw EPS is about $29$% for Innsbruck and $8$% for Hamburg (Fig. 5e, f). In terms of the LS (Fig. 5b, c), the skill scores are higher with an improvement of approximately $87$% and $33$% for Innsbruck and Hamburg, respectively. The predictive performance gain for the more flexible rotation-allowing setup RAM-0 compared to the specification BLM-0 is around $7$% for Innsbruck and $2$% for Hamburg in terms of the ES (Fig. 5e, f). The LS shows slightly less pronounced relative improvements for the more flexible setup ($4$% and $1$%; Fig. 5b, c). The distinct improvements in the scores for RAM-0 are as expected for Innsbruck due to a more flexible utilization of the ensemble information. For plain areas like Hamburg, we assume the better performance is based on an enhanced adjustment of the location parameters as both wind components are included in the linear predictors and not due to the smooth rotation (c.f. Fig. 4).

4.4 Rotation-allowing models with correlation

After investigating the two competing location or scale setups, we now focus on an extension of the model RAM-0 by explicitly estimating the underlying correlation structure. Different model specifications for the correlation parameter $\rho$ are tested employing the same linear predictors for $\mu_{\star}$ and $\sigma_{\star}$ (see Table 1).

Figure 6 shows correlation parameters predicted by different models for the forecast step $+12$ h for the full validation period. For comparison, the underlying correlation structure of the raw EPS is also shown. The latter is distributed similarly for Innsbruck and Hamburg and has almost the shape of a Gaussian distribution (Fig. 6a, e). The intercept-only model RAM-IC, with a varying intercept over the year, estimates correlations between $-0.27$ and $0.48$ for Innsbruck (Fig. 6b) and values near zero without clear seasonal variations for Hamburg (Fig. 6f). At both stations, the models with varying effects conditional on the wind direction have similarly distributed correlation parameters with a slightly larger range of predicted values for the model RAM-ADV (Fig. 6d, f) than for the setup RAM-DIR (Fig. 6c, g). The predicted correlation parameters are on average larger for Innsbruck than for Hamburg.

4.5 Predictive performance – models with correlation

Figure 7 shows the verification of bivariate wind speed predictions with an explicitly estimated correlation parameter for Innsbruck and Hamburg; the scores are aggregated over the forecast steps $+12$ h to $+72$ h on a $12$ hourly temporal resolution. As in Fig. 5, the predictive performance is validated in terms of the LS and the ES, based on the full predictive bivariate distributions. However, for comparing different predictive distributions with different correlation structures, the ES’ discriminatory ability is limited as it mainly focuses on the location part and hardly discriminates between different correlation structures (Pinson and Tastu, 2013). In Fig. 7, skill scores are shown for the different correlation models with the post-processed model RAM-0 as reference. The model RAM-EMP, employing the empirical correlation of the raw EPS, performs slightly worse than the reference model for both stations and both scores. This indicates that the raw dependence structure of the EPS has rather low skill. However, for all other models which explicitly model the correlation only little additional improvement in terms of the LS and the ES is found. At Innsbruck, the intercept-only model RAM-IC performs best in terms of the ES (Fig. 7c). Regarding the LS, minor benefits are present for the most flexible model setup RAM-ADV (Fig. 7a). For Hamburg, a similar picture is depicted in terms of the LS (Fig. 7b). For the ES (Fig. 7d), the model RAM-IC performs slightly worse than the reference model and the model setup RAM-ADV performs best.

To validate the calibration of the post-processed predictions, multivariate rank histograms (Gneiting, 2008) are exemplary shown for the model with no correlation RAM-0 and for the model with the most flexible regression splines in comparison to the raw EPS (Fig. 8). Although the latter is valid for the grid cell rather than for a single location, both model setups tested are clearly better calibrated than the highly under-dispersive raw ensemble. However, for Innsbruck the multivariate rank histograms of the post-processed forecasts are slightly over-dispersive (Fig. 8b, c) and for Hamburg slightly negatively skewed (Fig. 8e, f). The flexible model setup RAM-ADV shows no significant difference compared to the model with assumed zero correlation (RAM-0).

4.6 Evaluation for all sites

After the previous model comparison at two weather stations, Fig. 9 shows aggregated skill scores for groups of respective five stations classified as topographically plain, mountain foreland, and alpine valley sites (see Fig. 1). For the location or scale models, two comparisons are shown: The BLM-0 model is compared to the raw EPS as reference (Fig. 9a, d), and the more flexible rotation-allowing setup RAM-0 is compared to BLM-0 (Fig. 9b, e). For the correlation specification, the most flexible model RAM-ADV is compared to the correlation model RAM-0 employing a constant correlation of zero (Fig. 9c, f).

The post-processing employed by the simplest model BLM-0 already shows a distinct improvement over the raw EPS with largest values for alpine valley sites. In terms of the ES, the skill scores range between mean values of $10$% for the plain sites and $45$% for the alpine valley sites (Fig. 9d). A similar picture with an overall larger magnitude is shown for LS (Fig. 9a). In the comparison of the two different setups for the location or scale part (Fig. 9b, e), the more flexible setup is better regarding both scores for all station types; largest improvements are found for stations located in the foreland followed by stations within alpine valleys. The validation of the correlation models (Fig. 9c, f) shows that the flexible estimation of the correlation dependence structure is clearly superior only for station sites within an alpine valley.

5 Discussion and conclusion

In this study, we model the zonal and meridional wind components employing the bivariate Gaussian distribution in a GAM framework. In contrast to previous studies all distribution parameters, namely the means and variances for both wind components but also the correlation coefficient between them, are simultaneously estimated. The overall performance of the models is evaluated for three groups of station types classified as topographically plain, mountain foreland, and alpine valley sites.

Section 5.1 discusses the benefits of the rotation-allowing model setup RAM-0 over the baseline model BLM-0. In Sect. 5.2, the different correlation models are discussed regarding the potential reason for the low discrimination ability in terms of predictive performance. At the end, in Sect. 5.3, a recommendation is given which statistical model should be used in matter of simplicity and performance.

5.1 Rotation-allowing model setup

The rotation-allowing model (RAM-0) utilizes the ensemble information of both wind components for both location and both scale parameters in a GAM framework. This allows the statistical model to adjust for potential misspecification in the ensemble wind direction by a smooth rotation conditional on the day of the year and the forecasted wind direction. For stations in complex terrain, this may be particularly advantageous due to unresolved topographical features.

The estimated effects confirm a distinct wind rotation for the valley site (Innsbruck), while for the station in the plain (Hamburg) barely any adjustments of the forecasted wind direction is needed (see Fig. 4). In terms of predictive performance, the more flexible model RAM-0 outperforms the baseline model BLM-0 for almost all times and stations (see. Fig. 9b, e). However, the increase in predictive skill is similar for all three station types. This indicates that – even if no or only little rotation is needed – additional covariates usually yield a better adjustment of the distribution parameters and therefore an increased predictive skill. Furthermore, the results indicate that EPS wind forecasts in complex terrain are not solely tilted due to unresolved valley topographies, but show little skill on average. Thus, for alpine valley sites the rotation-allowing model mainly captures climatological properties conditional on the forecasted EPS wind direction. In accordance to this analysis, larger improvements can be found for stations located in the mountain foreland where the EPS has a higher information content and a certain rotation might be necessary.

5.2 Correlation specifications

Several different model specifications for the correlation parameter have been tested. Among others, a flexible setup employing wind direction and speed as potential covariates for the correlation parameter by non-linear smooth effects following the idea of Schuhen et al. (2012). The estimated correlation parameters seem to be reasonable, and show, on average, larger values for Innsbruck than for Hamburg (see Fig. 6 for forecast step $+12$ h). In terms of predictive skill, all models tested show only minor improvements compared to the models with zero correlation. The improvements are highest for stations located inside alpine valleys with a mean improvement of $1$% in terms of ES and $5$% in terms of LS skill scores (see. Fig. 9c, f).

As an illustration of the potential reasons for no more pronounced enhancements by an explicit estimation of the dependence structure, Fig. 10 shows exemplary forecasts for a station located within an alpine valley (Innsbruck; Fig. 10a) and in the plains (Hamburg; Fig. 10b). The figure shows the raw EPS members (gray points) plus the respective observations (black crosses), climatological estimates (gray dashed lines), and the corresponding post-processed bivariate distributions without (green lines) and with (blue lines) an explicitly estimated correlation parameter. For the valley station, the raw EPS has only little skill and the uncertainty of the post-processed bivariate distributions tends towards the climatological estimate. Although a distinct correlation is estimated by the RAM-ADV model, the variance is still in the same range as for the RAM-0 model. In contrast, for the station in the plain the uncertainty of the post-processed predictions is much smaller than the uncertainty of the climatological estimate due to a higher information content of the EPS. The estimated correlation is close to zero and the predictions of RAM-0 and RAM-ADV look almost identical with a similar elliptic shape as the raw EPS. This means that for locations where the ensemble provides only little information, the post-processed uncertainty is rather large and the statistical model tries to capture unexplained features by the correlation parameter. For stations where the predictive skill of the raw ensemble is already high, the statistical models get valuable information about the expected wind situation and are able to accurately specify the location and scale parameters. Thus, the correlation of the residuals becomes less important and typically smaller. This interpretation is supported by the probabilistic scores used in this study which show improvements in the RAM-ADV models mainly for alpine valley sites where the skill of the raw ensemble is rather low.

5.3 Proposed model specification

The study shows that the flexible rotation-allowing models bring significant performance benefits for stations located in complex terrain as well as for stations in the plain. Therefore, we propose using a similar setup employing both EPS wind components by a smooth rotation-allowing framework. For correlation, we have not found a clear distinction between the different correlation models tested for stations located in the plain and the foreland. For stations located within an alpine valley, minor improvements could be found. Despite these somewhat unexpected findings, this has clear advantages for operational usage: Estimating a single bivariate response distribution forcing the correlation dependence structure to zero is the same as post-processing each wind component separately in a univariate setup with marginal Gaussian response distributions. A univariate post-processing approach for each respective wind component simplifies the estimation process in terms of complexity of the required statistical models and reduces computational time with only little loss of predictive skill at least for the stations tested in this study.

Code availability

The bivariate Gaussian model estimation is performed in R 3.5.2 (R Core Team, 2018) based on the R package \pkgbamlss (Umlauf et al., 2018). The package provides a flexible toolbox for distribution regression models in a Bayesian framework. The computation of the ES is based on the R package \pkgscoringRules (Jordan et al., 2018).

Appendix A Model specification complements

A.1 Smooth functions

Generalized additive models (GAM, Hastie and Tibshirani, 1986) and generalized additive models for location, shape, and scale (GAMLSS, Rigby and Stasinopoulos, 2005) are generalizations of linear regression models which allow one to include potentially non-linear (and even multi-dimensional) effects in the linear predictors $\eta$. One frequently used form of non-linear terms are smooth functions, also often referred to as regression splines. These regression splines are directly linked to the model parameters as additive terms in $\eta$ and allow the statistical model to include non-linear transformations of a specific covariate, if needed. For further details a comprehensive introduction to GAMs is given in Wood (2017). An example for an additive predictor $\eta$ with a smooth function is:

[TABLE]

where $\alpha_{\bullet}$ are regression coefficients, $\text{x}_{\bullet}$ the covariates, and $\alpha_{1}\cdot\text{x}_{1}$ and $f_{1}(\text{x}_{2})$ are a linear and a potentially non-linear one-dimensional effect, respectively. Generally, $f_{1}$ can be any transformation of the covariate $\text{x}_{2}$ dependent on the specification of $f_{1}$. For periodic values smooth ‘cyclic’ splines are often applied, meaning that the function has the same value at its upper and lower boundaries. This is similar to applying a linear combination of (several) trigonometric functions, as, e.g., performed by Schuhen et al. (2012). In this study, we utilize ‘cubic’ smooth functions with cyclic constraints (Wood, 2017).

A.2 Time-adaptive training scheme

To account for seasonal variations of the intercept and the linear coefficients, seasonal cyclic splines are used. If the covariates provide sufficient information a time-adaptive training scheme might not be required. However, if the bias and/or the slope coefficient is not constant throughout the year or the covariate’s skill varies over the year these terms are mandatory to allow the statistical model to depict seasonal features.

We therefore fit one statistical model over a training data set including several years of data but allow the coefficient included in the linear predictor(s) $\eta$ to smoothly evolve over the year:

[TABLE]

As before, $\alpha_{\bullet}$ are the regression coefficients, x*∙* the covariates, and $f_{\bullet}(\text{doy})$ employ cyclic regression splines conditional on the day of the year (doy).

Appendix B Skill scores used for verification

To compare the different bivariate Gaussian models of this study, we employ skill scores. A skill score shows the improvements over a reference. For all measures with a perfect score of zero, the skill score simplifies to:

[TABLE]

In this study we use the logarithmic score (LS, Good 1952) and the energy score (ES, Gneiting and Raftery 2007) to validate the probabilistic performance of the bivariate Gaussian predictions of the statistical post-processing models. Both multivariate scores evaluate the full predictive distribution returned by the statistical models.

The calculation of the ES is based on the R package \pkgscoringRules (Jordan et al., 2018). For a predictive distribution $F$ on $\mathbb{R}^{d}$ given through $m$ discrete samples $\mathbf{X}_{1},\dots,\mathbf{X}_{m}$ from $F$ with $\mathbf{X}_{i}=(X_{i}^{(1)},\dots,X_{i}^{(d)})\in\mathbb{R}^{d},i=1,\dots,m$ the ES can be written as:

[TABLE]

where $\|\cdot\|$ denotes the Euclidean norm on $\mathbb{R}^{d}$, and $\mathbf{y}=(y^{(1)},\dots,y^{(d)})\in\mathbb{R}^{d}$ the multivariate observation.

The logarithmic score is defined as

[TABLE]

where the forecast $F$ on $\mathbb{R}^{d}$ admits a probability density function $f$, and $\mathbf{y}=(y^{(1)},\dots,y^{(d)})\in\mathbb{R}^{d}$ is, as before, the multivariate observation. For the bivariate Gaussian distribution, the density function $f$ is identical to the likelihood function $L$ in Eq. (1). To be able to calculate skill scores, we compute the LS by the cumulative distribution function for the neighborhood of the multivariate observation $\mathbf{y}$, defined by $\mathbf{y}\pm\epsilon$. Thus, the LS is limited upwards to zero (converges to zero the closer the forecasts are to perfect predictions), which in turn allows one to calculate skill scores as given by Eq. (B.1). In this study, we use $\epsilon=0.1$ ms*-1*; be aware that, depending on $\epsilon$, the quantitative values of the $LS$ might alter, but the qualitative results remain unchanged. The calculation of the bivariate cumulative distribution function is based on the R package \pkgmvtnorm (Genz and Bretz, 2009).

Acknowledgements

This project was funded by the Austrian Research Promotion Agency (FFG), grant no. 858537. We also thank the Zentralanstalt für Meteorologie und Geodynamik (ZAMG) for providing access to the data.