Space–time calibration of wind speed forecasts from regional climate models

Abstract

Numerical weather predictions (NWPs) are systematically subject to errors due to the deterministic solutions used by numerical models to simulate the atmosphere. Statistical postprocessing techniques are widely used nowadays for NWP calibration. However, time-varying bias is usually not accommodated by such models. The calibration performance is also sensitive to the temporal window used for training. This paper proposes space–time models that extend the main statistical postprocessing approaches to calibrate NWP model outputs. Trans-Gaussian random fields are considered to account for meteorological variables with asymmetric behavior. Data augmentation is used to account for the censoring of the response variable. The benefits of the proposed extensions are illustrated through the calibration of hourly 10-m height wind speed forecasts in Southeastern Brazil coming from the Eta model.

Appendices

Appendix A: Calculation of the roughness length

According to Hansen (1993), the vegetation present on a surface influences the aerodynamic roughness characteristics encountered by the mean wind flow over that surface, affecting both the mean wind speed and direction predicted by numerical models and various other atmospheric parameters. Therefore, the surface roughness length, $z_0$, defined as the height at which the wind speed equals zero, has an important role in the modeling of atmospheric processes.

The aerodynamic roughness length parameter, $z_0$, used as a covariate in the application, was estimated for each calibration site from key atmospheric variables, following some principles of the Monin-Obukhov similarity theory (Monin and Obukhov 1954). Particularly, assuming a logarithmic wind profile, the averaged wind speed $\overline{u_i}$ at height $z_i$ (10 m), given by:

$$\begin{aligned} \overline{u_i}= \frac{u_{*}}{k} \left[ \ln \left( \frac{z_i}{z_{0}}\right) -\Psi (\zeta _i) \right] \end{aligned}$$

was used to derive the roughness length $z_0$ as

$$\begin{aligned} z_{0}= z_{i}\exp \left( \frac{-\overline{u_{i}} k}{{u}_{*}} - \Psi (\zeta _{i})\right) , \end{aligned}$$

(12)

where k is the von Karman constant, $u_{*}$ is the friction velocity, $\Psi (\zeta _{i})$ is the stability correction function of the wind profile and $\zeta _i=\frac{z_i}{L}$ is the dimensionless stability parameter given by the height above ground, $z_i$, normalized by the Obukhov length, L.

Measurements of air temperature, air pressure, sensible heat flux, momentum flux, and wind stress can be used to derive the parameters $u_{*}$, L and $\Psi (\zeta _{i})$. In this application hourly air temperature and air pressure data were obtained from the available meteorological stations, while hourly reanalysis data from the CFSV2 model (Saha et al. 2014) were used as the source of heat and momentum fluxes, after interpolation from the CFSV2 regular grid to the calibration sites.

Hourly values of the roughness length, $z_0$, were first obtained from (12) for the 2 years of available data. Then, median values of $z_0$ by month and by hour within each month (288 values in total for each calibration site) were calculated, considering only those $z_0$ values that were estimated during neutral conditions of atmospheric stability (such condition is achieved when $|L|>500$).

Appendix B: Robustness of prior distributions

In this section, we verify the robustness of prior distributions assigned to the parameters of the proposed postprocessing models applied in Sect. 4. Based on these prior distributions in simulated experiments, the posterior and prior distributions associated with the parameter vectors from DGOP and STEMOS are presented in Fig. 10. For clarity of exposition, only the results of a subset of common parameter $\varvec{\theta }_0$ (vector with dimension r) are exhibited. Specifically, the intercept parameter at moment $t=0$ represented by $\theta _{0,0}$ is shown.

Even when using this specific set of prior distributions in a distinct application from that previously reported, Fig. 10 shows that the posterior distributions differ significantly from the prior distributions for all parameters. This is evidence that the gain of information is mostly provided by the data, ensuring that non-informative properties are preserved in general applications.

Fig. 10

Posterior and prior distributions associated with parameters of a DGOP and b STEMOS in simulated experiments. The vertical line represents the actual value assigned for parameters

Appendix C: Sensitivity analysis of discount factors

This sensitivity analysis was performed to choose the best combination of discount factors ($\tilde{\delta }_T$, $\tilde{\delta }_S$, $\tilde{\delta }_V)'$ for ($\delta _{T}$, $\delta _{S}$, $\delta _{V})'$. In particular, we also consider a subset of 59 weather stations and a random selection of 5 days per season from the temporal range of the available dataset. Following the same model comparison criteria described in Sect. 4, Table 4 reports the sensitivity analysis of discount factors through average MAE, RMSE, d, IS, DIC, and LPML values for 24-h wind speed forecasts at 10-m height. The setups of the dynamic models differ only in the discount factor values. The results indicate greater sensitivity for $\delta _T$, which significantly reduces the predictive performance when assuming lower values. In contrast, the variation of values for $\delta _{S}$ and $\delta _{V}$ only slightly alters the predictive performance. Thus, we determine ($\tilde{\delta }_T$, $\tilde{\delta }_S$, $\tilde{\delta }_V)' = (0.99, 0.95, 0.99)'$.

