A- A+
Alt. Display

# Estimating forecast model bias in coupled global and limited-area models

## Abstract

Idealised perfect model experiments suggest that performing data assimilation on a ‘composite’ state vector, constructed from global and limited-area model states, can be beneficial to both model states. Here, an illustrative scheme is implemented to account for systematic forecast errors attributed to the imperfect model dynamics. Results from numerical experiments suggest that even simple bias correction schemes can correct forecast errors in composite states.

Keywords:
How to Cite: Kretschmer, M., Hunt, B.R. and Ott, E., 2015. Estimating forecast model bias in coupled global and limited-area models. Tellus A: Dynamic Meteorology and Oceanography, 67(1), p.28040. DOI: http://doi.org/10.3402/tellusa.v67.28040
Published on 01 Dec 2015
Accepted on 2 Nov 2015            Submitted on 31 Mar 2015

## 1. Introduction

Limited-area models are commonly used to generate regional high-resolution atmospheric forecasts. Correcting systematic forecast model errors, known as forecast model bias, can be complicated in such systems by the limited-area model's need for lateral boundary conditions. Lateral boundary conditions are typically provided by a concurrently running, nested limited-area or global model, each of which frequently exhibits model errors of its own. In operational forecasting, data assimilation is used to correct global and limited-area forecasts towards recent observations, but this is typically done separately for each model. Performing data assimilation simultaneously on both global and limited-area model states has recently been shown to be beneficial to both global and limited-area model states (Yoon et al., 2012; Kretschmer et al., 2015). However, simultaneous assimilation introduces additional coupling between model states, further allowing errors in each model to affect the other.

The composite state framework forms a state estimate from a combination of limited-area and global model forecasts (Kretschmer et al., 2015). The purpose of this paper is to illustrate the possibility of applying bias correction to the composite state method, and further, to illustrate that doing this may have substantial benefit. Correcting composite state forecast model bias accounts for the overall effect of model biases present in coupled global and limited-area forecast models. For this purpose, we will use a bias correction scheme initially presented in Baek et al. (2006).

Data assimilation techniques have been successfully applied to estimate the effect of, and to correct for, forecast model biases using variational and ensemble techniques (Dee and Da Silva, 1998; Carton et al., 2000; Dee, 2005; Keppenne et al., 2005; Li et al., 2009). Often, this is accomplished through state-vector augmentation techniques (Jazwinski, 1970; Cohn, 1997; Anderson and Anderson, 1999). The technique presented here aims to estimate the cumulative effect of model error on short-term model forecasts. Note that this bias correction method only modifies the analysis procedure, and not the forecast model equations. While many studies explore and account for the effects of forecast model bias on global models, many fewer have investigated how to account and correct for forecast model bias when limited-area models are involved [e.g. see the review by Meng and Zhang (2011)].

## 2. Background

We test our bias correction scheme in numerical experiments with biased forecast models, similar to Baek et al. (2006). For this discussion, and in our experiments, the ‘truth’ model dynamics are denoted as

(1 )
$\frac{d{\mathbf{x}}^{T}}{dt}=\text{M}\left({\mathbf{x}}^{T}\right),$

and the vector xT represents the truth model state vector.1 We model forecast bias as an additive model error tendency; details of the specific biased forecast models used in our experiments are presented below.

### 2.1. Forecast model bias correction

The present work adapts the composite state method to correct for forecast model errors. Baek et al. (2006) presented a series of ‘bias models’ to account for forecast bias resulting from imperfect forecast model dynamics in ensemble forecasting systems. Data assimilation is performed on a chosen bias model to adaptively estimate both the model state and forecast bias. This is accomplished by augmenting the state vector x with a vector of ‘bias corrections’.

Here, we consider Bias Model I proposed by Baek et al. (2006). In particular, this bias model assumes that, after spin-up, forecasts are initialised at time tm−1 to the approximate truth state at tm−1, ${\mathbf{x}}_{m-1}^{T}$. Forecast errors then arise because of deviations between the truth and forecast model dynamics. Using data assimilation, the unknown vector b is estimated so that the truth at time tm is estimated by Bias Model I

