A dilemma of the uniqueness of weather and climate model closure parameters

2.1. On state estimation methods

One of the most common techniques for state estimation is Kalman filtering (Kalman, 1960). The basic idea in all filters is to obtain the posterior state estimate $x_k^{\text{est}}$ and some error statistics, typically the error covariance matrix $C_k^{\text{est}}$, given the prior information $x_k^{\text{p}}$ and the observations y _k at time step k.

The posterior state estimate can be transported to the next filtering time step to become the new prior state by using the forecast model. For the state vector, this is usually straightforward, but the fundamental question in filtering is how to transport the uncertainty in time. In EKF, this problem is solved by calculating the prior covariance matrix $C_k^p$ at time k as

where Q _k is the model error covariance matrix and M _k is the tangent-linear forecast model. If the dimension of the state space is high, such as in large-scale atmospheric models, eq. (1) becomes impossible to compute without approximations.

One approximative approach is ensemble filtering, where the uncertainty is transported to the next filtering time step by transporting an ensemble of states using the forecast model. The difficulty (and variety) of the ensemble methods lies in how to update the prior ensemble in the analysis step. There are many alternatives, for instance the ensemble adjustment Kalman filter (Anderson, 2001, 2010), which belongs to the family of the so-called ensemble square root filters (Tippett et al., 2003). In this study, the EAKF was chosen because it is a deterministic ensemble based filter, that is, it does not include randomisation of the ensembles. The non-randomness is important, since in Section 2.3 we want to create a likelihood function based on the filtering outputs and the randomness in the likelihood would complicate the parameter inference (Dowd, 2011; Hakkarainen et al., 2012). Also, EAKF was chosen because it is conveniently available in the data assimilation research testbed (Anderson et al., 2009, DART) and widely applied (e.g. Schwartz et al., 2011; Torn and Davis, 2012).

In EKF, the modelling errors of the forecast model are taken into account via the error covariance matrix Q _k . In EAKF (and in ensemble square root filtering in general), the neglect of unknown model errors can lead to over-confidence in the prior state and hence ignorance of the observations in the analysis, which can lead to filter divergence. This is effectively circumvented by introducing covariance inflation factors (Anderson and Anderson, 1999). In addition, the sampling errors related to the limited ensemble size (undersampling), are mostly removed by localisation (Hamill et al., 2001), where observations affect only nearby grid points.

2.2. On parameter estimation

Parameter estimation techniques developed for the atmospheric models can be divided in two categories: online and offline methods. In online methods, it is assumed that model parameters are not static quantities, but can evolve adaptively, for example, as a part of a state estimation process (e.g. Annan et al., 2005) and are thus computationally relatively cheap to implement. Typically, a statistical interpretation is lacking in online methods, since the rate of change of the parameters is controlled by the user. In this article, these online methods are not discussed further. Instead we use offline methodology, where a predetermined (training) set of observations y is used for tuning the parameter values.

In Bayesian methodology (e.g. Gelman et al., 2003), the knowledge about the unknown parameters is inferred from the posterior distribution:

which is evaluated using the prior p(θ) and the likelihood p(y|θ). The prior contains the information that we have about the parameters based on the accumulated information from the past. The likelihood function specifies how plausible the observed data are given model parameter values. Therefore, defining a proper likelihood function is the central question in parameter estimation.

The parameter estimation problem is often written in a (cost) functional form with parameters as arguments. Different numerical methods can be used for optimising the cost function and point estimates, such as the Maximum a Posteriori (MAP) estimate, can be obtained. In addition, Monte Carlo methods like MCMC can be used for producing samples from the posterior distribution instead of a single point estimate, and studying the uncertainty in the parameter values. For large-scale applications, applying MCMC is challenging, since the methods can involve thousands of repeated likelihood evaluations, but ways to improve the efficiency of MCMC for high-CPU models have been recently developed (Solonen et al., 2012).

2.3. Likelihood via filtering methods

In this section, the filter likelihood concept is discussed in the context of the extended Kalman filter and the ensemble adjustment Kalman filter. The approach is known in the parameter estimation of stochastic models (e.g. Singer, 2002), but less studied in connection with deterministic, chaotic systems (Hakkarainen et al., 2012). For the sake of completeness, the derivation of the likelihood computation is briefly reviewed in the Appendix.

When EKF is considered, the filter likelihood formula, i.e., the likelihood for observing y _1:n given the parameters θ, can be written as

