Comparison between the assimilation of IASI Level 2 ozone retrievals and Level 1 radiances in a chemical transport model

The prior information used for Level 2 (L2) retrievals in the thermal infrared can influence the quality of the retrievals themselves and, therefore, their further assimilation in atmospheric composition models. In this study we evaluate the differences between assimilating L2 ozone profiles and Level 1 (L1) radiances from the Infrared Atmospheric Sounding Interferometer (IASI). We minimized potential differences between the two approaches by employing the same radiative transfer code (Radiative Transfer for TOVS, RTTOV) and a very similar setup for both the L2 retrievals (1D-Var) and the L1 assimilation (3D-Var). We computed hourly 3D-Var analyses assimilating L1 and L2 data in the chemical transport model MOCAGE and compared the resulting O₃ fields among each other and against ozonesondes. We also evaluated the joint assimilation of limb measurements from the Microwave Limb Sounder (MLS) in combination with IASI to assess the impact of stratospheric O₃ on tropospheric analyses. Results indicate that significant differences can arise between L2 and L1 assimilation, especially in regions where the L2 prior information is strongly biased (at low latitudes in this study). In these regions the L1 assimilation provides a better variability of the free-troposphere ozone column. L1 and L2 assimilation instead give very similar results at high latitudes, especially when MLS measurements are used to constrain the stratospheric O₃ column. A critical analysis of the potential benefits and drawbacks of L1 assimilation is given in the conclusions. We also list remaining issues that are common to both the L1 and L2 approaches and that deserve further research.

1 Introduction

The global monitoring of atmospheric composition relies on a large number of dedicated satellite missions and on the sustained improvement of numerical forecast models. Today, research and operational centres provide both satellite-based reanalyses and forecasts of atmospheric composition for a large number of applications, spanning from stratospheric ozone monitoring (van der A et al., 2010) to climate change (Flemming et al., 2017) and air quality (Zhang et al., 2012; Marécal et al., 2015).

Satellite sensors measure the spectral signature of gases and aerosols on the radiation field that traverse the atmosphere. Retrieving the concentration of a given gas from the radiation measured at the satellite position represents an inverse problem that is in most cases ill-posed and underdetermined; i.e. finding the solution requires some type of mathematical regularization or prior information (Rodgers, 2000). The accuracy of the solution depends in general on the intensity of the spectral signature of the retrieved compound, the source of radiation (e.g. the Earth or the Sun), the observation geometry and the accuracy of the radiative transfer model (RTM). The last-mentioned factor also means correctly accounting for all the atmospheric constituents or surface properties that affect the radiation field but are not retrieved themselves (auxiliary RTM inputs).

When the retrieval is done within a Bayesian framework, like the optimal estimation method (Rodgers, 2000), the measurement errors, the RTM errors and the uncertainty in the prior information (also named background or a priori profile) are prescribed. The procedure then provides an estimation of the error covariance for the retrieved quantity and the averaging kernels (AKs), which quantify the sensitivity of the retrieval to the true state and are linked to the degrees of freedom (DOF) of the solution. The retrieval errors and the AKs (or the DOF) can be used first to diagnose the quality and the relevance of the atmospheric retrieval. They become even more important when retrievals are further assimilated in numerical forecast models because they weight the impact of the observations in the system.

Chemical transport models (CTMs) solve the chemical and physical processes within the atmosphere but are based on meteorological fields from a numerical weather prediction (NWP) model to advect the chemical species. Coupled chemistry–meteorology models (CCMMs) that simulate both meteorology and chemistry online became available later but are quite common today in operational centres (Zhang et al., 2008; Flemming et al., 2015). There are currently growing efforts to introduce even stronger coupling of the atmosphere with both ocean and surface models, which has given rise to so-called Earth system models (ESMs; Brown et al., 2012; Hurrell et al., 2013). ESMs provide a comprehensive tool for climate predictions and reanalyses, but they are also considered for state-of-the-art air-quality modelling (Neal et al., 2017).

Closely following the historical advances in modelling, the assimilation of satellite data was first introduced in CTMs (Geer et al., 2006; Lahoz et al., 2007), and it is now also well integrated in operational CCMMs (Flemming et al., 2017).

Today, numerous satellite retrievals of trace gases (e.g. O₃, CO, NO₂, CH₄, CO₂) and aerosols (aerosol optical depth, AOD) are assimilated daily within operational CTMs and CCMMs (Inness et al., 2015; Bocquet et al., 2015).

For a long time, meteorological variables such as temperature and water vapour profiles have been corrected by means of assimilating satellite radiances (Level 1 data) directly in NWP models. Therefore, the RTM became part of the observation operator of the assimilation system (Andersson et al., 1994). This avoided the introduction of biases in NWP that arose from poor prior information used in satellite retrievals at that time and neglect of the AKs (Eyre et al., 1993). On the other hand, chemical species and aerosols are mostly corrected by means of assimilating geophysical retrievals (Level 2 or L2 data) that are made available by satellite data providers. To remove the impact of the prior information when assimilating L2 retrievals, the AK of the retrieval must be multiplied by the modelled profiles before computing the innovation vectors (Rodgers, 2000; Eskes and Boersma, 2003; Migliorini, 2012). However, within standard methods based on the linearization of the RTM, like optimal estimation, issues might still arise when the prior information used in the retrieval is far from the true atmospheric state: this might challenge the linearization of the observation operator and result in sub-optimal retrievals. Since the AKs themselves are also a result of the retrieval (and depend upon its prior information), we expect that a perfect removal of the prior information within data assimilation (DA) cannot always be ensured.

The precise conditions that provide an equivalence between assimilating retrievals (using some kind of weighting function) and radiances have been formalized by Migliorini (2012) and further tested by Prates et al. (2016) on synthetic satellite observations. These authors conclude that the equivalence holds under the hypothesis of an almost linear radiative transfer (RT) regime and with careful selection of the prior error covariances in order to maximize the measurement information in the retrieval step. Nonetheless, testing the two approaches within an operational system and with real observations remains crucial to verify whether these conditions are met in practice. Moreover, the perfect equivalence only holds when all the auxiliary inputs of the RTM are exactly the same in both the retrieval and the radiance assimilation. It is clear that a climatological option for some RTM inputs will always be a more practical choice when computing L2 retrievals. On the other hand, the evolution towards strongly integrated ESMs will allow in principle to dispose of the most accurate prior information for all RTM inputs and favours the radiance assimilation approach. In this context, it appears important to introduce and evaluate the assimilation of radiances for chemical applications as well.

To the best of the authors' knowledge, the existent literature on this topic only concerns meteorological applications. Han and McNally (2010) explored the possibility of assimilating O₃ sensitive radiances within a NWP model but without comparing the two approaches. Similarly, Weaver et al. (2007) examined the assimilation of satellite radiances for aerosols, but the focus was on the impact of using modelled aerosol microphysical properties as auxiliary input for the RTM, and no comparison was provided. No other studies could be found concerning the assimilation of chemical compounds.

