## Abstract

Forecasts from numerical weather prediction models suffer from systematic and nonsystematic errors, which originate from various sources such as subgrid-scale variability affecting large scales. Statistical postprocessing techniques can partly remove such errors. Adaptive MOS techniques based on Kalman filters (here called AMOS), are used to sequentially postprocess the forecasts, by continuously updating the correction parameters as new ground observations become available. These techniques, originally proposed for deterministic forecasts, are valuable when long training datasets do not exist. Here, a new adaptive postprocessing technique for ensemble predictions (called AEMOS) is introduced. The proposed method implements a Kalman filtering approach that fully exploits the information content of the ensemble for updating the parameters of the postprocessing equation. A verification study for the region of Campania in southern Italy is performed. Two years (2014–15) of daily meteorological observations of 10-m wind speed and 2-m temperature from 18 ground-based automatic weather stations are used, comparing them with the corresponding COSMO-LEPS ensemble forecasts. It is shown that the proposed adaptive method outperforms the AMOS method, while it shows comparable results to the member-by-member batch postprocessing approach.

## 1. Introduction

Ensemble forecasts of meteorological variables provided by numerical weather prediction (NWP) systems are becoming increasingly popular in several decision-making problems concerning social and economic issues, such as water resource management, power generation, and protection from natural hazards. Ensemble forecasts are operationally employed for producing probabilistic forecasts of weather variables by accounting for the main sources of uncertainty in NWP model output (Buizza et al. 1999). The first main source of uncertainty is the error in the initial conditions, which is then amplified by the chaotic nature of the atmospheric processes. The second source of uncertainty is related to the errors introduced by the imperfection of the model formulation of the atmospheric processes, which also affects the forecasts at all lead times (Nicolis et al. 2009).

The ensemble spread is supposed to represent the uncertainty of the forecast. In the ideal case, the observed future meteorological state variable should fall within the predicted ensemble spread. However, all NWP models that are applied for operational weather forecasting tend to produce outputs that are systematically biased and underdispersive with respect to the local observations (e.g., those that are provided by ground weather stations). An NWP model can produce systematic or nonsystematic model errors in a generic point of the model grid, because of the limitations of the model ability to simulate the physical and dynamical properties of the system at the model resolution scale (Wilks 2011). Additional systematic or nonsystematic model errors are also noticed when the NWP model output is compared with point observations, due to the subgrid variability of the observations (Buzzi et al. 1994; Persson and Strauss 1995).

Reducing the impact of model errors is particularly important when weather forecasts are employed for site specific applications, such as irrigation management at farm level (Perera et al. 2014; Pelosi et al. 2016) and optimization of wind power farms (e.g., Cassola and Burlando 2012).

Statistical postprocessing techniques can be used to partly correct the impact of model errors in NWP models, and thereby improve the performance of the direct NWP model output (Vannitsem and Nicolis 2008; Vannitsem 2008).

Statistical postprocessing can be classified into two main categories:

Batch statistical methods, which are calibrated offline with a training set of numerical weather forecasts and ground observations. An overview of methods for deterministic and ensemble forecasts can be found in Wilks (2011) and in the recent comparison of member-by-member (MBM) and distribution-fitting approaches in Schefzik (2017).

Adaptive methods, that sequentially postprocess the forecasts by continuously updating the correction parameters as new ground observations become available.

An advantage of adaptive methods is that the updating procedure automatically accounts for NWP model changes (upgrades) and/or nonstationarity of the model prediction performance with respect to seasonal weather changes (Persson 1991). Adaptive methods are valuable for specific sites where ground weather stations have not been operating for a long period, and thus a dataset of ground observations is not available for a reliable training of a batch postprocessing method (Cheng and Steenburgh 2007; Delle Monache et al. 2006, 2011). Also, adaptive methods have the advantage of being computationally cheaper compared with more traditional batch estimation methods. This is because only the last available set of forecasts and observations is employed at a given time step for updating model correction parameters (Pinson 2012).

Several studies have implemented adaptive Kalman filtering procedures for correcting forecast errors of deterministic predictions, namely of temperature and wind speed (e.g., Persson 1991; Simonsen 1991; Homleid 1995; Galanis and Anadranistokis 2002; Roeger et al. 2003; Crochet 2004; Delle Monache et al. 2006; Galanis et al. 2006, 2009, 2016; Libonati et al. 2008; Louka et al. 2008; McCollor and Stull 2008; Cassola and Burlando 2012; Stathopoulos et al. 2013), entailing the sequential update of the parameters of regression equations employed for predicting the error or the corrected forecasts (predictand) as a function of the raw forecast variables used as predictors. These adaptive procedures are also referred to as adaptive model output statistics (AMOS), since they implement the same regression scheme adopted by MOS postprocessing batch procedures, which estimate the regression parameters by statistically matching the deterministic NWP raw forecasts against a long record of observations. One critical aspect emerging from the implementation of AMOS procedures is the identification of the two noise variances to be adopted in the state-space formulation employed within the Kalman filter: one for the state equation describing the random evolution of the regression parameters; and the other for the observation equation represented by the regression equation linking the predictand with the predictors.

