Adjoint inverse modeling of a CO emission inventory at the city scale : Santiago de Chile ’ s case

Adjoint inverse modeling of a CO emission inventory at the city scale: Santiago de Chile’s case P. Saide, A. Osses, L. Gallardo, and M. Osses Departamento de Ingenieŕıa Mecánica, Universidad de Chile, Santiago, Chile Departamento de Ingenieŕıa Matemática and Centro de Modelamiento Matemático (CNRS UMI 2807), Universidad de Chile, Santiago, Chile Departamento de Geof́ısica, Universidad de Chile, Santiago, Chile Received: 10 February 2009 – Accepted: 15 February 2009 – Published: 10 March 2009 Correspondence to: P. Saide (psaide@ing.uchile.cl) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
Emission inventories (EIs) are a compilation of emission estimates classified according to source types, processes, species, etc.They constitute key-tools for air quality management, allowing the assessment of responsibilities and identifying measures to curb air quality degradation.Air quality control policies generally start addressing acute problems (extreme pollution events), and aim at reducing emissions of primary pollutants Figures

Back Close Full Screen / Esc
Printer-friendly Version Interactive Discussion from stationary sources.As the attainment objectives become more ambitious (e.g., long-term health effects) and measures must deal with the more elusive control of secondary pollutants and mobile sources, the need for more reliable and accurate EIs becomes apparent.EIs have large uncertainties that hamper our ability to diagnose and forecast air quality and its impacts (Lindley et al., 2000).In the case of emissions from mobile sources uncertainties have multiple origins, ranging from random errors in measurements devices when establishing emission factors to systematic errors due to incomplete statistics (e.g., fleet distribution), and poorly constrained processes (e.g., gas-to-particle transformation).
A way to reduce the uncertainties referred to earlier is to use inverse modeling or data assimilation.There are in fact numerous applications of inverse modeling to enhance emission inventories, from global (P étron et al., 2002, e.g.,) to continental scales (e.g., Elbern et al., 2007) and local or urban scale (e.g., Qu élo et al., 2005).Many of them use a variational approach, often called 3-D-Var in meteorological data assimilation (Kalnay, 2002), which was also chosen in this work.However, modifications were made to apply this method for the city scale that includes solving the minimization function with constraints to assure positive solutions, add an error balancing factor and add a term to address the co-location of sources and observations.The latter is a major problem when applying inverse modeling at the city scale.This work aims at presenting and applying a methodology to improve an EI of carbon monoxide (CO) or any other linear atmospheric tracer over a city with observations from an air quality monitoring network with hourly observations.For that purpose we apply a state-of-the-art dispersion model previously used under diverse weather conditions (Polyphemus, Mallet et al., 2007).Regional and mesoscale weather simulations were conducted followed by direct CO dispersion runs feeding Polyphemus with the weather fields and with the EI.After comparing models outputs (meteorological and dispersion) against available observations, an inverse modeling approach was applied to correct the CO inventory, both in space and time.Finally, the optimized emissions were validated for a period out of the assimilation window.The case analyzed corresponds to Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
Santiago de Chile, a six million inhabitant city exposed to severe air pollution, but the method is general and it can be applied elsewhere.
In the next section, we describe observations and set up of the models used.Section 3 describes weather and direct dispersion simulations using the a priori inventory.The inverse methodology and simulations are presented in Sect. 4. We summarize our results and allude to future work in Sect. 5.

Meteorological observations
In our analysis we consider synoptic data from the reanalysis from the National Centers for Environmental Predictions (NCEP, Kalnay et al., 1996), and soundings (00:00, 12:00 UTC) collected by the Chilean Weather Office (DMC) at Santo Domingo (33.39 • S, 71.37 • W, 75 m.a.s.l) about 100 km south west of Santiago.In addition, we analyze hourly averaged meteorological data (wind speed and direction, temperature and relative humidity) sampled at the air quality and meteorological stations in Santiago (Fig. 1).Finally, we include observations collected at Santiago's international airport (DMC), in which in addition to wind, temperature and humidity, surface pressure and cloud cover are reported.
According to the reanalysis (not shown), January 2002 is characterized by the predominance of the Pacific high, particularly during 2-11 January.Between 12 and 20 January a trough disrupts the zonal flow in the mid troposphere.The zonal fl ow is recovered for the rest of the period.These changes are also apparent in the observations at the sounding station (Fig. 2) and at the airport (not shown).In fact, during 15 through 20 January there is a weakening of the subsidence inversion and a strengthening of the vertical mixing below 1000 m.a.s.l.Previously, around 10 through 15 January there is evidence of coastal low, which induces easterly winds over the Santiago basin as well as warmer temperatures and dryer air, generally resulting in a higher pollutant loading Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion (Gallardo et al. 2002;Garreaud et al., 2002).The latter is partially counteracted though by the strong radiatively driven mixing allowing a substantial development of the mixing layer.

