Variational bias correction for Mode-S aircraft derived winds

Abstract

A variational bias correction (VarBC) technique is proposed within the three-dimensional variational (3D-Var) assimilation system of the Météo-France convective scale model AROME to remove biases present in horizontal wind components derived from aircraft Mode-S data. A set of two predictors derived from the geometry between ground and air speeds has been chosen in order to correct the aircraft speed and heading. Data collected from a set of 8 ADS-B antennas located over France for two periods in 2018 (April and October) have been used to assess the VarBC. Thanks to the availability of anchoring data (radiosoundings, AMDAR aircraft reports), the adaptive bias correction scheme is very efficient for improving the quality of Mode-S derived winds to a level that is similar to those from AMDAR and Mode-S bias corrected by the KNMI. The heading bias correction is the most important one in order to reach background departure values around 2.5 m/s. The behaviour of the VarBC is assessed more precisely for specific aircraft where its capacity to adapt to rapid changes due for example to aircraft maintenance is demonstrated. Moreover, when examining separately predictors for both wind components they appear to converge most of the time towards similar speed and heading corrections. Additional quality controls have allowed to reach background departure values of 2.2 m/s while increasing by a factor of ten the number of aircraft observations to be assimilated in the AROME-France domain.

Keywords:

• aircraft winds
• variational bias correction
• data assimilation
• Mode-S derived data

1. Introduction

Data Assimilation (DA) systems are designed to provide accurate initial conditions (so-called analysis) to Numerical Weather Prediction (NWP) models that are run operationally by national meteorological centres. The analysis is obtained from an optimized combination of available observations with a numerical short-range forecast. During the last fifteen years considerable efforts have been undertaken to develop DA systems for regional NWP models at convective-scale (grid mesh of a few kilometers) as described in the review of Gustafsson et al. (2018). Most improvements have come from the development of advanced DA algorithms (variational and ensemble methods), from more realistic description of physical processes (explicit deep moist convection) and from the use of dedicated observing systems. Indeed, concerning this last point, whereas in global NWP systems the vast majority of observations comes from polar orbiting satellites, this is not the case for regional NWP models where such satellites overpass the domain of interest only twice a day. Therefore, more suitable observing systems are necessary, in particular those capable of sampling the small spatial and fast temporal scales characterizing mesoscale weather systems. This requirement explains recent developments in order to exploit ground-based weather radars in terms of Doppler wind velocity and reflectivity (e.g. Montmerle and Faccani, 2009; Caumont et al., 2010; Wattrelot et al., 2014), ground based GNSS receivers informative about tropospheric water vapour (e.g. Mahfouf et al., 2015; Hdidou et al., 2020) and geostationary satellite infrared radiances (e.g. Guedj et al., 2011) in regional DA systems. Despite these dedicated observing systems, the amount of observations available at a given time providing knowledge about the atmosphere at convective-scale remains at least two orders of magnitude lower than the number of atmospheric variables to be initialized in such NWP models. To reduce this gap a number of areas are being explored towards new observing systems that could sample the fast temporally and small spatially convective scales within regional DA systems. The dominance of winds for weather evolution at small scales makes the use of additional accurate aircraft derived data (such as Mode-S winds) particularly appealing. The preparation of the future hyperspectral infra-red sounder IRS (Infra-Red Sounder) on board the geostationary satellite MTG (Meteosat Third Generation) to be launched in 2024 (Guedj et al., 2014) is also relevant in that context since this innovative instrument will meet spatial and temporal requirements to describe convective scales mostly for temperature, water vapour and clouds.

The objective of the present study is to describe preliminary activities undertaken at Météo-France (MF) towards the use of Mode-S derived winds in the convective-scale operational model AROME (Seity et al., 2011) using a three-dimensional variational (3D-Var) DA system (Brousseau et al., 2016). Aircraft observations have always been used in global and regional DA systems since they provide valuable and accurate upper-air measurements of temperature and horizontal wind components (Benjamin et al., 1999; Cardinali et al., 2003; Petersen, 2016; James and Benjamin, 2017; Doerenbecher and Mahfouf, 2019). They can be considered as opportunity measurements since, contrary to meteorological satellites or radiosondes, they are not operated by national weather centres to probe the atmosphere. For that purpose, coordination between commercial airlines, civil aviation and Air Traffic Control (ATC) is necessary at various levels (national weather services, regional such as EUMETNET, global through WMO). Various programmes exist to allow the worldwide reception by weather centres of aircraft derived data in real time through the Global Telecommunication System (GTS). Besides classical AMDAR (Aircraft Meteorological Data Relay) reports, where on board measurements are processed to derive meteorological parameters and then downlinked to a ground station (Painting, 2003), a new approach has been considered over the last ten years where raw data¹ available from Mode-S (Mode Selective) reports can be transmitted to ground and processed to derive meteorological parameters (de Haan, 2011). The primary advantage is to potentially collect much more data than conventional reports that are only available from specific aircraft equipped with on-board dedicated software. Moreover, there is a financial cost for meteorological agencies associated with air-to-ground data transmissions, limiting the number of messages that can be used. The direct interrogation of aircraft transponders by ATC Tracking and Ranging (TAR) Secondary Surveillance Radars (SSR) leads to Mode-S/EHS (Enhanced Surveillance) reports where the magnetic heading and the true airspeed information is available, in addition to altitude, Mach number, ground speed and heading. These data can be collected from radars, but also from ADS-B (Automatic Dependent Surveillance – Broadcast) receivers. In this case, the aircraft speed (from Mode-S) and position (from ADS-B) are not synchronous, which may lead to less accurate derived winds. It has been shown that the accuracy of such derived meteorological parameters is comparable to that of AMDAR reports for wind components (around 1.5 m/s) but is worse for temperature (around 3 K) that is deduced from the Mach number not recorded with sufficient precision (de Haan, 2016; Mirza et al., 2016; Stone, 2018).

