Improving the handling of model bias in data assimilation

Patrick Laloyaux, Massimo Bonavita


Errors in numerical weather prediction arise from two main sources: incorrect initial conditions and deficiencies in the dynamics and the physical parametrizations of the forecast model. To correct initial errors, four-dimensional variational data assimilation (4D-Var) adjusts the initial state of the atmosphere to find the model trajectory that best fits the most recent meteorological observations. The new initial state from which forecasts start is called the analysis. In the standard 4D-Var formulation, known as strong-constraint 4D-Var (SC-4DVar), the model is assumed to be perfect and any systematic model errors (biases) which gradually accumulate in the short-range forecasts used in the data assimilation system (the first guess) are not taken into account. A version of 4D-Var which relaxes the assumption that the model is perfect, known as weak-constraint 4D-Var (WC-4DVar), has been run at ECMWF for some years, but without significant positive impact on analysis accuracy and forecast scores (Trémolet & Fisher, 2010). This has motivated the work reported in this article, which aimed to revise the operational WC-4DVar configuration to better mitigate the effect of model biases during the assimilation step.

Tests of the new WC-4DVar configuration have shown that analysis and first-guess temperature biases in the stratosphere are reduced by up to 50%. In view of these results, the revised WC-4DVar will be implemented in the next upgrade of the Integrated Forecasting System (IFS Cycle 47r1) planned for later this year. The initial implementation is restricted to the stratosphere, but work is under way to extend WC-4DVar to the troposphere using deep learning methods.

Identifying the problem

A standard approach to diagnosing the model error that develops during the data assimilation cycle is to compare the first-guess trajectory with accurate, unbiased observations. Figure 1 shows a time series of the difference between temperature retrievals from radio occultation satellite measurements and their model equivalents computed with the ECMWF operational first-guess trajectories. In the atmospheric layer between approximately 100 hPa and 10 hPa, the model is systematically colder than the observations (by about 0.6°C between 70 hPa and 10 hPa in 2019), while above 7 hPa the model is systematically warmer (up to about 0.6°C in 2019). This temperature bias of the model evolves over time due to changes in operational model cycles and partly due to seasonal variability. In March 2016, the bias increased when the horizontal resolution in high-resolution forecasts was increased from TL1279 (triangular–linear grid, corresponding to a grid spacing of about 16 km) to TCo1279 (triangular–cubic- octahedral grid, corresponding to a grid spacing of about 9 km) (IFS Cycle 41r2). It increased again when the radiative model was revised in the July 2017 upgrade (IFS Cycle 43r3).


FIGURE 1 Successive model changes have led to large temperature biases in parts of the stratosphere in the short-range forecasts used in ECMWF’s data assimilation system. The chart shows how such biases have evolved since 2014 as measured against temperature retrievals from radio occultation satellite data.


The radio occultation measurements provide near- homogeneous coverage that can be used to study the spatial structure of the model bias. Figure 2 shows the difference between temperature retrievals and first- guess trajectories from ECMWF operations between 70 hPa and 100 hPa from 1 October 2018 to 1 February 2019, when IFS Cycle 45r1 was operational. The model bias is predominantly cold with large-scale spatial variations. Colder biases are observed over areas of strong convection (e.g. Indonesia and Southern America). This could be linked to an insufficient representation of the effects of subgrid-scale gravity wave activity, which leads to missing momentum transfer from the troposphere to the stratosphere.


FIGURE 2 Difference between radio occultation temperature retrievals and first-guess temperatures from ECMWF operations between 70 hPa and 100 hPa over the period 31 August 2018 to 31 January 2019.


One of the main challenges in data assimilation is to attribute errors apparent from discrepancies between observations and the model to their proper sources. The diagnostic presented in Figure 2 suggests that model error which develops during the 4D-Var assimilation window contains identifiable large-scale structures, which are distinct from the small-scale errors in the first-guess initial conditions (the background) diagnosed by the Ensemble of Data Assimilations (EDA) system (Laloyaux et al., 2020a). Enforcing such scale separation in the errors that 4D-Var tries to correct opens a new perspective in the quest to disentangle model and background errors.