CO observations in Santiago
Since the late 1990s, the Metropolitan Sanitary Authority (www.asrm.cl) in Santiago maintains an air quality monitoring network, whose hourly averaged validated data are available over the internet.This network consists of eight stations that are distributed over Santiago in order to represent different exposure conditions in the city (Fig. 1), not necessarily dispersion patterns.CO, along with other criteria pollutants, is measured with nominal instrumental precision of ±0.1 ppm and detection limit of 0.04 ppm.The samplers are subject to zero-span validation procedures and a daily record sheet for each station is kept (V.Berríos and P. Villavicencio, personal communication).Station Providencia was excluded due to the fact that it is located just a few meters away from a highly loaded road, and it representativity is therefore doubtful.Our analysis refers to January 2002.We visually inspected CO time series of all stations to detect any misleading data.This led us to exclude part of the data from stations Cerrillos and El Bosque, which showed at that period values out of range suggesting a possible malfunctioning of the instruments.The inspection of the data showed similar diurnal cycles for all stations except for Parque O'Higgins and Las Condes.The majority of the stations showed a sharp morning maximum around 8 a.m.(LT) and a broad evening maximum around 21 p.m. (LT) associated to traffic rush hours (Fig. 3).In the case of Parque O'Higgins the afternoon maximum is absent.Las Condes showed a mid-morning maximum at 11 a.m.(LT) and an evening maximum.This behavior is present in all years of analyzed data (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006).
Another relevant aspect is that the variance of the observations is generally higher during night-time than during day-time.This is explained by the calm conditions (winds slower than 1 or 2 m/s, Fig. 4) and the ten times smaller emission rates that prevail Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc

Printer-friendly Version
Interactive Discussion during night compared to daytime (Fig. 1).

Emissions
Emissions used in this work correspond to the Santiago official inventory for year 2002 (www.conama.cl).Total emissions of CO estimated for year 2002 are 200 ktons CO per year, accounting for mobile and stationary sources.Mobile sources add up to 90% of total CO emissions within the Metropolitan Region.The inventory is temporally disaggregated in hourly emission rates for a typical working day, at a spatial resolution of 2×2 km 2 over a domain that covers the whole urban area (Fig. 1).
Mobile emissions are estimated according to a bottom-up methodology (MODEM model, Corval án and Osses, 2002) considering official traffic modeling results (ES-TRAUS model, De Cea et al., 2003), comprehensive traffic counts, analysis of databases for vehicle technology distribution and emission factors from COPERT III model (Ntziachristos and Samaras, 2000) and local measurements.
The emission rates show a typical diurnal cycle, strongly related to traffic activity, with a primary maximum around 8 a.m. and secondary but wider maximum between 17 and 20 p.m.There is a third peak which appears at noon with variable intensity depending on the area under analysis (Fig. 1).Night values of vehicle emissions are less reliable due to scarce and less accurate traffic information (Corval án and Osses, 2002).

Meteorological model
In this work, weather fields to drive the dispersion model are generated by the Pennsylvania State University/National Center for Atmospheric Research numerical model MM5 (Grell et al., 1995).We chose this model since it is used for producing operational forecasts by the DMC.The model is run in a nested mode at 54, 18, 6 and 2 km horizontal resolution, with the highest resolution domain of roughly 200×140 km 2 centered on Santiago.Vertically, the model covers from the surface up to 100 hPa in 31 levels.We use a similar set of parameterizations as used in the operational forecast at Introduction

Conclusions References
Tables Figures