where $r_k=y_k-\mathrm{\mathscr{H}}(x_k^p)$ is the innovation vector and $C_k^y=H_k{C}_k^p{H}_k^T+R_k$ is its error covariance matrix at time k. Operator $\mid ·\mid$ denotes the matrix determinant. In filtering, $\mathrm{\mathscr{H}}$ is the observation operator that maps from the state space to the observation space and H is its linearization. The prior state $x_k^p$ is given by the model, the prior covariance matrix $C_k^p$ is given by the Kalman filtering formulas [eq. (1)] and R _k is the observation error covariance matrix.

Essentially, the filter likelihood formula [eq. (3)] depends on the summation over the cost functions J _k , where the normalising term $log\mid C_k^y\mid$ has to be accounted for, because it implicitly depends on the parameters θ. The dependence comes from eq. (1) via the tangent-linear forecast model. It should be noted that eq. (3) itself does not require that $x_k^p$ and its covariance matrix $C_k^p$ are calculated using EKF. In principle, any state estimation method that produces these prior estimates could be used.

In EAKF, as implemented in the DART system, the discrepancy between the prior ensemble and the observations are calculated in the observation space for every single observation sequentially. The update of the ensemble and the weights are also calculated sequentially in the observation space. Thus, it is natural to adopt this approach also in the filter likelihood calculations. The cost function part of the filter likelihood formula [eq. (3)] can be calculated sequentially as a sum

where y _l is an individual observation and $\overline{z_l^p}$ is the individual prior ensemble mean in observation space and ${\sigma }_{y_l}^2$ and ${\sigma }_{z_l^p}^2$ are their variances, respectively. The prior ensemble mean $\overline{z_l^p}$ depends implicitly on y _l−1, because of the sequential update. The summation goes through all the individual observations y _1:L,k at time k. The DART system does not provide the above computation, but including these in the DART code can be easily implemented.

2.4. Model and data

To explore the filter likelihood calculation, a modified version of the Lorenz-95 system is used (Lorenz, 1995; Wilks, 2005). The full system consists of two interrelated variables X and Y, and it is written as

where i=1,…, I and j=1,…, JI. That is, each of the ‘slow’ state variables X _i are forced by a sum of the additional fast variables Y _j . The fast variables have dynamics similar to the slow variables, but they are also coupled with the slow variables. In the model, cyclic boundary conditions, i.e., X _I+1=X ₁ and Y _JI+1=Y ₁, are used. In this paper, values I=40, J=8, F=F _Y =10, h=1 and c=b=10, adopted from Leutbecher (2010), are used.

The parametrized forecast model for the full system reads

where $g(X_i,\theta )={\theta }_0+{\theta }_1{X}_i$ is the parameterization in which the effect of the missing ‘sub-grid’ scale fast variables Y _j are modeled using only the single, ‘resolved’, local variable X _i .

The aim of the experiment is to tune the parameters θ=(θ ₀, θ ₁) given a synthetic set of observations y _1:n , computed by the full set of equations in (5) and (6) with additive noise, using the filter likelihood technique. In the experiment, one filtering step (‘24 h’) is 8 model integration steps (‘3 h’). Here the number of assimilation steps n=100 (‘days’), and hence during each run the model is integrated 100×8 times. The experiment is adopted from Hakkarainen et al. (2012) and a more comprehensive introduction is given there.

2.5. Experiment setup

In the experiment, the parameter values are sampled from the posterior distribution {consisting of filter likelihood [eq. (3)] and uninformative prior} using an adaptive Metropolis algorithm (Haario et al., 2001). The chain length is selected to be 3000 and the first 500 samples (the burn-in period) are discarded. The filter likelihood MCMC experimentation is synthesised with the following pseudo-algorithm, where the superscript m on θ denotes MCMC chain index:

Step 0:

 Initialise the MCMC run and select θ ¹.

Step 1:

 Propose a new candidate $\hat{\theta }$ for the parameters.

Step 2:

 Run the state estimation system using $\hat{\theta }$ and evaluate formula [eq. (3)] during the run.

Step 3:

 Accept $\hat{\theta }$ with probability $min\left(1,\frac{p(y_{1:n}\mid \widehat{\theta })p(\widehat{\theta })}{p(y_{1:n}\mid {\theta }^m)p({\theta }^m)}\right)$.

Step 4:

 If $\hat{\theta }$ is accepted set ${\theta }^{m+1}=\hat{\theta }$, else set ${\theta }^{m+1}={\theta }^m$.