Table 4 Average MAE, RMSE, d, IS, DIC, and LPML of 24-h deterministic forecasts for wind speed at 10-m height (sensitivity analysis of discount factors)

Appendix D: Model comparison criteria

In this section, we briefly describe the model comparison criteria used to compare the prediction of the fitted models in Sect. 4. The first three criteria (RMSE, MAE and index of agreement) are appropriate to compare numerical predictions from the Eta model, which provides only deterministic estimates, with the proposed postprocessing models. The probabilistic forecasts are evaluated through IS, which takes into account the amplitude and coverage of the prediction intervals in a parsimonious way.

1.1 D.1  Mean absolute error and root-mean-square error

Standard measures of goodness of fit were also entertained in this study for comparison purposes. The root-mean-square error (RMSE) and the mean absolute error (MAE) are given by:

$$\begin{aligned} \text{ RMSE }=\frac{1}{nT}\sum _{i=1}^{n}\sum _{t=1}^{T}(y_{t}(s_i)-\hat{y}_{t}(s_i))^2 \ \ \ \text{ and } \ \ \ \text{ MAE }=\frac{1}{nT}\sum _{i=1}^{n}\sum _{t=1}^{T}|y_{t}(s_i)-\hat{y}_{t}(s_i)|, \end{aligned}$$

where $\hat{y}_{t}(s_i)$ is obtained through a Monte Carlo estimate of the posterior mean of the predictive distribution, that is, $E\left[ y_{t}(\mathbf{s}_i) \mid \mathbf{y}\right]$, across N draws. Smaller values of RMSE and MAE indicate better model fit.

1.2 D.2  Index of agreement

Willmott (1981) introduced a standard measure for assessing the quality of forecasts. The index of agreement (d) ranges between 0 (absence of agreement) and 1 (perfect agreement), and is given by:

$$\begin{aligned} d = 1 - \frac{\sum _{i=1}^n \sum _{t=1}^{T}\left( y_{t}(s_i) - \hat{y}_t(s_i)\right) ^2}{\sum _{i=1}^n \sum _{t=1}^{T}\left( |\hat{y}_t(s_i) - \bar{y}|+|y_t(s_i) - \bar{y}|\right) ^2}, \end{aligned}$$

where $\bar{y} = \frac{1}{n} \sum _{i=1}^n\sum _{t=1}^{T} y_t(s_i)$.

1.3 D.3  Interval score

The interval score (IS, Gneiting and Raftery 2007) is a scoring rule for interval predictions considering the symmetric prediction interval with level $(1-\alpha )\times$100%. The score is rewarded by accurate intervals and penalized when there is no coverage of the forecast. If actual values are contained in the prediction interval, this measure is reduced to the range amplitude. The average IS is given by:

$$\begin{aligned} \begin{aligned} \text {IS} =&\quad \frac{1}{nT}\sum _{i=1}^n\sum _{i=1}^T (\hat{u}_t(s_i)-\hat{l}_t(s_i)) \\&+\frac{2}{\alpha }(\hat{l}_t(s_i)-y_t(s_i))\mathbb {1}\left\{ y_t(s_i)<\hat{l}_t(s_i) \right\} \\&+\frac{2}{\alpha }(y_t(s_i)-\hat{u}_t(s_i))\mathbb {1}\left\{ y_t(s_i)>\hat{u}_t(s_i) \right\} \end{aligned} \end{aligned}$$

where $\hat{l}_t(s_i)$ and $\hat{u}_t(s_i)$ are, respectively, the lower bound obtained by the $\frac{\alpha }{2}$ quantile, and the upper bound, obtained by the $1-\frac{\alpha }{2}$ quantile based on the predictive distribution. The indicator function is represented by $\mathbb {1}$.

Smaller IS values indicate more efficient probabilistic forecasts.

1.4 D.4  DIC

Particularly useful in Bayesian model selection problems, the deviance information criterion (DIC, Spiegelhalter et al. 2002) is a hierarchical modeling generalization of the Akaike information criterion (AIC, Akaike 1974). This criterion is negatively oriented implying that models with smaller DIC should be preferred to models with larger DIC. The DIC is given by:

$$\begin{aligned} \text {DIC} = -2 \int \log \{p(\mathbf{y} |\Theta )\}p(\Theta |\mathbf{y} )d\Theta + \text {P}_\text {D}, \end{aligned}$$

