Estimation of the error covariance matrix for IASI radiances and its impact on ozone analyses

In atmospheric chemistry retrievals and data assimilation systems, observation errors associated with satellite radiances are chosen empirically and generally treated as uncorrelated. In this work, we estimate inter-channel error covariances for the Infrared Atmospheric Sounding Interferometer (IASI) and evaluate their impact on ozone assimilation with the chemical transport model MOCAGE (MOdèle de Chime Atmospheric à Grande Echelle). The method used to calculate observation errors is a diagnostic based on the observation and analysis residual statistics already adopted in numerical weather prediction 5 centers. We used a subset of 280 channels covering the spectral range between 980 and 1100 cm−1 to estimate the observation error covariance matrix. We computed hourly 3D-Var analyses and compared the resulting O3 fields against ozonesondes and the measurements provided by the Microwave Limb Sounder (MLS). The results show significant differences between using the estimated error covariance matrix with respect to the empirical diagonal matrix employed in previous studies. The validation of the analyses against independent data reports a significant im10 provement, especially in the tropical stratosphere. The computational cost has also been reduced when the estimated covariance is employed in the assimilation system.

For this study, L1c data 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 co-located land mask and cloud fraction values, ob-15 tained from the Advanced Very High-Resolution Radiometer (AVHRR) measurements, also on board Metop.