Step 5:

 Set m=m+1 and go to step 1 until m is equal to the chain length.

In our implementation, we use an uninformative (flat) prior, which means that the prior is a constant, p(θ) α1. Then, the prior term cancels out in Step 3 of the above algorithm, since $p(\widehat{\theta })/p({\theta }^m)=1$.

For readers not familiar with MCMC, a short introduction is in order. The MCMC algorithm works by generating candidate values from a proposal distribution, which in our case is Gaussian centred at the current parameter value. The candidates are then either accepted or rejected according to a simple rule that contains the ratio of the posterior densities at the proposed point and at the current point, see Step 3 of the above algorithm. One can see that moves ‘upward’ (to a point with higher posterior density) are always accepted. However, also moves downward can be accepted, and the resulting random walk ends up exploring the posterior distribution without converging to a single point. The algorithm is a standard tool in statistical analysis and can be shown to yield samples from the correct target density. Here, we apply an adaptive version of the algorithm Haario et al. (2001), where the proposal distribution is automatically tuned to better match the posterior distribution. For more details on MCMC, see Robert and Casella (2005).

Using the filter likelihood [eq. (3)] for MCMC, we are able to estimate both the model parameters θ and the filter-specific tuning parameters. In the experiments, for the EKF version, the model error covariance matrix is parameterised as follows

where I is an identity matrix. The parameter λ is tuned together with the model parameters θ ₁ and θ ₁. In the EAKF version of the filter likelihood, the prior covariance inflation factor and the cut-off radius related to localisation (Gaspari and Cohn, 1999), defined in the DART setup, are tuned together with the model parameters θ.

As a validation metric for the parameters, we use the average root mean squared error (RMSE) defined as

where $x_k^{\text{true}}$ is the true state and $x_k^{\text{est}}$ is the posterior estimate. We evaluate RMSE on a fixed grid for θ ₁ and θ ₂ keeping the filter-specific parameters fixed.

2.6. Results

The experiments are started by a step, where all parameters (i.e. filter parameters and the model parameters) are estimated together with MCMC. In Fig. 1, pairwise MCMC samples using EKF (left panel) and EAKF (right panel) are shown. It can be seen, that the marginal distributions are rather Gaussian and all parameters are identified. Based on the mean values of the MCMC chains, values −4.8824, 1.4310 and 0.1889 were chosen for λ, the prior inflation factor and the cut-off radius, respectively.

In the second step, to further underline the difference in the model parameter estimates in the usual situation where the filter parameters are fixed, we calculate the conditional θ distributions using MCMC by fixing the filter parameters to the mean values found in the previous step. For the validation purposes, the RMS error [eq. (9)] is evaluated on a fixed grid to explore the parameter space for ‘all’ possible pairs of θ ₀ and θ ₁.

The results are illustrated in Fig. 2, where the 50% and 95% probability contours calculated from the EKF and EAKF MCMC chains are plotted using solid and dashed lines, respectively. The colours on the background indicate the average RMSE values [eq. (9)] for EKF (left panel) and EAKF (right panel). We observe that the closure parameter values that correspond to the lowest average RMSE values for EKF and EAKF are not the same. This implies that the optimal parameter values for EKF are not optimal for EAKF and vice versa. This shows that there is a dilemma of the uniqueness of the model closure parameters: an identical model formulation can have more than one optimal parameter value combination. This can be explained by the fact that the two data assimilation components have different systematic approximation errors, and these errors are compensated, to some extent, by the different parameter values in the forecast model. In this sense, the forecast model parameter values should be considered as filter-dependent quantities. The consequences of this result are discussed below in Section 3.

Figure 2 also depicts that the filter likelihood method captures the optimal parameter values correctly: the lowest posterior RMSE values are close to the obtained posterior distributions. In addition, it can be seen that the results produced by the EKF and EAKF versions of the filter likelihood technique are similar, but not identical. The posterior distribution obtained using EKF is closer to the parameter values that yield low RMS error for EKF. And vice versa: the parameters estimated using EAKF are closer to the optimal values for EAKF than the ones obtained by EKF. That is, applying the filter likelihood technique seems to provide rather good filter-specific parameter estimates. In this study, the RMSE values produced by EKF are lower than those produced by EAKF (see the RMSE values at the background of Fig. 2), which was expected. In addition, it can be noted, that the choice of data assimilation component has a much larger effect on the RMSE than the exact values of the model parameters.