Assimilation of surface NO 2 and O 3 observations into the SILAM chemistry transport model

This paper describes the assimilation of trace gas observations into the chemistry transport model SILAM (System for Integrated modeLling of Atmospheric coMposition) using the 3D-Var method. Assimilation results for the year 2012 are presented for the prominent photochemical pollutants ozone (O3) and nitrogen dioxide (NO2). Both species are covered by the AirBase observation database, which provides the observational data set used in this study. Attention was paid to the background and observation error covariance matrices, which were obtained primarily by the iterative application of a posteriori diagnostics. The diagnostics were computed separately for 2 months representing summer and winter conditions, and further disaggregated by time of day. This enabled the derivation of background and observation error covariance definitions, which included both seasonal and diurnal variation. The consistency of the obtained covariance matrices was verified using χ diagnostics. The analysis scores were computed for a control set of observation stations withheld from assimilation. Compared to a free-running model simulation, the correlation coefficient for daily maximum values was improved from 0.8 to 0.9 for O3 and from 0.53 to 0.63 for NO2.


Introduction
During the past 10-15 years, assimilating observations into atmospheric chemistry transport models has been studied with a range of computational methods and observational data sets.The interest has been driven by the success of advanced data assimilation methods in numerical weather prediction (Rabier, 2005), as well as by the development of op-erational forecast systems for regional air quality (Kukkonen et al., 2012).Furthermore, the availability of remote sensing data on atmospheric composition has enabled construction of global analysis and forecasting systems, such as those described by Benedetti et al. (2009) and Zhang et al. (2008).Assimilation of satellite observations into stratospheric chemistry models has been demonstrated, e.g. by Errera et al. (2008).
Data assimilation is defined (e.g.Kalnay, 2003) as the numerical process of using model fields and observations to produce a physically and statistically consistent representation of the atmospheric state -often in order to initialise the subsequent forecast.The main techniques used in atmospheric models include the optimal interpolation (OI, Gandin 1963), variational methods (3D-Var and 4D-Var, Le Dimet and Talagrand, 1986;Lorenc, 1986), and the stochastic methods based on the ensemble Kalman filter (EnKF, Evensen, 2003Evensen, , 1994)).Each of the methods has been applied in air quality modelling.Statistical interpolation methods were used by Blond and Vautard (2004) for surface ozone analyses and by Tombette et al. (2009) for particulate matter.The EnKF method has been utilised by several authors (Constantinescu et al., 2007;Curier et al., 2012;Gaubert et al., 2014) especially for ozone modelling.The 3D-Var method has been applied in regional air quality models by Jaumouillé et al. (2012) and Schwartz et al. (2012), while the computationally more demanding 4D-Var method has been demonstrated by Elbern and Schmidt (2001) and Chai et al. (2007).Partly due to its significance in relation to health effects, the most commonly assimilated chemical component has been ozone.
The performance of most data assimilation methods depends on correctly prescribed background error covariance Published by Copernicus Publications on behalf of the European Geosciences Union.J. Vira and M. Sofiev: Assimilation of surface NO 2 and O 3 observations matrices (BECM).This is particularly important for 3D-Var, where the BECM is prescribed and fixed throughout the whole procedure, in contrast to the EnKF based assimilation methods, where the BECM is described by the ensemble of states, and to the 4D-Var method, where the BECM is prescribed but evolves implicitly within the assimilation window.
A range of methods of varying complexity have been employed to estimate the BECM in previous studies on chemical data assimilation.The "National Meteorological Centre" (NMC) method introduced by Parrish and Derber (1992) is based on using differences between forecasts, with differing lead times as a proxy for the background error.Kahnert (2008), as well as Schwartz et al. (2012), applied the NMC method for estimating the BECM for assimilation of aerosol observations.Chai et al. (2007) based the BECM on a combination of the NMC method and the observational method of Hollingsworth and Lönnberg (1986).This observational method was also used by Kumar et al. (2012) the in assimilation of NO 2 and O 3 data.
The BECM can also be estimated using ensemble modelling; this approach was taken by Massart et al. (2012) for global and by Jaumouillé et al. (2012) for regional ozone analyses.Finally, Desroziers et al. (2005) presented a set of diagnostics, which can be used to adjust the background and observation error covariances.This method has been previously applied in chemical data assimilation for example by Schwinger and Elbern (2010) and Gaubert et al. (2014).
In contrast to short and medium range weather prediction, the influence of initial condition on an air quality forecast has been found to diminish as the forecast length increases.For ozone, Blond and Vautard (2004) and Wu et al. (2008) found that the effect of the adjusted initial condition extended for up to 24 h.Among other reactive gases, NO 2 has been a subject for studies of Silver et al. (2013) and Wang et al. (2011).However, the shorter lifetime of NO 2 limits the timescale for forecast improvements especially in summer conditions.
An approach for improving the effectiveness of data assimilation for short-lived species is to extend the adjusted state vector with model parameters.Among the possible choices are emission and deposition rates (Bocquet, 2012;Curier et al., 2012;Elbern et al., 2007;Vira and Sofiev, 2012).
The aim of the current paper is to describe and evaluate a regional air quality analysis system based on assimilating hourly near-surface observations of NO 2 and O 3 into the SILAM chemistry transport model.The assimilation scheme was initially presented by Vira and Sofiev (2012); in the current study, the scheme is applied to photochemical pollutants and moreover, we discuss how its performance can be improved by introducing statistically consistent background and observation error matrices.The analysis fields are produced for the assimilated species at an hourly frequency using the standard 3D-Var assimilation method (Lorenc, 1986).The diagnostics of Desroziers et al. (2005) are applied in this work for estimating the background and observation error standard deviations, notably resolving their seasonal and diurnal variations.The evaluation is performed for the year 2012 using stations withheld from assimilation.In addition to assessing the analysis quality, the effectiveness of assimilation for initialising the model forecasts is evaluated.
The following Sect. 2 presents the model setup and briefly reviews the 3D-Var assimilation method.The procedure for estimating the background and observation error covariance matrices is discussed in Sect.3. The assimilation results for O 3 and NO 2 for the year 2012 are discussed in Sect. 4. Section 5 concludes the paper.