The objective of this study is to perform a first strict comparison between the assimilation of radiances and retrievals, with respect to O₃ estimation in the thermal infrared (TIR). To this end, systematic differences between the assimilation of retrieval and radiances have been minimized as much as possible, for example by means of employing the same RTM within the two approaches.

We consider the case of O₃ assimilation using the Infrared Atmospheric Sounding Interferometer (IASI) on board the European MetOp satellites (Clerbaux et al., 2009). Several IASI O₃ retrievals have already been well validated (Dufour et al., 2012) and used directly to provide multi-annual time series of the global O₃ budget (Wespes et al., 2016) or successfully assimilated within global (Peiro et al., 2018) and regional CTMs (Coman et al., 2012). However, an empirical correction of the retrievals has been found to be necessary to ensure globally unbiased reanalyses, and slightly degraded assimilation results are still found at mid-latitudes and high latitudes (Emili et al., 2014). Since the tropospheric O₃ signature in the selected IASI spectral window decreases over colder surfaces, the impact of the retrieval's prior information might become more relevant at high latitudes. In addition, the majority of IASI O₃ retrievals use a single a priori profile globally (Barret et al., 2011; Boynard et al., 2016), which might present very large local departures from the true O₃ profile. Hence, IASI O₃ assimilation represents a good benchmark to evaluate the differences between the assimilation of retrievals and radiances.

The IASI SOFRID (Software for a Fast Retrieval of IASI Data) O₃ product (Barret et al., 2011) and MOCAGE CTM have been used here to benefit from the experience of previous studies (Emili et al., 2014; Peiro et al., 2018). SOFRID and MOCAGE DA are based on a variational algorithm, and, since SOFRID employs RTTOV (Saunders et al., 1999), which is a community RTM developed originally for NWP applications, the same RTM has been implemented in the MOCAGE system. Global O₃ analyses are computed for July 2010, and the results are compared against all available radiosoundings to evaluate their accuracy. Since the sensitivity of IASI TIR measurements to O₃ is not uniform along the atmospheric column, we also investigate the impact of assimilating more accurate stratospheric profiles from the Microwave Limb Sounder (MLS) in combination with IASI radiances. This might reveal possible synergies when assimilating multiple instruments that sense different layers of the atmosphere.

The paper is organized as follows. The satellite measurements, the Level 2 retrievals and the validation measurements used for this study are described in Sect. 2, as well as main steps concerning the preprocessing for some of the datasets. The chemical transport model, the radiative transfer model, the assimilation algorithm and the setup of the experiments are described in Sect. 3. The assimilation of IASI retrievals and radiances is compared in Sect. 4.1, and the impact of MLS assimilation in combination with IASI is discussed in Sect. 4.3. The conclusions are summarized in the last section, where some recommendations are also given.

2 Measurements

2.1 IASI

IASI flies on board the series of polar-orbiting satellites MetOp operated by the EUropean organization for the exploitation of METeorological SATellites (EUMETSAT). It provides hyper-spectral measurements of the Earth's thermal radiation in the 3.62–15.5 µm (2760–645 cm⁻¹) window and serves meteorological and atmospheric chemistry applications (Clerbaux et al., 2009; Hilton et al., 2012). IASI is an operational mission meant to provide long-term (>20 years) time series of accurate TIR spectra at high spatial resolution. A total of three IASI instruments will be flying simultaneously at the end of 2019, providing nearly global coverage six times per day (morning and evening overpasses). Hence, they represent a great opportunity for both NWP and climate–chemistry reanalyses. Only MetOp-A data, available from 2008 to present, have been employed for this study.

2.1.1 L1 radiances