Developing the solution

Traditionally, 4D-Var accepts the model as a strong constraint without taking into account possible biases (Box A). This strong-constraint formulation of 4D-Var is designed to take into account random, zero-mean errors from the model forecast and the observations. However, significant systematic errors are generated by the forecast models used in global numerical weather prediction. The IFS is no exception as the forecasts it produces have large stratospheric temperature biases.


Strong-constraint 4D-Var

Strong-constraint 4D-Var compares the trajectory of a short-range forecast with observations. The assimilation system iteratively adjusts the initial conditions of the short-range forecast (the ‘background’, xb ) to compute a new analysis xa that achieves a better compromise between model forecasts and observational data within the assimilation window. In the strong-constraint formulation, the model is assumed to be perfect and the only way to adjust its trajectory is by changing the initial conditions at the beginning of the assimilation window.

To deal with this type of error, a modification of the standard 4D-Var algorithm, called weak-constraint 4D-Var, was proposed and was implemented in operations at ECMWF in 2009. In the original implementation, the model error covariance matrix (which describes the error statistics of model bias) was built using differences between 12-hour ECMWF ensemble forecasts (ENS). The forecasts were initialised from the same initial states but were run with model error perturbations based on the Stochastically Perturbed Parametrization Tendency scheme (SPPT) and the Stochastic Kinetic Energy Backscatter scheme (SKEB). These model error parametrizations produce perturbations whose energy is higher at synoptic and sub-synoptic scales than at larger scales (Shutts et al., 2011). Consequently, the resulting model error covariance matrix contains length scales similar to the ones that dominate the background error covariance matrix. This situation made it impossible for weak-constraint 4D-Var to attribute model and background error to the correct source. For this reason, the original version of WC-4DVar was only implemented in the upper stratosphere, as a full application of the algorithm would have led to an unacceptable degradation of forecast scores. Thus, the impact of weak-constraint 4D-Var on ECMWF analyses and forecasts has always been small.

To overcome this problem, a new estimate of the model error covariance matrix has been computed from a climatology of the model error vectors estimated by the current weak-constraint 4D-Var (Laloyaux et al., 2020a). The correlation length scale from this new model error covariance matrix is more consistent with the patterns observed in Figure 2. However, some long-range spurious correlations have been corrected through the application of horizontal and vertical localisation. Specifically, the horizontal correlations are tapered down linearly to zero in the 2,500–5,000 km range.

A similar approach is applied to vertical correlations by tapering them to zero in model space between 5 and 10 model levels. In this new implementation of weak-constraint 4D-Var, first-guess trajectory adjustments deal with large-scale temperature errors which vary on timescales longer than the assimilation window (i.e. longer than 12 hours), while state increments correct smaller-scale and more transient background errors (see Box B). This scale separation between the two types of errors is a necessary condition for weak-constraint 4D-Var to attribute model and background errors to their correct sources. This has been corroborated through theoretical considerations and experiments with a simplified quasi-geostrophic model (Laloyaux et al., 2020b). It has been shown that the forcing formulation of weak-constraint 4D-Var illustrated in Box B performs well if background and model errors vary over different spatial scales and if this is reflected in their respective error covariances.

A better analysis

We have carried out experiments to compare weak- constraint 4D-Var based on the new model error covariance matrix to the current weak-constraint 4D-Var configuration used in operations and to the strong- constraint 4D-Var configuration, in which no model error forcing is applied. All experiments were run at the operational high-resolution forecast (HRES) resolution (TCo1279) with 137 vertical levels and a model top pressure of 0.01 hPa. They were run from 1 October 2018 to 1 February 2019 using the standard 12-hour assimilation window and the background error covariance matrix estimates from the then operational 25-member EDA. Running these high-resolution experiments is computationally costly and the data throughput relatively slow. However, it is necessary to assess the impact of the new model error covariance matrix on HRES forecasts since the amplitude of the stratospheric bias is significantly bigger at increased resolution. All experiments were initialised on 1 October 2018 using the ECMWF operational analysis valid at that time. The model error term was cold started by setting it to zero in the first assimilation cycle.