(2 )
${\mathbf{x}}_{m}^{T}\approx \text{F}\left({\mathbf{x}}_{m-1}^{T}\right)+\mathbf{b},$

where F is obtained by integrating the (biased) forecast model from time tm−1 to time tm. We note that the operator F may operate on either a model state estimate or the true state. The analysis is performed on the augmented state vector x, consisting of the corrected forecasted state vector, ${\mathbf{x}}_{m}^{b}=\text{F}\left({\mathbf{x}}_{m-1}^{a}\right)+{\mathbf{b}}_{m}^{b}$ and the current bias correction estimate ${\mathbf{b}}_{m}^{b}$. Both model state and bias correction estimates can then be updated using observations and the covariances between model variables and bias correction estimates to yield updated analysis estimates ${\mathbf{x}}_{m}^{a}$ and ${\mathbf{b}}_{m}^{a}$, which may then be forecast to the next analysis time, tm+1. The vector of bias corrections is evolved to time tm+1 via ${\mathbf{b}}_{m+1}^{b}={\text{G}}_{b}\left({\mathbf{b}}_{m}^{a}\right)$. The exact form of the bias correction time evolution operator Gb depends upon the properties of the detected forecast errors and may be empirically tuned as necessary.

### 2.2. Composite state forecasting

We account for forecast model bias in a coupled system of global and limited-area models. Data assimilation is performed on a background composite state vector, assumed to provide an optimal state estimate (Kretschmer et al., 2015). This composite state is formed by combining information from the global and limited-area model forecasts (which include coupling of global model provided boundary conditions), and is essentially created from the highest spatial resolution forecast information available at a given geographical location. Upon completion of the analysis procedure, the analysis composite state estimate ${\mathbf{x}}_{m-1}^{a}$ is evolved to time tm using

(3 )
${\mathbf{x}}_{m}^{b}=\text{L}\left({\mathbf{x}}_{m-1}^{a}\right)+{\mathbf{b}}_{m}^{b}.$

The operator L initialises the global and limited-area model states to ${\mathbf{x}}_{m-1}^{a}$, evolves these states from time tm−1 to tm using the (biased) global and limited-area forecast models and from these forecasts forms the composite state ${\mathbf{x}}_{m}^{b}$. For additional details regarding the construction of the composite state vector and the form of L used below, see Kretschmer et al. (2015).

## 3. Composite state bias correction

We consider coupled data assimilation performed between limited-area and global model states. When forecast models and data assimilation are coupled, errors in either of the forecast models can affect both limited-area and global model forecasts. Instead of accounting for model error separately in each of these forecast models, we consider errors in the composite state forecast, which represent cumulative effects of limited-area and global forecast model biases.

We approximate the effect of errors in the composite state forecast model, L in eq. (3), with an additive correction, using the composite state vector when applying eq. (2). Here ${\mathbf{x}}_{m}^{T}$ represents the truth at time m, and both ${\mathbf{b}}_{m}^{b}$ and ${\mathbf{x}}_{m}^{T}$ have the same spatial resolution as the composite state. As we shall see, correcting the forecast model biases of the composite state system allows the composite state, the best estimate of the true system state, to be dramatically more accurate in the presence of forecast model error than it would otherwise be.

## 4. Experimental details

### 4.1. The Lorenz models

Our numerical experiments use the one-dimensional chaotic models of Lorenz (2005). The dynamics of Lorenz's Model II is given at grid point n by

(4 )
$\frac{d{Z}_{n}}{dt}={\left[Z,Z\right]}_{K,n}-{Z}_{n}+F.$

The bracket term in eq. (4) represents spatially smoothed nonlinear advection, and F is a constant forcing parameter. Specifically, the bracket term of eq. (4) depends upon a smoothing parameter K and can be calculated for two quantities X and Y at each grid point n using

(5 )
${\left[X,Y\right]}_{K,n}=\frac{1}{{K}^{2}}\sum _{j=-J}^{J}\sum _{l=-J}^{J}\left({X}_{n-K+j-l}{Y}_{n+K+j}-{X}_{n-2K-l}{Y}_{n-K-j}\right),$