IASI L1c data contain calibrated and geolocalized spectra at 0.5 cm⁻¹ spectral resolution (after apodization), i.e. 8461 radiance values for each ground pixel, with a footprint of 12 km for nadir observations. For this study, historical L1c data granules have been downloaded from the EUMETSAT Earth Observation data portal (https://eoportal.eumetsat.int, last access: 16 July 2019) in NETCDF format. Data files also contain the observation geometry (Sun and satellite angles) for each ground pixel and the co-located land mask and cloud fraction values, obtained from the Advanced Very High Resolution Radiometer (AVHRR) measurements, also on board MetOp.

Table 1 Number of validation profiles for July 2010.

2.1.2 SOFRID L2 retrievals

The Software for a Fast Retrieval of IASI Data (SOFRID) was developed at the Laboratoire d'Aérologie to retrieve O₃ (Barret et al., 2011) and CO (De Wachter et al., 2012) profiles from IASI. It is based on the Radiative Transfer for TOVS (RTTOV) RTM (Saunders et al., 1999) and the 1D-VAR scheme developed within the Numerical Weather Prediction Satellite Application Facilities (NWP SAF) programme. SOFRID retrieves the O₃ profile in volume mixing ratio (vmr) units at 43 pressure levels between the surface and 0.1 hPa using 469 spectral channels within the main IASI O₃ window (980–1100 cm⁻¹). The choices that are made in SOFRID and are relevant for this study are summarized in Table 2. Note that a single a priori profile and error covariance matrix are used globally and that the surface skin temperature (SST) is estimated within the retrieval.

The number of DOF of the SOFRID retrieval has been evaluated to be between 2 and 3 for the full atmospheric column, with about 1 DOF for the tropospheric column (Dufour et al., 2012). SOFRID's averaging kernels corresponding to the retrievals used within this study (see Sect. 2.4) are shown in Fig. 1. We remark that the largest sensitivities are found for the lower stratosphere (50–100 hPa) and upper troposphere (200–300 hPa) levels. The sensitivity in the free troposphere (400–600 hPa) is maximum at tropical latitudes and decreases towards the poles due to the decreasing thermal contrast. Very low sensitivities are in general found for levels below 700 hPa at all latitudes.

(Borbas and Ruston, 2010)(Borbas and Ruston, 2010)

Table 2 Summary of the configuration of SOFRID L2 retrievals and MOCAGE L1 assimilation.

Figure 1 SOFRID O₃ averaging kernels for the month of July 2010 averaged globally (first plot) and for five separate latitude bands (90–60^∘ S, 60–30^∘ S, 30^∘ S–30^∘ N, 30–60^∘ N, 60–90^∘ N). Each coloured line corresponds to a retrieval's level, and the corresponding pressure is indicated in the colour bar. Only SOFRID levels with a pressure >50 hPa are displayed for better clarity.

The accuracy of the retrieved O₃ depends on the latitude and the vertical level but is generally within 10 %–20 % of the corresponding radiosounding values, once the averaging kernels are applied. However, biases are found in the troposphere with SOFRID (+10 %), and positive biases of about 15 % are found in the upper troposphere–lower stratosphere (UTLS) region with all current IASI O₃ products (Dufour et al., 2012). The reasons for such biases are not yet fully understood and can impact data assimilation (Emili et al., 2014) or trend analysis (Gaudel et al., 2018) negatively. This study will provide further insights about the impact of the constant a priori profile on IASI retrievals.

The SOFRID V1.5 retrievals described in Barret et al. (2011) are available for the full MetOp-A period at http://thredds.sedoo.fr/iasi-sofrid-o3-co (last access: 16 July 2019). The V3.0 version of SOFRID retrievals has been used for this study and was obtained from Brice Barret (personal communication, 2019). The main difference to version 1.5 concerns the temperature and water vapour profiles employed in the radiative transfer computations, which are taken from the ECMWF NWP model instead of EUMETSAT L2 retrievals. The SOFRID v3.0 preprocessor retrieves the operational analysis (type “an”) at 00:00, 06:00, 12:00, and 18:00 UTC from the ECMWF NWP model and assimilation system (Integrated Forecast System, IFS), regridded to a resolution of $\mathrm{0.25}{}^{\circ }\times \mathrm{0.25}{}^{\circ }$. All the fields are then interpolated at the closest hour to the IASI pixel, and a nearest-neighbour interpolation is done to extract the corresponding profiles and surface properties. Since the CTM is also based on ECMWF NWP forcing fields (Sect. 3.1), this choice minimizes possible systematic differences between L2 retrievals and L1 assimilation. Also, SOFRID v3.0 is based on a more recent version of RTTOV and newer IASI coefficients (v11.1, coefficients on 101 levels) than the original L2 product (v9.0, coefficients on 43 levels). In addition to the O₃ retrieval and its error covariance, SOFRID files contain a number of auxiliary and diagnostic fields. The cloud fraction is based on a combination of the EUMETSAT L2 product (AVHRR) and a brightness temperature (BT) analysis at 11 and 12 µm to fill pixels with missing AVHRR data (Barret et al., 2011). An index based on the V-shaped sand signature computed as $\mathrm{\Delta }\mathrm{BT}=({\mathrm{BT}}_{\mathrm{829}{\mathrm{cm}}^{-\mathrm{1}}}-{\mathrm{BT}}_{\mathrm{972.5}{\mathrm{cm}}^{-\mathrm{1}}})+({\mathrm{BT}}_{\mathrm{1202.5}{\mathrm{cm}}^{-\mathrm{1}}}-{\mathrm{BT}}_{\mathrm{1096}{\mathrm{cm}}^{-\mathrm{1}}})$ is used to detect pixels affected by large aerosol load. Usage of these products will be detailed in Sect. 2.4, “Data preprocessing”.

Figure 2 Total number of IASI observations per model grid box ($\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$) retained after the selection procedure described in Sect. 2.4 and further assimilated in this study for the month of July 2010. The total number of assimilated observation for the entire globe (N) is given above the map.

2.2 MLS L2 retrievals

Since 2004 the Microwave Limb Sounder (MLS) has been flying on board the research mission AURA and measures thermal emission at the atmospheric limb (Waters et al., 2006). It provides about 3500 stratospheric profiles of multiple atmospheric constituents each day, including O₃ (Froidevaux et al., 2008). Since version 3 of MLS products has been in operation, O₃ profiles have been retrieved on 55 pressure levels, with a recommended range for scientific usage between 0.02 and 261 hPa for version 4.2 (Livesey, 2018). The biases of MLS O₃ profiles are typically within 5 % with respect to ozonesondes and lidar measurements (Hubert et al., 2016), with slightly higher values below 200 hPa. Given its good accuracy, MLS O₃ has been widely used both for trend analysis (Froidevaux et al., 2015) and assimilation experiments (Massart et al., 2010; Miyazaki et al., 2012; Inness et al., 2015). Similarly to previous studies (Emili et al., 2014), we retain only the most accurate data using MLS, i.e. above 170 hPa. The MLS V4.2 product used in this study has been downloaded from the Goddard Earth Sciences Data and Information Services Center (GES DISC) web portal (https://disc.gsfc.nasa.gov, last access: 16 July 2019).

2.3 Radiosoundings

Ozonesondes are launched on a weekly basis by meteorological services and provide accurate profiles of O₃ up to 10 hPa with a vertical resolution of 150–200 m. Electrochemical concentration cell (ECC) sondes, which represent the largest percentage of the global network, have a precision of about 5 % (Thompson et al., 2003). Radiosoundings are relatively sparse, and their geographical distribution is much more representative of the northern mid-latitudes. However, for several decades they have provided the most precise information on vertical ozone distribution in the troposphere. Therefore, they have been used to derive widely used tropospheric O₃ climatologies (McPeters et al., 2007) and validate both satellite products (Dufour et al., 2012) or models (Geer et al., 2006). They will be used in this study to validate all model simulations. Data are collected and distributed by the World Ozone and Ultraviolet Radiation Data Center (WOUDC; http://www.woudc.org, last access: 16 July 2019).

2.4 Data preprocessing

Some further preprocessing has been applied to the original L1c and SOFRID datasets to ease the interpretation of the assimilation experiments presented later in Sect. 4. The objective was to ensure that exactly the same spectra are used for both L1 and L2 assimilation.

Only the spectral channels that are used in SOFRID are extracted from IASI L1c granules, i.e. channel no. 1350 (980 cm⁻¹) to 1818 (1100 cm⁻¹). Some further screening is applied to remove channels that are affected by strong H₂O absorption, as also done in SOFRID.

The spatial resolution of the CTM ($\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$; Sect. 3.1) is much coarser than IASI pixel size. Since it is preferable to avoid all kind of spatial averaging of the observations, a significant reduction of ground pixels is needed. In return, we employ strict selection criteria to avoid contamination from clouds and bright surfaces as much as possible, which reduces the RT accuracy and increases retrieval or assimilation errors. The data selection is performed as follows.

First, only L1 pixels with both IASI and AVHRR highest-quality flags are kept.

Then, ground pixels from IASI L1 and SOFRID products are filtered using their respective cloud masks (Sect. 2.1.1 and 2.1.2) and only keeping pixels with a cloud fraction less than or equal to 1 %.

SOFRID pixels with a sand signature greater than 0.5 and with a number of retrieved levels lower than 35 (mountains) are also filtered out.

Resulting datasets are then matched; i.e. only common ground pixels that remained available after the previous L1 and SOFRID independent selections are kept. Finally, data thinning is performed using a regular grid of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ resolution and only keeping the first pixel that falls in every two grid boxes. This ensures a minimum distance of 1^∘ among assimilated observations. After the completion of the data selection procedure the final number of retained ground pixels for L1 and SOFRID is about 5000 per day, compared to about 10⁵ when only the cloud screening is applied. The total number of L1 and L2 observations resulting from the above selection and further assimilated in this study is displayed in Fig. 2.

3 Method

This section summarizes the main characteristics of the CTM (Sect. 3.1), the RTM (Sect. 3.2) and the assimilation algorithm (Sect. 3.3) used in this study. Further details on the particular selection of the main parameters of the assimilation experiments (e.g. the error covariances) are given in Sect. 3.4.

3.1 Chemical transport model

The chemical transport model (CTM) MOCAGE (Josse et al., 2004) is used in this study. A global configuration with a horizontal resolution of $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ and 60 hybrid sigma-pressure levels up to 0.1 hPa has been used. The vertical resolution varies from about 100 m in the planetary boundary layer to about 700 m in the upper troposphere, decreasing further to approximately 2 km in the upper stratosphere. Chemical mechanism, emissions and physical parameterizations follow the setup used for operational air-quality forecasts (Marécal et al., 2015), which includes about 100 species and 300 chemical reactions. A similar configuration has been employed by Barré et al. (2013) to assimilate IASI O₃ columns over Europe but with a lower model top at 5 hPa. Other authors favoured a simplified chemistry scheme but with a model top at 0.1 hPa to assimilate satellite O₃ products globally (Emili et al., 2014; Peiro et al., 2018).

For this study we considered the highest available model top because we need to simulate the full atmosphere to compute radiances. In addition, the 0.1 hPa top matches with the vertical grid used for SOFRID retrievals (Sect. 2.1.2), making the comparison of the two assimilation approaches (radiances versus L2) stricter. The full chemical scheme is chosen instead of a simplified chemistry to reduce biases of the modelled O₃ in the troposphere as much as possible. The main intent of this study is in fact to evaluate the impact of dynamical and accurate O₃ prior information on assimilation results.

The meteorological forcing comes from the ECMWF IFS, from which we retrieved the forecast (type “fc”) initialized with the analysis at 00:00 UTC each day. The NWP fields are interpolated on the horizontal grid of MOCAGE ($\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$) during the retrieval process and stored with a time step of 3 h. During the integration of MOCAGE, the meteorological forcing is linearly interpolated at the advection time step of MOCAGE (hourly) and on the CTM's vertical grid.

3.2 Radiative transfer model

RTTOV (Saunders et al., 1999) is a community RTM developed for operational NWP models. One of its main advantages is computational efficiency, which is achieved by running accurate but costly line-by-line RT simulations for a large number of satellite sensors, observation geometries and atmospheres and storing the corresponding coefficients in large lookup tables. RTTOV provides application programming interfaces (APIs) for the direct RT computations plus the tangent linear and adjoint model, which are needed in variational assimilation systems.

Version 11.3 of RTTOV (Saunders et al., 2013) has been used in this study for the L1 assimilation. This version includes coefficients for the IASI TIR channels computed using a fine atmospheric grid (101 vertical levels). The SST, 2 m temperature, 2 m pressure and 2 m wind vector are taken from high-resolution ($\mathrm{0.125}{}^{\circ }\times \mathrm{0.125}{}^{\circ }$) global IFS forecasts initialized with the analysis at 00:00 UTC each day and co-located (nearest neighbour) with satellite ground pixels prior to data assimilation. A linear interpolation from the 3 h forecast steps to the closest hour of the IASI observations is also performed for these fields. The surface emissivity is based on the RTTOV monthly TIR emissivity atlas (Borbas and Ruston, 2010). Only clear-sky RT computations are performed for this study, and no aerosols have been prescribed. The RTM configuration is summarized in Table 2. Due to the different processing chains, the auxiliary inputs of RTTOV could not be set exactly equal for L1 assimilation and L2 retrievals (Table 2). The potential impact of these residual differences is discussed in Sect. 4.1.

3.3 Assimilation algorithm

The assimilation suite for MOCAGE is based on a variational algorithm and was developed initially within the ASSET (Assimilation of Envisat data) project (Lahoz et al., 2007). The objective was to assimilate satellite products at a global scale, and a 3-D-FGAT implementation was chosen. It evolved later to provide air-quality reanalyses at the surface based on a 3D-Var implementation (Jaumouillé et al., 2012) and was extended to 4D-Var when employing linearized chemistry schemes (Massart et al., 2012; Emili et al., 2014). In all cases the minimization of the variational cost function is performed using the limited-memory Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm (Liu and Nocedal, 1989). In this study we used a 3D-Var algorithm with hourly assimilation windows and with O₃ as the control variable.

The 3-D background error covariance is modelled through a diffusion operator (Weaver and Courtier, 2001) and allows the specification of heterogeneous correlation length scales. Compared to previous studies using the MOCAGE assimilation suite, a new vertical correlation operator has been employed here: the vertical error correlation is now assigned by explicitly filling a positive definite matrix using the Gaussian formulation of Paciorek and Schervish (2006) and by numerically computing its square root. This avoids difficulties encountered with diffusion-based operators concerning the normalization in the presence of boundaries (e.g. the surface) and heterogeneity (Mirouze and Weaver, 2010). Since the vertical dimension of the model grid is relatively small, this choice does not impact the numerical cost and the memory requirements significantly with respect to the previous implementation based on diffusion.

The observation operator of MOCAGE allows a large number of measurements to be assimilated, spanning from columns of gases (Massart et al., 2009) to aerosol optical depth (Sič et al., 2016). Next, we give some details of the implementation used in this study to assimilate vertical profiles and radiances.

After the horizontal and temporal interpolation of the model fields at the satellite ground-pixel position, modelled profiles are linearly interpolated to the retrieval's vertical grid. When the averaging kernels are used (i.e. for SOFRID assimilation), the linear estimation equation (Barret et al., 2011) is used to remove the impact of the prior information from the innovation vector. The ensemble of these operations is stored as coefficients of a large sparse matrix and done through its multiplication by the model 3-D field. This approach is practical since numerous applications of the linearized and adjoint operator are needed during the minimization of the variational cost function. Differently from all previous studies involving IASI O₃ assimilation (Massart et al., 2009; Emili et al., 2014; Peiro et al., 2018), in which L2 profiles were first reduced to total or partial columns prior to assimilation, here we assimilate the full L2 profiles directly (43 levels). This avoids any loss of information and allows a fairer comparison between L2 and radiance assimilation. The error covariance matrix of the profile-type observations is diagonal in the latitude–longitude dimensions, but off-diagonal terms are allowed along the vertical dimension. Alternative approaches exist to optimally reduce the dimension of the L2 observation space based on the DOF of the retrievals (Migliorini et al., 2008; Mizzi et al., 2016), which are of interest to further reduce the numerical cost of SOFRID assimilation without loss of accuracy. However, this is left for future work.

To compute modelled radiances we employ the same horizontal and temporal interpolation as in the case of profile observations, except for the vertical interpolation. In fact, the RTTOV vertical interpolator is used for radiance computations instead of the MOCAGE one. All model levels (60) and corresponding pressure levels are given as input to RTTOV, which performs the vertical interpolation to the IASI coefficient levels internally. Since the model vertical resolution is lower than the one available in RTTOV for IASI coefficients (101 levels), we used the default option based on Rochon et al. (2007). Also, O₃ profiles above the CTM top (0.1 hPa) are completed using RTTOV climatological profiles. Auxiliary inputs for the radiance computation include the pressure, temperature and water vapour profiles, which are interpolated from the corresponding MOCAGE fields.

The MOCAGE control vector has been extended to include the SST, as in the SOFRID retrieval scheme. This proved to be important since small errors in the SST translate in significant differences between modelled and measured radiances. Not accounting for this would produce incorrect O₃ analyses. The SST does not belong to the MOCAGE prognostic fields, nor is it prescribed on the MOCAGE grid. Hence, the SST analysis is not propagated in time, and no spatial covariance model has been implemented so far. In this regard the treatment of the SST is equivalent to that done in L2 retrievals. In the context of MOCAGE 3D-Var, it can be interpreted as a variational bias correction term in the observation space (Dee and Uppala, 2009), with prior values given by the NWP model (IFS; see Sect. 3.2).

3.4 Setup of the experiments

We performed numerical experiments for the month of July 2010, which corresponds to the typical presence of summer O₃ maxima in the Northern Hemisphere linked to photochemical pollution. July 2010 is also interesting due to the development of a strong La Niña episode (Peiro et al., 2018). The main difference between assimilating L2 and L1 data consists in using a climatological (L2 assimilation) versus a dynamical a priori profile (L1 assimilation) for the inversion of the radiative transfer problem. The chosen period presents large local deviations of the O₃ field from climatological values. Therefore, it provides an interesting benchmark period with respect to the objective of this study.

The CTM was initialized on 1 June 2010 with a zonal climatology and run for a 1-month period (spin-up) to provide chemically balanced initial conditions on 1 July 2010 for all simulations.

The observation error covariance matrix (R) is prescribed according to the choices adopted in SOFRID V3.0. When the radiances are assimilated, a diagonal matrix (i.e. with no inter-channel correlation) is used, with a constant standard deviation of to 0.7 mW m⁻² sr⁻¹ cm for all channels. This is a simplified although common setting for most IASI O₃ retrievals (Barret et al., 2011; Boynard et al., 2016). The SST, which is controlled as well within radiance assimilation, has a prescribed standard deviation of 4 ^∘C for all ground pixels. When L2 profiles are assimilated we used the full non-diagonal error covariance matrix provided by SOFRID or MLS retrievals.

We considered a dynamical rejection of observations based on the relative differences between simulated and measured values with respect to simulated values. It avoids assimilating observations with departures from the corresponding model background that are too large. The threshold values are set to 12 % for L1 radiances and 2000 % for L2 profiles, and trespassing the threshold for any particular channel or profile level rejects the entire spectrum or profile. The strong difference between the two thresholds is a consequence of the very different nature of assimilated observations: the exponential shape of O₃ profiles can produce very large departures when the gradient is the steepest (tropopause), and a small rejection threshold would filter out most of the profile observations. This is not the case for radiances, which vary on a linear scale. Threshold values have been chosen based on misfit histograms in order to remove abnormal tails. As a consequence, L1 and L2 pixels that pass the selection and are further assimilated could differ. However, the relative number of rejected observations for the entire month of July is quite limited in both cases (3 % for L1, 6 % for L2), thus not affecting the results statistically.

The setup of the background error covariance (B) is a critical step both for L2 retrievals and data assimilation. We did benefit from past experiences using MOCAGE, IASI and MLS O₃ (Massart et al., 2012; Emili et al., 2014; Peiro et al., 2018) to define a first guess of B, and we tried to further derive an optimal parameterization for this study. Note that the B matrix (3-D) used in data assimilation is by definition different with respect to the one specified within SOFRID (1-D), but the same 3-D B is used for all data assimilation experiments (L1 and L2).

Figure 3 Relative root mean square error (RMSE) of the control simulation with respect to radiosoundings (solid lines) and MLS (dotted lines) averaged globally (first plot) and for five separate latitude bands (90–60^∘ S, 60–30^∘ S, 30^∘ S–30^∘ N, 30–60^∘ N, 60–90^∘ N). To compute the percentage, the RMSE statistics have been divided by the corresponding average profile of the observations (radiosoundings or MLS) for each band.

Concerning the standard deviation, Emili et al. (2014) and Peiro et al. (2018) employed vertically varying errors expressed as percentages of the background O₃ profile, with larger relative errors in the troposphere and smaller in the stratosphere. Since we use a more detailed chemistry model here (Sect. 3.1), we evaluated the root mean square error (RMSE) of the free model simulation (control) against ozonesondes and MLS profiles (Fig. 3). We remark that the model's RMSE reproduces the vertical features observed in previous studies, with smaller errors in the stratosphere (between 20 and 50 hPa), larger errors in the free troposphere and highest errors close to the tropopause and within the planetary boundary layer. Note also the zonal variability of the maxima, which appear to be linked to the variability of the tropopause height. Thanks to the detailed chemical mechanism, biases (Fig. 4) are generally smaller than in the studies cited previously but remain significant compared to standard deviation values (Fig. 5), especially around the tropopause. Interestingly, RMSE and standard deviation values computed against MLS are generally smaller than those evaluated against ozonesondes, whereas biases are more consistent between the two datasets. We attribute this effect to the larger number of MLS observations (Table 1), which provides more robust standard deviation statistics.

Figure 4 Relative bias of the control simulation with respect to radiosoundings (solid lines) and MLS (dotted lines). Same plots as in Fig 3.

The background standard deviation is prescribed through a smooth step function that takes values of 2 % above 50 hPa and 10 % below to reproduce roughly the patterns observed in Fig. 5. Values are smaller than those in Fig. 5 to account for the error reduction during the assimilation, which is particularly strong when MLS observations are used (see Sect. 4.3). Also, neglecting error correlations between IASI channels within R leads to a strong weight towards the observations: reducing the background standard deviation compensates in part for this effect. All the choices made to define B are a result of a large number of assimilation evaluations, in which different options were considered. For example, setting values of 5 % and 25 % leads to less accurate results for both L1 and L2 assimilation (not reported). The percent profile is multiplied by the hourly O₃ field of the control simulation once for the entire period and not at every forecast time step. Therefore, all assimilation experiments presented in this study are based on the same B matrix. This choice was made to permit a stricter comparison between L1 and L2 assimilation experiments.

Figure 5 Relative standard deviation of the control simulation with respect to radiosoundings (solid lines) and MLS (dotted lines). Same plots as in Fig. 3.

The vertical error correlation diffuses the assimilation increments between model levels and has been found to significantly impact the quality of O₃ analyses with current model vertical resolutions (not shown). In general, small values of vertical correlation are favoured in the stratosphere due to the stratification and to avoid injection of large stratospheric O₃ increments in the troposphere, whereas larger values are expected within the troposphere due to vertical mixing. In this study a constant value of one model level defines the length scale of the Gaussian correlation (Sect. 3.3). Different choices for the stratosphere and troposphere did not lead to particular improvements (not shown).

Finally, the exponential scale of the horizontal error correlation is set to be equal to 200 km, with the zonal component that is reduced towards the poles to account for the increasing resolution of the model's grid (Emili et al., 2014).

The choice of the background and observation errors is relatively simplistic in this study. Further improvements of the B parameterization could be achieved by diagnosing the forecast errors hourly (Desroziers et al., 2005) or using ensembles of model forecasts. However, more complex and costly estimations do not always improve the results of chemical assimilation systematically and significantly (Massart et al., 2012). Moreover, a good estimation of B cannot be done independently of that of R, which is kept fixed here on purpose. Additional research is needed in this regard, which is beyond the scope of this study.

4 Results

A total of six simulations for the month of July 2010 have been performed (Table 3), starting on 1 July: a free model simulation (control) and five 3D-Var analyses assimilating SOFRID L2 profiles (named L2a), IASI L1 radiances (L1a), MLS L2 profiles (MLSa), MLS plus SOFRID L2 profiles (MLS + L2a) and MLS plus L1 radiances (MLS + L1a). The first three simulations (control, L2a and L1a) are discussed in Sect. 4.1. The control simulation and the three analyses that include MLS are discussed in Sect. 4.3. All simulations have been validated against ozonesondes profiles to elucidate the differences of the resulting O₃ vertical distribution. A total of 219 radiosoundings are available globally for July 2010 (Table 1). The co-location of ozonesonde profiles with model fields in time and space is performed through the MOCAGE observation operator (Sect. 3).

Table 3 Names of experiments and assimilated data.

4.1 IASI assimilation

The average O₃ values of the control simulation are displayed in Fig. 6. O₃ fields have been first interpolated vertically from the model grid to a selection of pressure levels, covering both the stratosphere and the troposphere, and averaged afterwards. The maps show well known properties of the O₃ distribution such as the strong zonal gradients in the stratosphere and the presence of local minima in the tropical free troposphere due to deep convection. The average difference between the control simulation and the fixed a priori profile used in SOFRID retrievals (Sect. 2.1.2) is displayed in Fig. 7. Large differences (>100 %) are found at low latitudes both in the lower stratosphere and in the troposphere, with largest values close to the tropical tropopause (150 hPa). This is expected since the SOFRID a priori profile is based on a global ozonesonde climatology that is more representative of mid-latitude O₃ profiles (Sect. 2.1.2).

Figure 6 Average O₃ values of the control simulation in units of parts per billion (ppb) for July 2010. From left to right, different pressure levels are displayed covering the stratosphere (a, b, c) and the free troposphere (d, e, f). Average, maximum and minimum values of the displayed fields are given above each map.

Figure 7 Relative average differences between the control simulation and the SOFRID a priori profile on July 2010. Values are given as percentages (%) of the control simulation (Fig. 6). Same plots as in Fig. 6.

We discuss the geographical differences between L1a and L2a analyses by looking at the monthly bias between the two experiments, divided by the average O₃ of the control simulation (Fig. 6). Relative differences are displayed in Fig. 8. First, we remark that differences are generally significant both in the stratosphere and in the troposphere, with absolute values that can exceed 50 % of the O₃ field locally and global averages as high as to 20 %. Largest differences in the stratosphere are found at tropical latitudes, L1a showing larger O₃ values than L2a at 20 hPa and lower at 70 hPa. In the troposphere the strongest positive differences are still found in the tropics, especially over central Africa, eastern Asia, South America and Middle East regions. Differences become smaller when moving down to 750 hPa and tend to disappear at lower altitudes (not shown), which is normal considering the vertical sensitivity of IASI. At mid-latitudes and high latitudes, relative differences are smaller than at the tropics. This behaviour is consistent with the fact that the SOFRID prior information is much less accurate for tropical latitudes than for mid-latitudes and high latitudes (Fig. 7). Overall, these plots suggest that when the L2 a priori profile is strongly biased (>100 %), the equivalence between L1 and L2 assimilation in the thermal infrared is not verified for O₃, even when the averaging kernels are employed.

Figure 8 Relative average differences (%) between radiances and Level 2 assimilation (L1a minus L2a divided by the corresponding O₃ values of the control simulation in Fig. 6) for July 2010. Same plots as in Fig. 6.

To confirm that the observed differences are not a consequence of the slightly different NWP inputs used in L1 assimilation and L2 retrievals (Table 2) we rerun the L1a simulation using exactly the same SST a priori values used in SOFRID retrievals. Since we cannot use the same water vapour and temperature profiles of SOFRID within L1a due to the different vertical grids, a different approach has been used: we repeated all assimilation experiments but using ERA Interim (Dee et al., 2011) instead of the NWP forecasts as meteorological forcing for the CTM. This increases potential differences between the L1 and L2 assimilation, due to the different configurations of operational NWP and ERA Interim (model resolution, assimilated instruments, etc.). In all the above cases we obtained very similar results to those presented in Fig. 8 (not shown), which suggests that differences between L1a and L2a discussed previously do not depend on the auxiliary RTM inputs.

To further verify which one out of the L1a and L2a experiments reproduces the measured O₃ profiles better, we validated the three simulations against radiosoundings. Figure 9 reports the RMSE differences computed globally and for five different latitude bands. The displayed values are the differences between the RMSE of the assimilation experiment and the corresponding value for the control simulation (Fig. 3). Negative values in Fig. 9 indicate that the assimilation improved the O₃ field and decreased the relative RMSE with respect to ozonesondes by the amount displayed on the plot. Looking at the global averages we remark that below 70 hPa the gain is similar for both L1a and L2a experiments and quite significant at 200 hPa (20 %). Note, however, the strong similarity between the global and 30–60^∘ N statistics, due to the over-representation of ozonesondes for the Northern Hemisphere (NH) mid-latitudes (63 % of the total).

Figure 9 Relative difference of RMSE (ΔRMSE) with respect to radiosoundings for L1a (blue) and L2a (red). The difference is computed by subtracting the RMSE of L1a (L2a) from the RMSE of the control simulation (Fig. 3). Negative values mean that the assimilation improved (decreased) the RMSE of the control simulation, and positive values indicate degradation (increase) of the RMSE. The statistics are computed for the same latitudes as in Fig. 3.

In the NH the RMSE of the control simulation is effectively reduced between 70 and 300 hPa (up to 20 %). L1a shows a slightly better gain than L2a between 150 and 300 hPa. Interestingly, both L1a and L2a display increased RMSE between 300 and 400 hPa. This behaviour is also confirmed when the vertical error correlation is switched off in the 3D-Var B and with different choices for the vertical interpolation of O₃ optical coefficients within RTTOV (log-linear or Rochon, not shown). Since large negative biases were present in the control simulation (as low as −30 %; see Fig. 4), a possible explanation is that part of the strong positive correction of O₃ between 100 and 300 hPa is propagated downwards, where both absolute O₃ concentrations and relative biases are much lower. This can degrade the analysis accuracy below 300 hPa. Whether this propagation is carried out by the Jacobian matrix of the observation operator (either through the RTM or the retrieval's AK) or by vertical O₃ transport is not yet elucidated and would need further investigation. Also, other possible factors affecting the accuracy of the RTM exist, like inadequate vertical resolution close to the tropopause or uncertainties in meteorological profiles or in the impact of aerosols. Nonetheless, these errors impact both L1a and L2a in our study: further optimization of the L1 assimilation configuration with respect to the L2 retrievals is left for a future study. The RMSE is reduced again at about 500 hPa between 30 and 60^∘ N, although not very significantly. The assimilation increases the RMSE of the tropospheric profile (350–1000 hPa) at northern latitudes (60–90^∘ N). In general, the validation confirms that L1a and L2a have a very similar accuracy in NH at mid-latitudes and high latitudes, as also suggested previously by Fig. 8. However, the strongest positive corrections are confined to the UTLS.

At the tropics (30^∘ S–30^∘ N) the results differ more significantly. In the troposphere (below 100 hPa), both L1a and L2a reduced the RMSE of the control simulation, although by a smaller amount than in the NH (5 %). Note also that L1a RMSE reduction is larger than L2a between 400 and 600 hPa, whereas it is the other way around at about 250 and 800 hPa. Above 100 hPa we observe an increase of RMSE that peaks at 60 hPa with L2a and at 30 hPa with L1a but smaller in magnitude for L1a. This behaviour might be linked to the strong differences that exist between the SOFRID prior information and the modelled O₃ at the tropical tropopause, to some other factor affecting the RT computations, to overestimation of the background error covariances or to a complex combination of all previous causes. A full satisfactory explanation has not been found yet.

Results in the Southern Hemisphere (SH) (30–90^∘ S) are again similar: a lower RMSE than for the control simulation is found for both L1a and L2a in the upper and lower stratosphere (between 30 and 100 hPa and between 150 and 300 hPa). An improvement is also found in the troposphere (400–600 hPa) at mid-latitudes (30–60^∘ S), but the low number of ozonesondes available in this band (Table 1) requires a more careful interpretation.

Figure 10 Relative difference of RMSE (ΔRMSE) with respect to MLS profiles for L1a (blue) and L2a (red). Same plots as in Fig. 9.

Since radiosoundings do not provide a uniform global coverage, and vertical coverage is also lacking in the vicinity of the O₃ maximum, we validated the three simulations against MLS measurements. The RMSE differences for stratospheric profiles can be found in Fig. 10. These statistics are based on more than 10⁵ profiles for the global average and between 15 000 and 30 000 for zonal averages, depending on the latitude band (Table 1). The patterns observed in the stratosphere with respect to ozonesondes are also confirmed with MLS. The only exceptions are a smaller RMSE degradation at 50 hPa for L2a in the tropics and for both L1a and L2a at 150 hPa in the 30–60^∘ S band. Higher confidence should be given to the RMSE values provided by MLS than those obtained with radiosoundings (see also Fig. 3 and the relative discussion in Sect. 3.4). However, a similar RMSE behaviour is observed overall, and this bolsters the robustness of the conclusions derived with the radiosoundings in the troposphere.

Figure 11 Relative average differences between MLS + L1a and MLS + L2a (MLS + L1a minus MLS + L2a divided by the corresponding O₃ values of the control simulation in Fig. 6) for July 2010. Same plots as in Fig. 8.

Figure 12 Relative difference of RMSE (ΔRMSE) with respect to radiosoundings for MLSa (teal), MLS + L1a (dark blue) and MLS + L2a (red). Same plots as in Fig. 9.

4.2 Computational cost

The computational cost of L1 assimilation is necessarily higher than for L2 assimilation. Additional CPU time is due not only to online RTM computations but also to a higher number of iterations needed by the minimizer to converge. For a typical 24 h long simulation performed on Intel Xeon E5-2680 V3 CPU, the total CPU time is 3.9 CPU hours for L2a and 13.2 CPU hours for L1a. Note that the L2a time does not include the cost of the L1 to L2 processor but only the cost of the 3D-Var assimilation plus the model forecasts. Most of the CPU time for L1a is spent on the linearized and adjoint calls of the RTM (50 % of the total CPU time), whereas the corresponding time spent for the observation operator within the L2a experiment is about 1 %. However, the total CPU time can be significantly decreased by reducing the maximum number of iterations of the minimizer. A simulation with a halved number of iterations (75) showed very similar results to the ones that have been reported (150 iterations) and could be considered if computation time is a critical factor. Moreover, with standard high-performance computers, and thanks to the parallel nature of the observation operator and the RTM, we could obtain a speed-up of about 24 on the 24 CPU cores. This reduces the runtime of L1a to about 36 min for the 24 h long simulation versus 13 min for L2a. The extra cost of L1 assimilation therefore also seems acceptable for operational applications.

4.3 IASI and MLS assimilation

Some issues were identified in the previous section in the stratosphere, especially at tropical latitudes. Among possible reasons, one is that inversion of TIR measurements might be particularly sensitive to the vertical distribution of O₃ in the tropical stratosphere. We consider assimilating MLS L2 profiles in combination with IASI here to correct the model stratosphere and troposphere simultaneously, as also done in previous studies (Emili et al., 2014; Peiro et al., 2018). When the radiances are assimilated, the RT problem is solved for the entire atmospheric column within the iterations of the variational algorithm. Therefore, enhanced and better synergies could be observed than when only L2 products are assimilated.

We report in Fig. 11 the impact of assimilating MLS in combination with IASI L1 and L2 by computing the average differences between MLS + L1a and MLS + L2a. We remark that the differences in the stratosphere are highly reduced with respect to Fig. 8, which is expected due to the direct constraint of MLS observations. Significant differences (>10 %) remain below 150 hPa, with patterns and sign similar to those in Fig. 8. The amplitude of the differences is, however, also slightly reduced at 300 and 500 hPa.

Figure 13 Taylor diagrams of modelled tropospheric ozone columns (340–750 hPa) for the Control simulation (green), MLSa (violet), MLS + L1a (grey) and MLS + L2a (yellow) averaged globally and for five separate latitude bands. The Taylor statistics are computed against radiosoundings.

We compared the RMSE of MLSa, MLS + L1a and MLS + L2a computed against ozonesondes (Fig. 12) to evaluate if the joint assimilation improves the overall O₃ distribution. MLSa provides particularly accurate results down to 200 or 300 hPa, depending on the latitude, with a robust reduction of the RMSE with respect to the control simulation. The only exception is in the SH mid-latitudes below 250 hPa, where the MLSa RMSE increases. We suspect that this might be linked again to the combination of strong O₃ gradients at the tropopause height and the negative bias of the control simulation above the tropopause (see Sect. 4.1). Overall MLSa confirms results found in past studies (Massart et al., 2012; Emili et al., 2014) and represents much better prior information for assimilation of radiances or retrievals.

We remark that MLS + L1a and MLS + L2a now provide closer results in the NH and in the tropics compared to Fig. 9. The stratospheric O₃ gain is much more significant with MLS + L1a and MLS + L2a than with L1a and L2a and remains very close to MLSa, demonstrating that assimilating accurate stratospheric profiles remains essential for O₃ reanalyses. The only region where IASI further improves the UTLS profile with respect to MLSa is in the NH: a positive, albeit small, effect of assimilating IASI in combination with MLS is found between 150 and 300 hPa. On the other hand, below 300 hPa, the addition of MLS (Fig. 12) does not bring further improvements with respect to IASI alone (Fig. 9). We can conclude that MLS corrects most of the errors introduced by IASI assimilation in the stratosphere (Fig. 9), but no particular synergy is observed in the case of MLS and L1 assimilation in the troposphere.

In Fig. 13 we report the Taylor plots concerning the free troposphere O₃ column (340–750 hPa), to further evaluate the skills of the assimilation experiments in terms of variability. We examine the free troposphere here since it is where the direct impact of IASI assimilation is the largest and the impact of MLS the smallest (except for the 30–60^∘ S band). IASI assimilation improves the variability of the modelled O₃ field when looking at global averages, but this conclusion varies as a function of the latitude band. Robust and significant improvements are only found at the tropics and in the SH polar region; mixed results are obtained elsewhere. This confirms previous findings obtained with L2 assimilation (Emili et al., 2014) and adds the conclusion that better prior information does not necessarily solve all issues related to the assimilation of TIR measurements at mid-latitudes and high latitudes. Nevertheless, the assimilation of radiances provides slightly better results at all latitudes in general and permits more variability to be extracted from IASI spectra, especially at tropical latitudes.

5 Conclusions

In this study we addressed the following question: what are the differences between the direct assimilation of IASI radiances (Level 1) and the assimilation of Level 2 products for O₃ analyses and reanalyses? We used an experimental setup in which differences between the L2 retrieval and the L1 assimilation have been minimized as much as possible, for example by using the same RTM (RTTOV) and control vector (O₃ and SST) in both approaches. This allowed the impact of the O₃ prior information and its error covariance on the quality of the analysis to be delved into.

We performed twin assimilation experiments with the MOCAGE CTM and the SOFRID O₃ retrievals, using the same IASI ground pixels for both L1 and L2 assimilation, named L1a and L2a respectively. We compared the obtained analyses against each other and against ozonesondes and MLS profiles for the month of July 2010.

The results suggest that the accuracy of the O₃ prior information used in the L2 retrievals can influence the analysis, even when the averaging kernels are employed within the assimilation. When the O₃ prior information is strongly biased (at low latitudes in this study), L1a and L2a differ significantly (>10 %) and the analysis shows a better variability when assimilating directly L1 radiances instead of L2 profiles. L1a and L2a are otherwise very similar at mid-latitudes and high latitudes, where the SOFRID prior information is closer to the true O₃ profile.

We conclude that particular care should be taken before assimilating satellite retrievals with prior information that can, in some circumstances, differ significantly from the local ozone profile. Computing retrievals using an a priori profile issued from a model could be relevant in improving current IASI O₃ L2 products and might reduce the differences between L2a and L1a observed in our study. Preliminary results with SOFRID based on a modelled a priori profile also show significant differences with the original product (Brice Barret, personal communication, 2019), with patterns similar to those presented in this study. However, when the final purpose is data assimilation, the L1 approach is more practical and statistically consistent, especially in the case that the observations need to be assimilated within the same forecast model that was used to compute L2 retrievals.

A positive impact has been found when assimilating MLS profiles and IASI simultaneously (either L1 or L2), which corrected stratospheric biases due to IASI assimilation alone. Differences between L1 and L2 assimilation are globally reduced by MLS in the stratosphere but remain significant (>10 %) in the tropical troposphere. Also, MLS assimilation strongly improves the model's accuracy down to 200 hPa, and a clear added value of IASI assimilation (L1 or L2) can only be observed in the tropical troposphere. These results remind us that the information brought by limb sounders like MLS into the DA system remains essential to improve upper stratosphere O₃. Interesting perspectives for future work are to (i) verify whether the assimilation of O₃ retrievals from UV spectrometers like GOME-2 or TROPOMI also shows issues related to the a priori dependence and (ii) examine if UV assimilation could replace MLS when assimilated jointly with IASI and provide similar performances in the stratosphere. This will be important to ensure the capacity to carry out accurate O₃ reanalyses when the MLS instrument is phased out.

We reckon that L1 assimilation requires the full atmosphere to be modelled, which may not be available to some models, for example those conceived exclusively for tropospheric applications. Moreover, Level 2 products can be aggregated vertically to correct some model layers selectively and averaged spatially to fit models with coarser resolution than the satellite ground-pixel size. This cannot easily be done with radiances and should be addressed in future research.

In this study the observations, their error covariance and the RTM auxiliary inputs were kept almost identical between L1 and L2 assimilation on purpose. Further research is needed to address issues that are common to L1 and L2 assimilation, e.g. increased errors close to the tropopause in the NH or in the tropical stratosphere. Improvements are expected, for example, by increasing the vertical resolution of the model, including modelled aerosols within the RT or using more realistic observation error covariances. Including more modelled variables among the RTM inputs is in particular of interest in the context of the evolution towards ESMs, for which hyper-spectral sounders like IASI can provide very valuable constraints for multi-variate reanalyses (atmosphere plus surface). Including inter-channel and ground-pixel correlations in the observation error covariance matrix seems necessary to correctly weight very dense IASI observations within higher-resolution models than the one used in this study. All these aspects deserve further research.