Materials and methods
This section presents the SILAM dispersion model, the observation data sets used, and describes the assimilation procedure.

The SILAM dispersion model and experiment setup
This study employed the SILAM chemistry transport model (CTM) version 5.3.The model utilises the semi-Lagrangian advection scheme of Galperin (2000) combined with the vertical discretisation described by Sofiev (2002) and the boundary layer scheme of Sofiev et al. (2010).Wet and dry deposition were parameterised as described in Sofiev et al. (2006).
The chemistry of ozone and related reactive pollutants was simulated using the carbon bond 4 chemical mechanism (CB4, Gery et al., 1989).However, the NO 2 analyses were produced with separate simulations employing the DMAT chemical scheme of Sofiev (2000).This follows the setup used in operational air quality forecasts with the SILAM model, where the two model runs are necessary since the primary and secondary inorganic aerosols are only included in the DMAT scheme.The SILAM model has been previously applied in simulating regional ozone and NO 2 concentrations (Huijnen et al., 2010;Langner et al., 2012;Solazzo et al., 2012), for global-scale aerosol simulations (Sofiev et al., 2011), as well as for simulating emission and dispersion of allergenic pollen (Siljamo et al., 2013).The daily, European-scale air quality forecasts contributing to the MACC-II project are publicly available at http://macc-raq.gmes-atmosphere.eu.
In this study, the model was configured for a European domain covering the area between 35.2 and 70.0 • N and −14.5 and 35.0 • E with a regular long-lat grid.The vertical discretisation consisted of eight terrain-following levels reaching up to about 6.8 km.The vertical coordinate was geometric height.The model was driven by operational ECMWF IFS forecast fields, which were initially extracted in a 0.125 • long-lat grid and further interpolated to the CTM resolution.Chemical boundary conditions were provided by the MACC reanalysis (Inness et al., 2013), which uses the MOZART global chemistry-transport model.
The emissions of anthropogenic pollutants were provided by the MACC-II European emission inventory (Kuenen et al., 2014) for the reference year 2009.The biogenic isoprene emissions, required by the CB4 run, were simulated by the emission model of Poupkou et al. (2010).
Three sets of SILAM simulations were carried out in this study.First, the background and observation error covariance matrices were calibrated using 1-month simulations for June and December 2011.The calibration results were used in reanalysis simulations covering the year 2012.Finally, a set of 72 h hindcasts was generated for the period between 16 July and 5 August 2012, to evaluate the forecast impact of assimilation.The hindcasts were initialised from the 00:00 UTC analysis fields.The timespan included an ozone episode affecting parts of southern and western Europe (EEA, 2013).The reanalysis and hindcasts use identical meteorological and boundary input data, and hence, the hindcasts only assess the effect of chemical data assimilation.
The analysis and forecast runs were performed at a horizontal resolution of 0.2 • .The setup for calibrations runs (June and December 2011) was identical except that a coarser horizontal resolution of 0.5 • was chosen in order to reduce the computational burden.The model time step was 15 min for both setups.