According to the guidelines provided by the European Centre for Medium-Range Weather Forecasts (ECMWF 2015), in case of ensemble forecasts, AMOS technique should be applied to the control forecast or the ensemble mean, whereupon the updated regression equations should be applied to all members of the ensemble. The motivation is that if the filtering procedure would be applied to each ensemble member, the forecast anomalies represented by the ensemble spread would be erroneously treated as model errors, and thus the ensemble spread would be dampened after postprocessing.

In this paper, we present a Kalman filtering procedure for single-model ensemble forecasts (i.e., when the ensemble forecasts are produced by the same NWP model). A sequential correction of the impact of model errors of the ensemble members is obtained by exploiting the entire ensemble information content, instead of reducing the analysis to the ensemble mean or the control forecast. With the proposed method, one correction equation is retrieved and applied to all members; however, the parameters of the regression equations are retrieved by exploiting the second-order statistics of the raw forecast ensemble, while an explicit assessment of the noise variance of the observation equation is not required. The impact of the model errors from the ensemble NWP output are partly corrected without dampening the less predictable forecast anomalies represented by the ensemble spread. The procedure is formulated under the commonly held assumption that individual ensemble members represent possible solutions of the weather state for a given initial state (Sivillo et al. 1997). Each member bears a sample value of the impact of the model error conditioned to the forecasted weather state, beside the “random” variability, which can be considered as a sample value of the forecast anomalies represented by the ensemble spread.

In the following, we first provide an overview of the AMOS approach as it has been proposed in previous works. Then we illustrate the new Kalman filtering approach proposed for the ensemble forecasts (hereafter referred to as AEMOS), highlighting its advantage with respect to the traditional AMOS approach applied to the ensemble mean. The proposed AEMOS procedure is applied for postprocessing Consortium for Small-Scale Modeling Limited-Area Ensemble Prediction System (COSMO-LEPS) ensemble forecasts of 2-m temperature and 10-m wind speed in the Campania region (southern Italy). AEMOS performance is also evaluated in comparison with other postprocessing procedures used as benchmarks: AMOS applied to the ensemble mean as suggested by the ECMWF guidelines (ECMWF 2015) and a MBM batch postprocessing approach (Van Schaeybroeck and Vannitsem 2015).

## 2. Adaptive Kalman filtering for postprocessing

In this section, we illustrate the general structure of the adaptive Kalman filtering procedures (AMOS), as developed and applied for deterministic forecasts. For simplicity, we refer to the case that just one forecast variable is used as a predictor, as it has been usually applied. The procedure can be easily generalized to the case of multiple predictors.

We examine raw forecasts aggregated at time step , with lead times with *k* = 1, …, *r*, where is the forecast range.

Let be the raw forecast of a generic meteorological variable, issued at time for time *t* (i.e., with lead time ). The forecast error with respect to the observation is defined as follows:

Analogously to the MOS correction scheme, the forecast error is the predictand, which is assumed to be explained as a generic polynomial function of the forecast variable :

where is the unexplained noise and denotes the *j*th coefficient of the polynomial function of order *p*.

Instead of the forecast error, the actual corrected forecast can be used as predictand, as proposed both in adaptive (e.g., Crochet 2004) and batch (e.g., Van Schaeybroeck and Vannitsem 2015) postprocessing procedures. The two expressions are equivalent, since from a formal perspective the regression parameters with one approach are linear combinations of the regression parameters with the other approach. No significant differences are expected by using one approach or the other. In this section, we follow the structure originally implemented by Persson (1991) and that has been followed by most of the other studies implementing the AMOS approach.

The coefficients of the polynomial function are assumed to be time dependent, artificially evolving in time according to a random walk:

where represents a white Gaussian noise, affecting the dynamic evolution of the coefficients. The importance and usefulness of this simple Markovian model stem from the fact that the noise generated by the superposition of many small, independent, random effects is often Gaussian (Jazwinski 1970). It provides a common framework for introducing the stochastic forcing when using the Kalman filter, which ensures a constant increase of the state variance over the time unit (Evensen 2003). The noise component is independent and drawn from a zero-mean normal distribution with variance .

Equations (2) and (3) can be written in a matrix state-space form corresponding to the following system and observation equations, respectively:

In the equations above is the vector of the regression coefficients , denotes the parameter evolution noise that is assumed to be normally distributed with zero mean and covariance matrix , and is the random noise corrupting the predicted variables and is assumed to be normally distributed with zero mean and variance . Vector is the matrix representation of the explanatory variables in the polynomial observation equation [Eq. (2)]:

Given the regression coefficients estimated for the forecast issued at time valid for time , when the last analysis was performed, the prior estimate of the regression coefficients and their error covariance matrix of at time is defined as follows:

Given the raw forecasts available at time *t*, the corresponding forecast error is predicted as follows:

As an observation becomes available at the same time *t*, the observed error is computed as in Eq. (1). The innovation of the error with respect to the predicted error is defined as follows:

The current prior prediction is then combined with current observation information to provide an updated estimate (also known as posterior estimate) of the regression coefficients and of the corresponding error covariance . The posterior estimate is obtained by a weighted average of the prior prediction and the observation, with weights defined by a vector , known as the Kalman gain:

According to Eq. (11), a high gain implies that the filter places more weight on the last observation, while a low gain implies that the state updates are closer to the prior model predictions. The optimal weights are those providing updated state estimates with the minimum mean-square error. Under the Kalman filter theory (Kalman 1960), the optimal Kalman gain is estimated accounting for the relative certainty of the observation and the prior state prediction by means of the state and innovation covariances:

In Eq. (12), is the variance of the innovations and is estimated as follows:

Minimizing the mean-square error is equivalent to minimize the trace of the posterior covariance matrix , which is estimated as follows:

It can be shown that the Kalman filter algorithm outlined above is equivalent to an optimization algorithm that recursively determines the minimum state of a cost function. It can be shown that the cost function consists of a weighted prediction error and estimation error components given by (Wan and Nelson 2001):

Estimates of the initial state and error covariance are required before running the filter, although their impacts tend to rapidly disappear after the first few updates. The identification of the noise variances and is a more critical aspect, since they affect the filter performance. The ratio of the noise covariance to the observation variance affects the convergence rate and tracking performance. The larger this ratio, the more quickly older data are discarded. The covariance matrix is generally assumed to be diagonal (i.e., the noises of the regression coefficients are assumed to be mutually independent).

Homleid (1995) applied a Kalman filter procedure for correcting temperature forecasts and verified that automatic identification of the noise variances based on statistical procedures, such as the expectation maximization, resulted in a Kalman filter adapting too slowly to the new weather conditions. The best results could be obtained instead by manually tuning the noise covariances to achieve fast adaptation of the filter to new weather conditions.

This type of filter is generally not able to adapt to sudden changes in forecast deviation associated with fast changes in weather conditions, since the correction of future values is only based on previous differences between observations and forecasts and the model does not use any knowledge of the dependency of the forecast deviations on weather conditions (Homleid 1995). Persson (1991) suggested a way to improve the filter performances by introducing other prognostic variables as predictors, such as the temperature at 850 hPa for postprocessing 2-m temperature. However, the importance of different predictors is highly dependent on geographical position and it varies with time, with strong nonlinearity with respect to the predictand. Persson (1991) also found that the two noise covariances (or their ratio) should be manipulated through external interference, also automatically in accordance to weather changes identified in the NWP outputs. In this view, Crochet (2004) proposed some statistical adaptive procedures for estimating the noise variances, but the filter was found to be quite sensitive to instability emerging in the corrected forecasts. Galanis and Anadranistokis (2002) proposed to adaptively estimate the noise variances by simply computing them from the sample data of the last 7 days of predicted regression parameters and forecast errors. This approach was also implemented in more recent studies (e.g., Galanis et al. 2006; Cassola and Burlando 2012; Stathopoulos et al. 2013).

## 3. Adaptive Kalman filtering applied to ensemble forecasts

In this section, we illustrate a new adaptive postprocessing procedure for ensemble forecasts based on a Kalman filtering technique. We examine forecasts of an ensemble prediction system (EPS) aggregated at time step , with lead times with *k* = 1, …, *r*, where is the EPS forecast range.

Let be an ensemble of forecasts of a generic meteorological variable, issued at time for time *t* (i.e., with lead time ). As mentioned in the introduction, we assume that each ensemble member bears a sample value of the impact of model errors conditioned to the forecasted weather state. Thus, the forecast error can be computed by applying Eq. (1) to each *i*th member of the ensemble forecast with respect to the same observation . Analogously, the forecast error of the *i*th member of the ensemble forecast is assumed to be partially explained by a generic polynomial function of order *p* applied to the forecast variable [i.e., by applying Eq. (2) to each member]. The coefficients of the polynomial function are assumed to be time dependent, artificially evolving in time according to a random walk, as in Eq. (3), and thus the dynamic evolution of the regression parameters is expressed by Eq. (4). Instead, the observation equation is generalized with the following matrix form:

where is the vector of ensemble forecast errors , is the random noise vector corrupting the ensemble forecast errors , and is assumed to be normally distributed with zero mean and covariance matrix , which is assumed to be diagonal analogously to . Matrix is a Vandermonde matrix with dimensions *n* × (*p* + 1) and has the following structure:

The state-space structure outlined above is equivalent to that of a dynamic system with multiobservations or multisensors. This type of state-space formulation can be managed with the information fusion Kalman filtering (IFKF) theory, which yields the combination of information fusion theory and Kalman filtering theory (e.g., Sun 2004). IFKF has been widely applied in integrated navigation systems for maneuvering targets, for example airplanes, ships, automobiles, robots, etc. (Sun 2004).