Back Close
Full DMC.In particular, we use MRF scheme for characterizing the boundary layer (Hong and Pan, 1996).NCEP reanalyzes (Kalnay et al., 1996) are used for providing boundary and initial conditions.We extract three-dimensional fields for wind, temperature, humidity, etc., from the MM5 simulations for the inner most domain (2 km horizontal resolution).These fields are interpolated into a latitude, longitude and altitude domain with twelve vertical levels, which we use to drive the dispersion model (see Sect. 2.5).Further details can be found in Saide (2008).

Dispersion model
We use the Polyphemus platform (e.g., Mallet et al., 2007, and references therein).
The core of Polyphemus consists in the chemistry-transport-deposition model Polair3d (Boutahar et al., 2004), which is an Eulerian model that solves the continuity equation using operator splitting.Polair3d was configured to solve only transport for CO, given its relatively long turn-over time of 2-4 months (Seinfeld and Pandis, 2006).No scavenging is considered as summers are dry in Santiago and the solubility of CO is negligible at this scale.Dry deposition is not included since CO has a very low deposition velocity (<1 cm/s).To compute vertical diffusion coefficients we use the Louis parameterization (Louis, 1979).Boundary conditions are Dirichlet homogeneous and zero for the inflow, and Neumann homogeneous for the outflow.The zero inflow assumption is justified since emission rates in Santiago are very high and no relevant sources are found up-wind from the city.We start the simulations from null values, leaving 24 h of spin-up time.
The Polyphemus domain is centered at Santiago, with 70 (longitude)×63 (latitude) grid-cells of 2×2 km 2 , and twelve vertical levels: 13,38,188,463,750,1050,1400,1800,2250,2750,4000 and 6000 m above the surface.We use a fine resolution in the boundary layer (BL) since CO is mainly found near the surface.A sensitivity analysis was made increasing the vertical resolution from 12 to 24 vertical levels.The outputs show no significant difference, thus we kept the 12 levels configuration. 3 Direct simulations

Simulated meteorological fields versus observations
Figure 2 shows the temporal evolution of temperature and wind direction and speed below 5 km at Santo Domingo as observed and modeled.Overall, the model captures the observed synoptic variations, particularly the coastal low event described earlier.
However, increased vertical mixing between 15th and 20th January is underestimated.
Most weather models show an unsatisfactory representation of low level clouds (e.g., EOL, 2006, and references therein).Regarding wind, temperature and humidity conditions in the Santiago basin as described by the meteorological stations within the urban area of Santiago, the model is able to reproduce most observed features (Fig. 4).
Wind speed and temperature are generally reproduced within a couple of measurement units, showing a high correlation (>0.7).Nevertheless, modeled wind speed values are higher than observations, particularly during night-time.This might be related with shortcomings in the BL parameterization discussed later or an inadequate representation of the surface roughness.Relative humidity is generally represented within 25% of the observed values.Wind direction is modeled less accurately than the other variables, particularly during night-time when winds are weaker (<2 m/s).However, prevailing winds are reproduced within roughly 50 degrees or less.
The model is able to capture the marked diurnal cycle characteristic of the summer radiatively driven circulation.However, the amplitude of the diurnal cycle is somewhat underestimated.During night-time, the model tends to uphold warmer temperatures, more humid air and stronger winds than observed.This can be explained since BL parameterizations fail to represent very stable nocturnal conditions as those that possibly prevail in the Santiago basin (e.g., Zhong and Fast, 2003).Unfortunately, there are no observational data to support or disprove this hypothesis.However, by calculating the temperature gradient between station Independencia (downtown Santiago) and station Lo Prado (1060 m.a.s.l.over the mountain range west of the Santiago basin) we find that the model systematically underestimates nocturnal gradients, suggesting Introduction

Conclusions References
Tables Figures

Back Close
Full an exaggerated vertical mixing during night-time (not shown), which is coherent with the proposed hypothesis.
In summary, we evaluate the model performance as adequate for the air quality simulations we are interested in, particularly during daytime.