Weak-constraint 4D-Var

Weak-constraint 4D-Var introduces a forcing term η into the model to correct for the model bias which builds up in the model trajectory. Initially, η can be set to zero. In each data assimilation cycle, weak-constraint 4D-Var then optimises the fit of the short-range forecast trajectory to observations in two ways: by suitably changing the initial conditions (as in strong-constraint 4D-Var), and by suitably adjusting the forcing term η, which is applied at every model time step. In the case illustrated here, for a single parameter x (e.g. temperature), the forcing term η cools the trajectory at every time step to correct for the temperature warm bias in the model. In this formulation of 4D-Var, a model error covariance matrix which contains the error statistics of the model bias (i.e. standard deviations and spatial correlations) needs to be calculated offline. This matrix enters into the calculations through which the data assimilation system determines the optimal combination of initial state and forcing term adjustments.

Figure 3a shows vertical profiles of the mean analysis departures (observations minus analysis, O–A) and background departures (observations minus background, O–B) with respect to radiosonde data for strong-constraint 4D-Var (SC-4DVAR), the currently operational weak-constraint 4D-Var (WC-4DVAR-OPER) and the new weak-constraint 4D-Var (WC-4DVAR-NEW). In the weak-constraint experiments, the first-guess trajectory used to compute the background departures is corrected using the model error estimate from the previous assimilation. A smaller mean error in the background departures shows that the assimilation system is able to estimate and correct a significant amount of bias in the model. The weak-constraint 4D-Var currently used in operations corrects only a small fraction of the model bias over 40 hPa, while the new model error covariance matrix better corrects the diagnosed cold and warm biases of the model, reducing the mean error by up to 50%. Since radiosonde observations are available only up to 5 hPa, Figure 3b shows analysis and background departures with respect to radio occultation bending angles to assess the upper stratosphere. The results show that bias in the upper stratosphere located between 30 km and 45 km (11 hPa to 1.5 hPa) is also significantly reduced in the new system.


FIGURE 3 Vertical profiles of analysis departures (O–A) and background departures (O–B) with respect to (a) radiosonde observations and (b) radio occultation observations (GPS-RO) for strong-constraint 4D-Var (SC-4DVAR), the current operational weak-constraint 4D-Var (WC-4DVAR-OPER) and the new weak-constraint 4D-Var (WC-4DVAR-NEW). The statistics were computed for the period from 1 October 2018 to 1 April 2019.


Reducing biases in forecasts

In the experiments described thus far, the model error estimated by weak-constraint 4D-Var was only used to correct model integrations computed inside the assimilation cycle. This enables the assimilation system to cycle forward in time what it has learned about the model bias and produces a better first-guess trajectory.

Plotting the difference between radio occultation temperature retrievals and the first-guess trajectory over different periods reveals that the spatial structure of the model error varies little over time (not shown). This means that the same model correction could be applied to control the bias in longer-range forecasts. Building on this idea, experiments have been carried out at low resolution (TCo399, corresponding to a grid spacing of about 25 km) and over an earlier period (20 May 2018 to 6 September 2018) to test the impact of extending the model trajectory adjustments applied in the new weak-constraint 4D-Var to 10-day forecasts. Figure 4 shows the mean error of the 10-day forecast initialised by the current weak-constraint 4D-Var, the new weak-constraint 4D-Var and the new weak-constraint 4D-Var with model trajectory adjustments in the forecast for temperature at 50 hPa in the northern hemisphere. Scores were aggregated over the period 1 June 2018 to 15 September 2018 and radiosonde observations were used for verification. The impact of using the new weak-constraint 4D-Var configuration can be seen in the analysis (lead time 0), where the mean error is decreased from -0.55°C to -0.39°C. Figure 4 shows that the improvement is maintained throughout the forecast. Nevertheless, the mean error increases significantly over time as the model converges towards its own cold- biased climate. By contrast, when the model trajectory adjustment estimated in weak-constraint 4D-Var is applied as a constant forcing term throughout the 10-day model integration, the mean error in the forecast barely increases over time. This shows that weak-constraint 4D-Var potentially has the capability to significantly reduce biases not only in the analysis but also in forecasts.