When multiple observations are available for the same stochastic system, all measured sensor data can be processed into a single centralized filter, which yields optimal estimates with minimal information loss in linear cases. Centralized filters can be structured according to two different architectures, which are equivalent and optimum in the case of linear systems: a parallel and a sequential approach (Willner et al. 1976). The parallel approach considers all measurements together to retrieve the state update (Fig. 1). The sequential approach updates the state vector by exploiting the multiple measurements sequentially (Fig. 2). The two filters present the same prediction stage. While they are formally different in the update stage, they are mathematically equivalent in case of linear state-space systems.

In section 3a we present the structure of the prediction stage for a generic forecast lead time (*k* ≥ 1). In section 3b we illustrate the update stage of the parallel filter architecture for a generic forecast lead time (*k* ≥ 1). The same is illustrated for the sequential filter architecture in the appendix. In section 3c we suggest a generalization of the update stage for lead time (*l* > 1), by simultaneous analyses of the forecast issued for all lead times with *k =* 1, …, *l*. In section 3d we illustrate the strategy proposed for estimating the noise variances.

### a. Prediction stage

Given the regression coefficients estimated for the forecast issued at time valid for time , when the last analysis was performed, the prior estimates of the regression coefficients and error covariance matrix of are estimated according to Eqs. (7) and (8), respectively. Given the ensemble forecasts available at time *t*, an ensemble of predicted errors is computed as follows:

As an observation becomes available at time *t*, an ensemble of observed error is computed according to Eq. (1) and the ensemble innovation of the forecast error with respect to the predicted ensemble error is defined by applying Eq. (10) to each member.

In the following, we present the parallel filter for updating the regression coefficients and the error covariance matrix. The sequential one is discussed in the appendix.

### b. Updating with a parallel filter

With a parallel filter, all measurements in the ensemble are processed in parallel for updating and . The posterior estimate , representative of the regression coefficients valid for the forecasts issued at time for time *t* with a lead time , and the corresponding covariance are estimated according to the following update equations:

In Eq. (19), is the Kalman gain associated with the *i*th member and is computed as follows:

where is the noise relative to the *i*th member (i.e., the *i*th diagonal element of the noise variance matrix ).

The above algorithm entails the inversion of , which is exposed to errors when is close to be singular. The inversion of the covariance matrix can be avoided with a suboptimal solution (Sun 2004), as follows:

where the term at denominator is the innovation variance of the *i*th member.

### c. Simultaneous analyses of multiple forecast lead times

According to the algorithm presented above, the regression coefficients estimated for correcting the forecasts at are only based on the observed error at time of forecasts issued at time and, thus, do not account for the observed error at time for forecasts based on the weather state analysis at more recent times. This also implies that a set of regression coefficients is independently estimated for each forecast lead time, without accounting for the correlation of the errors observed at different lead times. Under the hypothesis of stationarity of the impact of the model error for different lead times, the regression coefficients estimated for the shortest lead time could be applied for greater lead times.

Here we explore the option of estimating the regression coefficients for correcting the impact of the model errors of the forecast with lead time , accounting for the errors observed at time with respect to all forecasts issued at times with . For a generic lead time , the same algorithm illustrated above is then implemented by expanding the number of observation equations, which simultaneously consider all members for :

where

having size (*n* × *l*, 1) and (*n* × *l*, *p*), respectively.

### d. Innovation variance and system noise identification

As commonly applied in Kalman filter applications (e.g., Entekhabi et al. 1994), the system noise variance is adaptively estimated to be proportional to the prior prediction of the state:

where denotes the identity matrix and *c* is a multiplicative coefficient that needs to be set.

With the procedure herein proposed, we avoid an explicit estimation of the system variance , by assessing directly the innovation variance from the second-order statistics of the ensemble. If we assume that the ensemble of innovations is homoscedastic, the innovation variance is equal for all members, and it is estimated by the ensemble variance of the innovations plus a term representative of the error variance of the ground measurements:

where *d* is the relative accuracy of the ground measurements. In case of simultaneous analyses of multiple forecast lead times, Eq. (28) is rewritten as follows:

In the approach outlined above, an explicit assessment of the noise variance is avoided and the Kalman gain is adaptively estimated in accordance with the weather changes identified in the ensemble NWP output.

### e. Forecast correction

After each update has been completed, a set of posterior estimates of the regression coefficients are available based on the analysis of the ground observations at time *t* and the forecasts issued at time and valid for time . These regression coefficients are employed for correcting the forecast issued at time *t* for the following lead time as follows:

The same correction equation is applied in case the regression coefficients are estimated by simultaneous analyses of multiple forecast lead times.

## 4. Assessment of forecast correction performances

### a. Study area and dataset