In eq. (5), J = K/2 if K is even and J=(K − 1)/2 if K is odd (Lorenz, 2005). Lorenz's Model III models the state variable Zn as the sum of long and short spatial scale components, with Zn=Xn+Yn. The long scale component Xn is defined as an average of Zn over nearby grid points, and the short scale component Yn is defined as the residual Yn=ZnXn. The dynamics of Model III are given by

(6 )
$\frac{d{Z}_{n}}{dt}={\left[X,X\right]}_{K,n}+{d}^{2}{\left[Y,Y\right]}_{1,n}+c{\left[Y,X\right]}_{1,n}-{X}_{n}-d{Y}_{n}+F.$

In the above model equations, d, c and K are parameters, and t represents time. See Lorenz (2005) for more details of these models.

For our numerical experiments, Model II describes the global model dynamics, and Model III describes both limited-area and truth model dynamics. The global and limited-area models are defined on different subsets of the N = 960 grid point lattice on which the truth model state is defined. The grid points of the truth lattice are indexed from 0, and all model state values are referenced by their index on the truth grid. The global model state is defined on every fourth point of the N=960 truth grid, and the limited-area model is defined on a subset of the truth model grid, over the grid point interval [240, 720].

The ‘true’ state is generated from a long free model run of Model III with parameters dT=10, cT=0.6, FT=15 and KT=32. In addition to the inherent bias of the global model due to its lower spatial resolution, we also introduce a bias by shifting the value of F used in the limited-area and global forecast models from the value used in the truth model, FT. The global model uses Kg=8 to match the length scale of K=32 used in Model III. Below, we use xm to represent the values of Z at time m at the composite state grid points, consisting of every grid point in [240, 720] and every fourth grid point elsewhere.

### 4.2. Data assimilation

In the numerical experiments, data assimilation is performed on a 32-member composite state ensemble using the Local Ensemble Transform Kalman Filter (Ott et al., 2004; Hunt et al., 2007), with a local analysis patch size of 81 grid points. The initial composite state ensemble is generated from samples of a free run of the global Model II. To prevent filter divergence, constant multiplicative covariance inflation is utilised (Anderson and Anderson, 1999). We achieve optimal results when composite state vector and bias correction parameter ensembles are inflated independently, by different amounts.

Observations are created by adding Gaussian white noise with mean 0 and standard deviation 1 to the true state at observation locations, and are directly compared to the value of the composite state at the observation location. The numerical experiments utilise a homogeneous observation network, with observations generated and assimilated every eight grid points.

### 4.3. Bias correction

We find that evolving the bias corrections in time using weak, numerical spatial diffusion allows for larger time steps and smaller ensemble sizes, in addition to speeding convergence, as previously reported in Baek et al. (2006). In our experiments, b is evolved in time via

(7 )

The parenthetical notation (n) in eq. (7) denotes the value of the given field at grid point n. We empirically tune the diffusion coefficient Db(n) over a range of magnitudes in order to minimise the root mean square (RMS) analysis error for each forecast model bias. We choose the diffusion coefficient Db(n) to be spatially constant where the composite state vector has constant spatial resolution.

In all experiments, the components of the forecast model bias are constant in time. When the global and limited-area forecast models are biased differently, their model biases are denoted with the subscripts g and r, respectively.

### 4.4. Error metric

We use the root mean square error (RMSE) of the analysis ensemble mean as a metric of analysis accuracy, and the root-mean square error of the forecast ensemble mean to measure forecast accuracy. Ensemble forecasts are initialised from the analysis ensemble. We define as the difference at location n and time m between the f-hour lead time forecast ensemble mean and the truth, respectively. The RMSE averaged over time is

(8 )
$RMSE\left(n,f\right)={\left\{\sum _{m=1}^{c}\left(\mathrm{\epsilon }\left(m,n,f\right){\right)}^{2}/c\right\}}^{1/2}.$

The RMSE of the analysis ensemble corresponds to f=0 in eq. (8). Unless noted, results shown use c=20000 in eq. (8).

## 5. Results and discussion

