1. Introduction
Intercomparison studies of global aerosol models [
1] have shown that considerable differences exist in spatial distributions and temporal evolutions of aerosol fields. Bates
et al.[
2] suggested, based on an analysis of two models, that uncertainties in emissions are the largest contribution to model differences, although attempts to homogenize aerosol emissions across models [
3] have only partly reduced those differences. Uncertainties in emissions are certainly large: Bates
et al.[
2] put regional emission differences at ∼30% for the precursor gas SO
2, ∼300% for carbonaceous aerosol and ∼500% for either dust or sea salt. Granier
et al.[
4] compared various emission databases for anthropogenic and biomass burning aerosol and found smaller differences in global emissions of ∼30% for either black carbon or precursor gas SO
2. Regionally, however, much larger differences of 100% to 300% existed. Textor
et al.[
1] showed that global emission variations among models (which is a likely underestimate of the uncertainties) amounted to 50% for dust, 200% for sea salt and ∼25% for other aerosols.
Emission inventories for fine aerosol (sulfate precursor gas, nitrate and various carbons) are at best monthly averages and often yearly averages (forest fire emissions being the exception), further increasing daily emission errors. In contrast, emissions for coarse mode aerosol (dust and sea salt) are based on parameterisations (in, e.g., surface wind speed), which may not be suited to the coarse spatial resolution of global models. Because of these uncertainties in the emission inventories and parameterisations, it is worthwhile to explore whether they may be improved through systematic comparison of simulations and observations.
A few attempts at top-down estimation of aerosol emissions have been made before. Zhang
et al.[
5] assimilated TOMS aerosol index in an optimal interpolation scheme to optimize monthly biomass smoke emissions. Hakami
et al.[
6] assimilated black carbon mass estimates in a 4D-VAR system to estimate daily emissions of black carbon mass measurements during ACE-Asia. Dubovik
et al.[
7] used a generalized least squares fitting algorithm to optimize fine and coarse mode aerosol emissions according to MODIS observations. Yumimoto
et al.[
8] employed a 4D-VAR system to assimilate Asian LIDAR observations in a regional model and optimized dust emission. Finally, Sekiyama
et al.[
9] assimilated CALIOP observations with an ensemble Kalman filter and optimized dust and sulfate emissions. Recently, Huneeus
et al.[
10] also used an ensemble Kalman filter to optimize monthly emissions of all major aerosol species.
Atmospheric aerosol mixing ratios at a time
t + Δ
t are the result of the initial mixing ratios at time
t and the emissions from
t to
t + Δ
t. The initial conditions are themselves dependent on earlier initial conditions and emissions. Disentangling the impact of initial conditions and emissions is a challenge, especially when those emissions are time-varying. If Δ
t is small the relative impact of initial conditions (and its errors) on estimated emissions will be very large while when Δ
t is large it is unlikely that the emissions are constant. In this paper, we propose to use an ensemble Kalman smoother to overcome these problems. This technique is similar to techniques devised and tested by Peters
et al.[
11] and Feng
et al.[
12] for the estimation of CO
2 emissions. Emissions during a period Δ
t will be estimated iteratively from later observations during a period Δ
T =
nΔ
t, n = 1, 2, 3, 4, . . ..
To our knowledge, this is the first paper that attempts a top-down estimation of aerosol emission using an ensemble Kalman smoother, which should solve the conflation of initial and boundary conditions discussed in the previous paragraph. In addition, we use (indirect) observations of both aerosol quantity (AOT) and aerosol size (AE from AERONET). Dubovik
et al.[
7] and Huneeus
et al.[
10] also used size information (MODIS fine mode fraction, in their case).
In the next section, we will explain the Kalman smoother technique. Section 3 shows the new emissions and the improvement in modelled AOT and AE. Section 4 discusses the new emission fluxes in the light of AEROCOM findings. Finally, we will summarize our work.
2. Method
Our emission estimation is very similar to techniques devised by Peters
et al.[
11] and Feng
et al.[
12] for CO
2 emission estimation. Three components are required to estimate emissions top-down: a model to forecast aerosol observations based on assumed emissions, actual observations and a data assimilation technique to estimate emissions from observations.
We use the global AGCM with aerosol module MIROC-SPRINTARS v3.54 [
13,
14] to establish a relationship between assumed emissions and aerosol mixing ratios. This model is run at a resolution of T42 with 20 σ–layers and has its meteorology nudged towards NCEP reanalysis values. It considers the emission, transport and deposition of the major aerosol species—dust sea salt, sulfate and carbons. Emissions of dust and sea salt are based on parameterisations of (primarily) wind speed (note: we use the Monahan parameterisation for sea salt, [
15]) while the emissions for carbons are given by inventories. The sulfate aerosol originates from the precursor gas SO
2, which is oxidized in the gas phase by OH and in the liquid phase by H
2O
2 and O
3 (for these gases, climatologies are used). SO
2 itself derives from direct emission (based on inventories) and conversion of DMS (whose emission over sea is based on a parameterisation).
The observations are AERONET level 1.5 direct sunlight retrieval AOT at 674 nm and AE at 870/440 nm, as well as MODIS Terra level 2 AOT at 550 nm over ocean. No MODIS observations over land were used as we wanted to use those for independent validation. Assumed errors in the observations are determined from an independent combination of measurement error and spatio-temporal representation error. The AERONET observations were averaged over 3 hours. Measurement errors in AERONET AOT are assumed 0.015, from which one can calculate associated AE measurement errors. AERONET representation errors are according to the analysis in [
16] (typically 10% of AOT). In the current study all AERONET sites are used. All MODIS observations within a three hour window are averaged over 0.5 by 0.5 degree boxes. We applied Zhang & Reid’s QA filtering to the MODIS data [
17], but not their empirical corrections for biases due to wind speed and fine mode fraction. Thus our MODIS observations are screened for cloudiness but still contain small biases. MODIS measurement errors are according to [
17] and the representation error is the standard error in AOT over 0.5 by 0.5 degree.
Figure 1 shows the distribution of AERONET sites and MODIS Aqua observations (to be used for validation) used in this study.
In this paper, we will use Bayesian inference to estimate emissions. The central idea is that emissions (and all derived variables like mixing ratios) are represented not only by a value but also by an (estimated) error in that value. The value will represent our best estimate of the emission for a certain species at a certain time and location, while the error represents the associated uncertainty. This allows us to use conditional probabilities to estimate emissions. In our case, we will use a fixed-lag Kalman smoother as a particular implementation of Bayesian inference (see [
18] for mathematical details). This fixed-lag Kalman smoother is in essence a Kalman filter that iteratively estimates emissions. Both Kalman smoother and filter assume a linear model and Gaussian errors (in model and observations). As we nudge our model to reanalysis data, these conditions are largely satisfied. For practical reasons, the emissions and their errors will be represented by an ensemble of emissions that will yield an ensemble of mixing ratios based on independent runs of MIROC-SPRINTARS. For the mathematical details of the ensemble Kalman filter, see [
19]. It uses the statistical information present in the ensemble, combined with observations, to improve upon our estimates for the emissions.
The emission ensemble is defined as in Schutgens
et al.[
16] by multiplying the standard emissions with perturbed factors of a specified mean and spread. Dust and sea salt emissions are perturbed around their parameterised values (dependent mainly on wind speed), emissions for the sulfate precursor gas and carbon aerosol are perturbed around their inventories. Here we do not use the original multi-year averaged monthly data as in [
13] but we will use a yearly average to ensure the largest number of sources is present in our initial emissions. The current Kalman smoother can estimate new emission strengths but leaves locations of sources intact. For the initial ensemble (32 members), we assume that the perturbed factors have a mean of 1 and a spread of 1 and follow a log-normal distribution. Perturbed factors for dust and sea salt are randomly generated, independently from each other and from carbons and SO
2. The emissions for the fine mode aerosol are perturbed in identical fashion, although their initial emission inventories differ.
New estimates for the perturbed factors are obtained from the solution to the Kalman equation
where
xf are the initially assumed perturbations (the mean of the so-called prior ensemble) and
xa the improved perturbations (the mean of the so-called posterior ensemble) for different species at different locations. The operator
H relates the prior emissions to forecasts of simulated observations. It consists of MIROC-SPRINTARS and an observation operator. This observation operator simulates aerosol observations based on the forecasted mixing ratios and their assumed scattering properties. The matrix
R contains estimates of the errors in the actual observations
y (including spatial representation errors), while
Pf describes the assumed errors in the perturbations
xf (the covariant of the prior ensemble). It is customary to formally include model transport errors in
R, but just like all researchers before us, we assume those errors to be negligible for reasons of simplicity. Note that the (propagated) errors in the improved perturbations are defined by
Pa.
So far, we have only estimated new perturbations based on a single ensemble forecast
H ·
xf and a single set of observations
y. Suppose that the latter contains all available observations during a period Δ
t. Those observations will contain some (!) information about the emissions before that time. But observations at a later time will also contain (some) information about those emissions. It then makes sense to estimate emission iteratively, see
Figure 2. We start by assuming initial perturbations and calculate an ensemble of forecasts for a longer period Δ
T > Δ
t (forecast 1). At the end of the forecast, we use observations to estimate new perturbations based on
Equations (1) and
(2) (analysis 1). Those new perturbations are then used for a new and extended forecast (forecast 2) that (temporally) partly overlaps the earlier forecast 1. To extend forecast 2, again perturbations have to be assumed for the most recent times. This new forecast is then combined with new observations to obtain a new analysis (2). Analysis 2 (temporally) partly overlaps with analysis 1 with two exceptions. Firstly the new observations allow estimation of the perturbations at the most recent times. Secondly, perturbations at the most distant time in analysis 1 are no longer estimated. Although this is done primarily for practical reasons of limitations in CPU time and memory, it is not a limitation in the technique. If Δ
T is smaller than the aerosol residence time, we still improve our perturbations but do not use all observations that contain relevant information. If Δ
T is much larger than the aerosol residence time, we have used all available information as much as we possibly can. In the current experiments, we will use Δ
t = 2 and Δ
T = 6 days. Partly this choice is based on limited computer resources. The typical aerosol residence timescale is less than 6 days, however.
A Kalman filter would use ΔT = Δt and would not recalculate the forecast for the new emissions. In the Kalman smoother, new emission estimates are made ΔT/Δt times using different observations (but subsequent in time). Note that the temporal resolution of the perturbed factors is Δt.
3. Results
The purpose of the assimilation system is to improve estimates of the emissions, here for January 2009. We will test the system by running a forecast with the new emissions and compare them against observations. Initially we will compare against the observations that we assimilated (AERONET AOT & AE as well as MODIS Terra AOT over ocean) to see how readily the system adjusts to those observations. Next we will validate the new emissions against independent observations by MODIS Aqua AOT over both ocean and land. Over land we will thus validate against a very different dataset, which implies that both observational issues (measurement errors, sampling) and data assimilation issues will influence the comparison. Over ocean, we effectively remove the impact of observational issues as MODIS Terra and Aqua use very similar sensors and algorithms. Note however that Terra and Aqua do have different temporal sampling. Comparison against observations will be made only for the last three weeks in January 2009 (the Kalman smoother requires a spin-up time of about ΔT).
Global AOT distributions for the standard and new emissions are shown in
Figure 3. Clearly, the new emissions have a large impact. Some statistics of the comparison with observations are given in
Table 1 where we present information on the overall bias, correlation and regression slope. We see that AOT and AE due to the new emissions yield smaller biases and better regression slopes than the standard emissions. Correlations improve but not a lot as the observations themselves are already quite noisy.
In
Figure 4, a density plot of the distribution of co-located MODIS Terra over ocean and experiment AOT is shown for forecast runs with either the standard or new emissions. We present density plots with both linear and logarithmic axis to show that results do not differ depending AOT magnitudes. For the standard emissions, AOT over ocean is generally estimated too low with a poor correlation with the observations. This bias is removed for the new emissions, just as correlation with the observations improves. MODIS Terra data were used to obtain the new emissions, so this is not a proper validation. The remaining random errors in AOT seem acceptable given typical random errors in the MODIS AOT product (black lines).
Figure 5 shows density plots of the distribution of co-located AERONET observations and experiment data (both AOT and AE). The picture is noisier than for MODIS Terra due to the lower observational count, but large biases in both AOT and AE are clearly reduced when using the new emissions. Again AOT for the standard emissions is negatively biased, but AE for the standard emissions is positively biased. Remaining random errors in AOT are relatively large compared with typical AERONET AOT errors but they appear acceptable (within 2 σ, this AERONET observational error includes both a retrieval error and a representation error). AERONET AOT and AE were used in the estimation of the new emissions and
Figure 5 cannot be considered a validation.
In
Figure 6 we show the results for the Irkutsk AERONET station where the new emissions yield far better AOT than the standard emissions. This is obvious not only in the mean value of AOT but also in its temporal evolution, which shows larger daily variations for the new than the standard emissions (note that the new emissions for fine mode aerosol have a 2 day time resolution). Irkutsk was purposefully selected to show the possibilities of our system. Again, this is not a proper validation (AERONET was used to estimate the emissions), but it serves to show how well the assimilation system can adapt the emissions to the observations.
Figure 7 shows biases
per site for either the standard or new emissions experiments. The bias is here defined as the median deviation from the observations during the last three weeks. When a site falls between the solid
y = ±
x lines, the new emissions reduce the bias at that site. The dotted lines indicate the amount of reduction in this bias. We see how the new emissions positively impact AOT and AE. In the case of AOT, there is a marked improvement, with a lot of sites showing a strong bias reduction. For AE, the figure also indicates an improvement but less clearly. Observational errors in AE are much larger than in AOT and consequently the Kalman smoother will allow less influence on the emission estimates from AE assimilation. There are some sites where the new emissions have no or an adverse effect even though these AERONET observations were used to estimate the new emissions. These sites are listed by name in the figures. When the model has inherent errors and multiple nearby observations are used to estimate emissions, irresolvable contradictions may arise. This is most obvious for Solar Village and Karachi, where the best the smoother can do is a large negative AE bias for Karachi and a large positive AE bias for Solar village in an attempt to reduce overall bias. Note that Camaguey, Lannion, Cordoba CETT, Avignon and Carpentras all have a very low AOT and consequently very large errors in observed AE.
For a proper validation, we will use MODIS Aqua data (
Figure 8). It is interesting to split this analysis into two parts: one over ocean and one over land. Over ocean, we previously assimilated MODIS Terra AOT so a comparison against MODIS Aqua is not entirely independent. After all, similar sensors and retrieval algorithms are used. This is not necessarily a limitation to the validation, as it allows us to study the emission estimation without concerning ourselves with observational biases (which are rather similar in Terra and Aqua). In
Figure 8, we see that the new emissions also yield significant improvement against Aqua AOT (Aqua and Terra make observations at different times).
Over ocean, we see a very similar picture as in
Figure 4. The standard emissions cause a negative bias in the forecast, which is effectively removed by the new emissions. Note that also the correlation and slope improve in a remarkably similar manner as for MODIS Terra. Observations by MODIS Aqua occur with a three hour delay from Terra and, due to cloudiness, will have different locations.
Over land, the comparison against Aqua also shows improvement in AOT due to new emissions. But as the assimilated data (AERONET) and validation data (Aqua) are quite different, there are also observational biases and sampling that affect the validation. The standard emissions yield, quite generally, an AOT that is too low (this was also seen for AERONET, see
Figure 5). This negative bias becomes more pronounced at high AOT (a feature not clearly seen in the AERONET data, probably due to different sampling). For the new emissions, the bias has mostly disappeared but now some observations show a positive bias against Aqua (around AOT ≈0.2). It can be shown that this is first of all a sampling issue (
Figure 1), with most of the positive bias observations from Aqua occurring in regions where there were either almost no AERONET sites at all (central Africa) or only few observations from the existing sites (South Asia).
Finally, we present the global fields of SO
2, sea salt and dust emission for both the standard and new emissions, see
Figure 9. As mentioned before, these new emissions result in a markedly different AOT distribution (
Figure 3). Global SO
2 emissions increase, mainly due to increases in Asia. Over Europe, there is reduction in emissions. Global dust emission reduces by only 26% in the new emissions but there are clear differences in location (and timing) of major dust outbreaks. Global sea salt emission, on the other hand, shows an almost 3-fold increase, while location (and timing) of the emission do not change much. It seems that the dependence of sea salt emission on wind speed is much stronger according to the new emission estimate (not shown). Similar maps of the new emission of carbons exist but are not shown.
4. Discussion
In
Table 2, we compare global annual emissions of the four major aerosols according to various sources. The first column is the AEROCOM mean emission according to [
1] while the second column is the diversity (standard deviation divided by mean) according to the same paper. These values serve as a reference but obviously cannot be used to validate any of our results. The third column lists emissions according to [
13], which used an older version of MIROC-SPRINTARS than we do. In particular, the sea salt emission scheme in [
13] uses the Eriksson parameterisation for sea salt concentration in the lowest model layer as opposed to the Monahan parameterisation for sea salt emission. The fourth and fifth column list the emissions used in our standard model and those estimated in this paper. It is important to note that here we estimated emissions for January 2009 but converted them to yearly emissions. Finally, the last column gives error estimates remaining in these emissions, as determined by the Kalman smoother (
Equation (2)). Values between brackets are global averages weighted by AOT.
Compared with the AEROCOM mean emissions, MIROC-SPRINTARS appears to underestimate sea salt emissions and overestimate dust emissions. The carbon emissions are, however, close to the mean (at least the Takemura (2000) results). No AEROCOM results for SO
2 gas are quoted as they are not available. Dust emissions have significant yearly variations so the difference between the emissions in [
13] and our standard run are not strange. Interestingly the new dust and sea salt emissions have moved closer to the AEROCOM mean values. In general, the newly estimated emissions have values that are not extraordinary compared with other data sources.
The increase of carbon and sea salt global fluxes is not surprising given the negative bias in AOT for the standard emissions against MODIS Terra, Aqua and AERONET (
Figures 4,
5 and
8). For the standard emissions, AE had a positive bias for low values which occurred for AERONET sites near deserts or coasts. The increase in sea salt emission corrects this.
We now turn to the remaining errors in global emissions, as estimated by the Kalman smoother. Note that these errors are determined by (1) a priori errors in the emissions, (2) transport information from MIROC-SPRINTARS, and (3) observational errors. The remaining errors can be interpreted as propagated errors after combining model and observational information. In contrast, the diversity among AEROCOM models cannot be seen as a proper error estimate. The a-priori error in our emissions was 100% but due to the assimilation of observations the posterior error is significantly lower. Interestingly, errors over land are significantly higher than over ocean. Two effects contribute to this: (1) the poorer observational sampling by the AERONET network; (2) the state of aerosol which over land is typically a combination of dust, carbons and sulfate making emission estimation more difficult than over the deep ocean where sea salt is the major contributing species. Note that there is little difference in straight global averages of the error or an error weighted by AOT. The latter of course favours regions with lots of aerosol.
Finally, in
Figure 10, we show a global graph of the relative spread (error normalised by emission) for the newly estimated sea salt emission. There is quite a bit of spatial variation in these errors, most of which can be understood in terms of spatial sampling of our observations and the nature of regional aerosol. For instance, near the north pole, the relative spread is close to 1 since we did not assimilate any observations in that region (polar night). Over the deep oceans, the spread is small as we only perturbed sea salt emissions over ocean (but not DMS) and also did not introduce any uncertainty in transport and deposition of aerosol. Consequently, it is relatively easy to link observed AOT to sea salt emissions. In comparison, large relative spread exists near the coasts of strongly emitting land areas (e.g., South America, Africa and India). Here the presence of an outflow makes it significantly harder for the Kalman smoother to estimate emissions. Interestingly, similar large spread in sea salt emission is not found round the outflows from the South-East Asian landmass or North America. We surmise that the local AERONET sites serve to constrain the emission of aerosols over land so that the outflow itself is rather well characterised. Consequently, Terra AOT over ocean can be used to good effect to estimate local sea salt emissions. To prove this idea would require additional experiments.
5. Summary
We have developed an algorithm for top-down estimates of aerosol emissions by combining information from remote sensing observations and a transport model. The new emissions are verified by using them in a one-month forecast run that shows significantly better agreement with independent remote sensing observations. In this paper, we used MODIS Terra AOT over ocean and AERONET AOT and AE for estimating emissions and MODIS Aqua (over land and ocean) for validation.
The methodology of estimating emissions is based on Bayesian inference and uses a fixed-lag ensemble Kalman smoother. Relations between emissions and observations are established through ensemble runs of the MIROC-SPRINTARS global aerosol model. In addition to practical advantages like reduced CPU time (the present system analyzes a month of data in ∼1.5 days wall-clock time on 8 processors of an SX-8) and memory requirements, the fixed-lag Kalman smoother solves the problem of disentangling the contributions of initial conditions and emissions on observations. Also, unlike 4D-VAR methods, the Kalman smoother provides error estimates on the new emissions.
Our algorithm is one of a few algorithms for top-down estimates of aerosol emission and the first to use an ensemble Kalman smoother for an aerosol transport model. Our paper is also one of a few papers that include information about aerosol size (AE in our case) in the assimilation. This also seems to be the first paper that shows that the new emissions substantially improve spatial and temporal AOT and AE simulation in a forecast run. Compared with the standard emissions, the AOT from the new emissions show significantly smaller bias (fivefold reduction), slightly higher correlation (0.6 instead of 0.5) and significantly better regression slopes (0.8 instead of 0.3) when compared with independent Aqua AOT observations.
For January 2009, the new global aerosol (pre-cursor) emissions are substantially different from the standard emissions: for sea salt, 9073 Tg/yr (was 3145 Tg/yr), for dust, 3244 Tg/yr (was 4470 Tg/yr), for carbons, 136 Tg/yr (was 83 Tg/yr), and finally for SO
2 gas, 219 Tg/yr (was 145 Tg/yr). These values are still within the range of emissions used in the AEROCOM model intercomparison [
1]. The remaining uncertainties in the new emissions are estimated to be 78% for sulfate and carbons, 62% for dust and 18% for sea salt. Sparsity of observations over land explains the larger errors for sulfate gas, carbons and dust when compared with sea salt.
The new sea salt emissions are not only characterised by significantly larger values but also a more pronounced latitudinal dependence. Inspection of time-series suggests that this is the result of an underestimation of the response to wind speed by the original emission parameterisation (or an underestimation of 10 m wind speeds by the model). The new dust emission shows significant differences in its spatial (and temporal) distribution from the standard emission although the global monthly flux is only reduced by 26%. The new emission inventory for fine mode aerosol is foremost characterised by its high temporal resolution (2 days as opposed to 1 month in the original emissions), although the limited sampling of assimilated observations over land imply that often no observations are available.
We feel it is important to point out that any emission estimate will only be as good as the transport model that is used. More research into the quality of simulation of aerosol is needed.
The system we present in this paper can be easily expanded, thanks to its ensemble approach. Points of particular interest include improved prescriptions for the emission ensemble, perturbation of aerosol sinks in the ensemble and increased observational coverage (MODIS or MISR observations over land, MODIS AE observations over ocean, CALIOP observations of the free troposphere). It may also be interesting to include direct observations of precursor gases like SO2.