The study area is the region of Campania, about 14 000 km^{2} of land in southern Italy, between the Tyrrhenian Sea and the Southern Apennines (Fig. 3). Weather forecasting is a challenging task in this region, as in other coastal regions of the central Mediterranean basin, where weather patterns are strongly influenced by the complex topography close to the coastline (e.g., Buzzi et al. 1994; Pelosi and Furcolo 2015).

According to the Köppen–Geiger climate classification, most of the region is characterized by dry-summer subtropical climates, known as the “Mediterranean climate.” The coastal zone presents warm summers, while the adjacent inland zones are subject to hot summers. The eastern border zone of the region, close to the Apennines range, has a continental climate, as frequently occurs at higher elevations adjacent to areas with a Mediterranean climate (Peel et al. 2007).

In the region, the Regional Meteorological Service manages a reference weather monitoring network, made of 18 ground-based automatic weather stations (AWS) distributed across the region as indicated by the triangles in Fig. 3 and listed in Table 1. This AWS network complies with the European National Meteorological Services (EUMETNET) technical specifications: each station is equipped with redundant sensors to provide measurements with high accuracy and precision standards.

The weather forecasts considered in the analysis are the numerical weather prediction outputs given by COSMO-LEPS, which is a limited-area ensemble prediction system, implemented by the Hydro-Meteo-Climate Regional Service of Emilia-Romagna (ARPA–SIMC), within European COSMO. COSMO-LEPS provides daily probabilistic forecasts, consisting of 16 members, for up to 5 days, with a spatial resolution of 7 km (Montani et al. 2011; Pelosi et al. 2016). In this study, we make use of the 0000 UTC runs. The locations of the COSMO-LEPS grid points overlaid with the reference AWS sites are shown in Fig. 3 (circles). Postprocessing procedures were applied for correcting daily averages of 10-m wind speed and air temperature at 2 m recorded in years 2014–15.

### b. Experimental protocol

The AEMOS procedure (both with the sequential and parallel approach) was compared with:

the AMOS procedure applied to the ensemble mean;

a state-of-the-art batch postprocessing method based on the MBM approach developed by Van Schaeybroeck and Vannitsem (2015).

AMOS was applied to the ensemble mean as illustrated in section 2, and then the updated model correction parameters were applied to all members of the ensemble with Eq. (30). For the sake of comparison, the noise parameters were estimated with an approach consistent with the one used for the proposed AEMOS procedure. More details are provided in the following section.

In the MBM procedure, regression parameters are optimized based on the minimization of a performance score. Each ensemble member is corrected individually by a linear mapping, resulting in an overall shifting and scaling of the ensemble mean and spread. MBM retains rank correlations and thus takes correlations between nearby locations, lead times, and multiple meteorological variables into account (Van Schaeybroeck and Vannitsem 2015). The MBM method has four variational parameters (the so-called ensemble mean nudging and scaling and ensemble spread nudging and scaling), which are here determined through minimization of the empirical continuous ranked probability score (CRPS) (Gneiting et al. 2005) for the observations and the corrected-forecast members, with no specific distribution assumed (Van Schaeybroeck and Vannitsem 2015). For our verification study, we calibrate the MBM method on the first year of data (training set to optimize the regression parameters). The second year is used for comparing the proposed AEMOS with AMOS and MBM postprocessing procedures.

Both AEMOS and AMOS postprocessing procedures are implemented with a linear correction equation [i.e., with a polynomial of order *p* = 1 in Eq. (2) and following equations]. In similar previous studies, higher-order polynomials have also been applied to improve the filtering performances (e.g., Cassola and Burlando 2012; Stathopoulos et al. 2013). In this study, we focused on the simple linear correction equation, since the aim was to test the feasibility of this new methodology as compared with previous postprocessing procedures, rather than to look for the optimal correction for the specific test case.

### c. Parameter specification for adaptive procedures

Both AEMOS and AMOS procedures require the choice of the initial estimates of and error covariance . Since the verification is performed for the second year, the initial estimates and are not relevant for the performance statistics illustrated in the following paragraph. The initial values of the two regression parameters in vector are set equal to zero (i.e., we assumed that the initial error of the forecasting model in use is nonsystematic). The initial covariance is assumed to be diagonal with values set after a few trials and errors to ensure the stability of the algorithm (Medina et al. 2014). However, the covariance matrix loses its diagonal structure after the first update according to Eq. (14). AEMOS also requires the identification of the parameter *c* in Eq. (27) for estimating the noise variance and the parameter *d* in Eq. (28) defining the relative accuracy of the ground measurements. Parameter *c* is set after a few trials, to ensure the stability and convergence of the filter. Parameter *d* is set equal to a value normally assumed for environmental sensors (e.g., Chirico et al. 2014). Table 2 summarizes the parameters that are used in this experiment.