Our first numerical experiment uses Lorenz Model III [eq. (6)] for the limited-area model dynamics with parameter values dr=10, cr=0.6, Kr=32 and forcing Fr=14, so that the limited-area model dynamics may be denoted as the truth dynamics with an additive bias given by ΔFr=FrFT=−1. We use Lorenz Model II [eq. (4)] with Kg=8 and forcing Fg=13 to govern the global model dynamics. The values of dr, cr, Kr and Kg are constant across all experiments described here. The RMSE of the analysis composite state ensemble mean, calculated for constant forecast model bias, is shown in Fig. 1. It is seen that the composite state method is able to achieve significantly improved results when applying bias correction (blue curve) compared with the non-bias corrected case (red curve).

Fig. 1

RMSE of the analysis ensemble mean [eq. (8)] when using the composite state method (CSM). Bias correction (blue curve) significantly increases the analysis accuracy compared with the analysis without bias correction (red curve), and approaches perfect (unbiased) global forecast model results (black curve). The inset figure plots the spatial dependence of the time-averaged estimated bias correction b (gold curve) and the estimate provided by eq. (9) (green curve).

The bias-corrected composite state analysis corrects for both the spatially independent global model bias, ΔFg=FgFT= − 2, and the lower resolution and imperfect global model dynamics (Model II). We note that, without bias correction, the limited-area model has substantial error near the left boundary of its domain, and this feature is not seen with bias correction. We interpret this as evidence that the bias-corrected composite state improves lateral boundary conditions to the limited-area model, as the Lorenz models exhibit rightward information propagation (Yoon et al., 2010).

When modelling forecast model bias with additive model error tendencies, as is the focus of this manuscript, we expect the bias corrections to converge to

(9 )
$\mathbf{b}\approx -\mathrm{\Delta }F\mathrm{\Delta }t,$

where Δt denotes the assimilation window (Baek et al., 2006), which for our experiments has the value Δt=0.05. For the constant model biases in this experiment, eq. (9) predicts b should converge to b≈0.1, where the global model is defined and b≈0.05, where the limited-area model is defined. The bias correction parameters estimated in this experiment are shown in the inset of Fig. 1 (gold curve), and closely approximate the predicted asymptotic values of b shown in green.

In practice, it is likely that model errors will have some spatial dependence, and to investigate how well such biases may be corrected for in composite state forecasts, we allow the forecast model biases to vary according to $\mathrm{\Delta }{F}_{g}\left(n\right)=\mathrm{\Delta }{F}_{r}\left(n\right)=\text{sin}\left(2\pi \frac{n}{960}\right)$. Figure 2 shows the RMS analysis error of the composite state when there is spatially varying forecast model bias. As in Fig. 1, the bias controlled composite state analysis (blue curve) again outperforms the uncorrected composite state analysis (red curve). Bias correction accounts very well for the imposed model error in ΔF as well as the improper global forecast model dynamics, as evidenced by the relatively flat behaviour of the analysis RMSE curve.

Fig. 2

Same as Fig. 1, but with spatially dependent bias ΔF. Bias correction leads to decreased analysis RMSE compared to the composite state analysis without bias correction (blue versus red curves, respectively). The inset figure shows averaged bias correction b (gold curve) versus the value predicted by eq. (9) (green curve).

The increased accuracy of the bias-corrected analysis implies that the effect of the model bias is being accurately estimated. According to eq. (9), the bias corrections are expected to converge to $\mathbf{b}=-0.05\text{sin}\left(2\pi \frac{n}{960}\right)$, and this estimate is plotted in the inset of Fig. 2 as a green curve, along with the time-averaged value of b (gold curve). Their close agreement illustrates how the bias correction scheme is effective when forecast model bias is spatially varying.

Bias correction can potentially improve forecast results as well. We consider forecast lead times that are multiples of the assimilation window, Δt. Every Δt hours after initialisation at time tm, each forecast ensemble member xF is adjusted according to ${\text{X}}^{\text{F}}\to {\text{X}}^{\text{F}}+{\text{b}}_{m}^{a}$. Figure 3 shows RMSE of 2-d forecasts resulting from integrating an ensemble of global model states, initialised from the composite state analysis ensemble, while applying this methodology. Correcting the forecast model bias, even in this rudimentary fashion, clearly improves forecast accuracy, even out as far as 2-d lead times.

Fig. 3