where $\mathbf{y}$ is a vector of observed values and $\Theta$ is the parameter vector. Thus, $p(\mathbf{y} |\Theta )$ represents the likelihood function and $p(\Theta | \mathbf{y} )$, the posterior distribution. Defined as a Bayesian measure of model complexity, $\text {P}_\text {D}$ is given by:

$$\begin{aligned} \text {P}_\text {D} = 2 \log \{p(\mathbf{y} |\tilde{\Theta })\} -2 \int \log \{p(\mathbf{y} |\Theta )\}p(\Theta |\mathbf{y} )d\Theta , \end{aligned}$$

with $\tilde{\Theta }$ denoting the Bayes estimator.

A Monte Carlo approximation of DIC is given by:

$$\begin{aligned} \widehat{\text {DIC}} = -2 \frac{1}{M}\sum _{m = 1}^{M} \log \{p(\mathbf{y} |\Theta ^{(m)})\} + \hat{\text {P}}_\text {D}, \end{aligned}$$

where

$$\begin{aligned} \hat{\text {P}}_\text {D} = 2 \log \{p(\mathbf{y} |\tilde{\Theta })\} -2 \frac{1}{M}\sum _{m = 1}^{M}\log \{p(\mathbf{y} |\Theta ^{(m)})\}, \end{aligned}$$

with ${\Theta ^{(m)}}$ denoting the m-th MCMC sample of $\Theta$ from posterior distribution $p(\Theta |\mathbf{y} )$, $m = 1, \dots , M$.

1.5 D.5  LPML

A component of the Bayes factor, the logarithm of the pseudo marginal likelihood (LPML, Gelfand 1996) is given by:

$$\begin{aligned} \text {LPML} = \frac{1}{n} \sum _{i=1}^{n}\log \{\text {CPO}_i\}, \end{aligned}$$

where $\text {CPO}_i$ represents the conditional predictive ordinate (CPO) for location i and is given by:

$$\begin{aligned} \begin{aligned} \text {CPO}_i =&\quad p\big (\mathbf{y} (s_i)| \mathbf{y} (s_{-i})\big ) \\ =&\quad \bigg (\int \frac{1}{p(\mathbf{y} (s_i)|\Theta )}p(\Theta |\mathbf{y} )d\Theta \bigg )^{-1}, \end{aligned} \end{aligned}$$

with $\mathbf{y} (s_{i})=\big ({y}_{1}(s_{i}), \ldots , {y}_{T}(s_{i})\big )'$ and $\mathbf{y} (s_{-i})=\big (\mathbf{y} (s_{1}),\ldots ,\mathbf{y} (s_{i-1}), \mathbf{y} (s_{i+1}),\ldots , \mathbf{y} (s_{n})\big )'$, $i = 1,\ldots ,n.$ Note that the CPO$_i$ is based on the leave-one-out-cross-validation process and estimates the probability of $\mathbf{y} (s_{i})$ given the observation of $\mathbf{y} (s_{-i})$. In particular, a Monte Carlo approximation of CPO$_i$ is given by:

$$\begin{aligned} \widehat{\text {CPO}}_i = \bigg (\frac{1}{M}\sum _{m = 1}^{M} \frac{1}{p(\mathbf{y} (s_i)|\Theta ^{(m)})}\bigg )^{-1}, \end{aligned}$$

with ${\Theta ^{(m)}}$ denoting the m-th MCMC sample of $\Theta$ from posterior distribution $p(\Theta |\mathbf{y} )$, $m = 1, \ldots , M$.

Finally, the Monte Carlo approximation of LPML is given by:

$$\begin{aligned} \widehat{\text {LPML}} = \frac{1}{n} \sum _{i=1}^{n}\log \{\widehat{\text {CPO}}_i\}. \end{aligned}$$

The preferred model maximizes this criterion.

Appendix E: Supplementary results

See Figs. 11, 12, 13, 14.

Fig. 11

Boxplots of average MAE of 24-h deterministic forecasts for wind speed at 10-m height during the seasons. The dense horizontal line and the $\mathbf {\times }$ symbol in each boxplot represent respectively, the median and mean values

Fig. 12

Boxplots of average RMSE of 24-h deterministic forecasts for wind speed at 10-m height during the seasons. The dense horizontal line and the $\mathbf {\times }$ symbol in each boxplot represent respectively, the median and mean values

Fig. 13

Boxplots of average index of agreement of 24-h deterministic forecasts for wind speed at 10-m height during the seasons. Absolute agreement between predictions and actual values occurs when $d=1$. The dense horizontal line and the $\mathbf {\times }$ symbol in each boxplot represent respectively, the median and mean values

Fig. 14

Normal quantile-quantile plots for residuals of each fitted model during the seasons