Simulated CO fields versus observations
Accordingly with the prevailing winds in the Santiago basin during summer, during daytime the large emissions of CO that take place in Central Santiago are transported towards and accumulated in Eastern Santiago.During night, when the winds turn to down-slope winds, CO is transported westwards.This is illustrated in Fig. 5.The daily CO values are generally captured by the model except at Parque O'Higgins, Las Condes and to a lesser extent at La Florida, where discrepancies are evident (Fig. 3).At Parque O'Higgins the morning maximum is somewhat overestimated and no evening maximum appears in the observations, differently from the simulation that shows a distinct secondary maximum.We attribute these mismatches largely to an inadequate characterization of the emissions, rather than to shortcomings in the representation of the transport processes.In fact, wind speed is overestimated at this station, which should result in higher and not lower concentrations.The observations at Las Condes and at La Florida show a mid-morning maximum (arround 11 a.m.) instead of (Las Condes) or in addition to (La Florida) the morning maximum associated to the morning rush-hour emissions prescribed in the model.Both the simulations and the observations show a sharp transition in wind direction between 8 a.m. and 10 a.m.denoting the transition from the nocturnal (easterly, downslope) and the diurnal (westerly, upslope) regimes (Fig. 4).Again, we suspect the emissions rather than the transport.Although the characterization of the diurnal cycle in emissions has uncertainties, we find no evidence of a missing mid-morning maximum, suggesting the presence of unknown stationary sources.
The nocturnal values are generally underestimated at all stations.We hypothesize that this may depend on too strong a mixing in the model.Since this is difficult to Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion assess and, as stated before, emission rates are less reliable for night-time, we will exclude the night-time values in the inverse simulations.
A normalized Taylor diagram (Taylor, 2001) in Fig. 6 shows the error statistics for the period between 2 and 29 January 2002.The model results are divided into three categories: diurnal, nocturnal and weekend's values.In addition to the disagreement between simulations and observations for nocturnal values, a mismatch is also apparent for weekend values.These values will also be excluded from the inverse modeling exercise.
Overall, the present simulations show better or comparable error statistics than previous works (e.g., Gallardo et al., 2000;Schmitz, 2005).It is worth noticing that these authors found similar systematic errors for the stations located on the easterly bound of Santiago.

Method
The basic idea of the 3-D-Var approach is to take the dispersion model as a function (H) that connects the emission (x) space to the observation (y) space.This function is linearly expanded around an a priori or first guess (x b ), in this case the official EI for CO.Since the chemistry is switched off and only advection and turbulent mixing is applied, the source-receptor relationship is linear afine, obtaining: where y 0 is the concentration coming from emissions not optimized, initial and boundary conditions and H is the Jacobian of the model or sensitivity matrix with respect to emissions and it is evaluated close to the a priori guess emissions, i.e., those that minimize the cost function: where R is the normalized error covariance matrix for the observations, B the normalized error covariance matrix for the emissions, α is a regularization parameter, diag() represents a diagonal matrix with the diagonal elements in brackets and F is a matrix of normalizing factors for B matrix to avoid spurious temporal variability for the emission areas co-located with observations.The terms in the expressions above are explained in the following subsections.It can be shown (e.g., Kalnay, 2002) that the optimal solution is given by: where x a is the analysis (a posteriori inventory), and P a is the covariance matrix for the analysis.
The key aspects of this algorithm are, on the one hand the calculation of the sensitivity matrix H, and on the other hand, the choice of the error covariance matrices.The choice of α and F, as we explain later, are useful for balancing the model and the parameter errors., 2007).In general, if the number of parameters is much larger than the number of observations, it is preferable to use the adjoint approach.In our case, H can be estimated either way.We use an approximate adjoint for a linear trace by using the direct model running backwards on time and with the opposite of the wind fields (Davoine and Bocquet, 2007).Figure 7 shows a comparison between both estimates showing a significant correspondence (correlation higher than 0.99).For computing sensitivities, we run the adjoint with a unitary volume emission at the observation sites and recover sensitivities at each emission time and site.