MLS
The Microwave Limb Sounder (MLS) provides vertical profiles of several chemical components, by measuring the microwave thermal emission from the limb of Earth atmosphere (Waters et al., 2006). More than 2500 vertical profiles are observed daily, including trace gases with a vertical resolution of approximately 3 km. Several studies benefited from MLS products, notably 20 5 https://doi.org/10.5194/amt-2020-179 Preprint. Discussion started: 14 September 2020 c Author(s) 2020. CC BY 4.0 License. the ozone profiles in assimilation experiments (Emili et al., 2014;Massart et al., 2009), thanks to its low bias in the stratosphere (<5%) (Froidevaux et al., 2008).
In our study, we use the ozone profiles retrieved from MLS (V4.2 Products) as independent data to validate our results. The data have been downloaded from the Goddard Earth Sciences Data and Information Services Center (GES DISC) web portal (https://disc.gsfc.nasa.gov). 5
It provides complete global maps of total column ozone on a daily basis. The OMI ozone data record starts in October 2004, shortly after the launch of Aura (McPeters et al., 2015). The total column averaged over the month of the study (July 2010), resulting from the OMI-TOMS version 8 algorithm (Bhartia, 2002), is used here to validate the results of the assimilation 10 experiments.

Ozonesondes
Ozonesondes are in situ instruments carried by a radiosonde transmitting continuously the measurements as the ascent in the atmosphere.

Setup of the numerical experiments
The main purpose of this study is to estimate the IASI observation error covariance and verify its impact on the quality of the ozone assimilation results. Thus, we kept the same setup of (Emili et al., 2019) in terms of the period of the study, the model 20 configuration, the choice of assimilated observations, and of the background error covariance matrix. The observation error covariance matrix will be discussed in the section of the results (Section 3).
The model was initialized with a zonal climatology and the spin-up time used is one month (June 2010). Then, our simulations were performed for the month of July 2010. The forecast error standard deviation was assumed to be proportional to the ozone concentration. In fact, Emili et al. (2019) have evaluated the standard deviation of the free model simulation against 25 independent data (profiles from ozonesondes and MLS), and found out a small free forecast error in the stratosphere, larger error in the free troposphere and highest error close to the tropopause. The background standard deviation was, thus, taken equal to 2% above 50 hPa and 10 % below to mimic the validation's behavior. Similar choices were employed in (Massart et al., 2012;Peiro et al., 2018).
The background error covariance matrix is split into a diagonal matrix filled with the standard deviation and correlation matrix modeled using a diffusion operator. The correlation, characterized by the length-scale, spreads the assimilation increments in space. In this study, we keep the same configurations of horizontal and vertical length scales as in (Emili et al., 2019) .
The same preprocessing described in Emili et al. (2019) has been applied to our data before their use in the assimilation system. In order to avoid any contamination from clouds, data were filtered using a cloud mask and only pixels with cloud 5 fraction less than or equal to 1 % were kept. The cloud fraction values are obtained from the AVHRR measurements onboard Metop. Since the spatial resolution of MOCAGE is coarser than the pixel size, the number of ground pixels was reduced by thinning the data using a grid of 1°x 1°of resolution and only keeping the first pixel that falls in every two grid boxes. A dynamical rejection of observations -with a threshold of 12 % -based on the relative differences between simulated and measured values with respect to simulated values was considered. Some channels affected by H 2 O absorption (1008-1019,1028-1030,

Desroziers diagnostics
The observations used in the assimilation system could have a margin of error. We can identify two types of errors, system-15 atic and random errors. The systematic error is ordinarily corrected before the data assimilation process. Nevertheless, this correction should not account for the model bias, which requires some independent data to prevent the drift of the analysis to unrealistic values. In our case, we are assimilating troposphere and stratosphere data, the missing of an anchor instrument in the top of stratosphere led us to assume that our observations are unbiased. This assumption, in fact, is adopted in most chemical analyses (Emili et al., 2019;Massart et al., 2012). Random errors can arise from the measurements (e.g. instrumental error), 20 forward model, representativeness error (e.g. difference between point measurements and model representation), or quality control error (e.g. error due to the cloud detection scheme missing some clouds within clear sky only assimilation). These types of errors should be accounted by the observation error covariances matrix R. According to (Weston et al., 2013), the instrument noise could be assumed to be uncorrelated. However, the IASI measurements are apodized, which may introduce correlations between neighboring channels, particularly in our case where we are assimilating a subset of adjacent channels.

25
The radiative transfer may also introduce correlations between channels. The statistics of error from the instruments noise are known, while the characteristics of other sources of error are not yet well understood.
In this paper, we estimate the total error using the statistical approach introduced by Desroziers et al. (2006) .
Where x a is the analysis state vector, x b is the background state vector, y is the vector of observations and H is the observation 30 operator that computes model counterpart in the observation space.
This method has been used to estimate observation errors and inter-channel error correlations (Stewart et al., 2009). It can potentially provide information on imperfectly known observation and background-error statistics with a nearly cost-free computation (Desroziers et al., 2006). However, this approach assumes that the R and B matrices used to produce the analysis are exactly correct, which may not be always the case. Furthermore, Desroziers diagnostics compute the total covariances, more efforts are needed to understand and distinguish the sources of the error.

Error results
The Desroziers statistics were computed on the output of a 3D-Var experiment using a diagonal matrix R (with a standard deviation of 0.7 mWm −2 sr −1 cm as in (Emili et al., 2019). The diagnosed R could not be used directly in the assimilation system. In fact, the estimated matrix was asymmetric and not positive definite. This might result from the violation of one of the assumptions of Desroziers statistics stipulating that the B and R matrices used to produce the analysis were correctly 10 defined. Similar unrealistic features in the diagnosed covariance matrices were encountered in (Stewart et al., 2014;Weston et al., 2013) where an artificial inflation of observation errors was assumed. R needs to be a valid covariance matrix before being used in the 3D-Var assimilation system. Therefore, we first symmetrize the estimated matrix by taking the mean of the original and its transpose. Then we impose the negative eigenvalues to be equal to the smallest positive eigenvalue (Charles F. Van Loan, 1996). An other method was tested here to recondition the estimated matrix, the one called ridge regression (Weston 15 et al., 2013;Tabeart et al., 2020b) which consists of increasing all eigenvalues of R by the same amount. We favoured the first method since the standard deviation and the correlation values remain closer to the initially estimated quantities.
In order to evaluate the impact of starting from an experiment based on a diagonal R matrix, we performed a second 3D-Var experiment using the diagnosed and modified matrix. Then, we used these new forecasts and analyses to estimate again R with the Desroziers diagnostics. The resulted standard deviation was significantly greater than the one estimated from the experiment 20 that use a diagonal matrix (not shown). The same goes for the correlations (not shown). It should be noted that the re-estimated matrix was positive definite, unlike the first estimation where some unrealistic features were encountered. We have followed the same process to reestimate two other matrices (introduce, each time, the estimated matrix into the assimilation system in order to have analyses that we use to re-estimate a new matrix). The differences between the estimations in terms of standard deviation and correlations became smaller as we reestimate the matrices (not shown), suggesting a sort of convergence of the 25 estimation. All matrices and simulations discussed in this paper are based on the second estimation (which uses the analyses derived from the experiment using the first conditioned estimation). The reason for this choice will be rediscussed later (section 5.2). common setting for most IASI O3 retrievals (Dufour et al., 2012). At first glance, we note that the standard deviation used in previous studies is highly underestimated for the SST sensitive channels and overestimated for some ozone sensitive channels (around 1040 and 1050 cm −1 ). The diagnosed standard deviation increases to reach 2 mW m −2 sr −1 cm for SST sensitive channels (the first and the last twenty channels of the band (980-1000 cm −1 and 1080-1100 cm −1 ) and the channels between the observations are greater for the SST channels than those of the ozone. The same goes for the corresponding background and the analysis values. Since these diagnostics are based on observation, background and analysis residuals, a larger standard deviation for the SST channels than for ozone channels might be expected. We remarked that the estimated standard deviation is proportional to the radiance value (either the observation, the background, or the analysis value), which gives a relative 5 standard deviation of about 2.5 % of the average radiances values for the entire spectral window (not shown).
The IASI instrumental error is provided by the CNES (Centre National d'Etudes Spatiales ), taking into account different known effects such as flight homogeneity and apodization effect (Le Barbier Laura, personal communication, September 18, 2019). The instrumental error covariance matrix is computed as described in (Serio et al., 2020). This error remains smaller (maximum of 1.06e −6 mW m −2 sr −1 cm) than that used in the previous studies (0.7 mW m −2 sr −1 cm). Then, the important 10 estimated standard deviation noticed in our estimation might be due to the radiative transfer inputs error.
To investigate the off-diagonal part of R we present the correlation matrix in figure 2. The results show high correlations between the majority of the channels (larger than 0.4 ). In particular, a very high correlation is observed among SST sensitive channels (around 0.9 to 1). The regions of, relatively, lower correlation (around 0.4 to 0.7) represent the ozone channels correlations and cross correlation between ozone and SST sensitive channels.

15
The high correlation found here was expected since previous studies have highlighted the same behaviour in this spectral region (Bormann et al., 2010;Stewart et al., 2014;Bormann et al., 2016). In fact, the use of the same radiative transfer model for all channels may introduce similar errors among these channels.
The diagnostic discussed above is based on a global estimation, without any distinction between the type of the surface (land or sea) nor the time of the observation (day or night). Since the emissivity varies according to the type of the surface, and the 20 skin temperature is strongly driven by the sun radiation, we evaluated R taking these differences into account. Figure 3 shows the difference between the standard deviation estimated using data localized over the sea and land separately.
The error over land reveals large values for the SST sensitive channels in comparison with that estimated over the sea which, in turn, reproduces a slightly different error in comparison with the global estimation. The two surfaces introduce also a slightly different error regarding the ozone band. Since the number of observations over the sea is greater than over the land (almost 25 65 %), the global estimation is dragged down to be close to the sea estimation. To explain the large difference in the region of SST sensitivity channels over the land, the surface properties need to be discussed.
The surface emissivity varies with vegetation, soil moisture, composition, and roughness (Nerry et al., 1988) with values between 0.65 and 1 in the thermal infrared range. Low values are found in deserts and high values over dense vegetation and water surface (Capelle et al., 2012). The variability of the soil type over land in comparison with the sea that is homogenous 30 can lead to a larger standard deviation over land than over sea. The second parameter that impacts the electromagnetic spectrum measured by the infrared sensor is the surface temperature. Its variability is larger over land than over sea, which also contributes to the difference of the standard deviation between land and sea. On the other hand, these differences are larger for the SST channels than for the ozone channels. This might be explained by the direct link between the radiance signal and the 9 https://doi.org/10.5194/amt-2020-179 Preprint. Discussion started: 14 September 2020 c Author(s) 2020. CC BY 4.0 License.  surface signal in the SST bands especially since the observation and the estimated standard deviation were found proportional (mentioned above).
The same behaviour as the global estimation is reproduced when the statistics were performed from the data measured separately from the day and from the night. Figure 4 shows small differences in the SST channels sensitivity band (between 1080 and 1100 cm −1 ), whereas the curves take approximately the same values for the rest of the channels.

5
The data assimilation process uses the full error covariance matrix in order to minimize the cost function. Hence, we discuss now the differences between the correlation matrices taking the time (day/night) and the surface type (sea/land) of the observations into account. Figure    In conclusion, the variability in terms of correlations is more pronounced when the surface type is considered than in the case of the observation time. For the latter, the smaller differences found for the standard deviation (Fig. 4) are also found for correlations. For the follow-up assimilation experiments, we adopted the global estimation and neglected the effect of the time and surface of the observations. The rationale for this choice will be given during the discussion of the validation results (section 5.2).  4 Assimilation results

Ozone fields
In this section, we discuss the impact of the observation error covariance estimated previously on the ozone analysis. To this end, three experiments for the month of July 2010 were carried out: i). a model run without data assimilation called hereafter the free run (or Control), and noted in the rest of this paper 5 ControlExp.
ii). A 3D-Var assimilation of IASI radiances using a diagonal observation error covariance matrix (as in Emili et al. (2019)).
It will be referred here by RefExp.
iii). A 3D-Var assimilation of IASI radiances using a full matrix estimated with the Desroziers diagnostic noted hereafter by RfullExp.
The first experiment (ControlExp) was run to evaluate the benefit of the assimilation experiments and to quantify the improvements of each of the two analyses when they are validated against independent data. The same setup of Emili et al. (2019) was adopted for RefExp, which was taken as a reference to characterize the impact of accounting for the estimated R in the 5 third simulation (RfullExp). Figure 6 shows the difference between the zonal average of the analysis from the two assimilation experiments in units of parts per billion volume (ppbv). The zonal values were averaged over the month of the study before performing the difference.
The impact of the estimated R varies with latitude. It varies also with the height, adding or reducing the amount of ozone.
Overall, the estimated R reduces the amount of ozone in the high latitudes of the free troposphere and the tropical high 10 stratosphere, whereas the amount is increased in the vicinity of the lower stratosphere. The maximum reduction of ozone is larger than the amount added. The amount of ozone reduction reaches 600 ppbv, whereas the increase does not exceed 300 ppbv. In high northern latitudes (30°N-90°N), a significant addition is found (300 ppbv) covering almost the whole stratosphere, in opposition to the other latitudes where the difference changes sign with altitude. On the other hand, an important reduction of ozone is observed in the tropics at 20 hPa (more than 600 ppbv). To better understand the impact of the estimated R we 15 validate the results with independent data in the section of validation (section 5). The emphasis, up to now, has been on the impact of the estimated observation error on the ozone analysis. We discuss in the following subsection the impact on the skin surface temperature analysis.

Surface skin temperature
The assimilated spectra include both ozone and surface skin temperature sensitive channels. The IFS skin temperature was taken as a background in the assimilation process. We have computed the difference between the SST analysis and the background at

Computational cost
In our assimilation setup, the cost function is minimized hourly. For each window, the minimizer needs to converge after a 15 certain number of iterations. The cost of each iteration is dominated by the cost of the radiative transfer operators (tangent linear, the adjoint model) and of the background error covariance operators. When the observation error was assumed to be uncorrelated (RefExp), the number of iterations needed for each hourly cycle is significantly higher than when the estimated observation error covariance matrix is used. In fact, the introduction of the estimated R reduces the number of iterations from 150 (a fixed value to stop iterations if the convergence criteria were not achieved to save computational time) to 90 iterations in average. This means that the CPU time is reduced by more than 150% for each assimilation cycle. This improvement is highly valuable when the study is extended to longer periods. This reduction is due to the fact that the diagonal matrix is 5 pulling much closer to the observations than the correlated matrix, which makes it harder to find a solution resulting in slower convergence. This increase of the convergence speed was encountered in the Met Office 1D-Var system (Tabeart et al., 2020a) where a correlated observation matrix was introduced in the system. Furthermore, in Tabeart et al. (2018) the matrix R and the observation error variance appear in the expression of the condition number of the Hessian of the variational assimilation problem, indicating that these terms are important for convergence of the minimization function. To understand which one 10 between correlations and variance has a higher impact on the convergence speed, we proceeded as follows: we considered an estimation that did not lead to a convergence (the first estimation) of the minimizer, we kept its correlations and replace its standard deviation with that of an estimation that leads to a convergence (the second estimation). The minimizer, consequently, converges allowing us to conclude that the variance has a higher impact on the convergence speed.
5 Validation of O3 analyses 15 5.1 Total column Figure 8 shows the difference of the ozone total column (in Dobson Unit (DU)) provided by OMI and that of the RefExp (a) and that of RfullExp (b). At first sight, we note smaller differences over the tropics between the OMI total column and the total column given by RfullExp in comparison with that given by the RefExp. This behaviour was expected since an important reduction of the amount of ozone was observed in these regions (see figure 6). In the high northern latitudes, the differences 20 were slightly increased when the estimated matrix was adopted. This is a consequence of the increase in the amount of ozone encountered in these regions the stratosphere, compared to the amount reduced in the same region in the troposphere ( figure   6). On the other hand, the global mean and the standard deviation of these differences are lower in the case of using the new estimated matrix.
Hence, we conclude that the estimated matrix R has slightly improved the results in terms of ozone total columns.

Vertical validation
In this section, we validate the two simulations against radiosoundings and MLS data. We use the root mean square error (RMSE) as the main statistical indicator to quantify the accuracy of the experiments.
We compute the relative (to the control simulation) difference of RMSE with respect to radiosoundings and MLS averages globally and for five different latitude bands. The difference is computed by subtracting the RMSE of each experiment from 5 that of the control simulation. Negative values indicate an improvement of the O3 profiles. It should be noted that the representativeness of the statistics given by the MLS in the stratosphere is better than that of the radiosoundings because the number of profiles provided by MLS is much higher compared to the radiosoundings ones. Consequently, higher confidence is given to the validation using the MLS data in the stratosphere. Figure 9 reports the RMSE differences with respect to the radiosoundings. Considering the global RMSE (ALL), we notice 10 that the experiment with the estimated matrix improves the results above 150 hPa, around 400 hPa and near the surface. All things considered, the introduction of the estimated R has globally improved the O3 profiles in the stratosphere and in the free troposphere, especially in the tropics. In spite of its degradation in the vicinity of UTLS, the improvement remains always advantageous with respect to the control run.
The matrix used for this study (see Sec. 3.2) will be now discussed in this section since the decision was also based on the 15 outcome of the assimilation experiments presented in this section. We performed sequentially three assimilation experiments using the first, second and the third estimation of R (Sec. 3.2). The results of validation against radiosoundings and MLS showed small differences (not shown). Therefore, to avoid the initial impact of using a diagonal matrix we have adopted the second estimation (which uses the analyses derived from the experiment using the first conditioned estimation). In an operational framework, we may estimate the matrix daily (weekly or monthly if the period of the study is considerably long) 20 using the analyses of the previous day (using the analysis of the previous week or month respectively ). In other words, we may use a diagonal matrix to produce analyses for the first day or spin-up period, these analyses will be used to estimate the matrix that will be used for the second day, and so on throughout the period of the study.
We have also discussed the type (sea/land) and the time (day/night) of the observations while estimating the matrices. To check the impact of these differences on the assimilation results, we ran an additional assimilation experiment using the matrix 25 estimated considering the type of the surface of each observation (since the differences were more important than if the time of the observation was considered). Only slight differences among the results have been noticed (not shown). Thus, for simplicity, it seems reasonable to adopt the global estimation of the matrix and neglect the effect of the time and the type of the surface of the observations.
The correct specification of the observation error becomes a critical issue to assimilate efficiently the increasing amount of satellite data available in the recent years. We have estimated the observation errors and their inter-channel correlations for clear sky radiances from IASI ozone sensitive channels. We have evaluated, then, the impact of the estimated R on the SST and ozone analysis within our 3D-Var assimilation system. The outcome has been compared with an assimilation experiment 5 where the observation error covariance matrix was assumed to be diagonal and the standard deviation assigned empirically like in previous studies. The results of the experiments were, then, validated against independent data: OMI, MLS and ozonesondes.
The Desroziers diagnostics were adopted here to estimate R. The diagnostics used the analyses derived from a variational data assimilation experiment. The results have shown high correlations between the majority of the IASI spectral channels, particularly among the SST sensitive channels.

10
The impact of the surface type (land/sea) and the time of the observation (day/night) on the estimation was evaluated.
Significant differences were noticed in terms of standard deviation when the surface type was considered. However, slight differences in the results have been noticed when the matrix estimated considering the time and the type of the surface of the observation was used in the assimilation experiment. A global estimation of the matrix could finally be adopted.
Significant differences between the results of the experiments were encountered. The introduction of the estimated R reduces 15 the amount of ozone in the free troposphere and in the high tropical stratosphere, whereas ozone is added in the vicinity of the lower stratosphere. The total column also has shown differences on average and in terms of geographical distribution. A validation against OMI has shown that the results were closer to the observations when the estimated matrix was adopted.
The validation against MLS and ozonesondes showed that the introduction of the estimated R has globally improved the results in the stratosphere and in the free stratosphere especially in the tropics. In spite of a slight reduction in accuracy in the 20 vicinity of UTLS, the improvement remains always advantageous with respect to the reference assimilation. Concerning the computational cost, the introduction of the estimated R significantly reduces the number of iterations needed for the minimizer to converge.
In summary, accounting for an estimated R improves significantly the ozone assimilation results. Furthermore, this approach might be adopted in the assimilation of other chemical species and also in a chemical retrieval process framework.

25
In this study, the estimation was computed without taking into account any distinction of the error sources and assuming that the observation error was unbiased. More efforts will be needed to tackle these points. It should also note that we kept the same experiment setup of Emili et al. (2019) in order to be able to quantify exclusively the impact of the R. The background error matrix was still defined using a relatively simple and empirical method. Further research might be needed to perform a better estimation of the background error. A new channel selection might also be performed to reduce the computational cost