Meanwhile other steps to reduce forecast biases in the stratosphere are being taken. One of these is quintic vertical interpolation (see separate article in this Newsletter), which will be implemented in IFS Cycle 47r1. Another is an increase in the number of vertical levels in ensemble forecasts from 91 to 137, which is planned to be implemented once ECMWF’s next high-performance computing facility has become operational in Bologna.


FIGURE 4 Mean error of the 10-day forecast initialised by the current weak-constraint 4D-Var, the new weak-constraint 4D-Var and the new weak-constraint 4D-Var with bias correction in the forecast for temperature at 50 hPa in the northern hemisphere extratropics (20–90°N). The errors were calculated with reference to radiosonde observations over the period 1 June 2018 to 15 September 2018. The vertical bars indicate 95% confidence intervals.



There are not only biases in the model, but many conventional and satellite observation departures show systematic errors due to instrument configuration, approximations in radiative transfer calculations and other causes. To take these observational biases into account, a variational bias correction scheme (VarBC) has long been used in the 4D-Var system to estimate observation biases by finding corrections that minimise the systematic analysis observation departures. The choice of bias predictor determines the type of corrections that can be made, and VarBC has shown some skill in distinguishing between observation and model biases.

There is a clear interaction between weak-constraint 4D-Var and VarBC. As the systematic model error is better represented and corrected in the new weak-constraint 4D-Var, the VarBC correction should become smaller. The selection of the predictors for the different channels for each instrument was made many years ago based on tests in the strong-constraint 4D-Var framework. It should now be revisited for the channels sensitive to stratospheric temperature. More specifically, the usefulness of including air-mass predictors should be reassessed. Reducing the number of predictors might also help to reduce the VarBC inertia and the resulting long spin-up time. There is also a longer-term need to understand how many anchoring observations are required by VarBC and weak-constraint 4D-Var to correctly identify the sources of biases.

In principle, weak-constraint 4D-Var is able to eliminate biases not only in the analysis but also in forecasts. Results show that stratospheric temperature model bias in the IFS has little temporal variability, at least in the medium range (up to 10 days ahead). This property has allowed us to apply the same model error forcing over the whole 10-day forecast to correct the temperature model drift in the stratosphere towards the model’s climatology. It is unclear to what extent such an approach will work for longer forecast ranges (monthly to seasonal). It is conceivable that some spatial and temporal variability will need to be introduced in the model error forcing. This could be done, for example, by computing the mean and standard deviation of the model error estimates over different seasons to create time-dependent forcings.

Weak-constraint 4D-Var also has potential applications for climate reanalyses produced at ECMWF. Spurious climate signals are often introduced in reanalyses when a new type of observation is assimilated over a region which was previously poorly observed. Reducing model bias is critical as it will reduce such jumps in climate data records. The EU-funded Copernicus Climate Change Service (C3S) implemented by ECMWF is planning the production of atmospheric model integrations over the 20th century. Model simulations of the past using prescribed observationally based forcings and boundary conditions provide an important tool for understanding and estimating climate change. The model bias estimated by weak-constraint 4D-Var over a recent well-observed period could be used as a forcing throughout the whole century to reduce the impact of model drift on the reanalysis of earlier periods.

Further reading

Laloyaux, P., M. Bonavita, M. Dahoui, J. Farnan, S. Healy, E. Holm & S. Lang, 2020a: Towards an unbiased stratospheric analysis. Accepted for publication by QJRMS .

Laloyaux, P., M. Bonavita, M. Chrust & S. Gürol, 2020b: Exploring the potential and limitations of weak-constraint 4D-Var. Under review by QJRMS .

Shutts, G., M. Leutbecher, A. Weisheimer, T. Stockdale, L. Isaksen & M. Bonavita, 2011: Representing model uncertainty: stochastic parametrizations at ECMWF, ECMWF Newsletter No. 129, 19–24.

Trémolet, Y. & M. Fisher, 2010: Weak constraint 4D-Var, ECMWF Newsletter No. 125, 12–16.