Covariance matrices and balancing term
The error covariance matrix for the emissions is assumed to be diagonal representing the error variance at each grid-cell and hour in the inventory.No covariance is assumed given the lack of information on the emission model errors.Since the whole EI is obtained with the same methodology (Corval án and Osses, 2002), we assume the error variance to be equal for every emission parameter.In the same way, the error covariance matrix for the observations is taken as diagonal with the same error variance for each station given that the monitoring stations are surveyed by the same institution and according to the same protocol.The use of diagonal matrices is a common assumption because the error covariances are very difficult to assess (e.g., Qu élo et al., 2005).These matrices can be decomposed into a correlation matrix and a diagonal matrix of the variances (Kalnay, 2002).This formula can be factorized by the error variance mean value, obtaining the normalized error covariance matrices that appear in Eq. ( 2): given observation).The former is easy to asses (around 5% for our case) but the latter is not clear and generally is bigger than the former.Also, error variances in emissions are not well constrained given the uncertainities in the emission inventories.For this reason, instead of taking uncertain values on the variances, a regularization parameter is defined: where α compensates the weights between the two errors (observations and parameters) in the cost function.There are many techniques to obtain this factor such as maximization-expectation algorithm, unbiased predictive risk estimator, generalized cross validation, generalized maximum likelihood, and others (Davoine and Bocquet, 2007).Here we choose the L-curve method because it is easy to implement, it has been used with good results and it is applicable when constraints on the solution are used (Hansen and O'Leary, 1993).The method consists in comparing the growth of the observation error term with the growth of the emissions error term in a log-scale plot when α is changed.The errors terms then would be: with oe and ee the errors in the observations and in the emissions respectively.An L-shaped curve is obtained by comparison.The chosen α corresponds to the point of maximum curvature of the L-curve, where there is a balance between both errors.In our application, given the form of the resulting curve, we tested several values of α finding a minor sensitivity to the choice of the exact value.We adopted a value of 1.4E+7 (Fig. 8).

Positive emissions
When using the 3-D-Var solution (Eqs. 4 and 5) the optimized emissions obtained are not constrained to be positive.To enforce this condition the L-BFGS-B code (Zhu et al., 1997) was used to minimize the cost function in Eq. ( 2).Results show that the emissions obtained with Eq. ( 4) are similar (when positive) with the ones obtained through minimization.

Co-localization of sources and observations
When retrieving emissions by assimilation of ground measurements, an unrealistic influence by the observation sites has been found when spatial resolution is increased (Bocquet, 2005).This problem raises because the sensitivity matrix has local maxima at the observation sites.In simple words, if the resolution is increased, in the limit one can only modify emissions at the grid-boxes containing the observation sites and leave the rest of the emission field unchanged (Issartel et al., 2007), which is the behavior found in our case when assimilating data without a co-localization factor.There are many ways to deal with this problem.For instance one can decrease the resolution, i.e., reduce the number of parameters by combining emission grid-boxes.This can be done by adding the columns in the sensitivity matrix (H) or by using a non-diagonal B matrix by means of defining a radius of influence (e.g., Peters et al., 2005).Both approaches reduce the local influence of the observations, however they do propagate the observed variability to the neighboring grid-boxes, which might not always be realis-Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion tic.Another way to solve the co-localization problem is to use a renormalized geometry by computing a renormalization weighting function that lessens the local influence of observations (e.g., Issartel et al., 2007).It can be demonstrated (using Eqs. 4 and 5) that when R=0, x b =0 (full confidence in observations and no prior) and B is taken as diagonal, equal to the identity the renormalization weighting function is equal to the diagonal of F. So, from another point of view, the resolution problem can be solved applying a factor over the B matrix that can be computed as the renormalization weighting function described in Issartel et al. (2007).This function has been developed for retrieving point sources and further work is needed to apply it to dense sources (like those associated to mobile sources in a city) co-localized with observations and when a guess or prior is used.
In this work, we developed an approach in which a F factor is chosen so that the reliability of the initial inventory is increased at the observation sites, reducing the local infuence of the observations.As stated before, the norm for the columns of the H matrix corresponding to emission spots co-located with observations is very high compared to that of the rest of the columns and the closer to the observation site, the higher the sensitivity.Also, there are differences between sensitivities at different hours of the day due to the changes in meteorological conditions affecting mixing.One would want to smooth these differences by applying a factor but, as seen in Eq. ( 5), in the solution B matrix multiplies H matrix, so applying a factor to the B matrix would finally mean to apply a factor to the H matrix columns.To obtain the factor, we first compute the module norm of all columns of the H matrix normalized by the maximum value obtaining a preliminary factor: where i and j are emissions spatial and the temporal indexes, k and l the observations spatial and temporal indexes and h stands for the elements of H.The anticorrelation between F and the atmospheric mixing is apparent with higher F values when advection (wind speed) and vertical mixing (BL) are low resulting in higher sensitivity.Thus, dividing matrix B by this factor will decrease the variances in places with higher sensitivities, producing the required effect.The final factor to apply to the B matrix has to be between the value of the preliminary factor and the square of the value of the preliminary factor (Eq. 5, B appears 2 times in the equation multiplied by H on both sides and from one side).We chose to use the value without modifications because better results are obtained when applying the methodology with synthetic observations (Sect. 4.3).This factor represents a way to include strongest a priori information over the system and avoids spurious propagation of the observations to neighboring grid-boxes.

System set-up
Regarding the selection of observations, in addition to the criteria discussed in Sect.3.2, we restrict our inversion to emissions over a threshold of 0.5 (µg/m 2 /s) (temporal mean) representing emissions occurring within the urban area.In this way we reduce the number of parameters to be improved, reducing the computing time.Now, the observations reflect the influence from all emissions whereas the H matrix only reflects the information for the reduced inventory.To account for the missing emissions, a run is made with all the removed emissions (night values and emissions below the threshold) to estimate their impact at the observation points.The estimated impact is subtracted from the observations.Altogether, the number of parameters to be determined is 3615, i.e., 241 cells, with 15 h of temporal resolution from 6 a.m. to 20 p.m.The total number of observations is 1344 given by hourly outputs for 8 days (total of 10 days subtracting the 2 days of the weekend) and 7 stations.Hence, matrix H has m=1344 rows and n=3615 columns.In sum, we use the same emissions for every week-day, considering a diurnal variation at each grid-cell, and excluding week-ends and nocturnal values.

Testing with synthetic observations
In order to test the methodology and see what can be expected from the modified 3-D-Var approach, a synthetic run is performed.We take the official inventory as accurate or "true" and produce synthetic or pseudo-observations.These pseudo-observations are perturbed with a Gaussian error of 5% to simulate the error in the observations.
Then an a priori inventory (x b ) is generated by adding to all grid cells a positive bias of 5 (µg/m 2 /s), and the inversion procedure is applied.Table 1 shows a comparison between different assumptions over the emission error covariance.Three tests were performed, considering no a priori information, covariance between nearby emissions given by a radius of influence (r.o.i.) of 2.5 grid cells (e.g., Peters et al., 2005) and using F as defined before.It can be seen that improved statistics are obtained using the F factor.
Figure 10 shows results obtained with the F factor.In both a) and b) it can be seen that the difference between simulated concentrations (using the corrected inventory, x a ) and pseudo-observations is diminished.An interesting aspect to note is that the lack of an arbitrary radius of influence given by the monitoring stations, allows to improve the emissions even far from the observation sites.Also, emissions at different hours of the day are modified in a similar way.These two features are the reason why this methodology shows better results than using a covariance matrix with a radius of influence (Table 1).Another aspect to take into account is that using the F factor the emission error covariance matrix is diagonal decreasing computing time with respect the use of a radius of influence, where this matrix is non-diagonal.