Observations
This study uses the hourly observations of NO 2 and O 3 at background stations available in the AirBase database (http: //acm.eionet.europa.eu/databases/airbase/)maintained by the European Environmental Agency.Separate subsets are employed for assimilation and evaluation.
Two sets of stations were withheld for evaluation.The first set, referred to here as the MACC set, had been used in the regional air quality assessments within the MACC and MACC-II projects (Rouïl, 2013, also Curier et al., 2012).The second set consisted of the stations reported as EMEP stations in the database.The MACC validation stations included about a third of the available background stations for each species, and were chosen with the requirement to cover the same area as the assimilation stations.The EMEP network is sparser and has no particular relation to the assimilation stations.It can be noted that the EMEP stations included in AirBase do not comprise the full EMEP monitoring network.
The in situ data are used for assimilation and evaluation under the assumption that they represent the pollutant levels in spatial scales resolved by the model.We expect this assumption to be violated, especially at many urban and suburban stations due to local variations in emission fluxes.For this reason, only rural stations were used for evaluation of the 2012 reanalysis.The NO 2 assimilation set also excluded both urban and suburban stations.For ozone, the data from suburban stations were assimilated, however, the observation errors were assessed separately for suburban and rural stations, as outlined in Sect.3. The station sets are presented on a map in Fig. 1.
The statistical indicators used for model evaluation were correlation, mean bias and root mean squared error (RMSE).Since air quality models are frequently used to evaluate daily maximum concentrations, the indicators were evaluated separately for the daily maximum values.

The 3D-Var assimilation
In the 3D-Var method, the analysis x a minimises the cost function: where x b is the background state, y is the vector of observations, and H is the possibly nonlinear observation operator.
The uncertainties of the background state x b and the observations y are described by the background and observation error covariance matrices B and R, respectively.In this study, the control variable x consisted of the three-dimensional airborne concentration for either NO 2 or ozone.The m1qn3 minimisation code (Gilbert and Lemaréchal, 1989) was used for solving the optimisation problem Eq. ( 1).
For the surface measurements, the operator H was linear and consisted of horizontal interpolation only, since the surface concentrations were considered to be represented by the lowest model level.Following the hourly observation frequency, the analysis was performed every hour followed by a 1 h forecast.The forecast provides the background field for the subsequent analysis.
In the current study, only a single chemical component was assimilated in each run.Since O 3 is not a prognostic variable in the DMAT scheme, it cannot be assimilated into the NO 2 simulation.Assimilating NO 2 observations into the CB4 simulation would be technically feasible; however, simultaneous assimilation of NO 2 and O 3 would require care due to the strong chemical coupling between the species.The background and observation error covariance matrices would also need to be jointly estimated.

Background and observation error covariance matrices
The numerical formulation of the BECM in the current work follows the assumptions made by Vira and Sofiev (2012).
We assume that the background error correlation is homogeneous in space, and its horizontal component is described by a Gaussian function of distance between the grid points.Furthermore, we assume that the background error standard deviation σ b is independent of location.This allows writing  the BECM as B = σ 2 b C, where C is the correlation matrix and σ b is the background error standard deviation.
For estimation of the parameters for the covariance matrices B and R, we combined the NMC method, which is used for determining the correlation matrix C, and the approach of Desroziers et al. (2005), which is used for diagnosing the observation and background error standard deviations.
In the NMC method, the difference between two forecasts valid at a given time is taken as a proxy of the forecast error.In this work, the proxy data set was extracted from 24 and 48 h regional air quality forecasts for the year 2010.The forecasts are generated with the SILAM model in a configuration similar to the one used in this study.Since no chemical data assimilation was used in the forecasts, the differences were due to changes in forecast meteorology and boundary conditions only.The lead times were chosen to allow sufficient spread to develop between the forecasts.The forecast data were segregated by hour, resulting in separate sets for hours 00:00, 06:00, 12:00 and 18:00 UTC, and the correlations interpolated for all other times of day.
The horizontal and vertical components of the correlation matrix C were estimated separately.The horizontal correlation was determined by the length scale L, which was obtained by fitting a Gaussian correlation function to the data set.First, the sample correlation matrix C of the forecast differences was calculated.Then, the Gaussian correlation function was fitted to the empirical correlations C ij by minimising where the fitted correlation function is ) and x and y are the Cartesian coordinates for each grid point.To reduce the effect of spurious long-distance correlations due to the limited sample size, the fitting was restricted to grid points r i closer than d = 1000 km to each other.The distances, shown in Table 1, were computed for the lowest model layer.
The vertical correlation function was obtained directly as the sample correlation across all vertical columns for each time of day.As an example, the correlation matrix obtained for NO 2 at 12:00 UTC is shown in Fig. 2.
Since the NMC data set includes only meteorological perturbations, it is expected to underestimate the total uncertainty of the CTM simulations.Hence, the standard deviations were not diagnosed from the NMC data set, but instead, an approach based on a posteriori diagnostics was taken.The where E denotes the expectation, y is the observation vector and y a = H(x a ) and y b = H(x b ) are evaluated from the analysis and background fields, respectively.The background error covariance matrix cannot be uniquely expressed in observation space.However, assuming that each observation only depends (linearly) on a single model grid cell (i.e.horizontal interpolation is neglected), then: E[(y (i)  a − y The identities (3) and (4) hold for an ideally defined analysis system, provided that the background and observation errors are normally distributed and assuming the observation operator is not strongly nonlinear.
Furthermore, Eqs. ( 3) and ( 4) can be used to tune the parameters σ obs and σ b by means of fixed point iteration.First, a set of analyses is produced using initial parameter values.Then, the left-hand sides of Eqs. ( 3) and ( 4) are evaluated as averages over the analyses, resulting in new parameter values.The procedure is then repeated using the updated σ b and σ obs to produce a new set of analyses.In this work, we stopped the iteration when the RMSE at validation stations was no longer improving.We chose this criterion to avoid overfitting the parameters to the calibration data.
In this work, the observation error covariance matrix R was assumed diagonal.The initial values for σ obs and σ b were set to 11.2 and 20.6 µg m −3 for O 3 , and 4.0 and 8.0 µg m −3 for NO 2 .The values corresponded to typical mean-squared errors for a free-running model, which were attributed to the model and observation error variances in the ratio of 80/20, respectively.The standard deviations, together with the correlation matrices obtained with the NMC procedure, were then employed in the iterations to calculate a set of hourly analyses for the two calibration periods spanning June and December 2011.
The choice of calibration periods representing both winter and summer conditions was motivated by the strong seasonal variations in both O 3 and NO 2 .Both σ obs and σ b were segregated by hour, while for O 3 σ obs was also evaluated separately for suburban stations.For the reanalysis of year 2012, the standard deviations, obtained separately for June and December, were interpolated linearly for all other months.
Finally, the overall consistency could be evaluated by checking the identity (Ménard et al., 2000) where χ 2 = 2J (x a ) is twice the value of cost function (1) at the minimum, and N is dimension of the observation vector y.The identity (5) tests the overall consistency of the analysis and is affected by both B and R.

Results and discussion
The SILAM model was run for year 2012 with and without assimilation.Since the 3D-Var analyses require no additional model integrations in form of iterations or ensemble simulations, the hourly analyses increased the simulation runtime by only 10-15 %.
The effect of assimilation to the yearly-mean concentrations on the lowest model level is shown in Fig. 3. On average, the ozone concentrations are increased by the assimilation especially around the Mediterranean Sea, which indicates corresponding low bias in the free model run.The main changes in NO 2 levels are confined to somewhat more limited areas; in particular, areas near major mountain ranges (Alps and Pyrenees) show enhanced NO 2 levels in the analysis run.

Background and observation error covariance matrices
Refining the background and observation standard deviations iteratively both improves the consistency of the assimilation setup as measured by the χ 2 indicator (Eq.5), and improves the model-measurement comparison on the validation stations over the calibration period.However, after five iterations (for both June and December), the changes in χ 2 become slow and the validation scores no longer improve.Hence, the values for σ obs and σ b in fifth iterations were taken as the final values for 2012 reanalysis.The changes in χ 2 and model-measurement RMSE are summarised in Table 2.
The diagnosed observation and background error standard deviations for O 3 and NO 2 are shown in Fig. 4. For June, the standard deviations for ozone range between 11 and 21 µg m −3 for rural stations.For December, the diurnal variation is flatter, but the absolute values are generally not reduced, in contrast to the overall seasonality of O 3 .Especially for summertime night conditions, the values are higher than the values adopted in most of the earlier studies (Chai et al., 2007;Curier et al., 2012;Jaumouillé et al., 2012).However, the errors are comparable to the observation errors diagnosed using the CHIMERE model by Gaubert et al. (2014).The main error component is likely to be due to lack of representativeness: using the AIRNOW observation network, Chai et al. (2007) found standard deviations between 5 and 13 ppb for observations inside a grid cell with 60 km resolution.The maximum values occurred during nighttime.
The diagnosed observation and background error parameters are subject to uncertainty, since they are not uniquely determined (Schwinger and Elbern, 2010).Also, the parameters depend on the assumptions made regarding the correlation function.Nevertheless, the relative magnitude of observation errors during nighttime is interesting for interpreting the model-to-measurement comparisons.
The diagnosed background errors for ozone are between 5 and 9 µg m −3 depending on and time of day.For June, the diagnosed errors are largest between 09:00-10:00 and 21:00-22:00 UTC, which coincides with transitions between stable and convective boundary layers in summertime conditions.For December, only minor diurnal variation is observed.
The observation error standard deviation for NO 2 varies between 2.8 and 5.2 µg m −3 for rural stations.Suburban or urban stations were not assimilated for NO 2 .Contrary to ozone, the diurnal variation of background and observation errors both positively correlate with the diurnal variation of the pollutant.
The BECM and OECM were adjusted to optimise selfconsistency for 2 months in 2011.To assess the robustness of the obtained formulations, the χ 2 indicator was also computed for all analysis steps for the 2012 reanalysis simulation.
As seen in Fig. 5, the analyses using the adjusted BECM and OECM generally satisfy the consistency relation better throughout the year, when compared to the first guess values.The yearly-mean values for χ 2 are 1.05 and 0.97 for ozone and NO 2 , respectively.
Overall, the assimilation system is based on rather simplistic assumptions regarding the background and observation error statistics.In addition to computational efficiency, this approach benefits from having few tuning parameters, and the remaining parameters (σ obs , σ b and L) can be estimated using an automated procedure.As shown in the following section, the refined background and observation error definitions provide a clear improvement on analysis scores at the control stations, despite the rather limited training data sets.

Evaluation against independent observations
Tables 3 and 4 present the analysis skill scores for runs with both first guess and final BECM and OECM, and for the freerunning model with no assimilation.
In terms of correlation and RMSE, both analysis and free model runs show better performance for predicting the daily maximum than hourly values.This applies to both O 3 and NO 2 , although the difference is more marked for ozone.The opposite holds true for bias, which tends to be higher when calculated for daily maxima.
The comparison reveals a number of contrasts between the "MACC" and "EMEP" validation stations.First, the freerunning model shows better performance for NO 2 on the EMEP stations, while for ozone, the performance is better on the MACC stations.On the other hand, the data assimilation has a stronger impact on the scores for the MACC validation stations.This is especially visible in the case for NO 2 , a result that is consistent with the shorter lifetime of NO 2 compared to O 3 .
The differences largely originate from the different representativeness and coverage of the MACC and EMEP station sets.As seen in Fig. 1, the EMEP network covers the computational domain more evenly than the MACC validation stations, which are concentrated in central Europe.Since the coverage of assimilation and MACC validation stations is similar, the average impact of assimilation is stronger on the MACC than EMEP stations.
For the free-running simulations, the better performance for O 3 at the MACC stations is consistent with the geographical variations in the model skill; the densest coverage of the MACC validation stations coincides with the parts of Europe where many regional air quality models perform best for ozone (e.g.Vautard et al., 2009).The scores for NO 2 also vary by region, however, due to the shorter chemical life-time, the forecasts of NO 2 are more sensitive to unresolved variations in local emissions.This probably explains the better scores for NO 2 on the EMEP stations, since the EMEP network is specifically aimed at monitoring the background levels of pollutants.
For ozone, the assimilation had a variable effect on the model bias.While the correlation and RMSE were always improved by assimilation, the analyses have slightly larger negative mean bias (−4.6 vs. −4.0µg m −3 on MACC stations) than the free model.This is confirmed by the average diurnal profile shown in Fig. 6.However, the diurnal variation of analysis errors is flatter, and the strongest bias no longer coincides with the afternoon hours, when the highest O 3 concentrations are typically observed.
For NO 2 , the analyses have only slight negative bias (−0.38 µg m −3 ) on the MACC stations, which turns positive (about 1 µg m −3 ) for the more remote EMEP sites.As seen in Table 4, the difference between the station sets is similar to that of the free-running model.Given the differences between the MACC and EMEP station sets, this suggests that the model overestimates the lifetime of NO 2 , which in turn results in the positive bias in the analyses.The long lifetime of NO 2 in the SILAM DMAT chemistry scheme was also noticed by Huijnen et al. (2010).
The analysis scheme assumes an unbiased model, and hence, the negative bias present in the free-running simulations is reduced but not removed in the analysis fields.The assimilation setup including tuned OECM and BECM produces more biased analyses compared to the first guess setup, as seen in Fig. 6.This is a consequence of the differences between the diagnosed and first guess background and observation error standard deviations.Contrary to the tuned setup, the first guess attributes most of the model-observation discrepancy to the background error, which results in stronger increments towards the observed values.Consequently, the analysis bias is smaller.However, the tuned assimilation setup has consistently better RMSE and correlation than the first guess assimilation setup.Since the analysis bias is mainly a consequence of a bias in the forecast model, the bias should be addressed primarily by improving the model.As shown by Dee (2005), model biases can, in principle, also be addressed by the assimilation system.However, a possible bias correction scheme should be implemented with care, since due to representativeness errors, also observational biases could arise.
In addition to computing the regular statistical indicators for daily maxima, we evaluated the hit rates (the number of correctly predicted exceedances divided by the number of observed exceedances) for the 180 µg m −3 threshold for O 3 , with and without assimilation.Assimilation also turns out to improve the hit rate, albeit only slightly: from 0.25 to 0.26 on average for rural MACC validation stations, and from 0.13 to 0.15 for EMEP stations.If the averaging is restricted to the stations with more than 10 exceedances during 2012, the values change from 0.32 to 0.36 for MACC and from 0.21 to 0.43 for the EMEP stations.Obviously, the hit rates are sensitive to the low bias in the daily maxima.
For NO 2 , a specific source of observational errors is due to the molybdenum converters used in the chemiluminescence technique, which is the most common measurement technique for monitoring NO 2 .As discussed by Dunlea et al. (2007) and Steinbacher et al. (2007), this technique is subject to positive interference by the NO z species such as PAN, HNO 3 and HONO.
The interference can lead to overestimation of NO 2 by up to factor of two, however, the error varies by location and time, and may depend on the features of the instrument (Steinbacher et al., 2007).We estimated the magnitude of this effect from the free-running CB4 simulation.On most continental EMEP sites, the contribution of the NO z species to the total NO z + NO 2 was about 10-20 % of the simulated yearly mean.However, for a few sites the contribution could reach 50 %.
The O 3 and NO 2 observations were assimilated into separate model runs.Assimilation of O 3 had only a minor influence on NO 2 in the CB4 simulation; however, the mean bias was reduced by about 5 % on average for the MACC validation stations.Because the DMAT simulation does not include ozone as a tracer, the impact of NO 2 assimilation on ozone fields was not evaluated in this study.

Forecast experiments
In order to quantify the usefulness of data assimilation in forecast applications, a set of simulations without data assimilation were generated using the analysis fields at 00:00 UTC as initial conditions.The forecast experiment covered the time period between 16 July and 5 August 012.
The effect of chemical data assimilation on forecast performance was assessed as a function of the forecast lead time.Figures 7 and 8 present the correlation and bias for the O 3  and NO 2 forecasts, respectively, and compare them to the corresponding indicators for the analyses and the control run.
For ozone, the forecast improvements due to data assimilation were largely limited to the first 24 h of forecast.Also, the forecast initialised at 00:00 UTC from the analysis shows a larger negative bias for the daytime than the free model run.This is a result of the corresponding nighttime positive bias of the free model run.The bias is effectively removed in the 00 analysis; however, the subsequent forecast is unable to recover the level observed during daytime.The correlation coefficient during daytime is nevertheless improved slightly (from 0.75 to 0.78) by initialising from the analysis.While the forecast shows somewhat reduced positive bias for the hours between 18 and 30 (i.e. the following night), the subsequent daytime scores are already almost unchanged by assimilation.The results in Fig. 7 are computed for the MACC station network; a similar impact is observed at the EMEP stations.
Due to the shorter chemical lifetime, the effect of initial condition on forecasts of NO 2 can be expected to fall away more quickly than for ozone.This has been confirmed in the previous works based on the assimilation of data from the ozone monitoring instrument (OMI).Under summer conditions, Wang et al. (2011) found assimilation to provide no improvement in RMSE with regard to surface observations, while Silver et al. (2013) reported the NO 2 concentration to relax to its background values within 3-4 h.
In the forecast experiments performed within this study, the effect of assimilation on NO 2 forecast scores was limited to the first 6 forecast hours, which coincides with the nighttime in most of the domain.Hence, at least under the photochemically active summertime conditions, the analyses are only marginally useful for improving forecasts of NO 2 .
The forecast for short-lived pollutants like NO 2 is poorly constrained by the initial condition, because the boundary layer concentrations become driven mainly by local emissions, chemical transformations and deposition.This limits effectiveness of any assimilation scheme based updating only the initial condition.A possible way to extend the forecast impact is to include more persistent parameters, such as emission rates, into the state vector.This has been demonstrated by Elbern et al. (2007) for forecasting an ozone episode.In general, such an approach requires that the obtained a posteriori emission rates can be extrapolated to the forecast window, and that the assimilation scheme is able to correctly attribute the observed discrepancies to the uncertain parameters.

Conclusions
An assimilation system coupled to the SILAM chemistry transport model has been described along with its application in reanalysis of ozone and NO 2 concentrations for year 2012.Furthermore, the impact of using the O 3 and NO 2 analyses to initialise forecasts has been assessed for an ozone episode occurring in July 2012.
The assimilation consistently improves the modelmeasurement comparison for stations not included in the assimilation.For daily maximum values, the correlation coefficient is improved over the free running model from 0.8 to 0.9 for O 3 and from 0.53 to 0.63 for NO 2 on rural validation stations.The respective biases are also decreased, however, a bias of −7.4 µg m −3 remains in the O 3 analyses due to a negative bias in the free-running model.
During a 3-week forecast experiment, initialising the forecasts from the analysis fields provided an improvement in ozone forecast skill for a maximum of 24 h.For NO 2 , the improvement was limited to a window of 6 h.The findings for NO 2 are similar to the results published in previous studies (Silver et al., 2013;Wang et al., 2011).
The diagnosed observation error standard deviations for ozone have a strong diurnal variation, and reach up to about 21 µg m −3 at nighttime.These values are higher than usually assumed in chemical data assimilation, but agree well with the results obtained by Gaubert et al. (2014) with similar diagnostics.
The 3D-Var based assimilation has a low computational overhead.This makes it especially suitable for reanalyses in yearly or longer time scales, as well as for high-resolution forecasting under operational time constraints.Future studies will include more accurate characterisation of station representativeness, as well as further investigation of model biases for O 3 .

Figure 1 .
Figure 1.The stations networks used for assimilation and validation for O 3 (left) and NO 2 (right).The assimilation stations for O 3 include rural and suburban stations, for NO 2 only rural stations.For validation, only rural stations are shown.The red and blue colours refer to the MACC validation and EMEP stations subsets.

Figure 3 .
Figure 3. Yearly mean concentration (µg m −3 , left-hand panels) on lowest model layer and difference (assimilated -not assimilated, righthand panels) due to assimilation of O 3 (top panels) and NO 2 (bottom panels).

Figure 4 .
Figure 4. Diagnosed background (dashed) and observation error (solid lines) standard deviations (µg m −3 ) on rural stations for O 3 (left) and NO 2 (right).Red lines correspond to the calibration made for June 2011, blue lines correspond to calibration for December 2011.

Figure 5 .
Figure 5.The χ 2 /N obs consistency indicator for hourly analyses of O 3 (left) and NO 2 (right).The values in blue and green are shown for the first guess and final assimilation setups, respectively.Note the different scales for O 3 and NO 2 .

Figure 6 .
Figure 6.Diurnal variation of model bias (µg m −3 ).The first guess assimilation setup is shown in red and the final setup in blue.The reference run with no assimilation is drawn in green.The values are shown for the rural MACC validation stations and averaged over each day of year 2012 and over the stations.

Figure 7 .
Figure 7.The model bias (µg m −3 ) and correlation for O 3 at the MACC validation stations as a function of forecast length (blue lines).The corresponding indicators, the analyses (black) and control run (green), are shown averaged by time of day and replicated over the forecast window.

Table 1 .
Correlation length scales L (km) diagnosed from the NMC data set.

Table 2 .
The χ 2 /N consistency indicator and RMSE on rural MACC validation stations during the first and fifth iterations for tuning the observation and background error standard deviations.

Table 3 .
Comparison of performance indicators for ozone in the 2012 reanalysis.The scores are given for station sets "MACC" and "EMEP" as defined in Sect.2.2.For the analysis runs, scores are shown for the different background error covariance matrices discussed in Sect.3. The unit of bias and RMSE is µg m −3 .

Table 4 .
Comparison of performance indicators for NO 2 in the 2012 reanalysis.The station sets MACC and EMEP and assimilation options are as in Table3.