A number of studies have been undertaken in the context of regional NWP (de Haan, 2011; de Haan and Stoffelen, 2012; Lange and Janjic, 2016) where the assimilation of Mode-S winds leads to positive impacts in terms of forecast skill scores since the amount of data that can be assimilated is about ten times larger than those provided by AMDAR. One current issue with Mode-S data is their near-real time availability to meteorological centres across countries. At national level, a network of ADS-B antennas can be set-up as it has been done over the Netherlands (de Jong et al., 2018), the United Kingdom (Stone and Pearce, 2016) and recently over France (B. Piguet, personal communication). Collaborations with ATC can be established for getting data acquired by their TAR SSR as done for by the KNMI over the MUAC (Maastrich Upper Area Control Centre) area (http://mode-s.knmi.nl/). The foreseen evolution of the EUMETNET ABO (Aircraft Based Observation) programme in the coming years should enable the near real time exchange of Mode-S data among European countries through a European Meteorological Aircraft Derived Data Centre (de Jong et al., 2018).

An important pre-processing step before using Mode-S reports in DA systems is to correct biases in the wind observations. Even though the speed and direction of the aircraft with respect to the ground are estimated accurately (radar signal, GPS positions), the on-board instruments leading to the true heading and airspeed have systematic errors. As a consequence, bias correction techniques have been proposed for reducing them. de Haan (2011, 2013) proposed to use heading during landing or AMDAR data to correct for heading biases, but in a NWP context, it is difficult to generalize to all aircraft. de Haan (2013) also proposed to consider a trial period (e.g. one month) where Mode-S winds are compared against NWP short-range forecasts (assumed to have small biases) to assess heading biases for all aircraft. However, such method is not able to consider new aircraft not available during the trial period nor can it adapt to changes if the bias evolves with time. We propose in this paper to use an adaptive Bias Correction (BC) scheme embedded in the 3D-Var DA system of the AROME Model (so-called variational bias correction (VarBC) scheme). This technique has been used for more than ten years for the assimilation of satellite radiances. Auligné et al. (2007) have shown that without observational constraints, the VarBC is contaminated, analysis after analysis, by the model bias. Nevertheless, due to the availability of numerous non-corrected observations that can be seen as additional constraints in a DA system, the VarBC has shown some good skill in distinguishing between NWP model errors and observational biases. This method uses implicitly the information redundancy between the various observing systems to identify the most likely source of bias. The VarBC method has been recently adapted to conventional observations such as AMDAR reports for the temperature bias correction (Isaksen et al., 2012; Gao et al., 2019) and ground-based GNSS signals (Sanchez-Arriola et al., 2016).

After an evaluation of the quality of the AROME model in terms of aircraft derived winds and the capability of VarBC to distinguish bias sources, we show here that it can also be applied to correct wind measurements from Mode-S reports, which has, to our knowledge, never been done before.

The outline of the paper is as follows. Section 2 summarizes the main differences between wind derived from various types of aircraft acquisition systems. Then, the 3D-Var system of the AROME-France model is described. In Section 3, the adaptation of the variational bias correction method to Mode-S winds is explained with an emphasis on the choice of relevant predictors. The efficiency of the method is preliminary assessed by examining, in Section 4, results from a 20-day observation monitoring experiment for AMDAR, raw and bias corrected Mode-S winds. In the same section, the efficiency of the method is explained and discussed in more detail by presenting the bias correction of Mode-S/ADS-B winds from specific aircraft. A summary of the results together with the main conclusions of the study are provided in Section 5.

2. Dataset description

2.1. Meteorological parameters from aircraft data

Horizontal wind components can be obtained from aircraft sensors that measure (or derive) the speed of the aircraft relative to the air, its position and ambient temperature and pressure. On one hand, a selection of such observations is transmitted to ground stations using the AMDAR system. On the other hand, Mode-S is a type of aircraft related meteorological information using both information from aircraft sensors and information derived from a tracking and ranging enhanced surveillance radar used for ATC. Raw data contain information on ground speed, track angle, magnetic heading, airspeed and Mach number. Every 4.2 s, ATC radars interrogate all aircraft in view (about 270 km around the radar; less near the ground, more around the flight levels due to the Earth curvature). Then, every aircraft transponder responds by sending Mode-S parameters within a specified radio frequency domain. These parameters are defined by the International Civil Aviation Organization (ICAO 2007). Within European designated Mode-S EHS airspace (like in France, United Kingdom and Germany), aircraft must disseminate messages when interrogated by a Mode-S EHS ATC radar. These Broadcast Dependent Surveillance (BSD) messages contain Mode-S parameters such as the registers BDS4.0, BDS5.0 and BDS6.0.

ADS–B is a surveillance technology where an aircraft determines its position via satellite navigation and periodically broadcasts it (in register BDS0.5 for instance), enabling it to be tracked.

In this study, we distinguish Mode-S/EHS and Mode-S/ADS-B data which differ by the way they are retrieved. Mode-S/EHS are gathered by the ATC means (where a precise aircraft position is given by the radar) whereas Mode-S/ADS-B are collected by ADS-B antennas which lead to less accurate aircraft locations from emitted ADS-B messages that are not synchronous with other BDS registers used to compute wind components (Table 1). A first dataset of Mode-S/ADS-B winds over 20 days is considered to describe and evaluate the bias correction method. Then a second period of one month is examined to tune more precisely the methodology. Such data have been collected in real time, by a semi-operational ADS-B antenna network set-up by Météo-France.

Table 1. BDS registers used to compute Mode-S/ADS-B winds.

Register	Message type	Content
BDS0.5	ADS-B	Aircraft position
BDS4.0	Mode-S	Pressure
BDS5.0	Mode-S	Air motion speed, ground motion speed, ground track angle
BDS6.0	Mode-S	Mach number, vertical speed, magnetic heading, air motion speed

The wind vector $\overrightarrow{V}$ can be derived from the difference between the aircraft motion relative to the ground and its motion relative to the air (Figure 1). Thus $\overrightarrow{V}$ is the difference in heading vector $\overrightarrow{V_a}$ defined by its module $\left|\overrightarrow{V_a}\right|$ and heading α, and the ground track vector $\overrightarrow{V_g},$ defined by ground speed $\left|\overrightarrow{V_g} \right|,$ and track angle α_g: $$\overrightarrow{V}=\overrightarrow{V_g}-\overrightarrow{V_a}.$$

Fig. 1. Wind vector composition. Aircraft air motion (red), aircraft ground motion (blue) and wind (green).

In this equation, both angles α and α_g are defined in the same reference framework. Mode-S heading α is given with respect to magnetic north (so-called magnetic heading α_m) whereas α_g is given relatively to the true north. Changing the reference introduces a wind error that has been assessed by de Haan (2011, 2013) with a number of statistical offline heading correction methods. An additional airspeed correction was also presented.

As ${\alpha }_m$ is referenced to magnetic North in Mode-S messages, $\alpha ={\alpha }_c\left(t\right)+{\alpha }_m+\Delta \left(t,\lambda ,\phi \right)$ where α_c(t) is a time (t) dependent correction in heading due to the reference change (from true to magnetic North) with onboard lookup tables (but that can be outdated) and Δ(t, λ, φ) is the magnetic declination. de Haan (2011) has shown that such formula is valid for most of aircraft in a regional domain. Figure 2 illustrates over France the difference between a possibly onboard magnetic declination computed with an old 2010 lookup table and the one obtained with the 2020 lookup table that could represent α_c(t) for a particular aircraft at a given time. The 0.3° difference over France in this case is lower than the resolution of ${\alpha }_m$ (1°). In this study, $\Delta \left(t,\lambda ,\phi \right)$ is based on the International Geomagnetic Reference Field IGRF-12. Generally, α_c(t) evolves slowly with time unless a major change is made onboard (maintenance, calibration …). de Haan (2013) showed that there is often a bias on $\left|\overrightarrow{V_a}\right|$ due to the quality of pitot probe measurements. He chose to correct $\left|\overrightarrow{V_a}\right|$ with a linear dependency on $\left|\overrightarrow{V_a}\right|$ itself. He showed that in most cases, the trend is close to one and that this correction leads to small corrections on horizontal wind components. So, here we assume that $\left|\overrightarrow{V_a}\right|=v_{\mathit{apitot}}+v_{ac}\left(t\right)$ with v_apitot being the true airspeed and v_ac(t) evolving slowly in time unless a major change made onboard (maintenance, calibration …).

Fig. 2. Difference in magnetic declination between 2020 and 2010 IGRF-13 lookup tables over France and surrounding countries.

These two corrections, corresponding to systematic errors, can be estimated from offline statistical comparisons between Mode-S parameters and their equivalent from other sources (NWP model outputs, AMDAR reports or direction of landing runway). They can be used, as at KNMI centre, to produce bias corrected winds from MUAC Mode-S/EHS data.

Table 2 displays the various sources of information necessary to derive wind and aircraft position from AMDAR, Mode-S/EHS and Mode-S/ADS-B.

Table 2. Summary of measured and derived meteorological parameters.

Parameter	AMDAR	Mode-S/EHS	Mode-S/ADS-B
Wind	Computed on board from va (airspeed), α (true heading) and ground track	Computed on ground from va (airspeed), αm (magnetic heading) and ground track	Computed on ground from va (airspeed), αm (magnetic heading) and ground track
Aircraft position	Computed on board from ground track	Computed on ground from ATC radar localization	Computed on ground from ADS-B location (computed on-board)

Figure 3 summarizes the various data streams from aircraft to bias correct wind files. Three distinct data streams are displayed: Mode-S/EHS data stream managed by KNMI over MUAC domain, Mode-S/EHS data stream currently under development managed by Météo-France with DSNA (French civil aviation: ‘Direction des Services de la Navigation Aérienne’) data and Mode-S/ADS-B data stream managed by Météo-France with data from its own antenna network. The main difference between KNMI and Météo-France data streams is that a preliminary step is necessary at KNMI to calculate heading and airspeed corrections with precomputed statistics whereas, for the Météo-France streams, corrections are dynamically estimated during the DA process to be described hereafter.

Fig. 3. Flowchart of Météo-France and KNMI Mode-S data fluxes and processing for data assimilation in AROME-France.

2.1.1. Mode-S data

Mode‐S/ADS-B data are collected using the Météo-France ADS-B network defined by 8 ADS-B antennas (Table 3). Each individual aircraft sends information when requested by ATC radars and can be in sight of several ADS-B antennas leading to possible redundancy. Figure 4 shows a map of the current network of antennas and their approximate ranges. Data are then gathered (10 s time window to make a wind observation with 4 BDS registers), a first quality control is performed (Annex 1) and redundant data are rejected. Data from all interrogated aircraft by ATC radars reach up to 8,000,000 non-redundant observations per day over France.

Fig. 4. Météo-France ADS-B antenna network (reception range ∼ 400 km). See Table 3 for exact positioning.

Table 3. Location of ADS-B antennas from the Météo-France network.

Antenna location	Longitude	Latitude
1.Nice	7.21 °E	43.66 °N
2.Clermont	3.17 °E	45.79 °N
3.Brest	4.41 °W	48.44 °N
4.Trappes	2.01 °E	48.77 °N
5.Caen	0.46 °W	49.18 °N
6.Lille	3.13 °E	50.61 °N
7.Strasbourg	7.63 °E	48.53 °N
8.Toulouse	1.37 °E	43.58 °N

Two periods, with different data samples have been considered in this study:

• from 5 to 24 April 2018 for Météo-France Mode-S/ADS-B data (the network was incomplete but sufficient for a preliminary study) and for KNMI Mode-S/EHS unbiased data.

• from 6 October to 6 November 2018 for Météo-France Mode-S/ADS-B data.

A comparison of daily density coverage of KNMI and Météo France maps is given in Figure 5. They encompass different spatial domains except for Belgium and Western Germany. The density of observations is the highest around major airports (Paris, Toulouse, Nice, Brussels, Schiphol, Frankfurt) and at flight level around 150 hPa. More data are available from the French ADS-B network (by a factor of two with respect to the KNMI).

Fig. 5. Daily density coverage of Météo-France Mode-S/ADS-B active data ² from 6 October to 6 November 2018 (a) and of KNMI Modes-S/EHS active data from 5 to 24 April 2018 (b) by pixels of 0.5°x0.5° (central panel) or 0.5°x50 hPa (right and top panels).

2.1.2. AMDAR data

AMDAR uses raw data available on board to monitor aircraft systems. Data rates typically vary from one sample per second to 16 samples per second. According to Painting (2003), temporal averages of raw data, depending on flight phase, are made to reduce the noise. This onboard computation is done before transmitting the data with typically averaging times for ascent and descent of 10 s and for flight level above an altitude of 20,000 ft (FL200) of 30 s.

The availability of AMDAR data also depends on the phase of flight (Petersen, 2016):

• Once every 7 min at flight level (3 min in option).

• In case of a ‘maximum wind event’ an extra observation is made.

• For an ascending aircraft, observations are made at intervals of 10 hPa (or 5 hPa in option) from ground level to 100 hPa below static pressure at take-off, followed by pressure observation interval of 50 hPa (or 25 hPa in option) until en-route height is reached.

• For a descending aircraft, observations are made at pressure intervals of 50 hPa (or optionally 25 hPa) until static pressure at touchdown minus 100 hPa. Below this level, the pressure observation interval is set to 10 hPa (or 5 hPa in option).

The availability also depends on contracts with equipped companies including budget issues associated with downlink costs.

Two periods with AMDAR data have been considered in this study:

• from 5 to 24 April 2018.

• from 6 October to 6 November 2018.

An example of data density of AMDAR data over Western Europe (October 2018) is displayed in Figure 6. The coverage is rather even over the domain with maximum density over the major European airports. Contrary to Mode-S data, the vertical density does not show maximum values at flight level. Indeed, the data selection policy of the EUMETNET E-AMDAR programme that manages AMDAR data over Europe favours profiles rather than en-route reports. The number of observations is about ten times lower than that of Mode-S data.

Fig. 6. Daily density coverage of AMDAR active data from 6 October to 6 November 2018 by pixels of 0.5°x0.5° (central panel) or 0.5°x50hPa (right and top panels).

2.2. The NWP model

The aircraft derived observations previously described are to be systematically compared with outputs from short-range forecasts of a convective-scale NWP model. Moreover, they will be included in its data assimilation system to examine the impact on resulting analyses and forecasts. The NWP model is the high-resolution limited area non-hydrostatic model AROME (Applications of Research to Operations at Mesoscale) run operationally at Météo-France since 2008 (Seity et al., 2011). The current configuration of the model has a horizontal grid size of 1.3 km (leading to an effective spatial resolution around 10 km as diagnosed from kinetic energy spectra by Ricard et al., 2013) with 90 vertical levels from the surface to 10 hPa. This model includes a comprehensive package of physical parametrizations for deep convection (explicitly resolved), turbulence and shallow convection, radiation and surface processes as described in Seity et al. (2011). A three-dimensional variational (3D-Var) DA system at 1.3 km resolution with an incremental formulation is run every hour to define the initial conditions (1 h window analysis). From the analysis, forecasts are run with various ranges from 3 h to 48 h. The 1 h forecast from the previous AROME run is used as background information for the next run in the analysis process. The 3D-Var system considers a wide variety of observing systems such as radiosoundings, synoptic ground stations and high density French surface network (RADOME), surface ocean ship and buoy reports, aircraft data, ground based GNSS Zenith Total Delays (ZTD), radar reflectivities and Doppler winds, Atmospheric Motion Vectors derived from satellite imagery, ocean wind vectors from scatterometers and satellite radiances from various polar orbiting platforms and from the geostationary satellite MSG. The lateral boundaries of the domain are provided by hourly interpolated forecasts from the Météo-France global NWP model ARPEGE (which does not assimilate any Mode-S data).

3. The variational bias correction scheme

3.1. Data assimilation experimental set-up

The least square estimation, so called BLUE (Best Linear Unbiased Estimator) is a classical linear technique to estimate the nearest state from the truth in the least square sense (with a minimal RMS error) assuming that both model data and observations are unbiased. In the 3D-Var context, the BLUE estimation of x_a can be obtained by the minimization of a cost function: $$J\left(x\right)=\frac{1}{2}{\left(y-H\left(x\right)\right)}^T{R}^{-1}\left(y-H\left(x\right)\right)+\frac{1}{2}{\left(x-x_b\right)}^T{B}^{-1}\left(x-x_b\right)=J_o\left(x\right)+J_b\left(x\right)$$ with: y: the observation vector, x_a: the analysis vector of the NWP model, x_b: the background vector of the NWP model (1-h forecast), B: the background error covariance matrix, R: the observation error covariance matrix and H: the observation operator (that projects the model space into the observation space).

The AROME 3D-Var system provides an analysis at a given time t₀ (the base time of the analysis). All observations from t₀ − 30 min to t₀ + 30 min and considered to be valid at time t₀ are gathered in a 1-h time slot.

All available observations go through a pre-processing task mainly composed of a quality control and a spatial and temporal thinning to insure a certain level of consistency between the observing system sampling and the model resolution. Data actually selected for the assimilation process are called ‘active observations’.

After redundancy checks, and a first level of quality controls during the creation of Mode-S winds (Annex 1), a second quality control (so-called ‘background check’) is undertaken using model counterparts (short-range forecast). This quality control definitely rejects data that are out of the bound $\left|y-H\left(x_b\right)\right|<5\sigma$ with $\sigma =\sqrt{{\sigma }_o^2+{\sigma }_b^2}$ given σ_o and σ_b, the standard deviations of observation and background errors.

For Mode-S wind data (u and v), σ_o is taken equal to the AMDAR values since de Haan (2011) has shown that errors for both sources are comparable. This observation error has a vertical dependency and σ_b is provided by the ARPEGE Ensemble Data Assimilation System.

The same spatial sampling used for AMDAR data has been chosen in a first stage for Mode-S data. For each aircraft only one active observation by 25 × 25 km² box model level is selected. When AMDAR and Mode-S data are available within the same selection box, priority is given to AMDAR. The so-called ‘passive observations’ have passed the quality controls but have been flagged by the spatial thinning. There is between 10 to 20 times more Mode-S passive data than active ones. In practice, passive observations are used in the 3D-Var system in such a way that they have no impact on the result of the minimization but contribute to the convergence of the VarBC coefficients.

3.2. Adaptive bias correction for Mode-S data

We assume that the observation biases b can be estimated with a simple linear parametric form: $$b=\sum _{i=1}^N\left({\beta }_i{P}_i\left(x,y_k\right)\right)$$ where

P_i(x,y_k) are suitable predictors that can depend on the control variable x and/or on metadata y_k varying with observation and β_i is the bias parameter for predictor i.

The variational bias correction (VarBC) is an adaptive method (Auligné et al., 2007) where the bias b (and so β_is) is computed at each analysis cycle, but the choice of predictors remains the same. In our case, VarBC tends to minimize systematic innovations for Mode-S wind components and, at the same time, tends to preserve (improve) differences between background and other observations in the analysis system.

Bias parameters are included in the control variable of the 3D-Var and can be considered as additional degrees of freedom in the J_o term: $$J_o\left(x,\beta \right)=\frac{1}{2}{\left(y-H\left(x\right)-\sum _{i=1}^N\left({\beta }_i{P}_i\left(x,y_k\right)\right)\right)}^T{R}^{-1}\left(y-H\left(x\right)-\sum _{i=1}^N\left({\beta }_i{P}_i\left(x,y_k\right)\right)\right)$$

Another term is added to the cost function to control the temporal inertia of the bias parameters: $$J_{\beta }\left(\beta \right)=\frac{1}{2}{\left(\beta -{\beta }_b\right)}^T{B}_{\beta }^{-1}\left(\beta -{\beta }_b\right)\mathrm{ }$$ where β is the predictor bias parameter vector, β_b is the predictor bias parameter vector from the previous analysis and B_β is the predictor bias parameter error covariance matrix.

In practice B_β is diagonal with an error variance set to σ²_o/N_b where N_b is a positive integer (stiffness parameter) which determines the adaptivity of β in the analysis cycle.

The bias parameters β are adjusted simultaneously with other analysis variables considering all available information. The adjustment is optimized by respecting uncertainties on the background, on the observations and the constraint imposed on the inertia of the Mode-S u and v bias parameters. In practice, a bias assessment of Mode-S winds that would lead to a degradation of the difference with other observing systems (a fortiori those not bias corrected) can hardly happen even if it reduces the differences between Mode-S wind and model counterpart. The analysis process uses redundant information given by other observing systems to decide what is the most probable error source. Indeed, in the 3D-Var system, data not bias corrected by VarBC, like radiosoundings and AMDAR reports (called ‘anchoring data’), still contribute to the cost function and act as a constraint on the assessment of the control variable and, in particular, of the VarBC bias parameters.

As explained above, it is possible to distinguish between so-called passive and active observations in the AROME 3D-Var. The active data are those that are used in the minimization to find the analysis state x_a. Passive ones are only considered in addition to the active ones, to estimate the bias parameters. So, data that have not been selected during the thinning process are nevertheless used to compute the bias (increasing the data sampling and therefore improving the convergence of the predictor bias parameters).

We now explain how relevant predictors have been defined in order to bias correct Mode-S/ADS-B winds. As shown in Figure 1, the two horizontal components can be written as: (1) $$u=u_g-u_a=u_g-\left|\overrightarrow{V_a}\right|\hspace{0.17em}\text{sin}\hspace{0.17em}\left(\alpha \right)$$(1) (2) $$v=v_g-v_a=v_g-\left|\overrightarrow{V_a}\right|\hspace{0.17em}\text{cos}\hspace{0.17em}\left(\alpha \right)$$(2)

The uncertainties du, dv, dv_a and dα, linked to observation errors of the zonal wind (u), the meridional wind (v), the aircraft air speed ($\left|\overrightarrow{V_a}\right|$) and its true North heading (α), can be written by a first-order derivation of Equations (1)(1) $$u=u_g-u_a=u_g-\left|\overrightarrow{V_a}\right|\hspace{0.17em}\text{sin}\hspace{0.17em}\left(\alpha \right)$$(1) and (2)(2) $$v=v_g-v_a=v_g-\left|\overrightarrow{V_a}\right|\hspace{0.17em}\text{cos}\hspace{0.17em}\left(\alpha \right)$$(2) and assuming no error on aircraft ground speed: (3) $$du=\underset{du1}{\underbrace{\stackrel{P_{u1}}{\overbrace{-\left|\overrightarrow{V_a}\right|\hspace{0.17em}\text{cos}\hspace{0.17em}\left(\alpha \right)}}\stackrel{{\beta }_{u1}}{\overbrace{d\alpha }}}}\underset{du2}{\underbrace{\stackrel{P_{u2}}{\overbrace{-\hspace{0.17em}\text{sin}\hspace{0.17em}\left(\alpha \right)}}\stackrel{{\beta }_{u2}}{\overbrace{{dv}_a}}}}$$(3) (4) $$dv=\underset{dv1}{\underbrace{\stackrel{P_{v1}}{\overbrace{\left|\overrightarrow{V_a}\right|\hspace{0.17em}\text{sin}\hspace{0.17em}\left(\alpha \right)}}\stackrel{{\beta }_{v1}}{\overbrace{d\alpha }}}}\underset{dv2}{\underbrace{\stackrel{P_{v2}}{\overbrace{-\hspace{0.17em}\text{cos}\hspace{0.17em}\left(\alpha \right)}}\stackrel{{\beta }_{v2}}{\overbrace{{dv}_a}}}}$$(4)

By choosing the bias predictors P_u1, P_u2, P_v1 and P_v2 as shown in the above equations and assuming a constant bias, the VarBC shall tend to make β_u1 and β_v1 converging towards dα and β_u2 and β_v2 converging towards dv_a. Obviously, β_u and β_v will be correctly assessed by the VarBC method only if predictors are not always zero (North to South flights or West to East flights leading either to null values of cos(α) or sin(α)). Note that the predictors depend only on observation metadata. Nevertheless, since α_m is available in Mode-S/ADS-B data instead of α, we assume here that ${\alpha }_m+\Delta \left(t,\lambda ,\phi \right)$ is a good approximation of α, and that v_apitot is a good approximation for $\left|\overrightarrow{V_a}\right|.$

Finally, the set of four chosen predictors (two for each wind component) is: (5) $$P_{u1}=-v_{\mathit{apitot}}\hspace{0.17em}\text{cos}\hspace{0.17em}\left({\alpha }_m+\Delta \left(t,\lambda ,\phi \right)\right)$$(5) (6) $$P_{u2}=-\hspace{0.17em}\text{sin}\hspace{0.17em}\left({\alpha }_m+\Delta \left(t,\lambda ,\phi \right)\right)$$(6) (7) $$P_{v1}=v_{\mathit{apitot}}\hspace{0.17em}\text{sin}\hspace{0.17em}\left({\alpha }_m+\Delta \left(t,\lambda ,\phi \right)\right)$$(7) (8) $$P_{v2}=-\hspace{0.17em}\text{cos}\hspace{0.17em}\left({\alpha }_m+\Delta \left(t,\lambda ,\phi \right)\right)$$(8)

β_u1 and β_v1 are named ‘heading bias parameters’ and β_u2 and β_v2 ‘airspeed bias parameters’.

P_u1 and P_v1 are named ‘heading bias predictors’ and P_u2 and P_v2 ‘airspeed bias predictors’.

It is important to notice that heading ${\alpha }_m$ is given with 1° resolution. The more observations are available for a given aircraft in an assimilation window, the more accurate the bias assessment will be.

The VarBC method uses reduced centred predictors. All predictors are by definition centred. Standard deviations are the same for P_u1 and P_v1 (given by the dataset statistics from 5 to 24 April 2018, that is 150 m/s), on one hand, and for P_u2 and P_v2 (known sin and cos standard deviations, that is 0.707), on the other hand.

Therefore, the VarBC method has been implemented in AROME using the above predictors where airspeed (resp. heading) predictors have been divided by 0.707 (resp. 150 m/s).

4. Results

A number of statistics (mean and standard deviation) are presented over various time periods in terms of BackGround Departure (BGD) that is a useful parameter to compare observed data with model counterparts (in terms of short-range forecasts): $$\mathit{BGD}=y-H\left(x_b\right)\text{ at basetime }{t}_0$$

In the following, we compare raw Mode-S BGD of u and v (BGD_raw) with BGD obtained with bias corrected Mode-S data using the VarBC (BGD_varbc): $${\mathit{BGD}}_{\mathit{raw}}=y_{\mathit{raw}}-H\left(x_b\right)$$ $${\mathit{BGD}}_{\mathit{varbc}}=y_{\mathit{raw}}-H\left(x_b\right)-\sum _{i=1}^2\left({\beta }_i{P}_i\left(y_k\right)\right){\text{ with β}}_i\text{ computed at }{t}_0$$

The BGD statistics, computed only on active data, can be used to assess the main characteristics of the observation errors knowing the background error σ_b and assuming that short-range forecasts are mostly unbiased.

4.1. Monitoring statistics from the first period

A first experiment has been undertaken from 5 to 24 April 2018 to monitor Météo-France (MF) Mode-S/ADS-B and KNMI Mode-S/EHS winds for a preliminary assessment of the VarBC method. Quality controls and sampling choices applied to the KNMI Mode-S/EHS and MF Mode-S/ADS-B VarBC winds are the same as for AMDAR winds before entering the AROME 3D-Var. For MF Mode-S/ADS-B wind data, VarBC has been configured with a very low inertia (small value of Β_β) in this first experiment to allow the bias predictor parameters to change significantly from one analysis to another (J_β∼0). Such experiment was not cycled (except for β_i) which means that each background is taken from a 1 h forecast of the operational suite AROME such that the NWP model fields for Mode-S comparisons do not contain any Mode-S information (by contrast to AMDAR winds that are assimilated).

During this period, we compare different sources of aircraft wind data: AMDAR (assimilated, cycled), MF Mode-S/ADS-B (assimilated to compute bias, not cycled) and KNMI Mode-S/EHS (monitored). It is important to stress that VarBC Mode-S and MUAC data have been processed independently. Potentially, data from a given aircraft can be quasi collocated for each of the three datasets. If all of them were assimilated, such thinning would not be optimal, but in a VarBC context, AMDAR data will be used as anchoring observations for the VarBC MF Mode-S/ADS-B data.

4.1.1. Global behaviour

Table 4 presents global BGD mean and SD values of u and v wind components for raw MF Mode-S/ADS-B, VarBC MF Mode-S/ADS-B, KNMI unbiased Mode-S/EHS and AMDAR data over the 20-day period. Note that the comparison is done with different samples over the period of interest (see Section 2.1). One can notice the large positive impact of the VarBC on the overall quality of MF Mode-S/ADS-B winds by comparing results with other similar observing systems, with a reduction of the SD values by a factor of 1.6-1.7. By construction, constant heading and airspeed biases lead to variable bias values in terms of wind components (see Equations (3)(3) $$du=\underset{du1}{\underbrace{\stackrel{P_{u1}}{\overbrace{-\left|\overrightarrow{V_a}\right|\hspace{0.17em}\text{cos}\hspace{0.17em}\left(\alpha \right)}}\stackrel{{\beta }_{u1}}{\overbrace{d\alpha }}}}\underset{du2}{\underbrace{\stackrel{P_{u2}}{\overbrace{-\hspace{0.17em}\text{sin}\hspace{0.17em}\left(\alpha \right)}}\stackrel{{\beta }_{u2}}{\overbrace{{dv}_a}}}}$$(3) and (4)(4) $$dv=\underset{dv1}{\underbrace{\stackrel{P_{v1}}{\overbrace{\left|\overrightarrow{V_a}\right|\hspace{0.17em}\text{sin}\hspace{0.17em}\left(\alpha \right)}}\stackrel{{\beta }_{v1}}{\overbrace{d\alpha }}}}\underset{dv2}{\underbrace{\stackrel{P_{v2}}{\overbrace{-\hspace{0.17em}\text{cos}\hspace{0.17em}\left(\alpha \right)}}\stackrel{{\beta }_{v2}}{\overbrace{{dv}_a}}}}$$(4) ). Therefore, the VarBC method is designed to reduce u and v bias variability, which explains why it reduces BGD SD values. No changes are noticed on mean values that remain close to zero. Indeed, even if for individual aircraft wind component biases can be rather large, given the number of aircraft providing Mode-S data, when averaging all aircraft winds, the mean value is close to zero even for biased data. It is important to notice that close to zero mean differences between AROME and KNMI unbiased Mode-S/EHS show that AROME model is not biased compared to unbiased non-assimilated data. AMDAR and KNMI unbiased Mode-S/EHS winds have similar values with the lowest BGD mean and SD values. Results with VarBC Mode-S/ADS-B winds are very similar. Due to on-board computations AMDAR winds have a better quality than Mode-S ones. Mode-S/EHS data have better statistics than those from Mode-S/ADS-B that is explained by a better precision of aircraft locations. Other differences can be explained by different quality controls applied by each data producer.

Table 4. Global BGD mean and SD of u and v wind components for raw MF Mode-S/ADS-B, VarBC MF Mode-S/ADS-B, KNMI unbiased Mode-S/EHS and AMDAR data from 5 to 24 April 2018.

u BGD mean (m/s)	u BGD SD (m/s)	v BGD mean (m/s)	v BGD SD (m/s)	Number of active data	Source
0.2	4.5	−0.1	4.3	∼5,300,000	Raw MF Mode-S/ADS-B
0.2	2.7	−0.1	2.7	∼5,300,000	VarBC MF Mode-S/ADS-B
0	2.4	−0.1	2.4	∼2,200,000	KNMI unbiased Mode-S/EHS
−0.1	2.4	−0.1	2.4	∼400 000	AMDAR

4.1.2. Behaviour of an individual aircraft (aircraft A1)

The behaviour of the VarBC is assessed more precisely by examining an individual aircraft. We have selected an aircraft having one of the largest amount of Mode-S data during the period (4583 active data) with high BGD values on raw wind components.

Figure 7 shows the improvement on BGD statistics by the VarBC method. The bias correction narrows the sample distribution and centers it around its median with mean values close to zero. Moreover, the u (resp. v) component BGD standard deviation decreases from 10 m/s to 2.5 m/s (resp. 6.9 m/s to 2.6 m/s).

Fig. 7. Boxplots of u (blue) and v (red) BGD for aircraft A1 for raw (boxplots 1 and 3) and VarBC (boxplots 2 and 4) Mode-S/ADS-B data from 5 to 24 April 2018. Dots: standard deviation, circles: mean.

Figure 8 (top panel) displays the time evolution of the zonal wind component BGD values. The large positive and negative periodic biases present in raw data are almost suppressed by the VarBC method. The diagnosed biases are shown in Figure 8 (bottom panel). The flights of this aircraft contain a good sample of all possible directions (displayed by the green circles). Obviously, the sign of the computed heading bias depends on the heading of the aircraft (see Equation (3)(3) $$du=\underset{du1}{\underbrace{\stackrel{P_{u1}}{\overbrace{-\left|\overrightarrow{V_a}\right|\hspace{0.17em}\text{cos}\hspace{0.17em}\left(\alpha \right)}}\stackrel{{\beta }_{u1}}{\overbrace{d\alpha }}}}\underset{du2}{\underbrace{\stackrel{P_{u2}}{\overbrace{-\hspace{0.17em}\text{sin}\hspace{0.17em}\left(\alpha \right)}}\stackrel{{\beta }_{u2}}{\overbrace{{dv}_a}}}}$$(3) ). For this particular aircraft the heading bias mainly explains the total bias values (as blue and black lines are quasi superimposed) that oscillate between −10 to 10 m/s.

Fig. 8. Evolution for each observation of a aircraft A1 from 5 to 24 April 2018 of Mode-S/ADS-B u component raw (black) and VarBC (blue) BGD in top panel and VarBC computed biases (total bias black, heading bias blue, airspeed bias red) with aircraft heading (green dots) in bottom panel.

Figure 9 shows the evolution, analysis after analysis, of the u (blue curve) and v (red curve) bias parameters for this aircraft (that can evolve freely since J_β∼0). For both components, heading (top panel) and airspeed (middle panel) bias parameters fluctuate around their mean values. The hypothesis that dα and dv_a can be considered as quasi-constants approaching α_c(τ) and v_ac(τ) values and that the VarBC method is able to assess dα and dv_a appears to be correct for this example. Indeed, bias parameters fluctuate closely around a constant value supposed to be dα and dv_a and both independent computations of bias parameters (through u and v variables) lead to very similar asymptotic values. Fluctuations around the mean of the bias parameters can also be seen as a limitation of the VarBC method that cannot separate model and observation errors in cases where there is no nearby anchoring data during analysis process. Increasing the temporal inertia on the bias parameters could reduce such behaviour. It is interesting to note that, for u and v heading parameters, it takes only four analyses at the beginning of the period to converge towards the asymptotic value from a cold start (zero values).

Fig. 9. Evolution, for a aircraft A1, by analysis where data are assimilated from 5 to 24 April 2018 of u (blue) and v (red) diagnosed heading bias (top panel), airspeed bias (middle panel) and passive³ data number by analysis (bottom panel).

4.2. Monitoring statistics from the second period

During the period from 6 October to 6 November 2018, two VarBC Mode-S experiments were conducted. These experiments were cycled (i.e. the analysis uses its own previous short-range forecast for the background). For both of them, the stiffness parameter N_b of the VarBC has been increased with respect to the first period from a value of 10 to a value of 250 corresponding the median of all aircraft number of observations (when available) per assimilation window over the period. This means that, for a given aircraft, only 1-h time slot with at least 250 observations is necessary to significantly change the bias parameter values during analysis (instead of 10 for the first experiment). The number of observations per aircraft and per assimilation window (when available) varies between 1 and 1,500 over the period.

The goal of the first experiment was to compute the bias parameters for all aircraft present during the period. By considering values at the end of the period, aircraft with no conclusive bias corrected results (BGD SD above 3 m/s and absolute mean larger than 0.5 m/s) are blacklisted. This leads to a 6% reduction of active data (in the AROME 3D-Var) and to a 29% reduction of the number of used aircraft.

For the second experiment, aircraft with suitable bias parameters at the end of the first experiment have been selected and assimilated. This represents around 467 000 active and 7 866 000 passive Mode-S/ADS-B data per day. By comparison, during the same period, only 24 000 daily active AMDAR data have been assimilated (5% of Mode-S data). Results shown hereafter describe this second experiment.

4.2.1. Global behaviour

Figure 10 shows the evolution of main BGD statistics after the use of VarBC with two predictors for both u and v parameters to unbiase Mode-S data as described above. These results are compared to similar statistics for AMDAR winds. Aircraft with more than 200 active Mode-S and AMDAR data over the period have been considered to provide reliable BGD statistics, representing 6025 (resp. 738) aircraft over the 7603 (resp. 1110) available after quality controls during the period for Mode-S (resp. AMDAR) data. It is important to stress that u and v raw Mode-S/ADS-B BGD mean and SD values can be very different from an aircraft to another. Therefore, even if half of the aircraft have BGD means close to zero for raw data, the remaining ones exhibit mean differences reaching up to 10 m/s. For the same sample, after bias correction, the spread of BGD mean values around zero has been reduced, with maximum values around 2.4 m/s. Among the raw Mode-S/ADS-B dataset, 75% of both wind components BGD SD values are above 3 m/s and can reach values above 15 m/s. The VarBC allows to reduce significantly BGD SD values. Results are better than for the first period mainly because a number of aircraft have been blacklisted. The global BGD SD value is around 2.2 m/s for u and v Mode-S/ADS-B data. These results can be directly compared to those of AMDAR data (2.4 m/s for u and 2.45 m/s for v). This comparison confirms the good behaviour of the VarBC.

Fig. 10. Boxplots of aircraft u (blue) and v (red) BGD mean (light colors) and SD (dark colors) for raw MF Mode-S/ADS-B (left), VarBC MF Mode-S/ADS-B (middle) and AMDAR (right) datasets from 6 October to 6 November 2018.

Complementary information on random and systematic errors of Mode-S winds can be obtained using the triple collocation method (Stoffelen, 1998) with additional datasets from AMDAR aircrafts and AROME 1 h forecasts. This method allows to estimate simultaneously modelling and calibration errors in order to avoid pseudo biases that could be introduced by other intercomparisons. Results are shown in Appendix A. Collocated AMDAR and VarBC Mode-S u and v data differences are much lower than those between AMDAR and AROME and those between VarBC Mode-S and AROME. This result indicates that the variational scheme, constrained by all other non-corrected assimilated data, shows significant skill in distinguishing between the sources of biases. It also allows to assess random observation errors (without representativeness error) for both AMDAR and VarBC Mode-S (resp. ∼1 and ∼0.77 m/s) and 1 h AROME forecast error (∼1.74 m/s).

From Equation (1)(1) $$u=u_g-u_a=u_g-\left|\overrightarrow{V_a}\right|\hspace{0.17em}\text{sin}\hspace{0.17em}\left(\alpha \right)$$(1) , we define the following quantities: $$W_{ui}^p=\overline{\left\{\frac{100\times \left|{\beta }_{ui}^p{P}_{ui}^p\right|}{{\sum }_{i=1}^2\left|{\beta }_{ui}^p{P}_{ui}^p\right|}\right\}}\text{ and }{W}_{vi}^p=\overline{\left\{\frac{100\times \left|{\beta }_{vi}^p{P}_{vi}^p\right|}{{\sum }_{i=1}^2\left|{\beta }_{vi}^p{P}_{vi}^p\right|}\right\}}$$ representing the u and v mean absolute weights in percent of the absolute bias amplitude explained by predictor i for all Mode-S observations from aircraft p.

We also define: $$B_u^p=\overline{\left\{{\sum }_{i=1}^2\left|{\beta }_{ui}^p{P}_{ui}^p\right|\right\}}\text{ and }{B}_v^p=\overline{\left\{{\sum }_{i=1}^2\left|{\beta }_{vi}^p{P}_{vi}^p\right|\right\}}$$ that are the mean absolute amplitudes of the u and v VarBC biases for all observations from aircraft p.

Similarly, we define: $$B_{Ru}^p=\frac{100\times B_u^p}{{\mathit{max}}_p\left(B_u^p\right)}\text{ and }{B}_{Rv}^p=\frac{100\times B_v^p}{{\mathit{max}}_p\left(B_v^p\right)}$$ as being the u and v relative amplitudes with respect to the maximum amplitude value of all observations from aircraft p. We found over the one month period ${\mathit{max}}_p\left(B_u^p\right)=13.6m/s$ and ${\mathit{max}}_p\left(B_v^p\right)=16.8m/s.$

Figure 11 shows values of W^p_i and B_R^p for u and v by individual aircraft ordered by increasing values of W^p_i. For 80% of the aircraft, the heading bias predictor explains most of the VarBC bias absolute amplitude. For both u and v components after correction, the larger (resp. smaller) the amplitude of the absolute bias is, the more important the weight of the first (resp. second) predictor is. This is coherent with the results of de Haan (2013) where he has shown that corrections on $\left|\overrightarrow{V_a}\right|$ have a much lower impact on u and v parameters quality than the one on α.

Fig. 11. Fractional mean absolute bias amplitude by predictor W^p_i (solid: heading, dashed: airspeed) and relative total mean absolute bias amplitude B_r^p (grey dots with linear regression) for u (top panel) and v (bottom panel) by individual aircraft ordered by increasing values of W^p₁ from 6 October to 6 November 2018.

To assess globally the predictor choice hypotheses, Figure 12a presents, for aircraft with at least 30 analyses, the u and v heading and airspeed biases. For these aircraft, the u heading bias remains close to its mean value (percentile 5 and 95% close to the mean), and both heading biases are close to each other and vary between −2° and 5°. Airspeed biases exhibit the same kind of behaviour but with more spread around their mean values. It is mainly due to the fact that the VarBC technique uses normalized predictors. They are presented here without normalization in order to assess directly the order of magnitude of heading and airspeed biases.

Fig. 12. (a) u (blue) and v (red) computed heading (top panel) and airspeed (bottom panel) biases means by aircraft ordered by increasing u bias parameter mean from 6 October to 6 November 2018. Grey: percentile 5 and 95 of u and v bias parameters means by aircraft. (b) Distribution of estimated heading biases from u (blue) and v (red) wind components (6 October to 6 November 2018).

Distribution of heading biases shown in Figure 12a is displayed in Figure 12b. There are three main modes with values around 0.1 (25% of cases close to zero), 1.15 and 1.7 degrees. Similarly, de Haan (2013) has also found three classes but with other numerical values mainly because his study used on ground declination tables older than ours, considered another set of aircraft over a different domain with different on board declination tables for some of them.

In addition, u and v airspeed biases can, in particular cases, be very different. The study of three particular aircraft will highlight a number of contrasted behaviours.

4.2.2. Behaviour of individual aircraft

Three contrasted aircraft have been chosen. The first one (aircraft A1) has a large heading bias (top right in Figure 12a, same as studied previously in Section 4.1). The second one (A2) presents a particular behaviour with good general results in terms of BGD statistics but with large changes in bias parameter values during the period. The third one (A3) can be spotted in Figure 12a on airspeed bias (right hand side plot): it exhibits very different u and v values.

First example: aircraft A1. For the same aircraft studied in Section 4.1 and mentioned above, the VarBC method leads to similar results for this second period: BGD u (resp. v) SD values goes down from 9.6 m/s to 2.07 m/s (resp. 8.4 m/s to 2.2 m/s). Figure 13 shows the evolution, analysis after analysis, of the u and v bias predictor parameters for this aircraft. Since the bias parameter values do not start with zero values, the VarBC method is rather efficient from the beginning of the experiment. By comparison with the first period, bias parameter values fluctuate less around their mean due to the larger inertia. Values of u and v heading bias means are around $\overline{d\alpha }\sim 3.68°$ and values of u and v airspeed bias parameters means around$\overline{{ dv}_a}\sim 0.8m/s.$ Note that $\overline{d\alpha }$ has evolved since the first period (5 months ago) where it was $\overline{d\alpha }\sim 3.42°.$ For this aircraft, a statistical study of the differences between ${\alpha }_m+\Delta \left(t,\lambda ,\phi \right)$ for take-off and landing and axes of the Toulouse/Blagnac airport runway during the same period (11 values retained) led independently to a value of α_c(t)∼3.45° (knowing the 1° precision of Mode-S/ADS-B heading data).

Fig. 13. Evolution of u (blue) and v (red) diagnosed heading bias (top panel), airspeed bias (middle panel) and passive data number by analysis (bottom panel) for aircraft A1 by available analysis from 6 October to 6 November 2018.

Second example: aircraft A2. Results for this aircraft are particularly interesting because, during the period, a maintenance took place around the 17th of October (one day without data) to improve the quality of its heading. Until that date, zonal wind heading bias correction (du1 in Equation (3)(3) $$du=\underset{du1}{\underbrace{\stackrel{P_{u1}}{\overbrace{-\left|\overrightarrow{V_a}\right|\hspace{0.17em}\text{cos}\hspace{0.17em}\left(\alpha \right)}}\stackrel{{\beta }_{u1}}{\overbrace{d\alpha }}}}\underset{du2}{\underbrace{\stackrel{P_{u2}}{\overbrace{-\hspace{0.17em}\text{sin}\hspace{0.17em}\left(\alpha \right)}}\stackrel{{\beta }_{u2}}{\overbrace{{dv}_a}}}}$$(3) ) fluctuates with amplitude of ∼7 m/s (Figure 14). Later, the data quality is improved and heading bias correction amplitude goes down to ∼2 m/s. From Figure 15 it appears that the VarBC method only takes four analyses to adapt the correction and to reach a new equilibrium. Note that, at the beginning of the period, the VarBC predictor bias parameters where set to the values at the end of the period from the first experiment. This explains why four analyses are necessary to stabilize them.

Fig. 14. Evolution of Mode-S/ADS-B u component raw (black) and VarBC (blue) BGD in top panel and VarBC computed biases (total bias black, heading bias blue, airspeed bias red) with aircraft heading pointed in green in bottom panel (from 6 October to 6 November 2018 for each observation of aircraft A2). The vertical dotted line corresponds to the 17 October 2018.

Fig. 15. Evolution, of u (blue) and v (red) diagnosed heading bias (top panel), airspeed bias (middle panel) and passive data number by analysis (bottom panel) for aircraft A2, by available analysis from 6 October to 6 November 2018. The vertical dotted line corresponds to the 17 October 2018.

Third example: aircraft A3. The characteristics of this particular aircraft are shown in Figure 16. It appears that the VarBC method impacts mainly the u BGD mean value. Figure 17 shows that the u heading predictor explains most of corrected bias. Even if the VarBC method seems to be efficient for this aircraft (see Figure 16), two very different and stable airspeed bias parameters are exhibited in Figure 18: the first one is around 0 m/s and the second one is around 10 m/s. However, this large value has no impact on the computed bias because, as displayed in Figure 17, headings are always close to zero (flights through the AROME domain are very often from South to North). This means that P_u2 and P_v1 are always close to zero (see Equations (6)(6) $$P_{u2}=-\hspace{0.17em}\text{sin}\hspace{0.17em}\left({\alpha }_m+\Delta \left(t,\lambda ,\phi \right)\right)$$(6) and (7)(7) $$P_{v1}=v_{\mathit{apitot}}\hspace{0.17em}\text{sin}\hspace{0.17em}\left({\alpha }_m+\Delta \left(t,\lambda ,\phi \right)\right)$$(7) ). Whatever the values of β_u2 and β_v1 found by the VarBC method, only β_u1 and β_v2 impact the computed bias. They appear to be realistic in terms of airspeed and heading biases. The hypothesis that u and v predictor parameters converge toward airspeed and heading biases is only valid with non-zero predictors (see Section 3.2). Therefore, the VarBC method is efficient even in cases where two predictors are null. Note that, in this particular situation, the u and v bias predictor parameters difference is not a relevant criterium to reject such aircraft.

Fig. 16. Boxplots of u (blue) and v (red) BGD for aircraft A3 for raw (boxplots 1 and 3) and VarBC (boxplots 2 and 4) Mode-S/ADS-B data from 6 October to 6 November 2018. Dots: standard deviation, circles: mean.

Fig. 17. Evolution of Mode-S/ADS-B u component raw (black) and VarBC (blue) BGD in top panel and VarBC computed biases (total bias black, heading bias blue, airspeed bias red and aircraft heading pointed in green in bottom panel) from 6 October to 6 November 2018 for each observation of aircraft A3.

Fig. 18. Evolution, of u (blue) and v (red) diagnosed heading bias (top panel), airspeed bias (middle panel) and passive data number by analysis (bottom panel) for aircraft A3 by available analysis from 6 October to 6 November 2018.

4.3. Discussion

Model errors (linked to unresolved scales, parametrisations and initial conditions) can be persistent and/or of large scale. Depending on the available aircraft data, they can have an impact on both static and dynamic unbiasing methods. The VarBC method, being dynamic, could, if it was not able to distinguish the origin of the bias, convert any model error into an observational bias as shown by Auligné et al. (2007). The VarBC method applied to Mode-S data avoids this pitfall for the following reasons:

• It relies on AMDAR data (in addition to other wind measurements: AMVs, radiosoundings, surface weather stations, radar Doppler winds, etc.) as anchoring data. The AMDAR data are of the same type, at the same levels as Mode-S data and are also measured by the same aircraft.

• It uses a large amount of Mode-S data per hourly analysis window. Given the typical speed of an aircraft (∼200 m/s), it can travel up to 720 km within an assimilation window thereby providing a lot of data along its route (every 4 s). The small-scale local phenomena it will encounter, and the related model problems, will be largely mitigated by the averaging effect implied by the VarBC method when minimizing the cost-function of the 3D-Var. In addition, a large amount of anchoring data will be available that can potentially be used by the VarBC method to constrain the bias correction along a specific route.

• It can stabilize the evolution of the bias parameters using a temporal inertia. Without inertia, if an aircraft is scarcely present in the assimilation window (reduced amount of data) without nearby anchoring data, the bias correction could convert model error into observational bias. This can be limited, in this case, by adjusting the inertia parameter. Thus, even if the VarBC method has sufficient constraints in terms of anchoring data to manage cases with reduced amount of data for a given aircraft, tuning the inertia value allows us to give priority to aircraft travelling long distances within an assimilation window (ie. providing large amounts of data). In summary, within an assimilation window, a well set-up inertia value will insure bias parameters to vary only for well sampled aircraft and to stay quasi unchanged for scarcely ones.

• It operates mostly in free atmosphere regions where small-scale processes are less important than close to the surface.

Finally, aircraft with only few data can be problematic, but it is true for both VarBC and any other static method. Based on the BLUE principle of data assimilation, the VarBC method allows for a rapid convergence towards realistic results.

5. Conclusions and perspectives

The main objective of the study described in this paper was to apply the VarBC correction method to reduce, within a NWP DA system, the biases on horizontal wind components derived from Mode-S/ADS-B messages. This is, to our knowledge, the first time that such method is applied to Mode-S derived winds. For that purpose, we have considered the 3D-Var DA system of the Météo-France convective scale model AROME and Mode-S/ADS-B data gathered by a network of 8 antennas installed by Météo-France over the French metropole. To set-up the VarBC method, two trivial predictors related to known biases of Mode-S/ADS-B wind data have been chosen. For a given aircraft, the first one is linked to the heading whereas the second one is linked to the airspeed.

When sufficient wind Mode-S/ADS-B data are available for a specific aircraft, the airspeed and heading biases can be efficiently estimated by the VarBC method for both wind components. Therefore, for these aircraft, there is a significant improvement in Background Departure (BGD) (observation minus NWP short-range forecast) statistics for mean and Standard Deviation (SD) values.

From a first 20-day monitoring experiment, a preliminary assessment has shown the capacity of the VarBC method to compute heading and airspeed biases, and thus to improve the derived wind components from Mode-S wind data. First results exhibit BGD SD and mean values of bias corrected Mode-S/ADS-B winds comparable to the ones obtained with AMDAR or Mode-S/EHS (as processed and unbiased by the KNMI) data. It also shows that the AROME model is not biased at aircraft flight levels.

After removing aircraft having bias corrected wind BGD SD values larger than 3 m/s and absolute mean values larger than 0.5 m/s, 94% of the initial sample from the first monitoring period has been kept for a second period (6 October to 6 November 2018) leading to a set of 467,000 Mode-S/ADS-B wind observations suitable for assimilation over France. For this revised sample, the bias corrected BGD means are close to zero and SD values reach 2.2 m/s (comparable to statistics shown by de Haan, 2013, for Mode-S/EHS winds).

For 80% of aircraft with more than 200 active data, the heading bias predictor explains most of the bias correction. For both u and v components, the larger (resp. smaller) the bias correction is, the more important the weight of the heading (resp. airspeed) predictor is.

By definition within the VarBC method, heading and airspeed bias parameters tend towards their asymptotic values (i.e. actual heading and airspeed biases), if they are constant during the period of interest and if the predictors are not always zero. To insure that the VarBC method does not absorb significant NWP model errors, nearby ‘anchoring’ observations (Eyre, 2016) are necessary. This is the case for Mode-S data that often sample the same surrounding atmosphere as other good quality aircraft data (AMDAR, AIREP). To prevent bias parameters from changing significantly from one analysis to another, their associated inertia could be increased. We have verified with triple collocation method (Stoffelen, 1998) using 1844 AMDAR, Mode-S and AROME data that VarBC Mode-S u and v data are much closer to AMDAR data than to AROME ones. It tends to prove that the variational scheme, constrained by all other non-corrected assimilated data, shows significant skill in distinguishing between the bias sources.

However, we have also demonstrated the capacity of the VarBC method to rapidly modify the bias parameters (few hours) to changes in the raw data due for example to aircraft maintenance. Such behaviour gives a significant advantage to the VarBC method compared to offline statistical methods for bias corrections. Therefore, a compromise has to be found for the control of the inertia between these two opposite goals, that should be guided by the size of the data sample. Nevertheless, this is not straightforward since aircraft samples can be very different (infrequent long-haul flights, frequent medium-haul flights, aircraft crossing the domain in the middle or on the sides). Since the same inertia is used for all aircraft samples, the VarBC behaviour can be rather different from one aircraft to another. A pre-sampling step before assimilation could help to filter the samples in order to lead to a similar impact of the inertia on each of them.

It has been shown that for specific aircraft, flights patterns can be such that some predictors are always (or over a number of assimilation windows) close to zero. In these cases, the bias predictor parameters can tend to unrealistic large values in terms of airspeed or heading biases or can fluctuate significantly from one analysis to another. Fortunately, by construction of the predictors, this behaviour does not have any impact on the efficiency of the VarBC method.

The impacts of assimilating bias corrected Mode-S/ADS-B wind components on AROME forecasts will be presented in a forthcoming study. Quality controls have to be made by aircraft, before considering the assimilation of its data given the known weaknesses of the VarBC method. Indeed, the method is not necessarily stable when there are only few aircraft data (few data, few analyses). A number of underlying hypotheses are not always fulfilled (constant aircraft heading and speed bias). So, over a given time period, for a specific aircraft, a set of criteria based on the number of analyses where the data is used, the total number of data, the global behaviour of the VarBC method (by examining the FGD SD values for instance), etc. has to be defined. Preliminary results tend to prove that rejected aircraft represent a small fraction of the total amount of Mode-S data. An important aspect to examine in the context of DA studies is the sampling of Mode-S data given their availability at high temporal and spatial densities, compared to AMDAR data. A number of encouraging preliminary studies have been already undertaken (data selection within boxes of various sizes both on the horizontal and on the vertical, superobbing with revised observation errors, …), but work still remains to sample en-route phase of flights, to better handle nearby Mode-S data and to specify the inertia of the bias predictor parameters to reach the best compromise between noise reduction and freedom to change.

The VarBC method is very useful and efficient on its own to bias correct data. Nevertheless it could also be used as an interesting complement to the static unbiasing method (de Haan, 2013). Indeed, the use of both methods in DA context means that if the static method is efficient for a given aircraft, the VarBC method will not add any correction on top of it. Alternatively, when a drastic bias change takes place due to a maintenance or a sensor replacement, the VarBC method will be able to rapidly correct wind biases until a revised static bias correction method accounting for this change is available.

A large number of observations from Mode-S data can be provided to convective scale NWP models, with an overall good quality when an efficient method, such as the one proposed in this study, removes most of the biases. As a consequence, it would be extremely valuable that raw data collected at national level (as done in France with 8 ADS-B antennas) could be exchanged in real-time time for operational purposes at the continental and global levels. The importance of raw data is vital for being able to set-up a VarBC scheme (definition of relevant predictors) within NWP DA systems.

Acknowledgments

The technical support from F. Guillaume (Météo-France/CNRS) has been very useful during the course of the study. Olivier Henry (Météo-France) provided the plots of magnetic declination difference (Figure 2). S. de Haan (KNMI) is acknowledged for his pionnering work on Mode-S data in the NWP context and for providing a MUAC dataset for our own evaluations. Finally, the UKMO software has been essential to compute derived winds from BDS registers collected by Météo-France antennas with the help of Olivier Traullé, Guillaume Gamelin, and Jean-Marc Lefèvre (Météo-France) to insure data production.

Notes

1 Raw data: non-meteorological quantities recorded onboard the aircraft.

2 See Section 3.1 for active data definition.

3 See Section 3.1 for passive data definition.

Appendix A:

error estimate from triple collocation

A.1 Triple collocation

The method was carried out over a one-month period (6/10/2018 to 6/11/2018). The collocation parameters in terms of altitude, horizontal distance and time where set-up to 1 hPa, 2 km and 5 min respectively. AROME 1 h forecasts data are projected to Mode-S locations.

These criteria led to a sample of 1,844 of collocated values. Even though it is possible that collocated AMDAR and Mode-S data could correspond to the same aircraft, their observation errors are assumed to be uncorrelated. Indeed such measurements, obtained from different acquisition systems, are not taken at the same place and time.

A.2 Error modelling

We consider three observing systems measuring values ${\mathit{x}}_{\mathit{i}},\mathit{ }\mathit{i}=1,2,3$ of the same quantity t (true value). Following Stoffelen (1998), each measurement ${\mathit{x}}_{\mathit{i}}$ is decomposed as follows: (9) $${\mathit{x}}_{\mathit{i}}={\mathit{a}}_{\mathit{i}}(\mathit{t}+{\mathit{\epsilon }}_{\mathit{i}})+{\mathit{b}}_{\mathit{i}}$$(9) where ${\mathit{a}}_{\mathit{i}}$ is the trend, ${\mathit{b}}_{\mathit{i}}$ the bias and ${\mathit{\epsilon }}_{\mathit{i}}$ a random error. The mean $<{\mathrm{\epsilon }}_{\mathrm{i}}>$ is supposed to be zero and the variance $<{{\mathrm{\epsilon }}_{\mathrm{i}}}^2>$ noted ${{\mathrm{\sigma }}_{\mathrm{i}}}^2\mathrm{ }$ is supposed not to depend on t.

Four different systems are examined: AMDAR (1), Mode-S (2), Mode-S without bias correction (2b) and AROME (3). Afterwards, the triple collocation involving systems 1, 2 and 3 is noted TC, whereas the one involving 1, 2b and 3 is noted TC_b.

AMDAR data are chosen to be the reference for calibration ($a_1=1$ and $b_1=0$). AMDAR and Mode-S data can represent smaller scales that are not described by the AROME model. By assuming that the true signal t is at the AROME model scales, the variances resolved by AMDAR and Mode-S data and not by AROME are integrated into ${\epsilon }_1$ and ${\epsilon }_2$ as random errors. AMDAR and scaled Mode-S data resolving the same scales can have correlated errors: $<{\epsilon }_1{\epsilon }_2>=r^2.$

We can notice than ${\epsilon }_1$ also contains collocation errors due to the fact that AROME 1 h forecasts data are projected to Mode-S locations.

For this study, we have considered AMDAR and Mode-S data that have been assimilated in the AROME model. Moreover, we assume that AROME forecast errors are uncorrelated with observation errors ($<{\epsilon }_3{\epsilon }_2>=0$ and $<{\epsilon }_1{\epsilon }_3>=0$).

With the above hypotheses, the first and second order statistical moments can be deduced using Equation (9)(9) $${\mathit{x}}_{\mathit{i}}={\mathit{a}}_{\mathit{i}}(\mathit{t}+{\mathit{\epsilon }}_{\mathit{i}})+{\mathit{b}}_{\mathit{i}}$$(9) : (10) $$a_2=\frac{c_{23}}{c_{31}}\text{ and }{a}_3=\frac{c_{23}}{c_{12}-a_2{r}^2}$$(10) (11) $$b_2=M_2-a_2{M}_1\text{ and }{b}_3=M_3-a_3{M}_1$$(11) (12) $${{\sigma }_1}^2=c_{11}-\frac{c_{31}{c}_{12}}{c_{23}}+r^2$$(12) (13) $${{\sigma }_2}^2=c_{22}/a_2^2-\frac{c_{31}{c}_{12}}{c_{23}}+r^2$$(13) (14) $${{\sigma }_3}^2=(c_{33}-a_3{c}_{31})/a_3^2$$(14) (15) $$\mathrm{where}\mathrm{ }{M}_i=<x_i>,\mathrm{ }{M}_{ij}=<x_i{x}_j>\text{ and }{c}_{ij}=M_{ij}-M_i{M}_j$$(15)

Knowing the representativeness error r, all above other parameters can be estimated from the statistics of collocated data.

Observations errors can be also estimated using the a posteriori diagnostics proposed by Desroziers et al. (2005).

These estimates ${\sigma }_{id}$ also contain a representativeness such as: ${{\sigma }_{id}}^2=a_i^2{{\sigma }_i}^2,$ i = 1,2 (16) where ${{\sigma }_i}^2$ are given from Equations (12)(12) $${{\sigma }_1}^2=c_{11}-\frac{c_{31}{c}_{12}}{c_{23}}+r^2$$(12) and (13)(13) $${{\sigma }_2}^2=c_{22}/a_2^2-\frac{c_{31}{c}_{12}}{c_{23}}+r^2$$(13) . It is therefore possible to estimate $r^2.$ The Desroziers method also provides an estimate of the forecast error σ_3d that can be compared to $a_3$σ₃ deduced from Equation (14)(14) $${{\sigma }_3}^2=(c_{33}-a_3{c}_{31})/a_3^2$$(14) .

Quality of the collocation can be estimated using differences of AROME data at each AMDAR (a) and Mode-S (b) collocated different points. $$<{(x_{3a}-x_{3b})}^2>={a_3}^2(<{\left(t_a-t_b\right)}^2>+2{{\sigma }_3}^2-2<({\epsilon }_{3a}{\epsilon }_{3b})>).$$

As (a) and (b) are close, we suppose ${\epsilon }_{3a}$ and ${\epsilon }_{3b}$ very correlated. We have then: (17) $$\begin{array}{cc}{{\sigma }_c}^2\mathrm{ }\le \mathrm{ }<{(x_{3a}-x_{3b})}^2>/{a_3}^2& \mathrm{with}\mathrm{ }{{\sigma }_c}^2=<{(t_a-t_b)}^2>\end{array}$$(17)

A.3 Results

A.3.1 Global statistics

Figure A1 shows that the VarBC method allows Mode-S data to be more consistent with AMDAR data by reducing significantly the differences in terms of standard deviation values for u (resp. v) from 3.79 m/s (resp. 3.13 m/s) to 1.21 m/s (resp. 1.25 m/s). In addition, Mode-S VarBC and AMDAR data have very close statistics when compared to AROME.

Moreover, Mode-S VarBC data are closer to AMDAR data than to AROME data. AMDAR data selected for the triple collocation have been used to anchor the VarBC method. However, due to their reduced number, they represent only a fraction of the total observations (other distant aircraft, AMV, soundings …) that steered the bias correction for these aircraft.

A.3.2 Error estimations

To avoid problems caused by outliers (order 2 moment sensitivity) that can be seen in Figure A1, they have been filtered using an iterative scheme. At each step all collocated triplets are tested and triplet j $(\overline{{\mathit{x}}_1^{\mathrm{j}}},\overline{{\mathit{x}}_2^{\mathrm{j}}},\overline{{\mathit{x}}_3^{\mathrm{j}}})$ is rejected as an outlier using the ‘4 σ’ rule if ${\left(\overline{{\mathit{x}}_1^{\mathrm{j}}}-\overline{{\mathit{x}}_3^{\mathrm{j}}}\right)}^2\mathrm{ }\ge 16\mathrm{ }({\overline{{\sigma }_1}}^2+{\overline{{\sigma }_3}}^2)$ or ${\left(\overline{{\mathit{x}}_2^{\mathrm{j}}}-\overline{{\mathit{x}}_3^{\mathrm{j}}}\right)}^2\mathrm{ }\ge 16\mathrm{ }({\overline{{\sigma }_2}}^2+{\overline{{\sigma }_3}}^2)$ or ${\left(\overline{{\mathit{x}}_1^{\mathrm{j}}}-\overline{{\mathit{x}}_2^{\mathrm{j}}}\right)}^2\mathrm{ }\ge 16\mathrm{ }({\overline{{\sigma }_1}}^2+{\overline{{\sigma }_2}}^2-2r^2)$ where ${\overline{x}}_{\mathrm{i}}=\frac{x_i- b_i}{a_i}$ (99.99% of data kept in a normal distribution). $\overline{{\sigma }_i}$ is defined in the same way as ${\mathrm{\sigma }}_{\mathrm{i}}$ but for ${\overline{x}}_{\mathrm{i}}$ and is calculated with Equations (12)–(15) applied to ${\overline{x}}_{\mathrm{i}}.$ Then ${\mathrm{x}}_{\mathrm{i}}$ model error parameters are calculated with remaining triplets before going to the next step.

After filtering outliers for u and v parameters, it remains 1,825, resp. 1,830 triplets in TC, and 1,832, resp. 1,836 triplets in TC_b.

Table A1 shows the different possibilities of estimating $r^2$ from Equations (12)(12) $${{\sigma }_1}^2=c_{11}-\frac{c_{31}{c}_{12}}{c_{23}}+r^2$$(12) , (13)(13) $${{\sigma }_2}^2=c_{22}/a_2^2-\frac{c_{31}{c}_{12}}{c_{23}}+r^2$$(13) , and (16) from the filtered triple collocation TC for u and v.

$r^2$ is not so dependent on the parameter or on the data sample, thus we assume $r^2$ = 0.827 m²/s².

Table A2 shows that the trends are generally close to one and biases close to zero, except for the unbiased Mode-S data (residual biases between 0.1 and 0.3 m/s). AROME is not biased. The results for u and v are very similar. Scaled forecast errors estimated by the Desroziers method are close to those estimated with triple collocations since the same scales are treated in both cases, with values around 1.7–1.8 m/s. Observation errors for AMDAR and Mode-S VarBC data are close to each other: 1.3 m/s respectively 1.1 m/s.

The small difference in favour of Mode-S data can partly be explained by their number. Indeed, the same thinning is applied in the analysis to both data sources (1 aircraft data per 25 km box size, per model level, per analysis, by choosing the time closest to the analysis date). Since there are many more Mode-S data per analysis than AMDAR data, those selected for assimilation will be on average closer to the analysis date and potentially of better quality (by an improved time, and therefore position consistency). Moreover, the temporal accuracy of AMDAR and Mode-S data acquisition is 1 min and 1 s, respectively. This leads to a less accurate collocation of AMDAR data inducing an additional positioning error.

Finally, we can estimate the quality of the collocation using AROME data at AMDAR (a) and Mode-S (b) collocations from Equation (17)(17) $$\begin{array}{cc}{{\sigma }_c}^2\mathrm{ }\le \mathrm{ }<{(x_{3a}-x_{3b})}^2>/{a_3}^2& \mathrm{with}\mathrm{ }{{\sigma }_c}^2=<{(t_a-t_b)}^2>\end{array}$$(17) , and for both u and v parameters: $${\sigma }_c\le 0.165\mathrm{ }m/s \text{from }u \text{and }{\sigma }_c\le 0.152\mathrm{ }m/s\text{ from }v.$$

Despite the very different behaviour of Mode-S and Mode-S VarBC statistics, the error model parameters of AMDAR and AROME data are rather similar using either TC or TC_b triple collocation (a difference of about 0.1–0.2 m/s).

Annex 1:

Mode-S/ADSB- first level of quality controls

These quality controls are performed with the UK MetOffice wind computation software (flwindcalc) used at Météo-France to process data from the network of ADS-B antennas

Atmospheric winds are computed when the following criteria are satisfied on raw data:

• True airspeed (TAS): 75 kt < TAS < 550 kt

• Indicated airspeed (IAS): 100 kt < IAS <550 kt

• Difference between the calculated and reported ground track angle less than 5°

• Difference between the calculated and reported ground speed less than10 kt

• Roll angle less than 2°

• Time between all four required BDS registers (BDS 4.0, BDS 5.0, BDS 6.0 and BDS 0.5): 10 s

• Difference between the last and current reported heading less than 20°

• Difference between the aircraft heading and reported ground track angle less than 20°

Fig. A1. Boxplots of aircraft u (blue) and v (red) for Mode-S minus AMDAR (first from the left), Mode-S VarBC minus AMDAR (second), Mode-S minus AROME (third), Mode-S VarBC minus AROME (fourth) and AMDAR minus AROME (right) collocated datasets from 6 October to 6 November 2018. Star: standard deviation, square: mean.

Table A1. Representativeness error $r^2$ estimation.

${{\sigma }_{id}}^2$ / $r^2$	AMDAR (Equations (12)(12) $${{\sigma }_1}^2=c_{11}-\frac{c_{31}{c}_{12}}{c_{23}}+r^2$$(12) and (16))	Mode-S VarBC (Equations (13)(13) $${{\sigma }_2}^2=c_{22}/a_2^2-\frac{c_{31}{c}_{12}}{c_{23}}+r^2$$(13) and (16))
u (m^2/s^2)	1.586/0.769	1.234/0.785
v (m^2/s^2)	1.778/0.851	1.406/0.901

Table A2. Summary of the parameters of the error model decomposition (Equation (9)(9) $${\mathit{x}}_{\mathit{i}}={\mathit{a}}_{\mathit{i}}(\mathit{t}+{\mathit{\epsilon }}_{\mathit{i}})+{\mathit{b}}_{\mathit{i}}$$(9) ) derived from the TC and TC_b collocations of AMDAR, Mode-S and AROME 1 h forecasts data for u (top) and v (bottom) wind components.

u collocation TC/TCb	ai	bi	${\sigma }_i$	${\sigma }_{id}$/ ai
TC: AMDAR	1	0	1.283	1.259
TC: Mode-S VarBC	1.011	0.082	1.117	1.098
TC: AROME1	0.987	−0.029	1.685	1.732
TCb: AMDAR	1	0	1.455	1.265
TCb: Mode-S no BC	1.02	0.258	3.645	NA
TCb: AROME1	0.989	−0.037	1.553	1.771
v collocation TC/TCb	ai	bi	${\sigma }_i$	${\sigma }_{id}$/ ai
TC: AMDAR	1	0	1.325	1.334
TC: Mode-S VarBC	1.006	−0.107	1.147	1.178
TC: AROME1	0.985	−0.042	1.788	1.803
TCb: AMDAR	1	0	1.46	1.347
TCb: Mode-S no BC	1.009	0.109	3	NA
TCb: AROME1	0.986	−0.051	1.706	1.812