Using real observations
When applying the method with actual observations, we split the analysis into two parts.First we analyze the spatial changes and then the temporal changes.Given that we only consider 10 days of data, we concentrate on the average diurnal cycle.
Figure 11 shows the difference between the a priori inventory and the a posteriori Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion inventory for the time average of morning hours (left) and afternoon hours (right).The assimilation of observations suggests that the official inventory overestimates the emissions in the center-west part of Santiago, and it underestimates the emissions over the eastern bound of the city.When looking at the diurnal cycles of the retrieved emissions at different areas of Santiago (Fig. 12), we find that overall the changes are larger in the morning hours, between 6 and 12 a.m.(LT), than in the afternoon.During the rest of the day, the a posteriori inventory shows a similar behavior to that of the a priori inventory.Thus, the assimilation of CO concentrations does not bring major changes in emissions for any of the stations in the afternoon.Changes occur mainly in the morning coinciding with the morning traffic rush-hour.These features are illustrated in Fig. 12.In summary, the a posteriori inventory shows a decrease in total emissions of 8% with respect to the a priori inventory.Nevertheless, locally over 100% changes are found which is the case for emissions located near Las Condes (north-east of the city) during morning hours (Fig. 12).

Validation of the improved inventory
To address the question of the robustness of the changes to the inventory, we performed simulations for a ten-day period starting on 3 January, when no observations were assimilated.We run the model with both the a priori and a posteriori inventories.Overall, the new inventory results in improved error statistics (see Table 2), particularly in the morning.For instance, Fig. 3 shows the comparison between the simulated diurnal cycles with the old and the new inventories and the observations.Overall, and as explained earlier (Sect.4.3), major changes occur only during the morning.Thus, the changes to the a priori inventory seem robust, and adequate for both periods simulated.Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version Interactive Discussion