For the sake of comparison with AEMOS, the AMOS procedure is implemented by estimating the variance of the innovations as for AEMOS [i.e., by applying Eqs. (29) and (30) and thus by avoiding and explicit assessment of the observation noise variance]. The noise variance is instead set as in Eq. (27) with a parameter *c* equal to *n* × *l* times the values adopted for AEMOS, to ensure stability and convergence. The inflation of parameter *c* estimated for AEMOS is required to compensate for the fact that in AMOS the update is obtained by applying one Kalman gain to the mean ensemble, while in AEMOS the update is obtained by summing the gains applied to each of the *n* × *l* members.

### d. Statistical scores

The performances of the postprocessing procedures are assessed with the following verification tools: mean absolute error of the ensemble mean (MAE), CRPS, and binned spread–skill plot. Additionally, an uncertainty estimate is performed for each verification scores. We use a bootstrapping procedure based on 1000 random samples with replacements to assess the 90% confidence interval for each score.

The MAE is calculated as the mean of the absolute error of the ensemble mean over the verification period and it is used to evaluate the accuracy of the forecast in a deterministic fashion. The CRPS considers instead the entire distribution of the ensemble. It measures the integrated square difference between the cumulative distribution function (CDF) of the forecast variable, , and the corresponding CDF of the observed variable, 1{*x* ≥ *O*_{i}} as follows:

where is the step function that assumes the value 1, if the condition is met and 0 otherwise. Here, the CRPS is averaged across *n* pairs of forecasts and observations (with *n* = 365 days):

For computation of the CRPS, we use the expression of Hersbach (2000). The CRPS is equal to the MAE in case of deterministic forecasts.

Finally, to better examine the properties of the corrections and their statistical consistency (Anderson 1996), the binned spread–skill plot is employed. It compares the average ensemble variance to the mean squared error (MSE) of the ensemble mean over small class of intervals of the ensemble variance (e.g., Delle Monache et al. 2013). The edges of the bins are defined based on equal quantile intervals of the ensemble variance. Good performances, which here means statistical consistency, require that the average ensemble variance matches the MSE (i.e., the results are aligned along the one to one line).

### e. Results

Raw and postprocessed outputs of COSMO-LEPS are verified with the corresponding ground-based observations at the 18 sites shown in Fig. 3. The results refer to the second year (i.e., 2015) as specified above in the experimental protocol. Both postprocessing approaches (i.e., parallel and sequential filter) are evaluated. The two approaches give results that are undistinguishable and, therefore, only the results from the application of the parallel filter are shown.

#### 1) Wind speed at 10 m above ground

