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

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.


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)].

Background
We test our bias correction scheme in numerical experiments with biased forecast models, similar to Baek et al. (2006 , allowing third parties to copy and redistribute the material in any medium or format and to remix, transform, and build upon the material for any purpose, even commercially, provided the original work is properly cited and states its license.

SERIES A DYNAMIC METEOROLOGY AND OCEANOGRAPHY
For this discussion, and in our experiments, the 'truth' model dynamics are denoted as (1) and the vector x T 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.

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 spinup, forecasts are initialised at time t m (1 to the approximate truth state at t m(1 , x T mÀ1 . 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 t m is estimated by Bias Model I where F is obtained by integrating the (biased) forecast model from time t m(1 to time t m . 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, m and the current bias correction estimate b b m . 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 x a m and b a m , which may then be forecast to the next analysis time, t m'1 . The vector of bias corrections is evolved to time t m' The exact form of the bias correction time evolution operator G b depends upon the properties of the detected forecast errors and may be empirically tuned as necessary.

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 x a mÀ1 is evolved to time t m using The operator L initialises the global and limited-area model states to x a mÀ1 , evolves these states from time t m (1 to t m using the (biased) global and limited-area forecast models and from these forecasts forms the composite state x b m . For additional details regarding the construction of the composite state vector and the form of L used below, see Kretschmer et al. (2015).

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 x T m represents the truth at time m, and both b b m and x T m 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.

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 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 In eq. (5), J 0 K/2 if K is even and J 0(K ( 1)/2 if K is odd (Lorenz, 2005). Lorenz's Model III models the state variable Z n as the sum of long and short spatial scale components, with Z n 0X n 'Y n . The long scale component X n is defined as an average of Z n over nearby grid points, and the short scale component Y n is defined as the residual Y n 0Z n (X n . The dynamics of Model III are given by dZ n dt ¼½X ; X K;n þ d 2 ½Y ; Y 1;n þ c½Y ; X 1;n À X n À dY n þ F : (6) In the above model equations, d, c and K are parameters, and t represents time. See Lorenz (2005) for more details of these models. The 'true' state is generated from a long free model run of Model III with parameters d T 010, c T 00.6, F T 015 and K T 032. 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, F T . The global model uses K g 08 to match the length scale of K032 used in Model III. Below, we use x m to represent the values of Z at time m at the composite state

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 observa-tion network, with observations generated and assimilated every eight grid points.

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). (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 D b (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 D b (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.

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 Eðm; n; f Þ ¼ x f m ðnÞ À x T m ðnÞ 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 The RMSE of the analysis ensemble corresponds to f00 in eq. (8). Unless noted, results shown use c 020000 in eq. (8).

Results and discussion
Our first numerical experiment uses Lorenz Model III [eq. (6)] for the limited-area model dynamics with parameter values d r 010, c r 00.6, K r 032 and forcing F r 014, so that the limited-area model dynamics may be denoted as the truth dynamics with an additive bias given by DF r 0F r (F T 0(1. We use Lorenz Model II [eq. (4)] with K g 08 and forcing F g 013 to govern the global model dynamics. The values of d r , c r , K r and K g 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). The bias-corrected composite state analysis corrects for both the spatially independent global model bias, DF g 0 F g (F T 0 (2, and the lower resolution and imperfect global model dynamics (Model II). We note that, without bias correction, the limited-area model has substantial 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. 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 b % ÀDF Dt; (9) where Dt denotes the assimilation window (Baek et al., 2006), which for our experiments has the value Dt 00.05. 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 DF g ðnÞ ¼ DF r ðnÞ ¼ sinð2p n 960 Þ. 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 DF as well as the improper global forecast model dynamics, as evidenced by the relatively flat behaviour of the analysis RMSE 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 b ¼ (0:05 sinð2p n 960 Þ, 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, Dt. Every Dt hours after initialisation at time t m , each forecast ensemble member x F is adjusted according to x F ! x F þ b a m . 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.
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.

Acknowledgements
This work was supported by ONR grant N000141210785.