Conclusions
We have presented a methodology to improve emission inventories at a city scale.This methodology takes into account the question of having positive solutions, and more importantly, it deals with the co-localization of sources and observations.This is particularly important when dealing with relatively high resolution inventories.We applied this inverse modeling procedure over Santiago de Chile for hourly observations for a ten day period in summer 2002, optimizing the emission temporal and spatial resolution running the Polair3d dispersion model and its approximate adjoint.Overall, the assimilation indicates that total amounts of emissions in the official inventory are slightly overestimated (around 8%).Nevertheless, the assimilation induces significant changes at several areas (e.g., near Las Condes station in morning hours).Also, an underestimate of the emissions in eastern Santiago is apparent.In particular, we found a mid-morning maximum at stations Las Condes and La Florida that is difficult to reconcile with observed traffic patterns, suggesting the presence of undetermined stationary sources.
The results described above are robust for the analyzed diurnal summer conditions, when circulation in the Santiago basin has a strong radiatively driven diurnal cycle, and no significant synoptic changes occur.In fact, agreement between simulated and observed CO values is also improved outside the assimilation window and when using the new corrected inventory instead of the official one.In future work, we will evaluate the assimilation approach during winter conditions that are by far more synoptically variable than summer conditions.
We tackled the problem of co-location of observations and emissions by weighting the emission error covariance matrix (B) so that the reliability of emissions is increased at the observation sites.Also, the weighting is useful to avoid a radius of influence around monitoring stations and differences in the estimates in different hours of the day.The choice of appropriate weighting factors, including the balancing effect introduced by the L-curve approach used here, remains somewhat arbitrary.It appears necessary Introduction

Conclusions References
Tables Figures