Figures 4 and 5 show results for lead times of 1 day and 5 days, respectively, pertaining to AWS No. 5, taken here as sample case. Figure 4a shows the forecast error of the raw and AEMOS corrected forecasts with 1-day lead time, during the second year (i.e., 2015). At this station, raw forecasts of wind speed underestimate the variable (i.e., exhibit a negative error according to Eq. (1). The raw ensemble spread is narrow, with an order of magnitude smaller than the actual forecast error. AEMOS postprocessing reduces the error both at a low and high speed values. In some instances, AEMOS overcorrects the forecasts, especially around mid-March when there is an evident spike in the wind speed error related to AEMOS postprocessed data. The reason is that Kalman filter attributes errors from antecedent states to current states and this may lead to an overcorrection of the states if there is a sudden variation of the forecast error, which results in spurious peaks and oscillations in the predicted state (Ridler et al. 2014; Li et al. 2013, 2014). The algorithm is more vulnerable to overestimation when the predicted data are far from the observations at a particular data assimilation time step. Figures 4b and 4c show the time evolution of the corresponding regression coefficients estimated with the AEMOS procedure. The intercept () coefficient varies between −1.06 and −1.43 m s^{−1}, and the angular coefficient () varies between −0.36 and −0.55. The time pattern of the two coefficients exhibits a positive correlation. Following Eq. (30), the forecast correction was obtained by multiplying the raw forecast by a factor equal to (1 − ). Thus, a negative value of the angular coefficient implies that the correction factor is larger than 1 and the spread of the corrected ensemble increases.

Figure 5 illustrates the analogous results but for a forecast lead time of 5 days. In this case, the spread of the ensemble has the same order of magnitude of the forecast error. The coefficient varies between −1.43 and −1.82 m s^{−1}, while varies between −0.13 and −0.16. The two retrieved coefficients appear to be negatively correlated. As for the case of a lead time of 1 day, a negative value of implies that the spread of the corrected ensemble increases.

Global performance statistics averaged among the examined 18 AWS are depicted in Figs. 6–8, as a function of the lead time. The two statistics MAE and CRPS of the raw forecast are almost independent of the lead time. Note also that the results are statistically significant for 12 stations over the 18 studied.

Figure 6 presents the MAE statistics along with their 90% confidence intervals. The average MAE of the raw forecast is around 1.0 m s^{−1}. The average MAE after AEMOS postprocessing is equal to 0.62 m s^{−1} at a lead time of 1 day, and it increases with lead time, reaching 0.80 m s^{−1} with a lead time of 5 days. The improvement of the performances is statistically significant for AEMOS and then, for MBM. With MBM, MAE has a trend like AEMOS for increasing lead times but it is smaller than AEMOS on average by 0.1 m s^{−1}, with statistically significant differences between methods after a lead time of 1 day. AMOS performance is more sensitive to the forecast lead time: AMOS provides MAE of 0.68 m s^{−1} at a lead time of 1 day, when there is no statistically significant difference with AEMOS; from a lead time of 3 days AMOS shows no reduction of MAE.

Figure 7 presents CRPS statistics along with their 90% confidence intervals. The raw forecasts display small sensitivity to the lead time, with values around 0.75 m s^{−1}. AEMOS reduces the CRPS from 0.51 m s^{−1} at lead time of 1 day to 0.61 m s^{−1} at lead time of 5 days, while CRPS with MBM ranges between 0.38 m s^{−1} and 0.49 m s^{−1}. Significant differences in CRPS are evident comparing the raw forecasts with all correction methods’ outputs for all lead times, except for AMOS after lead time of 3 days when no CRPS reduction is observed.

Figure 8 shows the binned spread–skill plot, which compares the MSE with the mean ensemble variance for the different correction methods, for a lead time of 1 day (Fig. 8a) and a lead time of 5 days (Fig. 8b). The raw ensemble is highly underdispersed (i.e., the ensemble variance is smaller than the MSE) at all ensemble variance values both for all lead times. The postprocessed ensembles are still underdispersed for all methods and lead times. However, the improvement of the statistical consistency of the ensemble after postprocessing is satisfying, for both MBM and AEMOS. Figure 8a shows that, for lead time of 1 day, AEMOS and AMOS slightly correct the underdispersion of the raw ensemble while MBM performs better. For longer lead times (Fig. 8b shows the case of lead time of 5 days), AMOS fails in correcting the underdispersion of the ensemble but the AEMOS performances are similar to the MBM ones.

#### 2) Air temperature at 2 m above ground

Figures 9 and 10 show results, respectively, for lead times of 1 and 5 days, for the observed 2-m air temperature recorded at AWS No. 1. Figure 9a shows the forecast error of the raw and AEMOS corrected forecasts with a 1-day lead time, during the second year of observations. At this station, raw forecasts of air temperature are overestimated [i.e., exhibit a positive forecast error according to Eq. (1)]. The raw ensemble spread is very narrow, with two orders of magnitude smaller than the actual forecast error. Figures 9b and 9c show the time evolution of the corresponding regression coefficients estimated with the AEMOS procedure. The largest forecast errors occur in January–February and September–October. The correction is obtained with the coefficient varying between 0.51° and 2.92°C, with the highest values in the month of the largest forecast error, while the impact of the angular coefficient is negligible, with values varying between −0.0045° and 0.011°C. A negligible angular coefficient also implies that the postprocessing has no impact on the spread.

Figure 10 illustrates the analogous results but for a forecast lead time of 5 days. In this case, the spread of the ensemble has an order of magnitude smaller than the forecast error. The coefficient varies between 0.010° and 0.048°C, while varies between 0.025 and 0.29. Most of the correction is achieved with angular coefficient and the highest correction is obtained in February–March.

Global performance statistics averaged among the 18 AWS are depicted in Figs. 11–13, as function of the forecast lead time. The three examined postprocessing procedures present similar performance statistics.

Figure 11 presents the MAE statistics along with their 90% confidence intervals. The MAE of the raw forecast increases from 1.35° to 1.53°C. AEMOS produces the largest MAE reduction, with values ranging from 0.89° to 1.27°C. MBM provides the worst result at a lead time of 1 day (0.98°C), but similar MAE values produced by AEMOS at a lead time of 5 days. AMOS produces intermediate results at a lead time of 1 day (0.92°C) and slightly worse results at a lead time of 5 days (1.31°C). However, there are no statistically significant differences among methods.

Figure 12 presents the CRPS statistics along with their 90% confidence intervals. For all the three methods of postprocessing, the reduction of the CRPS of the corrected forecasts compared with the raw forecasts is significant but, similarly to what is observed with the previous statistics, there are no statistically significant differences among methods. Raw forecasts present a CRPS ranging between 1.0° and 1.1°C. AEMOS and MBM present very close CRPS statistics, between 0.67° and 0.87°C. The CRPS for the forecasts after the correction with AMOS varies from 0.68° to 0.92°C, with the worst results for increasing lead times.

Figure 13 depicts the binned spread–skill plot by comparing the MSE with the ensemble variance for the different correction methods, for lead times of 1 (Fig. 13a) and 5 days (Fig. 13b). All the three methods provide better performances in terms of statistical consistency if compared with the raw forecasts, which means that the high underdispersion of the raw ensemble is significantly corrected at all range of variation of the ensemble variance by all three postprocessing methods. For a lead time of 1 day (Fig. 13a), AEMOS and AMOS show the same plots while MBM performs better. For a lead time of 5 days (Fig. 13b), the three methods correct the underdispersion of the raw forecasts in a similar way, except for high ensemble variance values where AEMOS and AMOS show a slight overdispersion.

## 5. Summary and conclusions

A new Kalman filtering approach for the adaptive correction of the forecast error for single model ensemble forecasts is proposed. Our method has two main innovative aspects:

the second-order statistics of the ensemble forecast are employed for assessing the innovation variance, which is in turn employed for computing the Kalman gain without an explicit estimation of the observation noises;

the parameters of the correction equation are estimated by combining the ensemble forecast error within a multiobservation equations according to the information fusion Kalman filtering theory.

The first point allows one to simplify one of the most critical aspects in the implementation of the traditional Kalman filtering adaptive techniques, which is the simultaneous estimation of the noise and observation covariances. The exploitation of the second-order statistics for the assessment of the innovation variance is applied to the traditional adaptive MOS (AMOS) procedures, which entail the sequential correction of the ensemble mean.

The second innovation is proposed as an improvement with respect to the traditional Kalman filtering approaches that are applied to the ensemble mean or the control forecast. The rationale is that the correction performance can be improved by exploiting the fact that each ensemble member bears a sample value of the impact of model errors conditioned to the forecasted weather state.

Both the proposed innovations have been implemented in a new adaptive procedure that we refer to as AEMOS.

Both AMOS and AEMOS methods were applied for postprocessing COSMO-LEPS forecasts of the 2-m temperature and the 10-m wind speed in the Campania region (southern Italy) for forecast lead times up to 5 days. The forecasts were verified at 18 sites across the region. The member-by-member approach was used as a benchmark for the evaluation of the proposed methods.

AMOS and AEMOS present comparable performances in the correction of the 2-m temperature over the entire range of lead times. The multiobservation equation approach implemented in AEMOS does not lead to significant improvements, since the coefficient of variation of the forecast ensemble of the 2-m temperature is small. AMOS and AEMOS prediction performances are also comparable with the MBM approach, with AEMOS performing slightly better than MBM.

In case of the wind speed, AEMOS performs slightly worse than MBM; however, the correction performances of the AEMOS procedure are significantly higher than AMOS at lead times greater than 1 day. At a greater lead time the coefficient of variation of the forecast ensemble is much higher, and thus the implementation of the multiobservation approach has a greater impact on the correction performances than in the case of the 2-m temperature. This result also suggests that adaptive procedures applied to ensemble mean require additional tuning of the system noise for increasing lead times, in order to keep stable correction performances.

As seen in the present work, the Kalman filtering approach can be competitive even for ensemble forecasts as compared with statistical methods that are now flourishing in the recent literature. The latter approaches necessitate the use of a sufficiently large number of past forecasts on which the statistical approach has to be trained. This is only possible provided the model has not been substantially changed for at least a year (provided we do not need postprocessing during this first year), or the user has access to reforecasts that require considerable computer resources. These resources are usually not available, except at specific world-class meteorological centers, making Kalman filtering a very appealing approach.

## Acknowledgments

This study was partly supported by the PIRAM Project, funded by the European Union and Campania region, within the rural development program 2007–13 (PSR 2007–2013 Campania-Mis. 124HC). COSMO-LEPS forecasts were provided by ARPA Emilia Romagna–Servizio Idro-Meteo-Clima. The Regional Meteorological Service provided ground-based data. MATLAB code of the AEMOS procedure is available upon request. The authors acknowledge the editor and two anonymous reviewers for useful comments that helped to improve the manuscript.

### APPENDIX

#### Sequential Architecture of the AEMOS Update Stage

In the sequential architecture, the state vector and state noise covariance matrix are updated applying the filter sequentially to each member as a sequence of independent blocks. The updates in state estimate and covariance of the previous member in the sequence are regarded as a priori estimations for the following member updates.

The update algorithm is initialized as follows:

Then for *i* = 1, …, *n*, the regression coefficients and corresponding error covariance are recursively updated:

The final updated regression coefficients and error covariance are

## REFERENCES

*Stochastic Processes and Filtering Theory*. 1st ed. Mathematics in Science and Engineering Series, Vol. 64, Academic Press, 376 pp.

*Lectures and Papers Presented at the WMO Training Workshop on the Interpretation of NWP Products in Terms of Local Weather Phenomena and Their Verification*, H. R. Glahn, Ed., PSMP Report Series, Vol. 34, World Meteorological Organization, 27–32 .

*Lectures and Papers Presented at the WMO Training Workshop on the Interpretation of NWP Products in Terms of Local Weather Phenomena and Their Verification*, H. R. Glahn, Ed., PSMP Report Series, Vol. 34, World Meteorological Organization, 33–37.

*Kalman Filtering and Neural Networks*, S. Haykin, Ed., John Wiley & Sons, 123–174.

*Proc. IEEE Conf. on Decision and Control*, Clearwater, FL, Institute of Electrical and Electronics Engineers,

## Footnotes

© 2017 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).