RMSE of 2-d forecast ensemble mean, for spatially dependent bias as in Fig. 2. Blue and red curves compare forecasts with and without bias correction, respectively. The black curve shows ensemble forecast RMSE when forecasting with the truth model.

Current atmospheric forecast models contain myriad possible sources of error that can contaminate forecasts. When forecasting with coupled limited-area and global forecast models, these errors can affect the output forecasts in complicated and non-trivial ways, especially when performing coupled data assimilation on limited-area and global model states. Using a simple illustrative setup, we present here evidence that an appropriate method can account for cumulative forecast model errors in a coupled forecasting system, and that our results suggest that the composite state method with bias correction may be useful for producing forecasts with imperfect model dynamics.

## Notes

11In general, the dynamics represented by the operator $\text{M}$ may depend explicitly on time.

## 6. Acknowledgements

This work was supported by ONR grant N000141210785.

## References

1. Anderson J. L. , Anderson S. L . A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble assimilations and forecasts . Mon. Weather Rev . 1999 ; 127 : 2741 – 2758 .

2. Baek S.-J. , Hunt B. R. , Kalnay E. , Ott E. , Szunyogh I . Local ensemble Kalman filtering in the presence of model bias . Tellus A . 2006 ; 58 : 293 – 306 .

3. Carton J. A. , Chepurin G. , Cao X. , Giese B . A simple ocean data assimilation analysis of the global upper ocean 1950–95. Part I: methodology . J. Phys. Ocean . 2000 ; 30 : 294 – 309 .

4. Cohn S. E . An introduction to estimation theory . J. Meteorol. Soc. Jpn . 1997 ; 75 : 147 – 178 .

5. Dee D. P . Bias and data assimilation . Q. J. Roy. Meteorol. Soc . 2005 ; 131 : 3323 – 3343 .

6. Dee D. P. , Da Silva A. M . Data assimilation in the presence of forecast bias . Q. J. Roy. Meteorol. Soc . 1998 ; 124 : 269 – 295 .

7. Hunt B. R. , Kostelich E. J. , Szunyogh I . Efficient data assimilation for spatiotemporal chaos: a local ensemble transform Kalman filter . Physica D . 2007 ; 230 : 112 – 126 .

8. Jazwinski A. H . Stochastic Processes and Filtering Theory . 1970 ; Academic Press, San Diego, CA .

9. Keppenne C. , Rienecker M. , Kurkowski N. , Adamec D . Ensemble Kalman filter assimilation of temperature and altimeter data with bias correction and application to seasonal prediction . Non. Proc. Geophys . 2005 ; 12 : 491 – 503 .

10. Kretschmer M. , Hunt B. , Ott E. , Bishop C. H. , Rainwater S , co-authors . A composite state method for ensemble data assimilation with multiple limited-area models . Tellus A . 2015 ; 67 26495, doi: http://dx.doi.org/10.3402/tellusa.v67.26495 .

11. Li H. , Kalnay E. , Miyoshi T. , Danforth C. M . Accounting for model errors in ensemble data assimilation . Mon. Weather Rev . 2009 ; 137 : 3407 – 3419 .

12. Lorenz E. N . Designing chaotic models . J. Atmos. Sci . 2005 ; 62 : 1574 – 1587 .

13. Meng Z. , Zhang F . Limited-area ensemble-based data assimilation . Mon. Weather Rev . 2011 ; 139 : 2025 – 2045 .

14. Ott E. , Hunt B. R. , Szunyogh I. , Zimin A. V. , Kostelich E. J. , co-authors . A local ensemble Kalman filter for atmospheric data assimilation . Tellus A . 2004 ; 56 : 415 – 428 .

15. Yoon Y.-N. , Hunt B. R. , Ott E. , Szunyogh I . Simultaneous global and limited-area ensemble data assimilation using joint states . Tellus A . 2012 ; 64 : 18407 . doi: http://dx.doi.org/10.3402/tellusa.v67.18407 .

16. Yoon Y.-N. , Ott E. , Szunyogh I . On the propagation of information and the use of localisation in ensemble Kalman filtering . J. Atmos. Sci . 2010 ; 67 : 3823 – 3834 .