Back Close
Full All in all, we have shown a method that provides an inexpensive way to improve and update emission inventories.The key-constraint for the use of this method, and any other inverse method, is the availability of reliable and representative observations enabling the accurate evaluation of models and meaningful data assimilation.In the case of Santiago the network is reliable but its representativity still requires amelioration.In this respect, inverse modeling is again an useful tool that we will apply in our future work.Full x b .Then a cost function (J(x)) that balances the errors in the simulated values versus the observations (model error) and in the emissions (parameter error) is minimized to obtain an estimate of the optimal Introduction and C p are the correlation matrices for R and B respectively, D o and D p the diagonal matrices of error variances and σ 2 obs and σ 2 par the mean value of the error variances.The error variance in the observations must take into account instrumental error and representativeness error (errors due inability of the model to reproduce one Figure9displays the spatial and temporal distribution of this factor, showing the features described above.
the use of more physically based corrections.

Fig. 2 .Fig. 4 .
Fig. 2. Observation and model vertical profiles of temperature and wind for Santo Domingo station, values are obtained at 00:00 and 12:00 UTC daily for January 2002.

Fig. 3 .Fig. 4 .
Fig. 3. Diurnal cycle for ten days outside the inversion period for observations and results of the dispersion model using optimized and guess EI.Vertical bars represent standard deviation from observations.Units in µg/m 3 .

Fig. 5 .Fig. 6 .
Fig. 5. Concentration (µg/m 3 ) and wind (m/s) fields at the first level of the dispersion model.Panels (a) and (b) are mean fields for 8 a.m. and 16 p.m. LT respectively.Point filling colors represent observation mean values in each station at the same hour.For topography values and station number refer to Fig. 1.Concentration scale is logarithmic.

Fig. 7 .
Fig. 7. Scatter-plot comparing direct and adj Direct obs.are obtained running the model to be optimized, and adjoint obs.are obtaine trix trough the adjoint and multiplying by the µgr/m 3 .

Fig. 6 .Fig. 7 .
Fig. 6.Taylor diagram of CO simulated and observed, distinguishing weekend, night and day values for 2 to 29 January using the official inventory.

Fig. 7 .
Fig. 7. Scatter-plot comparing direct and adjoint concentrations.Direct obs.are obtained running the model with the emissions to be optimized, and adjoint obs.are obtained computing H matrix trough the adjoint and multiplying by the emissions.Units in µg/m 3 .

Fig. 8 .
Fig.8.L-curve to obtain the optimal value for the alpha parameter.Numbers at the right of each spot are the alpha value for that point.α * is the optimal value chosen.

Fig. 8 .
Fig.8.L-curve to obtain the optimal value for the alpha parameter.Numbers at the right of each spot are the alpha value for that point.α * is the optimal value chosen.

Fig. 9 .
Fig. 9. Spatial and temporal distribution of the F factor.Left panel: Time average of the F factor in a log scale.Note the maxima at the observation sites and the decay when increasing the distance from them.Right panel: Normalized values of diurnal cycle of the inverse of F, wind speed and boundary layer height (BL).

Fig. 10 .
Fig. 10.Synthetic simulation results.a) Mean in time of the difference between real and optimized EI.Note that values near 0 (green) show good agreement thus good inversion and values near -5 mean that the analyzed emissions are close to the guess emissions.b) Temporal profiles of the real, optimized and guess emissions.Units in µgr/m 2 /s.

Fig. 10 .
Fig. 10.Synthetic simulation results.(a) Mean in time of the difference between real and optimized EI.Note that values near 0 (green) show good agreement thus good inversion and values near −5 mean that the analyzed emissions are close to the guess emissions.(b) Temporal profiles of the real, optimized and guess emissions.Units in µg/m 2 /s.

Fig. 11 .Fig. 11 .
Fig. 11.Time average difference between the a priori (official) inventory and the optimized inventory assimilating real observations.Left panel: morning, right panel: afternoon.Note that values over 0 (red) represent a decrease in the emissions and negative values (blue) represent increase in emissions with respect to the official inventory.Units in µg/m 2 /s.

Fig. 12 .
Fig. 12. Daily cycle for reference and optimized emissions for a central zone of the city (left panels) and for an eastern zone of the city (right panels).Units in µgr/m 2 /s.Station names represent the grid where they are located.

Fig. 12 .
Fig. 12. Daily cycle for reference and optimized emissions for a central zone of the city (left panels) and for an eastern zone of the city (right panels).Units in µg/m 2 /s.Station names represent the grid where they are located.
To obtain H there are at least two methods: either use the direct model or use the adjoint of the model (e.g., Davoine and Bocquet

Table 1 .
. and Pandis, S.: Atmospheric Chemistry and Physics, Wiley,USA, 2006.6331Taylor,K.: Summarizing multiple aspects of model performance in a single diagram, J.Geophys.Res., 106, 7183-7192, 2001.6334 Zhong, S. and Fast, J.: An evaluation of the MM5, RAMS, and Meso-Eta models at subkilometer resolution using VTMX field campaign data in the Salt Lake valley, Mon.Weather Rev., 131, Introduction Statistics for synthetic case (see Sect. 4.3 for details).Three inversion methods are compared for different assumptions regarding the emission error covariance matrix.

Table 2 .
Error statistics simulations using the optimized and the guess inventories for a period out of the assimilation window.Correlation coefficient (R), root mean square error (RMSE, units in µg/m 3 ) and bias (units in µg/m 3 ) were calculated for each station.

Table 1 .
Statistics for synthetic case (See Section 4.3 for details).Three inversion methods are compared for different assumptions regarding the emission error covariance matrix.R RMSE [µgr/m 2 s] Bias [µgr/m 2 s] F factor 0.99 0.9 0.4