# Inverse modeling of nitrogen oxides emissions in urbanized regions of the Northern Hemisphere: Application of the CIF-CHIMERE inversion system to OMI and TROPOMI observations

Dilek Savas

► To cite this version:

Dilek Savas. Inverse modeling of nitrogen oxides emissions in urbanized regions of the Northern Hemisphere: Application of the CIF-CHIMERE inversion system to OMI and TROPOMI observations. Environmental Sciences. Université Paris Cité, 2023. English. NNT: tel-04349344

HAL Id: tel-04349344

https://hal.science/tel-04349344

Submitted on 17 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Université Paris Cité

Ecole doctorale des Sciences de l'Environnement d'Ile de France (ED129)

Laboratoire Interuniversitaire des Systèmes Atmosphériques (LISA)
UMR7583)

# Inverse modeling of nitrogen oxides emissions in urbanized regions of the Northern Hemisphere : Application of the CIF-CHIMERE inversion system to OMI and TROPOMI observations

Par Dilek Savas

Thèse de doctorat en Météorologie et Qualité de l'Air

Dirigée par Gaëlle Dufour

Présentée et soutenue publiquement le 8 Juin 2023

Devant un jury composé de :

Matthias Beekmann, Directeur de recherche, CNRS & Université Paris Cité & Université Paris Est Créteil, Président du jury

Claire Granier, Directrice de recherche, CNRS & Université de Toulouse III, Rapporteure

Vincent Guidard, Directeur de recherche, Météo-France, Rapporteur

Trissevgeni Stavrakou, Directrice de recherche, BIRA-IASB, Examinatrice

Solène Turquety, Professeur, Sorbonne Université, Examinatrice

Gaëlle Dufour, Directrice de recherche, CNRS & Université Paris Cité & Université Paris Est Créteil, Directrice de thèse
<figure><img src="image_2.png" /></figure>
# Titre : Modélisation inverse des émissions d'oxydes d'azote dans les zones urbanisées de l'hémisphère Nord : Application du système d'inversion CIF-CHIMERE aux observations OMI et TROPOMI

Résumé : Selon l'Organisation mondiale de la santé, l'exposition à la pollution de l'air ambiant serait à l'origine de 4.2 millions de décès prématurés par an dans le monde. Ce problème environnemental majeur est apparu au début du XXe siècle en raison de l'industrialisation et de l'urbanisation intenses des sociétés, obligeant les décideurs et les gouvernements à contrôler les émissions anthropiques de polluants atmosphériques. Depuis 1990, des initiatives importantes ont été prises dans les pays industrialisés, comme l'Europe, pour réduire les émissions de polluants en mettant en œuvre des stratégies d'atténuation. Plus récemment, la mise en œuvre de politiques d'atténuation a commencé à être observée dans des pays en développement comme la Chine.

L'évaluation des politiques de réduction des émissions repose souvent sur des inventaires officiels, établis à partir des déclarations d'activités des émetteurs et de données statistiques. Cette approche dite "bottom up" estime les émissions en extrapolant les mesures des facteurs d'émission qui ne sont disponibles qu'à des résolutions spatiotemporelles réduites. Cette méthode conduit des incertitudes élevées, car elle ne tient pas compte de la grande variabilité spatio-temporelle des flux d'émission. En outre, les incertitudes liées aux inventaires des émissions entraînent des incohérences dans l'évaluation des stratégies les plus appropriées pour prévenir les épisodes de pollution atmosphérique. Par conséquent, il est important de pouvoir suivre avec précision l'évolution des émissions et la tendance des polluants pour évaluer les politiques de réduction. L'un des moyens prometteurs pour cela est l'approche par inversion atmosphérique, qui utilise les observations atmosphériques disponibles pour contraindre les émissions par le biais de la modélisation atmosphérique et de méthodes d'inversion. Cette approche dite "top down" apporte des informations complémentaires aux inventaires "bottom up". Elle estime les émissions tout en exploitant la grande variabilité spatio-temporelle des observations satellitaires et la puissance de calcul.

Les oxydes d'azote ($NO_x = NO+NO_2$) sont parmi les polluants les plus réglementés en tant que précurseurs d'autres polluants atmosphériques, tels que l'ozone et les aérosols secondaires. Dans le cadre de cette thèse, nous avons tout d'abord mis en place et testé la capacité d'inversion des émissions quotidiennes de $NO_x$ du système de modélisation inverse CIF, couplé au modèle de chimie-transport CHIMERE et son adjoint, à une résolution modérée, en utilisant les observations du satellite OMI. Les résultats nous ont permis de déterminer les paramètres et les sensibilités du système CIF pour les inversions de $NO_x$.
Le système CIF-CHIMERE a d'abord été appliqué pour évaluer l'impact des fortes réglementations sur les émissions de NOₓ mises en place par la Chine depuis 2011. Nous avons assimilé les observations de NO₂ des instruments du satellite OMI et estimé les émissions de NOₓ pour 2015 et 2019 avec une résolution de 50x50 km² sur la Chine. L'année 2010 a été choisie comme base de référence ou inventaire a priori, car elle se situe juste avant la mise en œuvre des réductions d'émissions. Les résultats montrent que la réduction des émissions de NOₓ est limitée aux zones urbanisées et industrialisées tout en restant dans les limites des objectifs d'atténuation (10-15 %). Les émissions de NOₓ estimées ont également été utilisées pour simuler les concentrations de NO₂ en surface, et permettent d'améliorer légèrement la comparaison avec les mesures au sol.

Le système CIF-CHIMERE a également été appliqué pour évaluer les changements abrupts des émissions de NOₓ en Europe causés par la pandémie de COVID-19. En particulier, nous avons exploré le potentiel d'assimilation des observations à haute résolution de TROPOMI NO₂ pendant la période de confinement de 2020.

Mots clefs : Pollution Atmosphérique, Qualité de l'Air, Émission de NOx, Assimilation de Données, Modélisation Inverse, Oxydes d'Azote
# Title: Inverse modeling of nitrogen oxides emissions in urbanized areas of the Northern Hemisphere: Application of the CIF-CHIMERE inversion system to OMI and TROPOMI observations

**Abstract:** According to the World Health Organization, exposure to ambient air pollution is estimated to cause 4.2 premature deaths per year in urban and rural areas worldwide. Poor air quality also leads to damage to infrastructure, soil and water resources. This major environmental problem emerged at the beginning of the 20th century due to the intense industrialization and urbanization of societies, forcing policy makers and governments to control anthropogenic emissions of air pollutants. Since 1990, important initiatives have been taken in industrialized countries, such as Europe, to reduce emissions of pollutants by implementing mitigation strategies. More recently, the implementation of mitigation policies has started to be seen in developing countries such as China.

The assessment of emission reduction policies is often based on official emission inventories derived from emitters' statements of activities and statistical data. This so-called bottom-up approach estimates emissions by extrapolating emission factor measurements that are only available in a sparse spatial and temporal network. This involves high uncertainties, as it does not incorporate the high spatiotemporal variability of emission fluxes. Besides, uncertainties in emission inventories lead to inconsistencies in the assessment of appropriate strategies to prevent air pollution episodes. Therefore, being able to accurately monitor the development of emissions and the trend of pollutants is an important issue for evaluating reduction policies. One of the promising ways to overcome this problem is the atmospheric inversion approach, which uses available atmospheric observations to constrain emissions through atmospheric modeling and the inversion system. This so-called top down approach brings complementary information to bottom-up inventories. It estimates emissions while exploiting the high spatiotemporal variability of the satellite observations and computational power.

Nitrogen oxides ($NO_x = NO+NO_2$) are among the most regulated pollutants as precursors of other air pollutants, such as ozone and secondary aerosols. In the framework of this thesis, first, we set up and tested the daily $NO_x$ emission inversion capability of the state-of-art inverse modeling system CIF, embedded with the CHIMERE Chemistry Transport Model and its adjoint at moderate resolution using OMI satellite observations. The results lead us to determine the settings and sensitivities of the CIF system for the $NO_x$ inversions.

The CIF-CHIMERE system was applied first to evaluate the impact of strong $NO_x$ emission regulations implemented by China since 2011 within its 5-Year Plans. We assimilated $NO_2$ observations from OMI satellite instruments and
estimated NOₓ emissions for 2015 and 2019 with a resolution of 50x50 km² over Eastern China. The year 2010 was chosen as a baseline or a priori bottom-up inventory, as it was just before the mitigation implementation. The results show that the reduction in NOₓ emissions is limited to urbanized and industrialized areas but remains within the mitigation targets (10-15 %). The estimated NOₓ emissions were also used to simulate NO₂ surface concentrations, and the inversion approach was found to improve the comparison with ground-based measurements slightly.

The CIF-CHIMERE system was also applied to assess the abrupt changes in NOₓ emissions in Europe caused by the COVID-19 pandemic. In particular, we explored the potential of assimilating high-resolution TROPOMI NO₂ observations during the 2020 lock-down period.

Keywords: Atmospheric Pollution, Air Quality, NOx Emission, Data Assimilation, Inverse Modeling, Nitrogen Oxides
**Acknowledgements**

I would like to express my deep gratitude to everyone who has made meaningful contributions to my academic and personal development during this time. Above all, I sincerely thank my supervisor Gaëlle Dufour for her valuable guidance, profound expertise and feedback throughout the whole process. Her support and encouragement were of critical importance to my success, and I am deeply grateful for her mentorship.

I am also grateful to the committee members; Vincent Guidard, Claire Granier, Trissevgeni Stavrakou, Matthias Beekmann, and Solène Turquety for their insightful feedback and constructive criticism, which have contributed additional perspectives to my research.

I would like to thank the agencies that supported this research; the French National Research Agency (ANR) and Carole Daniel from the French National Centre for Space Studies (CNES) for the doctoral grant. Their financial assistance has made this research possible, and I am deeply grateful for their support.

I also would like to thank Adriana Coman for her continuous support and generous knowledge sharing. It has always been a great pleasure for me to talk to her and learn from her expertise. Without her support, much of this work would not have taken the shape it has. I am honored to have had the privilege of working with her.

I would like to acknowledge the support and encouragement of my colleagues at the LISA, and our collaborators at the LSCE within the ANR-ARGONAUT project for their help and scientific insights.

Additionally, I am grateful to my family and friends who have been a constant source of love, support, and encouragement throughout my academic journey. Their unwavering support, understanding, and patience have been critical in helping me stay focused and motivated especially during difficult times such as COVID-19 pandemic.

I would like to give special thanks to my mom and dad, Vildan-Metin Savas, for their unconditional love and tremendous support over the years to help me pursue my dreams.

Last but not least, I would like to express my sincere gratitude to Atahan Atadil, my dear life partner who has always made me feel his support, even from thousands of kilometers away. I cannot wait to explore the universe on our starship.

Thank you all! This journey would not have been possible without your support and encouragement, and I am deeply grateful for your contributions.
"On ne découvre pas de terre nouvelle sans consentir à perdre de vue, d'abord et longtemps, tout rivage." – André Gide

Dedicated to those brave enough to challenge themselves beyond their comfort zones.
Abbreviations

<table><tr><td>CH<sub>3</sub>CHO</td><td>acetaldehyde</td></tr><tr><td>AMF</td><td>air mass factor</td></tr><tr><td>AN</td><td>alkyl nitrates</td></tr><tr><td>ARGONAUT</td><td>pollutAnts and gReenhouse Gases emissiOns monitoring from spAce at high resolUTion</td></tr><tr><td>AKs</td><td>averaging kernels</td></tr><tr><td>BLUE</td><td>best linear unbiased estimator</td></tr><tr><td>CO</td><td>carbon monoxide</td></tr><tr><td>CO<sub>2</sub></td><td>carbon dioxide</td></tr><tr><td>CNES</td><td>Centre national d'études spatiales</td></tr><tr><td>CCD</td><td>charge coupled device</td></tr><tr><td>CTMs</td><td>chemistry transport models</td></tr><tr><td>CNEMC</td><td>China National Environmental Monitoring Centre</td></tr><tr><td>CAAQS</td><td>Chinese ambient air quality standards</td></tr><tr><td>UTC</td><td>Coordinated Universal Time</td></tr><tr><td>CAMS</td><td>Copernicus Atmosphere Monitoring Service</td></tr><tr><td>CFCs</td><td>chlorofluorocarbons</td></tr><tr><td>CLAES</td><td>Cryogenic Limb Array Etalon Spectrometer</td></tr><tr><td>°E</td><td>degrees east</td></tr><tr><td>°N</td><td>degrees north</td></tr><tr><td>DOAS</td><td>Differential Optical Absorption Spectroscopy</td></tr><tr><td>EOS</td><td>Earth Observing System</td></tr><tr><td>EEA</td><td>European Environmental Agency</td></tr><tr><td>ENIAC</td><td>Electronic Numerical Integrator and Computer</td></tr><tr><td>ECCAD</td><td>Emissions of atmospheric compounds and compilation of ancillary data</td></tr><tr><td>EDGAR</td><td>Emissions Database for Global Atmospheric Research</td></tr><tr><td>enKF</td><td>ensemble Kalman filter</td></tr><tr><td>4D-EnVar</td><td>Ensemble 4D Variational</td></tr><tr><td>EEA</td><td>European Environment Agency</td></tr><tr><td>ESA</td><td>European Space Agency</td></tr></table>
<table>
   <tr>
    <td>e.g.</td>
    <td>exempli gratia</td>
   </tr>
   <tr>
    <td>EPA</td>
    <td>Environmental Protection Agency</td>
   </tr>
   <tr>
    <td>Eq</td>
    <td>equation</td>
   </tr>
   <tr>
    <td>EU</td>
    <td>European Union</td>
   </tr>
   <tr>
    <td>ECMWF</td>
    <td>European Centre for Medium-Range Weather Forecasts</td>
   </tr>
   <tr>
    <td>EKF</td>
    <td>extended Kalman filter</td>
   </tr>
   <tr>
    <td>FMI</td>
    <td>Finnish Meteorological Institute</td>
   </tr>
   <tr>
    <td>4D-Var</td>
    <td>four dimensional variational</td>
   </tr>
   <tr>
    <td>GLCF</td>
    <td>Global Land Cover Facility</td>
   </tr>
   <tr>
    <td>GOME</td>
    <td>Global Ozone Monitoring Experiment</td>
   </tr>
   <tr>
    <td>GHG</td>
    <td>greenhouse gas</td>
   </tr>
   <tr>
    <td>hPa</td>
    <td>hectopascal</td>
   </tr>
   <tr>
    <td>HTAP</td>
    <td>Hemispheric Transport of Air Pollution</td>
   </tr>
   <tr>
    <td>HO<sub>2</sub></td>
    <td>hydroperoxyl radical</td>
   </tr>
   <tr>
    <td>OH</td>
    <td>hydroxyl radical</td>
   </tr>
   <tr>
    <td>INERIS</td>
    <td>Institut National de l'Environnement Industriel et des Risques</td>
   </tr>
   <tr>
    <td>INCA</td>
    <td>Interaction with Chemistry and Aerosols</td>
   </tr>
   <tr>
    <td>AR4</td>
    <td>IPCC Fourth Assessment Report</td>
   </tr>
   <tr>
    <td>TAR</td>
    <td>IPCC Third Assessment Report</td>
   </tr>
   <tr>
    <td>C<sub>5</sub>H<sub>8</sub></td>
    <td>isoprene</td>
   </tr>
   <tr>
    <td>KF</td>
    <td>Kalman filter</td>
   </tr>
   <tr>
    <td>LMD</td>
    <td>Laboratoire de Météorologie Dynamique</td>
   </tr>
   <tr>
    <td>LMDz</td>
    <td>Laboratoire de Météorologie Dynamique Zoom</td>
   </tr>
   <tr>
    <td>LISA</td>
    <td>Laboratoire Interuniversitaire des Systèmes Atmosphériques</td>
   </tr>
   <tr>
    <td>lat</td>
    <td>latitude</td>
   </tr>
   <tr>
    <td>LT</td>
    <td>local time</td>
   </tr>
   <tr>
    <td>lon</td>
    <td>longitude</td>
   </tr>
   <tr>
    <td>CH<sub>4</sub></td>
    <td>methane</td>
   </tr>
   <tr>
    <td>MEE</td>
    <td>Ministry of Ecology and Environment</td>
   </tr>
   <tr>
    <td>MEGAN</td>
    <td>Model of Emissions of Gases and Aerosols from Nature</td>
   </tr>
   <tr>
    <td>MEIC</td>
    <td>Multiresolution Emission Inventory for China</td>
   </tr>
   <tr>
    <td>NASA</td>
    <td>National Aeronautics and Space Administration</td>
   </tr>
   <tr>
    <td>NAAQS</td>
    <td>national ambient air quality standards</td>
   </tr>
  </table>
<table>
   <tr>
    <td>NIR</td>
    <td>near infrared</td>
   </tr>
   <tr>
    <td>NIVR</td>
    <td>Netherlands' Institute for Air and Space Development</td>
   </tr>
   <tr>
    <td>NO<sub>3</sub></td>
    <td>nitrate radical</td>
   </tr>
   <tr>
    <td>HNO<sub>3</sub></td>
    <td>nitric acid</td>
   </tr>
   <tr>
    <td>NO</td>
    <td>nitric oxide</td>
   </tr>
   <tr>
    <td>N<sub>2</sub></td>
    <td>nitrogen</td>
   </tr>
   <tr>
    <td>NO<sub>2</sub></td>
    <td>nitrogen dioxide</td>
   </tr>
   <tr>
    <td>NO<sub>x</sub></td>
    <td>nitrogen oxides</td>
   </tr>
   <tr>
    <td>N<sub>2</sub>O</td>
    <td>nitrous oxide</td>
   </tr>
   <tr>
    <td>NMVOCs</td>
    <td>non methane volatile organic compounds</td>
   </tr>
   <tr>
    <td>NMB</td>
    <td>normalized mean bias</td>
   </tr>
   <tr>
    <td>NGST</td>
    <td>Northrop Grumman Space Technology</td>
   </tr>
   <tr>
    <td>NWP</td>
    <td>numerical weather prediction</td>
   </tr>
   <tr>
    <td>OFFL</td>
    <td>offline</td>
   </tr>
   <tr>
    <td>OI</td>
    <td>optimal interpolation</td>
   </tr>
   <tr>
    <td>RO<sub>2</sub></td>
    <td>organic peroxy radical</td>
   </tr>
   <tr>
    <td>OECD</td>
    <td>Organization for Economic Co-operation and Development</td>
   </tr>
   <tr>
    <td>O<sub>2</sub></td>
    <td>oxygen</td>
   </tr>
   <tr>
    <td>O<sub>3</sub></td>
    <td>ozone</td>
   </tr>
   <tr>
    <td>OMI</td>
    <td>Ozone Monitoring Instrument</td>
   </tr>
   <tr>
    <td>PDEs</td>
    <td>partial differential equations</td>
   </tr>
   <tr>
    <td>CH<sub>3</sub>C(O)OO</td>
    <td>peroxyacetyl radical</td>
   </tr>
   <tr>
    <td>PANs</td>
    <td>peroxyacetyl nitrates</td>
   </tr>
   <tr>
    <td>parts per billion</td>
    <td>ppb</td>
   </tr>
   <tr>
    <td>PDFs</td>
    <td>probability density functions</td>
   </tr>
   <tr>
    <td>QA4ECV</td>
    <td>Quality Assurance for Essential Climate Variables</td>
   </tr>
   <tr>
    <td>REAS</td>
    <td>Regional Emission Inventory in Asia</td>
   </tr>
   <tr>
    <td>RMSE</td>
    <td>root mean square error</td>
   </tr>
   <tr>
    <td>SCIAMACHY</td>
    <td>Scanning Imaging Absorption Spectrometer for Atmospheric Cartography</td>
   </tr>
   <tr>
    <td>SOAs</td>
    <td>secondary organic aerosols</td>
   </tr>
   <tr>
    <td>S5P</td>
    <td>Sentinel-5 Precursor</td>
   </tr>
   <tr>
    <td>S5P-PAL</td>
    <td>Sentinel-5P Product Algorithm Laboratory</td>
   </tr>
  </table>
<table><tr><td>SCD</td><td>slant column density</td></tr><tr><td>SNR</td><td>signal to noise ratio</td></tr><tr><td>SD</td><td>standard deviation</td></tr><tr><td>SNAP</td><td>Standardized Nomenclature for Air Pollutants</td></tr><tr><td>SO<sub>2</sub></td><td>sulfur dioxide</td></tr><tr><td>SAOZ</td><td>Système d'Analyse par Observation Zénitale</td></tr><tr><td>TLH</td><td>tangent linear hypothesis</td></tr><tr><td>Tg N yr<sup>-1</sup></td><td>teragram nitrogen per year</td></tr><tr><td>3D</td><td>three dimension</td></tr><tr><td>TOMS</td><td>Total Ozone Mapping Spectrometer</td></tr><tr><td>TGCC</td><td>Très Grand Centre de Calcul</td></tr><tr><td>TEMIS</td><td>Tropospheric Emission Monitoring Internet Service</td></tr><tr><td>TROPOMI</td><td>TROPOspheric Monitoring Instrument</td></tr><tr><td>USA</td><td>United States of America</td></tr><tr><td>UARS</td><td>Upper Atmosphere Research Satellite</td></tr><tr><td>ultraviolet-visible</td><td>UV-Vis</td></tr><tr><td>Var</td><td>variational</td></tr><tr><td>v</td><td>version</td></tr><tr><td>VCD</td><td>vertical column density</td></tr><tr><td>VOCs</td><td>volatile organic carbons</td></tr><tr><td>WHO</td><td>World Health Organization</td></tr><tr><td>WMO</td><td>World Meteorological Organization</td></tr></table>
# Contents

<table><tr><td>ACKNOWLEDGEMENTS</td><td>7</td></tr><tr><td>ABBREVIATIONS</td><td>11</td></tr><tr><td>INTRODUCTION</td><td>17</td></tr><tr><td>FIRST PART : CONTEXT</td><td>21</td></tr><tr><td>CHAPTER 1 : NITROGEN OXIDES IN THE ATMOSPHERE</td><td>21</td></tr><tr><td>SECOND PART : METHODOLOGY AND INSTRUMENTS</td><td>43</td></tr><tr><td>CHAPTER 2 : INVERSE MODELING AND ITS TOOLS</td><td>43</td></tr><tr><td>THIRD PART : SCIENTIFIC APPLICATIONS</td><td>96</td></tr><tr><td>CHAPTER 3. INVERSE MODELING OF NO<sub>x</sub> EMISSIONS OVER CHINA</td><td>96</td></tr><tr><td>CHAPTER 4. NO<sub>x</sub> EMISSIONS IN EUROPE DERIVED FROM SATELLITE OBSERVATIONS</td><td>123</td></tr><tr><td>GENERAL CONCLUSIONS AND PERSPECTIVES</td><td>142</td></tr><tr><td>ANNEXES</td><td>147</td></tr><tr><td>ANNEX A : STATISTICAL FORMULAS</td><td>147</td></tr><tr><td>ANNEX B : 4D-VAR COST AND THE GRADIENT EVOLUTION</td><td>148</td></tr><tr><td>ANNEX C : NON-LINEAR OPERATORS</td><td>150</td></tr><tr><td>ANNEX D : STATISTICAL INDICATORS</td><td>151</td></tr><tr><td>ANNEX E : RATIO FORMULAS</td><td>152</td></tr><tr><td>ANNEX F : ANNUAL CONCENTRATIONS 2015</td><td>153</td></tr><tr><td>ANNEX G : STATISTICAL INDICATORS' RATIOS 2015</td><td>154</td></tr><tr><td>ANNEX H : IN SITU COMPARISON</td><td>155</td></tr><tr><td>ANNEX I : INNOVATION & RESIDUE MAP</td><td>156</td></tr><tr><td>ANNEX J : RÉSUMÉ SUBSTANTIEL EN FRANÇAIS</td><td>157</td></tr><tr><td>BIBLIOGRAPHY</td><td>163</td></tr><tr><td>TABLE OF CONTENTS</td><td>185</td></tr><tr><td>TABLE OF FIGURES</td><td>189</td></tr><tr><td>LIST OF TABLES</td><td>193</td></tr></table>
Introduction

According to the World Health Organization, exposure to ambient air pollution is estimated to cause 4.2 million premature deaths per year in urban and rural areas worldwide. Poor air quality also leads to damage to infrastructures, soil, and water resources. This major environmental problem emerged at the beginning of the 20th century due to the intense industrialization and urbanization of societies, forcing policymakers and governments to control anthropogenic emissions of air pollutants. Since 1990, important initiatives have been taken in industrialized countries, such as Europe, to reduce emissions of pollutants by implementing mitigation strategies. More recently, the implementation of mitigation policies has started to be seen in developing countries such as China.

Among the various air pollutants, nitrogen oxides (NOx = NO + NO2) are considered one of the major pollutants. Besides their direct impact on health, such as respiratory infections, asthma, and chronic lung diseases, and on ecosystems and cultures as damaging vegetation, NOx are precursors of other air pollutants, such as ozone and secondary aerosols, which also have adverse effects on human and animal health, vegetation, and climate by being short-lived climate forcers. Their dissolution in water leads to acids, and their transformation into nitrate particulate also damages vegetation, buildings, and materials and contributes to acidification and eutrophication of aquatic and terrestrial ecosystems. Then, NOx, and especially NO2, are among the most regulated pollutants. In Europe, road transport and energy (production, distribution and industrial use) are the main emission sectors followed by tertiary and residential sectors, agriculture and non-road transport (in 2019). According to the European Environmental Agency, European NOx emissions have been reduced by almost 60 % between 1990 and 2019. In China, following a significant increase in NOx emissions between the 90s and late 2000s due to rapid economic development, the 12th and 13th 5-Year government plans in 2011 and 2016, respectively, decided to implement NOx emission reductions and set targets of 10 % reduction in 2015 and 15 % reduction in 2020. In China, electricity generation and industrialization, which depend heavily on coal consumption, are the main contributor and responsible for 63 % of total NOx emissions, while road vehicles are the second largest contributor to NOx emissions by 31 %.
In parallel with the mitigation strategies, observational capabilities of pollutants and especially NO₂ concentrations have largely increased during the last three decades with the development of air quality measurement networks in Europe and more recently in China sampling various types of situations (background, rural, mountains, suburban, urban, traffic). In the meantime, satellite observations of NO₂ tropospheric columns have shown unprecedented improvements from the GOME satellite instrument with a spatial resolution of a few hundreds of kilometers in the mid-90s to OMI and GOME-2 (resolution of a few tens of kilometers) in the mid-2000s and TROPOMI (resolution of a few kilometers) since 2018. The satellite observations provide a global description of NO₂ distribution, now at high resolution, all over the globe once a day. These developments brought new insights for pollution and emission monitoring from the regional to the global scales. In addition, atmospheric modeling capabilities have also largely developed following the increase in computational power capabilities. Nowadays, chemical transport models are used routinely and operationally to simulate the tridimensional distributions of pollutant concentrations from local to global with a high spatial and temporal resolution. Operational systems, such as the Copernicus Atmospheric Monitoring Service (CAMS) in Europe or the Prev'Air national platform in France, manage to provide estimates of air pollutants at various scales by exploiting the complementarity of modeling and observations and the development of powerful data assimilation systems.

Once emission mitigations are implemented, one key issue is to evaluate their effectiveness and their efficiency. The assessment of emission reduction policies is often based on official emission inventories derived from emitters' statements of activities and statistical data. This so-called bottom-up approach estimates emissions by extrapolating emission factor measurements that are only available in a sparse spatial and temporal network. This involves high uncertainties, as it does not incorporate the high spatiotemporal variability of emission fluxes. Besides, uncertainties in emission inventories lead to inconsistencies in the assessment of appropriate strategies to prevent air pollution episodes. Therefore, being able to monitor accurately the development of emissions and the trend of pollutants is an important issue for evaluating reduction policies. One of the promising ways to overcome this problem is the atmospheric inversion approach, which uses available atmospheric observations, such as satellite observations, to constrain emissions through atmospheric modeling and an inversion system (data assimilation). This so-called top-down approach brings complementary information to bottom-up inventories. It estimates emissions while exploiting the high spatiotemporal variability of the satellite observations and computational power with a fast
update capability. Inverse modeling approaches have acquired a good level of maturity for the last decade. At the French level, the Community Inversion Framework (CIF), firstly dedicated to the inversion and monitoring of greenhouse gases at the global scale, has been in development during the last few years. Developing and testing its capabilities to perform atmospheric inversions of highly reactive species at the regional scale and a high resolution was one of the objectives of the ARGONAUT ANR project and the framework of this PhD work.

In this context, the first chapter of the thesis describes the current state of knowledge about NO₂ air pollutants, their sources, emission inventories, and the purpose of this thesis. The second chapter explains all the necessary tools and instruments to be used in inversions. It is where we set up and tested the daily NOₓ emission inversion capability of the CIF inverse modeling system, embedded with the CHIMERE Regional Chemistry Transport Model and its adjoint at moderate and high-resolution using OMI and TROPOMI satellite observations. In the third chapter, the CIF-CHIMERE system was applied first to evaluate the impact of the NOₓ emission regulations implemented by China since 2011 within its 5-Year Plans. We assimilated NO₂ observations from OMI satellite instruments and estimated NOₓ emissions for 2015 and 2019 with a resolution of 50 × 50 km² over Eastern China. The year 2010 was chosen as a baseline or a priori bottom-up inventory, as it was just before the mitigation implementation. In the fourth chapter, the CIF-CHIMERE system was applied to assess the abrupt changes in NOₓ emissions in Europe caused by the COVID-19 pandemic. In particular, we explored the potential of assimilating high-resolution TROPOMI NO₂ observations during the 2020 lockdown period. The conclusion summarizes the main results of the thesis and discusses the possible perspectives of this work.
# First part : Context

## Chapter 1 : Nitrogen oxides in the atmosphere

<figure><img src="image_3.png" /><figcaption>Figure 1 Parisian Sunset at the top of the Montparnasse tower (December 2021).</figcaption></figure>

### 1.1 Earth's atmospheric composition and structure

The atmosphere is the air mass surrounding a planet or other celestial body. The name originates from two ancient Greek words, “atmós” and “sphaîra,” meaning vapor and sphere. Earth’s atmosphere (air) is one of the four Earth subsystems, as well as the lithosphere (land), the hydrosphere (water), and the biosphere (living beings). It is invisible, odorless, and essential to the survival of every life form on the Earth. The components of the Earth’s atmosphere (hereinafter only atmosphere) have evolved since the formation of the Earth, and today, nitrogen (78.1 %), oxygen (20.9 %), and a small amount of argon (0.9 %) account for over 99.9 % of the air we breathe. The remaining 0.1 % consists of trace constituents relevant to climate
and human health, e.g., water vapor, carbon dioxide, methane, ozone, and nitrogen oxides. Besides gaseous species, solid or liquid particles of different sizes, known as aerosols, are also found in the atmosphere e.g., dust, sea salt, fog, and fume.

The atmosphere is composed of several interconnected layers, which have different characteristics mainly due to vertical variations in temperature and pressure (Figure 2).

<figure><img src="image_4.png" /><figcaption>Figure 2 Atmospheric layers (Image taken from Seinfeld and Pandis, 2016)</figcaption></figure>

The troposphere is the lowest layer from the Earth's surface to the tropopause at an altitude of 10-15 km. The upper-boundary of the troposphere (so-called tropopause) varies with the season (lower in winter, higher in summer) and the latitude (lowest over the poles, highest at the equator). It contains 80 % of the total atmospheric mass and almost all atmospheric water vapor, including clouds. It is, therefore, the layer in which most of the Earth's meteorology occurs. It is also, where air pollutants are emitted, vegetation grows, and interactions with the oceans occur. In the troposphere, the warmest part is the lowest area due to the incoming solar radiation absorbed by the Earth's surface and transferred into the heat. However, the temperature decreases by 6.5°C per km with an increasing altitude until the tropopause. In addition, this vertical temperature gradient allows for vertical mixing and meteorological disturbances.

Above the tropopause, the stratosphere stretches from about 11 km to the stratopause at 50 km. The low temperature of the tropopause acts like a barrier to water being transported into the
stratosphere, which makes the stratosphere quite dry compared to the troposphere. The stratosphere contains approximately 10-20 % of the total atmospheric mass. Within the stratosphere, there is the ozone ($O_3$) layer produced and destroyed naturally (also as a result of human activities) by incoming ultraviolet solar radiation interacting with molecular and atomic oxygen while releasing heat. This phenomenon determines the thermal structure of the stratosphere as the temperature progressively increases from 20 to 50 km until the stratopause. One of the biggest discoveries of the last century and still a hot topic is the degradation of the ozone layer in the Polar Regions due to the emissions of chlorofluorocarbons (CFCs) because of human activities (Farman, Gardiner, and Shanklin 1985; Solomon et al. 1986). Indeed, CFCs, which are entirely human-made, have caused a 'hole in the ozone layer' by breaking the natural cycle of ozone production in the stratosphere. The implementation of agreements such as the Montreal Protocol, signed in 1987, which limits the use of these ozone depletion gases, has resulted with the recovery of the ozone hole, although full recovery is not expected until the end of this century according to the projections (Osani and Haysmith 2023; World Meteorological Organization (WMO) 2018).

Above the stratopause, the mesosphere lies until the mesopause around 80 km, which is the coldest region in the atmosphere with a temperature of less than -80°C.

Above the mesopause, the thermosphere lies and is bounded by the thermopause between the altitudes 500 to 1000 km. The temperature increases at the lower levels of the thermosphere because of the absorption of the short wave solar radiation by nitrogen and oxygen at these levels. International Space Station (at 400 km) and low orbit artificial satellites (e.g., NASA Aura at 705 km, Sentinel-5P at 824 km) reside at this layer.

Atop of the thermopause, the exosphere lies, extending into the near vacuum of outer space.

## 1.2 Atmospheric lifetime

The concentration of the chemical gases and particles in the atmosphere is either directly determined by emissions (through anthropogenic or natural processes), depositions, and atmospheric transport, or indirectly by atmospheric chemical reactions (via production and loss). Atmospheric lifetime, that is the time for a molecule to be converted into another substance by atmospheric chemical reactions (or the time the molecule remains in a steady state), is defined as an indicator of the characteristic of a molecule. Depending on their
atmospheric lifetime, molecules show spatial and temporal variabilities. Some of them are chemically inert (or stable) and identified as long-lived species such as methane (CH₄), nitrous oxide (N₂O) and CFCs, whose atmospheric lifetimes are around 10 years for CH₄ and 100 years for N₂O and CFCs (Figure 3). This means that they have enough time to be transported away from their sources over long distances (on a global scale) and at higher altitudes before they chemically interact with other chemical substances and break down (e.g., stratospheric O₃ depletion due to N₂O, Portmann, Daniel, and Ravishankara 2012).

<figure><img src="image_5.png" /><figcaption>Figure 3 Atmospheric lifetime of different species from short-lived to long-lived species (Image taken from Seinfeld and Pandis, 2016)</figcaption></figure>

Therefore, they are much more mixed in the atmosphere and spread over the entire Earth. On the other hand, substances with a relatively shorter lifetime do not find a chance to be transported far from their emission sources. These short-lived species can be detected close to their emission source compared to the long-lived species before interacting with other substances. Hydroxyl radical (OH) is the freest reactive radical and a prominent short-lived species with a lifetime only of seconds. It reacts with many species and controls the oxidizing capacity of the atmosphere. Other species such as Isoprene (C₅H₈), Nitrogen Oxides (NOx) and Sulfur Dioxide (SO₂) have lifetimes from hours to days. They are mainly found at the lower levels of the troposphere in urban and regional scales and considered as air pollutants. They are important for the assessment of air quality.
Atmospheric lifetime of the species depends on various chemical and physical factors such that the presence of oxidizers (e.g. O₃, OH), the amount of incoming solar radiation and meteorological conditions (e.g., humidity, temperature). For instance, the amount of incoming solar energy (so-called insolation) differs from tropics to higher altitudes due to the angle of the Sun above the horizon (so-called solar elevation angle). This also yields the differences in the duration of the surface to the solar energy exposure. In conclusion, these phenomena are linked to the seasonality and diurnally. Therefore, atmospheric lifetime of the species, for example in the case of NOₓ, shows differences according to the seasons and the time of the day. NOₓ has a short lifetime but studies have been shown a longer NOₓ lifetime in winter than in summer (e.g., 21h vs 5.9 h in winter and summer, respectively (Shah et al. 2020) and even a longer daytime lifetime (29h) than night-time lifetime (6.3h) for a winter (Kenagy et al. 2018).

## 1.3 NOₓ chemistry (mainly tropospheric)

Nitrogen dioxide (NO₂) and nitric oxide (NO) are together designated the NOₓ family. Most of the NOₓ is emitted in the form of NO through high-temperature oxidation of the nitrogen present in combustion air. At high temperatures, molecular oxygen (O₂) dissociates into its atomic states (R1) and contributes to a series of reactions with nitrogen (N₂) in order to generate NO (R2-3). This process is shown in the following chemical equilibriums:

$$O_2 \leftrightarrow O + O \qquad (R1)$$

$$O + N_2 \leftrightarrow NO + N \qquad (R2)$$

$$N + O_2 \leftrightarrow NO + O \qquad (R3)$$

After the production of NO, the rapid cycling between NO and NO₂, also known as the Leighton cycle occurs by following the reactions:

$$NO + O_3 \rightarrow NO_2 + O_2 \qquad (R4)$$

$$NO_2 + h\nu \rightarrow NO + O(^3P) \qquad (R5)$$

$$O(^3P) + O_2 + M \rightarrow O_3 + M \qquad (R6)$$

Noting that M is a generic term stands for a non-reacting molecule but absorbing the energy of the reaction and stabilizing the production molecule(s) of the reaction. During the daytime, NO is oxidized by O₃ (by the titration of ozone), forming NO₂ (R4). Later, NO₂ is photo-dissociated back by radiation with a wavelength between 290 and 420 nanometers (visible and near ultra-violet) to NO, leading to the production of O₃ (R5-6). It is worth noting that this NO-NO₂
diurnal cycle is in a steady photochemical state and does not contribute to the net O$_3$ production. However, with the existence of different radicals (e.g., hydroperoxyl radical (HO$_2$) and organic peroxy radical (RO$_2$)), the NO-NO$_2$ balance in the Leighton cycle can be broken down as the radicals contribute to the transformation of NO into NO$_2$ by the sequence of reactions. The sequence is usually triggered by the reaction of VOC and carbon monoxide (CO) with the OH radical, resulting in the production of HO$_2$ and RO$_2$ (R7-8).

$$ VOC + OH \xrightarrow{O_2} RO_2 + H_2O \qquad (R7) $$

$$ CO + OH \xrightarrow{O_2} HO_2 + CO_2 \qquad (R8) $$

Noting that R in the reaction is the generic chemical formula of any alkyl organic group members e.g., organic peroxy radicals for RO$_2$ and alkoxy radical for RO. They are closely linked to the formation of secondary organic aerosols (SOAs) (Peng et al. 2022; Molteni et al. 2018; Orlando and Tyndall 2012).

In the presence of NO, this is followed by the conversion of NO to NO$_2$ through the reactions with RO$_2$ and HO$_2$ radicals, which also generates OH (R9-10). Thus, in particular, R10 also increases the oxidizing capacity of the atmosphere.

$$ RO_2 + NO \rightarrow NO_2 + RO \qquad (R9) $$

$$ NO + HO_2 \rightarrow NO_2 + OH \qquad (R10) $$

Again, the photolysis of NO$_2$ (from R9 and R10) leads to O$_3$ production (R5-6). O$_3$ production, therefore, depends on the amount of NO available and the presence of radicals (RO$_2$ and HO$_2$).

NO affects the production of RO (R9) as it reacts with the RO$_2$ radical formed through VOC oxidation (R7). In other words, the interactions between VOC and NO$_x$ contribute to SOA production. Therefore, NO$_x$ has a role in dispersion and the formation of SOA (Lane et al.,2008). During the daytime, the main sink of NO$_x$ from the atmosphere is via the oxidation of NO$_2$ with OH to nitric acid (HNO$_3$) (R11). NO$_2$ also reacts with O$_3$ and forms the nitrate radical (NO$_3$) predominantly through the reaction of (R12).

$$ NO_2 + OH + M \rightarrow HNO_3 + M \qquad (R11) $$

$$ NO_2 + O_3 \rightarrow NO_3 + O_2 \qquad (R12) $$
$NO_3$ produced during the daytime is rapidly converted back to $NO_2$ or NO by photolysis or by the reaction with NO.

During the nighttime, there is no photolysis of $NO_2$ (as in R5) so all $NO_x$ converted to $NO_2$ due to R4, reacts with $O_3$ to produce $NO_3$ which is also considered to be an important night-time oxidant (R12). The $NO_3$ radical reacts with $NO_2$ and produces gaseous dinitrogen pentoxide, $N_2O_5$ (R13, also an equilibrium) which is a nocturnal reservoir of $NO_x$. It can photo chemically (only during the daytime) or thermally decompose back to $NO_2$ and $NO_3$. However, during the nighttime and in the relatively low temperature (enough to prevent thermal decomposition) it can form $HNO_3$ by hydrolysis (R14). R12-14 are also considered one of the major paths of the removal of $NO_x$ from the atmosphere.

$$NO_2 + NO_3 + M \leftrightarrow N_2O_5 + M \qquad (R13)$$

$$N_2O_5 + H_2O \xrightarrow{\text{aerosol}} 2HNO_3 \qquad (R14)$$

$HNO_3$ is extremely water-soluble and removed by precipitation (wet deposition) causing acid rain or taken up by plants or soil (dry deposition) (R11 and R14). $NO_x$ is then a precursor of atmospheric $NO_3$ and $HNO_3$, which contributes to soil acidification, eutrophication and aerosol radiative forcing (Galloway et al. 2004; Liao 2005).

$NO_x$ has a short lifetime (Section 1.2) which varies from a few hours to about several days. However, peroxyacetyl radical ($CH_3C(O)OO$) also reacts with $NO_2$ to form more stable peroxyacetyl nitrates (PANs) which act as a carrier and a reservoir for $NO_x$ by the following reactions (R15-17).

$$CH_3CHO + OH \rightarrow CH_3CO + H_2O \qquad (R15)$$

$$CH_3CO + O_2 + M \rightarrow CH_3C(O)OO + M \qquad (R16)$$

$$CH_3C(O)OO + NO_2 + M \rightarrow PAN + M \qquad (R17)$$

where $CH_3CHO$ (also called acetaldehyde) is a widespread gas naturally produced by plants and widely used in industries. PAN is not soluble in water and is not removed by deposition. It allows long-range transport of $NO_x$ at low temperatures in the higher altitudes of the troposphere for weeks. Once PAN begins to descent into the warmer lower troposphere, it releases the $NO_x$ due to the thermal decomposition (R18).

$$PAN (+\text{heat}) \rightarrow CH_3C(O)OO + NO_2 \qquad (R18)$$
This means that, even though NOₓ has a short lifetime it can be transported far from its sources or from polluted regions to remote regions (Figure 4).

<figure><img src="image_6.png" /><figcaption>Figure 4 Long-range transportation of NOₓ from PAN. (Taken from Solène Turquety's "Atmospheric chemistry and air quality" course slides, 05/10/2020)</figcaption></figure>

## 1.4 NOₓ as an air pollutant : discovery and impacts

Air pollution is defined as an atmospheric condition in which the presence of certain substances or air pollutants is higher than normal ambient levels because of anthropogenic or natural activities and has direct or indirect harmful impacts on habitats.

The family of NOₓ first introduced as an air pollutant in the early 1950s with the discovery of the chemist Dr. Arie Haagen-Smit in California. During this time, a smog event occurred in the city of Los Angeles, causing reduced visibility, eye irritation, heavy smell, and damage to vegetation. He discovered the smog event, which is caused by the increased number of vehicles in the traffic and the emissions of NOₓ species from these vehicles. In fact, the emissions of NOₓ and VOCs (e.g., isoprene, ethane) from vehicle exhausts reacted with the incoming solar radiation (or UV) to produce O₃ (which is the key component of the smog) and atmospheric aerosols, and caused so-called brownish-grey photochemical smog. This finding brought attention to the role of NOₓ in the photochemical mechanisms related to air pollution in the troposphere (Section 1.3).
Therefore, not only NOₓ as a primary pollutant harm the ecosystems and human health, but also secondary pollutants, such as O₃ and HNO₃, and SOA, which are a product of various chemical reactions through NOₓ, cause serious environmental and health problems.

In terrestrial ecosystems, exposure to high concentrations of NOₓ and O₃ has direct negative impacts on plant growth and the productivity of agricultural crops. Leaves, mainly through their stomata, absorb these pollutants directly, resulting in reduced seed production and crop loss (de Vries 2021). Moreover, when NOₓ reservoir gas HNO₃ dissolves in water, acid rain damages plants, aquatic ecosystems, and cultural heritage. In aquatic ecosystems, elevated NOₓ and NO₃⁻ concentrations can cause eutrophication and acidification in freshwater, affecting nutrient balances and oxygen levels in the water.

Both short and long-term exposures to high NOₓ and O₃ concentration levels affect various human body systems and organs, such as in the respiratory system from respiratory infections, asthma to lung cancer and in the cardiovascular system, from angina to heart diseases (Kampa and Castanas 2008).

NOₓ emissions have sophisticated impacts on climate change, both speeding up and slowing down the effects (Lasek and Lajnert 2022). Although NOₓ is a reactive trace gas, it is responsible for the formation of tropospheric O₃, which in addition to be a pollutant is also a greenhouse gas (GHG) (Mitchell 1989; Szopa et al. 2021). Ground level O₃ absorbs terrestrial and solar radiation, causing positive radiative forcing which results trapping of the radiation in the troposphere (Stordal et al. 2003). By this way, NOₓ indirectly affects climate warming through tropospheric O₃ production. On the other hand, NOₓ is also responsible for the formation of OH radicals or cleaning agents. On a global scale, OH reacts with CO to form carbon dioxide (CO₂). Nevertheless, it also reacts with CH₄, which is also one of the most abundant anthropogenic GHG after CO₂, and reduces the CH₄'s lifetime.

## 1.5 NO₂ monitoring tools

It was only in the late 19th century the first monitoring of air pollutants began with O₃ at the Pic du Midi Observatory in southwestern France with the intention of using O₃ as a medicine (Marenco et al. 1994; Palmer 2017).

After the discovery of the NOₓ family of being the most significant precursor of photochemical O₃ production and HNO₃, and their possible impacts, monitoring of NO₂, its vertical profiles,
and the global trend in NO₂ concentrations aroused interest for use in studies. Vertical profiles
can be derived using balloons (Berthet et al. 2017), airborne (Merlaud et al. 2011), and ground-
based instruments (e.g., Système d'Analyse par Observation Zénitale (SAOZ), (Kreher et al.
2020)) but their deployment is limited spatially and temporally. On the other hand, in situ
surface measurements are usually more accurate and provide measurements continuously over
longer period. However, an important concern of the surface in situ instruments is the
representativeness of the measurement as they can be influenced easily with the surrounding
environment (e.g., instruments situated at the stations nearby the traffic). Even if there are some
techniques to evaluate the measurement type (Flemming, Stern, and Yamartino 2005; Lachatre
et al. 2020), quantifying the representativeness remains difficult. Moreover, it can be costly to
set up a network of ground-based instruments to cover the area of interest and operate them
continuously. Therefore, observations from space are complementary to the other measurement
techniques for NO₂ monitoring especially thanks to their increased spatial and temporal
coverage. They are able to provide global view of NO₂. The analysis of satellite measurements
requires precise comprehension of the radiative properties of the atmosphere and sophisticated
methods to retrieve information on NO₂ amounts. Satellites do not measure surface
concentrations; they measure vertically integrated concentrations as columns. One of the
limitations of these observations are the clouds, which prevent observations. Long-term
continuous monitoring of NO₂ concentrations using satellites have been possible since the
1990s (Anenberg et al. 2022).

1.6 NO₂ concentration levels in the atmosphere

Official authorities (e.g., governments) or non-profit organizations (e.g., World Health
Organization (WHO)) determine thresholds and establish guidelines to keep ambient air
pollutant levels below hazardous levels. As NO₂ is an air pollutant, these regulations usually
include limits on the amount of NO₂ that can be emitted to keep ambient NO₂ concentrations
below certain thresholds. According to the WHO's latest report in 2021, air quality standards
for the protection of public health are defined as follows: ambient NO₂ concentrations should
not exceed an annual average of 10 µg/m³ and a 24-hour average of 25 µg/m³ (WHO Global
Air Quality Guidelines 2021). However, despite the WHO's relatively restrictive limits, the
European Environment Agency's (EEA) "Air quality in Europe 2022" report shows that, NO₂
concentrations are recorded above WHO levels in all Organization for Economic Co-operation
and Development (OECD) member countries (at total 38 countries) (European Environment
Agency 2022). The EEA's less restrictive threshold is set at 40 µg/m³ on an annual average. Similar to the EEA, China's ambient air quality standard limits (CAAQS) are 40 µg/m³ annual average and 80 µg/m³ 24-hour average for all areas (including urban and industrial areas) nationwide. And according to the National Ambient Air Quality Standards (NAAQS) of the United States Environmental Protection Agency (EPA), ambient NO₂ concentrations should not exceed 53 parts per billion (ppb) (or equivalent to 99.64 µg/m³ using the WHO's conversion factor) (EPA).

Figure 5 shows the global spatial distribution of tropospheric NO₂ columns from the satellite instrument, Ozone Monitoring Instrument (OMI) averaged over the years from 2005 to 2018 with a spatial resolution of 0.25° × 0.25°. It is quite apparent that the red spots with higher NO₂ concentrations are mostly located in the northern hemisphere, such as Eastern Asia (Eastern China, Japan), Western Europe (e.g., The Netherlands, Belgium, France), and the United States of America (USA) (e.g., New York City, Philadelphia). These places indeed have been observed the highest levels of NO₂ concentrations over the last two decades (Richter et al. 2005; R. J. van der A et al. 2006; Paraschiv et al. 2017; Jamali, Klingmyr, and Tagesson 2020).

<figure><img src="image_7.png" /><figcaption>Figure 5 Global map of mean tropospheric NO₂ concentrations from 2005 to 2018 as observed by the OMI satellite instrument. Figure is taken from Jamali, Klingmyr, and Tagesson 2020.</figcaption></figure>

Figure 6 represents the spatial distribution of the surface annual means of NO₂ concentrations (in µg/m³) in 2019 for Europe taken from the surface stations. The measurement data is based on the European air quality database. There are at total 3,480 stations available for this year. Of these, 230 exceed the EEA's limit standards (corresponds to 7 %), while 2820 (81 %) exceed the WHO's limits on an annual base.
<figure><img src="image_8.png" /><figcaption>Figure 6 Ambient annual mean of surface NO₂ Concentrations for 2019 (in µg/m³), extracted from https://www.eea.europa.eu/data-and-maps/dashboards/air-quality-statistics</figcaption></figure>

Similar to Figure 6, Figure 7 shows the spatial distribution of surface annual means of NO₂ concentrations (in µg/m³) in 2019 for China. There are 626 stations, of which 43 exceed the CAAQS recommended level (corresponds to 7 %) and 597 exceed the WHO recommended level (corresponds to 95 %).

Accordingly, a significant number of stations do not meet the requirements of WHO levels, while a significant number of stations fall below the CAAQS levels. Ultimately, the aim of policy makers or authorities is to implement mitigation policies to keep pollutant concentrations below these harmful levels. Where measurements cannot be made, for example for forecasting, or to evaluate the best mitigation policies, air quality models are used to simulate NO₂ concentrations in order to predict future ambient levels and impacts of the chosen policies.
<figure><img src="image_9.png" /><figcaption>Figure 7 Spatial distribution of annual mean surface NO₂ concentration in 2019 for China (Visualized by using data of the CNEMC network).</figcaption></figure>

It is also important to note that NO₂ concentrations show daily, seasonal and inter-annual trends. For instance, Figure 8 represents the trends of the monthly mean of the tropospheric NO₂ vertical columns over East Central China from 1996 to 2004 (Richter et al. 2005). NO₂ measurements were taken from the early phase of the satellite era, from the Global Ozone Monitoring Experiment (GOME) (Burrows et al. 1999; European Space Agency 1995) satellite instrument between 1996 and 2002, and then from the SCIAMACHY (Bovensmann et al. 1999) instrument for the years 2003 and 2004. The time series shows a slight increase in the NO₂ vertical columns from 1996 to 2004 for East Central China. Especially due to seasonality, higher values are observed in NO₂ measurements in winter months and these values have been increasing year by year (in particular after the 2000s).

<figure><img src="image_10.png" /><figcaption>Figure 8 Inter-annual monthly averages of NO₂ tropospheric vertical columns over East Central China measured by GOME and SCIAMACHY satellite instruments. Figure taken from Richter et al., 2005 (Shaded areas show the standard deviation)</figcaption></figure>
Looking at the global trends, NO₂ concentrations show a slight decreasing trend between 2004 and 2008, followed by an increase from 2008 to 2018, with a significant increase in the first 3 years (Figure 9, from Jamali, Klingmyr, and Tagesson 2020). However on a continental scale, this situation differs with a sharp decreasing trend in Europe from 2006 to 2008 and from 2016 to 2018. Still between 2008 and 2016 a slight increase is observed. China has increasing trends of NO₂ concentrations until 2011, followed by a sharp decrease until 2016 and a slight increase (or insignificant change) between 2016 and 2018. These figures show that despite the increasing global trends in the NO₂ concentrations the contribution of different regions to this global trend shows substantial regional differences. This may be a result of the implementation of different mitigation policies in different time periods.

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

Figure 9 Time series of annual mean of tropospheric NO₂ concentration using DBEST (Jamali et al., 2015) trend analysis method and OMI satellite measurements for the Globe, Europe and China. Figures are taken from Jamali, Klingmyr, and Tagesson 2020
For example, in China, the 12th Five-Year Plan, implemented by the Chinese government in 2011, aims to reduce annual NOₓ emissions by 10% within five years (de Foy, Lu, and Streets 2016). This is also consistent with the decrease in NO₂ satellite concentrations starting in 2011 as seen in Figure 9. In European Union (EU), the Gothenburg protocol has been conducted to improve combustion processes, thus reducing NOₓ emissions since 2005 (United Nations Economic and Social Council 2015). The task group of the Gothenburg protocol have developed policy scenarios in the framework of Long-range Transboundary Air Pollution, and set long-term goals to reduce NOₓ emissions by about 50 to 80% in the EU, North America, and the Western Balkans, and about 20% in the Eastern Europe, the Caucasus and Central Asia between 2005 and 2030. In addition, for instance, there are also global conventions such as the Marpol Convention, which aims to reduce global NOₓ ship emissions since 2005 (www.imo.org, Last Access: 10/09/2022).

These regulations commonly include limits on the amount of NOₓ that can be emitted by vehicles and industrial facilities, as well as requirements for the use of emissions control technologies.

## 1.7 NOₓ emissions

### 1.7.1 NOₓ emission sources

Humankind first began to influence atmospheric composition by emitting pollutants through agricultural activities (e.g., developing land for crops and transportation near water sources or along trade routes, and raising livestock for meat and dairy products). However, with the Industrial Revolution that began in the late 18th century, increasing anthropogenic activities began to increase air pollutant emissions significantly in industrialized regions. These activities became a dominant source of changes in the atmospheric compositions of air pollutants (Abram et al. 2016).

Increase in NOₓ emissions due to urbanization, industrial and economic development and consequently augmented energy consumption and transportation has a direct effect on the increase of NO₂ concentrations, although there is a non-linear relationship between NO₂ observations and NOₓ emissions (Vestreng et al. 2009; Konovalov et al. 2010; Huszar et al. 2021).
NOₓ emissions originate from both biogenic (or natural) and anthropogenic sources. In the nature, NOₓ can be produced by lightning with a considerable amount in the upper troposphere (more than 1.2 billion flashes occur annually (NASA/Goddard Space Flight Center 2009), globally at total 5 Tg N yr⁻¹ as in Table 1 or estimated between 2 and 8 Tg N yr⁻¹ by Schumann and Huntrieser 2007 (not mentioned in the table). Because during lightning, release of very high temperatures can split the molecular oxygen, which can latter react with abundant nitrogen (R1-3). Soil is also another natural source of NOₓ (Almaraz et al. 2018; Xiao Lu et al. 2021). NO is the predominant forms of soil NOₓ emissions produced by biogenic sources, such as microorganism activities in the soil via nitrification and denitrification processes. These processes can be affected by nitrogen and ammonia availability in the soil, and moisture, temperature, structure, type and pH of the soil. Therefore, human also influence soil NOₓ emissions through different agricultural activities such as using nitrogen fertilizer, irrigation, crop harvesting, and land clearance. Biomass burning (e.g., deforestation and agriculture) can also be both biogenic and anthropogenic sources of NOₓ emissions by burning of living or dead vegetation. Nevertheless, the significant contribution of NOₓ emissions comes from anthropogenic sources, particularly from combustion processes (Table 1). The combustion engines are included in the sectors like; energy sector (e.g., thermal power plants especially using fossil fuels such as coal as a raw material), industry, transportation (e.g., traffic, shipping, aircraft), residential (e.g., home heaters, gas stoves). Figure 10 visualizes the main NOₓ emission sectors of road transport, electricity and heat generation, and manufacturing and construction for EU countries in 2019 (European Union Emission Inventory Report 2021).

<table><caption>Table 1 Overview of global NOₓ emission budgets and their uncertainties in parenthesis for selected sources and sectors</caption><thead><tr><th rowspan="2">Source</th><th colspan="3">Estimated NOₓ emissions for different sources and sectors (Tg N yr⁻¹)</th></tr><tr><th>Lee et al. 1997<sup>a</sup></th><th>TAR<sup>b</sup></th><th>AR4<sup>c</sup></th></tr></thead><tbody><tr><td><strong>Anthropogenic Sources</strong></td><td></td><td></td><td></td></tr><tr><td>Fossil fuel combustion</td><td>22 (13-31)</td><td>33 (20-24)</td><td>25.6 (21-28)</td></tr><tr><td>Biomass Burning</td><td>7.9 (3-15)</td><td>7.1 (2-12)</td><td>5.9 (6-12)</td></tr><tr><td>Aircraft</td><td>0.9</td><td>0.7 (0.2-0.9)</td><td></td></tr><tr><td>Agriculture</td><td></td><td>2.3 (0-4)</td><td>1.6</td></tr><tr><td>Atmospheric Deposition</td><td></td><td></td><td>0.3</td></tr><tr><td><strong>Anthropogenic Total</strong></td><td><strong>30.8</strong></td><td><strong>43.1</strong></td><td><strong>33.4</strong></td></tr><tr><td><strong>Biogenic Sources</strong></td><td></td><td></td><td></td></tr><tr><td>Soils</td><td>7 (4-12)</td><td>3.3 (3-8)</td><td>7.3 (5-8)</td></tr><tr><td>Lightning</td><td>5 (2-20)</td><td>5 (2-12)</td><td>1.1-6.4 (3-7)</td></tr><tr><td><strong>Biogenic Total</strong></td><td><strong>12</strong></td><td><strong>8.3</strong></td><td><strong>8.4-13.7</strong></td></tr><tr><td><strong>Total</strong></td><td><strong>42.8</strong></td><td><strong>51.4</strong></td><td><strong>41.8-47.1</strong></td></tr></tbody></table>

<footnote>a Best estimate source and uncertainties adapted from Table 3 in the paper.</footnote>

<footnote>b The IPCC Third Assessment Report (TAR) (Houghton et al. 2001), Table 4.8.</footnote>

<footnote>c The IPCC Fourth Assessment Report (AR4) (Core Writing Team, Pashauri, and Reisinger 2008), Table 7.7.</footnote>
<figure><img src="image_14.png" /><figcaption>Figure 10 EU Anthropogenic NOₓ emission sources by sectors for 2019 (Taken from EU emission inventory report, 2021)</figcaption></figure>

## 1.7.1 NOₓ emission inventories

A series of emission sources or inventories are compiled into a comprehensive report accounting the amount of various pollutants released into the atmosphere from different sources (or emitters) over a period, usually annually. This so-called emission inventories are then used to track changes in pollutant levels over time, assess the effectiveness of pollution reduction policies and regulations, and provide data for scientific research on the impacts of pollution on the environment and human health. They are generally compiled by government agencies, but can also be compiled by private companies, non-governmental organizations, and research institutions.

### 1.7.1.1 Bottom-up approach

There are different ways to build emission inventories. The so-called bottom-up approach is one of the main and widely used approaches. Inventories using bottom-up method is based on the individual emitters. Using this approach, emissions are calculated, sector by sector, by multiplying the amount of pollutants emitted per unit of activity (such as fuel consumption or electricity generation) by the number of activities that occurred during a given period. This requires detailed information on the characteristics of each source, such as the type of fuel used, emission activity rate, emission factor, the efficiency of the technology, and the operating conditions. Due to the high number of emission calculations, and input data collected from one-to-one sources, it takes time, effort and becomes complicated to build a complete bottom-up inventory of emissions. And even more complicated for a larger region. For instance, estimating the amount of NOₓ emissions from an individual car is a challenging task, since it depends on
various factors such as, the car's model, the speed during the trip, the length of the route, and the maintenance of the vehicle. This becomes even more complicated when one tries to calculate the total emissions from every car in a country. Then combining this information from each activity and accurately reporting the true total amount of anthropogenic emissions involves high uncertainties. These uncertainties come from the activity data, the emission factor, spatial allocation and timing of emissions. For example, Zhao et al. 2011 estimated the uncertainty of national, bottom-up inventory for NOₓ emissions using the Monte Carlo approach to 37% for China mainly due to the oversimplified source classifications and roughly estimated emission factors. However, despite the uncertainties, bottom-up inventories provide highly detailed information. Therefore, this bottom-up approach is often used in national and regional emissions inventories.

### 1.7.1.2 Top-down approach

In contrast, the top-down approach estimates emissions by measuring pollutant concentrations in the atmosphere and using models to inverse the emissions that would be necessary to produce those concentrations. In particular, recent advances in satellite technology have enabled the use of observations from these sources, and the increase in computational power has made it possible to solve the highly complex problems. As a result, top-down approach is conducted using mathematical atmospheric inversion techniques with modeling and observations (widely used satellite observations) to constrain emissions on high performance computer clusters. Despite the fact that, it relies on assumptions in the models and uncertainties come from the model and instruments, the use of satellite observations brings some advantages. For instance, the observations cover a larger area at the same time (spatial) and provide continuous measurements over the area of interest on a daily base (temporal). Also, the increase in the resolution of satellite observations helps to better capture emission sources (higher resolution). All of these advantages helps to improve uncertainties contained in bottom-up inventories. For this reason, establishing top-down emission inventories by exploiting the potential of satellite observations is promising and still developing research area of interest.

### 1.7.1.3 Examples and comparison of two approaches

National or regional inventories are built using local pieces of information, while global inventories are created using the mosaic method by combining these inventories from different countries or regions into a uniform emissions inventory over a larger domain (Figures 11 and
12). For example, the Multiresolution Emission Inventory for China (MEIC) (Liu et al. 2017; Zheng et al. 2018) v1.0 bottom-up inventory, with a 0.25° × 0.25° spatial resolution, covers China as in Figure 11 for the period 2007 to 2017, and is used together with other local inventories to compile the regional MIX inventory for the whole of Asia (Li, Zhang, et al. 2017a). Furthermore, used in the mosaic of regional and global HTAP v2.2 inventory as in the Figure 12 (Janssens-Maenhout et al. 2015).

<figure><img src="image_15.png" /><figcaption>Figure 11 Compilation of regional MIX Asia bottom-up emission inventory with local bottom-up inventories (Figure is taken from Li M. et al. 2017)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 12 Compilation of several regional bottom-up emission inventories to build a global HTAP v2.2 emission inventory (Figure is taken from Janssens-Maenhout et al., 2015)</figcaption></figure>

The Emissions Database for Global Atmospheric Research (EDGAR) (Janssens-Maenhout et al. 2015) v4.3.1 (Crippa et al. 2016) is another bottom-up inventory with a resolution of 0.1° × 0.1° and covers a relatively large domain for the period from 1970 to 2010 as in Figure 12 compared to other inventories. CAMS (Copernicus Atmosphere Monitoring Service) operates
the ECCAD (Emissions of atmospheric Compounds and Compilation of Ancillary Data) database, which contains various emission inventory datasets at the global (1° × 1°) and regional scales (up to 0.05° × 0.1°). It includes national EMEP, TNO bottom-up inventories for Europe as well as, inventories obtained by using top-down approach.

Figure 13 shows an inter-annual comparison between different emission inventories: bottom-up inventories (REAS v1.2, REAS v2.2, MEIC and EDGAR) and top-down inventories (DECSO-OMI, DECSO-GOME2a, EnKF-MIROC and EnKF-CHASER) for the mainland China (taken from Ding et al. 2017).

<figure><img src="image_17.png" /><figcaption>Figure 13 Annual time series of NOₓ emissions from different emission inventories over mainland China. (Figure taken from Ding et al., 2017)</figcaption></figure>

According to the figure, top-down inventories (dotted lines) provide recent information on the emissions compared to the bottom-up inventories (plain lines), which is related to the difficulties of the real-time updates of bottom-up inventories. Because time is required for collecting and processing information from each emitter. Moreover, bottom-up inventories such as MEIC and Regional Emission Inventory in Asia (REAS) v2.2 contain relatively short-term information, such as a couple of years. Furthermore, all emission inventories show differences in the amount of annual NOₓ emissions compared to each other. The intersection of the values for some years only exists in top-down inventories (e.g., DECSO-GOME2a and EnKF-MIROC in 2011.)

To conclude, bottom-up emission inventories are subject to uncertainties and limitations (e.g., 50.7 % in NOₓ emissions for 2012 in EU28 (Crippa et al. 2018). However, accurate and up-to-date emission inventories are crucial for understanding the contribution of various human activities, modeling and projecting the associated changes in the atmospheric composition, and
developing the most feasible strategies to mitigate emissions. To this end, benchmarking studies
have shown that although top-down inventories also consist of uncertainties, they can provide
complementary information to improve bottom-up emission inventories by using their
advantages (N. Elguindi et al. 2020).

**1.8 Objectives of the thesis**

As shown previously, top-down approaches are promising to complement bottom-up
inventories providing regional to global information with a fast update capability. These
methodologies have reached a good level of maturity in the last decade with the development
of several tools (Müller and Stavrakou 2005; Miyazaki et al. 2017; Qu et al. 2017; Fortems-
Cheiney et al. 2021) and new satellite observations. At the beginning of this PhD work, the new
satellite instrument TROPOMI was starting to provide NO₂ tropospheric columns at an
unprecedented spatial resolution of a few kilometers. At the French level, a new inversion
system was in development, the Community Inversion Framework (CIF) (Berchet et al. 2021)
based on a previous less versatile version of PYVAR (Fortems-Cheiney et al. 2021). This tool
was first dedicated to the inversion and monitoring of greenhouse gases (CO₂, CH₄) at the global
scale. Capabilities to apply top-down approaches to highly reactive species at the regional scale
and the high resolution were limited. In this context, the main objective of the ANR
ARGONAUT project was to address the question of the French's capabilities to develop the
regional inverse modeling at the high resolution of reactive species, taking advantages of the
CIF development, the national regional model CHIMERE, and the newly TROPOMI
observations. The PhD work has been developed in this framework and co-funded a grant from
CNES. The main objective of my work was to test the capabilities of the CIF tool coupled with
the CHIMERE model. To do so, I focused my work mainly on NOₓ emissions because NO₂
satellite products are the most mature in terms of product for primary pollutants with reasonable
uncertainties and good temporal resolution (daily). Considering atmospheric chemistry, NOₓ
chemistry is also better known and represented in the model than other species, such as VOCs.
NOₓ inversion was then a good candidate to start with the evaluation of a new inversion system.
The objective was to test the CIF-CHIMERE in various configurations (domains, resolutions,
and observations) to evaluate its performances and conditions of work. Three domains have
been considered Europe, France, and China, two at medium resolution (50 km for Europe and
China) and one at high resolution (10 km for France). The objective was also to evaluate,
somehow, the gain of using observations at different resolutions, comparing inversion with
OMI and TROPOMI.

Scientifically, I focus my work first on Chinese NOₓ emissions. Indeed, the Chinese government developed new mitigation strategies for these emissions during the last decades. I evaluated the performances of the new CIF-CHIMERE system and the NO₂ observations from OMI to follow and quantify the emission reductions for two years (2015 and 2019) at the end of each 5 Year-Plan of the Chinese government. It was not possible to use the new observations from TROPOMI for this study as they started only in 2018. Secondly, I made preliminary tests to apply the CIF-CHIMERE in Europe and France using TROPOMI and OMI to evaluate improvements given by the high resolution of TROPOMI and also to analyze the potential limitations. I focus mainly on the first lockdown period in April 2020 due to the COVID-19 pandemic.

In the second part of the thesis (Chapter 2), the “Theoretical” Tools used in this PhD work are described: the historical development and theoretical bases of inverse modeling and the state-of-the-art CIF inverse modeling system are discussed, as well as the CHIMERE CTM, including the model inputs, and the different satellite and surface observations used. The various tests done to characterize the system and set up the most appropriate configuration are also presented.

Chapter 3 is devoted to the first case study, the inverse modeling of NOₓ emissions over China using OMI satellite observations, followed by Chapter 4, which presents the preliminary tests in Europe and France applied to the COVID-19 shutdown period, and discusses expected improvements from TROPOMI compared to OMI, difficulties and possible limitations. A general discussion of the results and findings and suggestions for the extension of this work concludes the manuscript.
# Second part : Methodology and
# Instruments

**Chapter 2 : Inverse modeling and its
tools**

In this thesis, the variational data assimilation technique is applied to inverse NOₓ emissions using NO₂ satellite observations. Therefore, this chapter describes the development of the data assimilation theory under the framework of inverse modeling with more emphasis on the variational approach. Moreover, it also includes descriptions of the models and tools used during the inversion and the data sources used to assess the inversion results.

**2.1 The concept and theory**

**2.1.1 What is an inverse problem?**

To predict the state or future of an atmospheric system, it is necessary to build physical or
mathematical models. A direct or forward problem finds the effects of a given cause (e.g., used
as inputs in a model) using these models, while an inverse problem does the reverse and aims
to find the cause of a given effect. In 1976, the American mathematician J.B.Keller called two
problems inverse if “the formulation of each problem contains all or part of the solution of the
other” (J. B. Keller 1976). However, the main difference between these problems is that direct
problems are well-posed, while the inverse problems are ill-posed (e.g., there is no solution of
the problem, or no unique solution, or the solution is unstable). This concept was introduced by
the French mathematician Hadamard (Hadamard 1923). On the other hand, regularization or
penalization methods (e.g., Tikhonov regularization (Tikhonov 1963)) are designed to solve
these ill-posed inverse problems and aim to find an existing, unique, and stable (continuous
dependence on the parameters) solution of the problem (Hansen 2010). General information
about inverse problems and their uses, including other fields of science, can be found briefly in
Kabanikhin's survey paper (Kabanikhin 2008).
## 2.1.2 Overview of data assimilation: past to present

### 2.1.2.1 Assimilation as an inverse problem

Inverse modeling or state-estimation theory is often referred to as data assimilation in geosciences (Carrassi et al., 2017), and data assimilation methods are used to solve inverse problems. For example, an inverse modeling that aims to find a set of model parameters (e.g., emissions) that best fits the observed data (e.g., concentrations) requires the assimilation of observations to estimate the system's state. Asch, Bocquet, and Nodet 2016 describes in detail the strong mathematical connections between techniques in inverse problems and data assimilation methods. The mathematical theory behind data assimilation comes from the 18th-century English statistician Thomas Bayes. His idea simply consists of modifying our prior belief with objective information; therefore, the posterior updated belief combines both prior and that new information (Bayes 1763). He described a theorem (after named the Bayes' theorem), which deals with the conditional probability of an event (y) based on the values of known probabilities of prior event (x) or in other words, to deal with the probabilities of the combination of states of information following the *Bayes' rule* (Equation 1):

$$p(x|y) = \frac{p(y|x) p(x)}{p(y)} \qquad (\text{Eq.1})$$

Tarantola has extensively elaborated this probabilistic approach to inverse modeling problems, from Bayes up to the day of his book's publication in 2005 (a rewritten and the latest version of the book from 1987 by the same author) (Tarantola 2005).

Data assimilation aims to combine optimally all the information from model simulations and observations, taking into account their uncertainties, to produce a better estimate of systems representing reality. This estimation process can aim to find the best estimate of the parameters of a model or the state of the system based on available observations. However, during the process, there are issues that make data assimilation challenging and add some difficulties in the calculations of the estimations: (i) systems are complex (e.g., the non-linearity of atmospheric systems or poorly known interactions). (ii) observations are noisy, sparse, and irregularly distributed in space and time (heterogeneous observations). (iii) the problem is ill-posed (regularizations modifies the main problem; thus the result would be the approximation of the true solution).
The invention and development of the least squares method for studying the position of stars by Gauss and Legendre in the early 1800s happened to be a pivotal event for today's data assimilation techniques (Stigler 1981) and also for the probability theory. The idea of the method was based on the minimization of the sum of the squares of the residuals (or errors) that is the departure between an observed value and the fitted value provided by a simulation model. In addition to that, the Gauss-Markov theorem guarantees that this minimization is optimal under some assumptions (e.g., the errors are uncorrelated, have zero mean and have equal variances) in a linear model and the best linear unbiased estimator (BLUE, described in detail later) of the coefficients is the least-squares estimator. While these advances in data assimilation theory were groundbreaking, atmospheric models and applications of atmospheric data assimilation using these models were still in their infancy.

### 2.1.2.2 First usages of data assimilation in atmospheric science

Data assimilation was first used in numerical weather prediction (NWP) models with the aim of short-term forecasting of meteorological conditions in 1960s, while the history of NWP models goes back much further. At the beginning of the 20th century, Cleveland Abbe described the physical laws of the atmosphere (e.g., principles of thermodynamics and hydrodynamics) in the need of forecasting the meteorological conditions (Abbe 1901). Vilhelm Bjerknes suggested an estimation or prediction of weather could be as a deterministic initial value problem subject to its initial conditions such as temperature, humidity, pressure, density (Bjerknes 1904). He developed a graphical method for solving a system of non-linear partial differential equations. He also found that solving these equations numerically was challenging, and no analytical solutions were possible (Lynch 2008). Then in 1922, Lewis Fry Richardson, an English meteorologist, influenced by Bjerknes, brought a novelty to the weather prediction by applying numerical methods (NWP models) and solving the equations analytically (Richardson 1922). He proposed the first data assimilation technique by interpolating the observations. However, in his analysis, the observations had not been properly assimilated because he observed significant errors due to the inaccuracies of the fields used as the initial conditions. Nevertheless, in early 1950, under the leadership of an American meteorologist Jule Charney, the development of the first generation electronic computer (called the Electronic Numerical Integrator and Computer, ENIAC), the necessary mathematical analysis and design
of a numerical algorithm made it possible to have a first successful 24-h NWP forecasts (Lynch
2008), despite the fact that, in their experiment they used Richardson's simplified equation set.

In the late 1950s, George Cressman's successive correction algorithm based on the method of Bergthörsson and DÖÖ brought the idea of weighting the background by observations, and adding linear combinations of the observational increments (observation minus background) to the backgrounds in three dimensions (3D) (Cressman 1959; Bergthörsson and DÖÖs 1955). These weights were determined statistically. They used previous forecasts for the initialization (first guess) of the model to provide complete information (background) over the domain. The algorithm combined the observations with models in a statistical sense via the least square method and with the intention to find a 'best fit', which is called optimal interpolation (OI).

2.1.2.3 Sequential approach: the Kalman filter and
its variants

By the 1960s, OI methods began to be globally operational in many meteorological centers to
make weather predictions. The optimal interpolation algorithm is a pioneering and simpler
version of the Kalman Filter (KF) method (Kalman 1960). KF and its variants are one of the
two main data assimilation approaches called sequential (statistical) data assimilation. KF is an
optimal linear estimator of the parameters, updating the previous optimal estimates as new
observations (noisy) available (therefore, sequential) and providing the uncertainty of the
analysis at each step. The detailed descriptions and the analysis of the method can be found in
Maybeck and Siouris 1980. However, KF is only optimal in the case of Gaussian statistics and
linear systems, so it could not be directly applied to nonlinear systems. As the real life problems
involve non-linear functions, the extension of KF that includes Taylor series expansion, to help
to get a linear approximation of the non-linear functions was necessary. This is through an
approximation of the non-linear function using first derivative of Taylor series called Jacobian
matrix. This version of the filter is called the extended Kalman Filter (EKF) and was developed
by Jazwinski 1970. However, the numerical cost for large atmospheric systems is very high and
the dimensionality of the computational grid and storage cost of the matrix led to the subsequent
development of other KF methods, such as the ensemble Kalman filter (EnKF) developed by
Evensen 1994. The method generates an ensemble of initializations to be used to ensemble
forecast in order to reduce the dimensionality. In the end, this method is a simplification and
close approximation to the KF and can provide similar results if the ensemble size is large
enough and if the model errors are well represented by the ensemble. Houtekamer and Mitchell
implemented EnKF for atmospheric data assimilation for the medium-range ensemble
forecasting systems in various of their work (Houtekamer and Mitchell 1998; 2001;
Houtekamer et al. 2005).

Assimilating data to predict the chemical state of the troposphere is even much more
computationally costly than weather forecasting. This is because atmospheric chemistry
involves chemical reactions and the number of parameters taken into account in the calculations
is greater than in meteorology (even higher if aerosol chemistry is considered too). The
application of the data assimilation techniques to the atmospheric trace gases only began in
1993 with Hartley and Prinn. They investigated the feasibility of using KF method to determine
regional surface fluxes for CFCs by comparing observations and models in 3D and found that
the higher differences between the model and observations lead higher biases in the estimations
(Hartley and Prinn 1993). Haas-Laursen, Hartley, and Prinn 1996 used KF to constrain trace
gases and to investigate their temporal source/sink cycle. Napelenok et al. 2008 developed a
method using KF, to constrain current NOₓ emission inventories by NO₂ column observations
from GOME satellite instrument. Mijling and van der A 2012 developed an algorithm (so-called
DESCO) using KF method for short-lived species. Miyazaki et al. 2017 estimated global
surface NOₓ emissions for a long term from 2005 to 2014 using EnKF and various observation
sources (e.g., OMI, SCIAMACHY satellite instruments).

2.1.2.4 Variational approach

Meanwhile, in addition to the sequential approaches, in the 80-90s, so-called variational
approaches using optimal control theory were introduced into the 3D variational (3D-Var)
method and later into the 4D variational (4D-Var) method, which considers the spatial and
temporal distributions of observations (Phillips 1982; Dimet and Talagrand 1986; Talagrand
and Courtier 1987; Courtier et al. 1998; Rabier et al. 1998). The 3D-Var and 4D-Var approaches
aim to combine observations and background information in an optimal way (indeed using
optimization techniques e.g., descent algorithm) to produce the best possible estimate of the
state of the system (their formulation is explained in the following sub-sections). It should be
noted that, under linearity assumptions, Kalman filter method, its variants and the 4D-Var
provide the same results at the end of the assimilation period (Lorenc 1986), since they all
satisfy BLUE (explained in the section 4D-Var and its pioneering 3D-Var).
First, in 1995 Fisher and Lary, performed 4D-Var approach for the stratosphere to chemically active trace gases using a reduced photochemical box model and satellite measurements using NO₂ observations from the Cryogenic Limb Array Etalon Spectrometer (CLAES) Instrument onboard the Upper Atmosphere Research Satellite (UARS) (Fisher and Lary 1995). They developed the adjoint model of reduced stratospheric gas phase mechanism combined with the box model. Later, Eskes et al. 1999 implemented 4D-Var approach using 2D model and total ozone column data from GOME to estimate global ozone field, in particular over the poles.

For the troposphere, first in 1997, Elbern et al. developed and studied the feasibility of the 4D-Var with the adjoint model of a tropospheric gas phase mechanism by assimilating the O₃ ground-based measurements with the assimilation window of hours only (Hendrik Elbern, Schmidt, and Ebel 1997). They mentioned the necessity of a 3D meteorological modeling and its adjoint to include diffusion and convection operators, and cloud parameters. Few years later, in 1999, Elbern and Schmidt presented the first application of the 4D-Var method for a 6h data assimilation window with a full comprehensive regional 3D Eulerian CTM and its adjoint, including chemical reactions, transport, and diffusion by using surface observations such as NO₂ and O₃ to obtain a chemical state of the troposphere (Hendrik Elbern and Schmidt 1999). They found the 4D-Var method to be feasible and useful but emphasized the importance of having appropriate prior knowledge on the concentration levels to be used as a background. So far, the parameters to be optimized were initial values, such as concentration levels. However, H. Elbern et al. 2000 optimized the emissions within an assimilation window of 6h. They investigated the performance of the 4D-Var to estimate NOₓ and VOCs' (O₃ precursors) emission rates by considering SO₂ and O₃ surface observations. Hendrik Elbern and Schmidt 2001 focused on the optimal analysis of chemical evolution of the ozone episode occurred over central Europe. This was the first 4D-Var application on a real case event by assimilating the surface observations of NOₓ and O₃ as well as O₃ radiosonde records with an extended data assimilation window to 9h. However, a longer assimilation window, such as 24 hours representing a full day, can be a strong constraint as it can reflect the diurnal profile of emissions. Thus, H. Elbern et al. 2007 performed 4D-Var within 24h assimilation window to estimate emission rates of all gaseous species emitted at each surface grid point using surface in-situ observations and a full comprehensive regional CTM over Europe.

In addition to the ground-based networks used in the assimilations, dedicated flight missions and various satellite instruments also used as a source of independent information. For example,
Talbot 2003 performed the first tropospheric airborne assimilation using 4D-Var to determine reactive nitrogen species over the Western Pacific. Müller and Stavrakou 2005 performed the first co-assimilation of annual global emissions (at continental scales) of different ozone precursor chemical species (CO and NOₓ) within 6 hours of assimilation window, using observations from the GOME satellite instrument. The simultaneous assimilation of CO and NOₓ can account better the chemical response of tropospheric ozone production cycle as they are both taking a role in the cycle. However, the reliability of the optimized emissions are not only depend on the independent observations but also chemical feedback between other species.

Kurokawa et al. 2009 developed an assimilation system with a parameterized NOₓ chemistry scheme in order to optimize NOₓ emissions over eastern China using GOME instrument. They run an assimilation system over a daily assimilation window containing observations from different times. Since then, there have been many studies that have widely used the 4D-Var assimilation approach in air quality assessments using different satellite observations and CTMs at a spatial resolution of a few kilometers to several hundred square kilometers (from a point source to a global) according to the interests (Pison, Menut, and Bergametti 2007; Boisgontier et al. 2008; Chai et al. 2009; Fortems-Cheiney et al. 2021; Qu et al. 2022).

In the last decade, hybrid approaches have received increasing attention and have become a popular research area (e.g., Ensemble 4D Variational (4D-EnVar)). These approaches (called hybrid) combine the two main methods, variational and statistical, by aiming to exploit the advantages of both methods, and they are still evolving (Bannister 2017). Moreover, there are studies in which these approaches are supported by deep learning models with the aim of providing more robust results (C. A. Keller and Evans 2019; Petetin et al. 2020; He et al. 2022).

### 2.1.3 General notations of data assimilation

This section defines standardized notations and describes the main elements of data assimilation that form the basis of the model theory used in this study. Ide et al. 1997 proposed these standardized notations to help the development of more advanced data assimilation methods.
### 2.1.3.1 The state vector, control variable and observations

#### 2.1.3.1.1 The state vector

Let us consider a discrete dynamical model (e.g. CTMs) for a physical system (e.g., atmospheric, oceanic) that evolves over time. The state of this model is represented by a collection of numbers equivalent to a column matrix called the **state vector** $x$ of dimension $n$ in the state space $X$. If we define $t_i$ as a time index at the $i$-th time, $x(t_i)$ and $x(t_{i+1})$ become two state vectors at the consecutive times $t_i$ and $t_{i+1}$. This system from time $t_i$ to time $t_{i+1}$ can be described as follows:

$$x(t_{i+1}) = \mathcal{M}_{t_i \to t_{i+1}} x(t_i) \qquad (\text{Eq.2})$$

where $\mathcal{M}$ is the dynamics operator of a computer model. It can be non-linear and consists of a discretization of non-linear partial differential equations (PDEs). It is assumed that the model is deterministic and perfect, so more or less it represents the reality to avoid considering the model error. However, a model is limited by its resolution, so it cannot simulate the real or **true state** ($x^t$) of the system. Therefore, the state vectors represent an approximation of the state of the atmosphere. Before the assimilation experiment, the true state can be estimated by the so called first guess of the true state or **background state** ($x^b$). This is often called initial or prior state and mostly derived from a CTM model, a previous prediction or a climatology before the **analysis** ($x^a$) is carried out.

#### 2.1.3.1.2 Control variable

In the complex systems, it is not easy to solve the analysis problem for all components of the state vector. It is because either the interactions between components are not completely well known (e.g., chemical interactions in an atmospheric model) or we are limited with the computational power. For this reason, the analysis problem corresponds to the search for a correction or **increment** ($\delta x$) of the background which belongs to a subspace called **control variable** space as follows:

$$x^a = x^b + \delta x \qquad (\text{Eq.3})$$

where $x^a$ is as close as possible to $x^t$. By this transition, essentially, our problem is to search for $x^a - x^b$ instead of searching for $x^a$.
### 2.1.3.1.3 Observations

The collection of all available observations (e.g., from spectrometers, radars, in-situ measurements, satellite instruments) are defined by a variable **observation vector y** in the observation space **Y**. All the parameters of the control space may not be observable and are not necessarily observed. So, there are fewer observations than parameters in the model. Moreover, to use the information contained in the observation vector, they need to be comparable with the state vector. However, observations are not regularly disposed and they do not correspond to a point of the grid of the model. Therefore, in order to compare the observations to the model, a function of **H** or an **observation operator** is defined. It maps the state space onto observation space. **H** is supposed to be known, most of the time it is assumed to be linear. If it is linear, then a matrix **H**, called tangent linear observation operator, exists. However, it can also be complex and non-linear.

Noting that, if the dimension of the observation vector **y** at *i*-th time is **p**ᵢ, usually **p**ᵢ ≪ **n**, (*n* being the dimension of the state space) because all the parameters of the state space may not be observed. There are not enough observations to estimate the state of space accurately. This fact is causing the inverse modeling being an ill-posed problem. Nevertheless, it is important to note that for some applications the dimensions of both state vector and observation vectors are quite large (≈ e.g., from 10⁶ to 10⁸).

### 2.1.3.2 Analysis, innovation and residuals

Analysis (xª) is an estimation process of the true state of a system for a given time using a combination of M computational model, y independent observations from instruments and xᵇ background information as well as initial and boundary conditions. It is the result of the so-called **data assimilation** process. During 4D-Var assimilation process, the **assimilation window** defines the time window of all observations considered at once. Data assimilation uses all this information to estimate the best possible representation of the true state. An important element of an analysis is the use of discrepancies between the observations and the state vector mapped by the operator as follows:

$$y - \mathcal{H}(x) \qquad (\text{Eq.4})$$

When it is calculated using xᵇ, the difference is called **innovations** and when xª is used the term is called **analysis residuals**.
### 2.1.3.3 Errors

One of the crucial elements in data assimilation is the representation of the unknown errors associated with the uncertainties found in the background, observations and analysis. Even if it is difficult to determine the exact value of the errors, it is possible to characterize them using probability density functions (PDFs) and using statistics associated to the PDFs. Therefore, it is not expected to know what the exact error is but statistics are calculated for the estimation of this error. The statistical formulas are recalled in Annex A: Statistical formulas.

**Background errors** represent the difference between the background state and true state of the system as follows:

$$ \varepsilon^b = x^b - x^T \qquad (\text{Eq.5}) $$

Its covariance **B** (background error covariance matrix) is represented by

$$ \mathbf{B} = \mathbb{E}(\boldsymbol{\varepsilon}^b \boldsymbol{\varepsilon}^{b^T}) = (\overline{\boldsymbol{\varepsilon}^b - \overline{\boldsymbol{\varepsilon}^b}})(\boldsymbol{\varepsilon}^b - \overline{\boldsymbol{\varepsilon}^b})^T \qquad (\text{Eq.6}) $$

where $\mathbb{E}$ represents the mathematical expectation, $T$ represents the transposition and the average is $\overline{\boldsymbol{\varepsilon}^b}$. The diagonal terms of the **B** matrix are the variances and give the confidence that we have in the background information. The off-diagonal terms of the matrix give how this information propagates spatially from the background knowledge locations. Therefore, they critically influence the outcome of the inversion. Background error consists of misrepresentation of the physical processes of the model or systematic and random errors. However, constructing a **B** matrix is a challenging task. Sometimes there are experts who determine and make available these values using different methods (e.g. Monte Carlo simulations) to determine the effects of uncertainties in the background information.

**Observation errors** represent the difference between the observations and the corresponding state in the reality as follows:

$$ \boldsymbol{\varepsilon}^o = \mathbf{y} - \mathcal{H}(\mathbf{x}^T) \qquad (\text{Eq.7}) $$

with the average of $\overline{\boldsymbol{\varepsilon}^o}$. Its covariance **R** (observation error covariance matrix) is represented by
$$ \mathbf{R} = \mathbb{E}(\boldsymbol{\varepsilon}^o \boldsymbol{\varepsilon}^{oT}) = \overline{(\boldsymbol{\varepsilon}^o - \overline{\boldsymbol{\varepsilon}^o})(\boldsymbol{\varepsilon}^o - \overline{\boldsymbol{\varepsilon}^o})^T} \quad (\text{Eq.8}) $$

R is assumed to be diagonal because observations are considered mutually uncorrelated from each other. It contains instrumental errors (measurement error and errors raising during the conversion of satellite measurement into concentrations or columns), the modeling errors comes from the definition of the observation operator (H or H), and representativeness errors due to the discretization of PDEs.

Analysis errors represent the difference between the analysis and the true state of the system as follows:

$$ \boldsymbol{\varepsilon}^a = \mathbf{x}^a - \mathbf{x}^t \quad (\text{Eq.9}) $$

with the average of $\overline{\boldsymbol{\varepsilon}^a}$. A measure of amplitude $||\boldsymbol{\varepsilon}^a - \overline{\boldsymbol{\varepsilon}^a}||$ is provided by the trace of the analysis error covariance matrix A,

$$ \mathrm{Tr}(\mathbf{A}) = \overline{||\boldsymbol{\varepsilon}^a - \overline{\boldsymbol{\varepsilon}^a}||^2} \quad (\text{Eq.10}) $$

The optimization results with in an analysis error smaller than the background error. The averages calculated over the errors are called biases.

### 2.1.3.4 Hypotheses

For simplification, the concept and the techniques of data assimilation rely on some assumptions as follows (Lorenc 1986; Daley 1991; Ghil 1989):

* Tangent linear hypothesis (TLH) allows the observation operator $\mathcal{H}$ to be linearized (denoted by H), e.g., the variables ($\mathbf{x}$) of the observation operator in the neighborhood of the background ($\mathbf{x}^b$) are assumed to be linear:

$$ \mathcal{H}(\mathbf{x}) - \mathcal{H}(\mathbf{x}^b) = \mathbf{H}(\mathbf{x} - \mathbf{x}^b) \quad (\text{Eq.11}) $$

* Errors $\boldsymbol{\varepsilon}^b$ and $\boldsymbol{\varepsilon}^o$ are non-trivial: The covariance matrices B and R are positive definite.

* Errors $\boldsymbol{\varepsilon}^b$ and $\boldsymbol{\varepsilon}^o$ are unbiased: $\mathbb{E}(\boldsymbol{\varepsilon}^b) = \mathbb{E}(\boldsymbol{\varepsilon}^o) = \overline{\mathbf{x}^b - \mathbf{x}^t} = \overline{\mathbf{y} - \mathcal{H}(\mathbf{x}^t)} = \mathbf{0}$

* Errors $\boldsymbol{\varepsilon}^b$ and $\boldsymbol{\varepsilon}^o$ are uncorrelated: $(\mathbf{x}^b - \mathbf{x}^t)(\mathbf{y} - \mathcal{H}(\mathbf{x}^t))^T = \mathbf{0}$
* The analysis $\mathbf{x}^a$ is linear: The analysis, defined by corrections to the background state, depends linearly on the differences between the background and the observations.

* The analysis $\mathbf{x}^a$ is optimal: The analyzed state is as close as possible to the true state in accordance with the root mean square, e.g., an estimate of the minimum variance.

### 2.1.4 Assimilation techniques

All data assimilation techniques rely on the classical Bayesian inversion framework with Gaussian assumptions (in other words, PDF of error is fully described by their mean and covariance). Having a priori information about the state of the atmosphere, the search for analysis is the search for the most probable state of the atmosphere knowing the observations. Therefore, the Gaussian Bayesian formulation consists in computing the following PDF:

$$p(\mathbf{x}|\mathbf{y}^o, \mathbf{x}^b) \sim N(\mathbf{x}^a, \mathbf{A})$$

where the solution of an inverse problem $\mathbf{x}^a$ is a random vector with a Gaussian or normal distribution assigned with $\mathcal{N}$.

#### 2.1.4.1 Statistical data assimilation

The aim is to produce an estimate $\mathbf{x}^a$ of the true input parameters $\mathbf{x}^t$ of the system. To this end, statistical interpolation (also called OI, which is a simpler version of KF) used to combine two independent sources of information. This method allows writing them as a linear combination between the background ($\mathbf{x}^b$) vector and observations ($\mathbf{y}$) vector as follows:

$$\mathbf{x}^a = \mathbf{x}^b + \mathbf{K}(\mathbf{y} - \mathcal{H}(\mathbf{x}^b)) \quad (\text{Eq.12})$$

in which an innovation $\mathbf{y} - \mathcal{H}(\mathbf{x}^b)$ is weighted by the linear operator $\mathbf{K} = \mathbf{BH}^T(\mathbf{HBH}^T + \mathbf{R})^{-1}$ called gain or weight matrix. Noting that, $\mathbf{H}$ is the tangent linear version of the operator $\mathcal{H}$. Under the assumptions above, the best estimate of $\mathbf{x}^a$ in Eq.12 is called **BLUE-Best Linear Unbiased Estimator** via least squares analysis. In other words, BLUE predicts an unbiased, optimal weighting for a linear correction of the state of the model from the observations. Thus, the problem seeks to find the best unbiased estimation of $\mathbf{x}$ that minimizes the variance of the analysis errors (Eq.9) (optimality criterion). The covariance matrix of the analysis error $\mathbf{A}$ is
calculated using the following error equivalent in Eq.13 (This equation is obtained using Eq.5, Eq.7, and Eq.9 in Eq.12):

$$ \boldsymbol{\epsilon}^a = \boldsymbol{\epsilon}^b + \mathbf{K}(\boldsymbol{\epsilon}^o - \mathbf{H}(\boldsymbol{\epsilon}^b)) \quad (\text{Eq.13}) $$

Therefore, for all **K** (not necessarily to be an optimal),

$$ \begin{aligned} \mathbf{A} = \mathbb{E}(\boldsymbol{\epsilon}^a \boldsymbol{\epsilon}^{aT}) &= \mathbb{E}\left[\left(\boldsymbol{\epsilon}^b + \mathbf{K}(\boldsymbol{\epsilon}^o - \mathbf{H}(\boldsymbol{\epsilon}^b))\right)\left(\boldsymbol{\epsilon}^b + \mathbf{K}(\boldsymbol{\epsilon}^o - \mathbf{H}(\boldsymbol{\epsilon}^b))\right)^T\right] \\ &= \mathbb{E}\left[\left((\mathbf{I} - \mathbf{K}\mathbf{H})\boldsymbol{\epsilon}^b + \mathbf{K}\boldsymbol{\epsilon}^o\right)\left((\mathbf{I} - \mathbf{K}\mathbf{H})\boldsymbol{\epsilon}^b + \mathbf{K}\boldsymbol{\epsilon}^o\right)^T\right] \\ &= \mathbb{E}\left[\left((\mathbf{I} - \mathbf{K}\mathbf{H})\boldsymbol{\epsilon}^b(\boldsymbol{\epsilon}^b)^T(\mathbf{I} - \mathbf{K}\mathbf{H})^T\right)\right] + \mathbb{E}[\mathbf{K}\boldsymbol{\epsilon}^o(\boldsymbol{\epsilon}^o)^T\mathbf{K}^T] \\ &= (\mathbf{I} - \mathbf{K}\mathbf{H})\mathbf{B}(\mathbf{I} - \mathbf{K}\mathbf{H})^T + \mathbf{K}\mathbf{R}\mathbf{K}^T \end{aligned} \quad (\text{Eq.14}) $$

where **I** is the identity matrix. If **K** is the optimum gain (so it minimizes the scalar error Eq.10 or the matrix entries along the principal diagonal of **A**), then the analysis error covariance matrix becomes **A** = (I - **K**H)B. Full details of the computation of **K** and **A** can be found in Burgers, Jan van Leeuwen, and Evensen 1997.

### 2.1.4.2 Variational data assimilation

As seen in the previous section, BLUE requires the direct computation of an optimal gain matrix **K**. However, in the complex systems, this matrix is very large and therefore expensive in computing time. Nevertheless, BLUE can be obtained without computation of the gain **K** using variational calculus, e.g., through the minimization of a functional. Therefore, variational methods, based on the optimal control theory, are promising in data assimilation. They consist of minimizing the least squares functional. The proof of their equivalence to the optimal interpolation problems can be found in (Lorenc 1986; Talagrand 1997).

#### 2.1.4.2.1 Formulation of variational data assimilation

Let a state vector $\boldsymbol{x}(t) \in \mathcal{X}$ and time $t \in \mathcal{T}$, which contains any physical variable(s) of the state of the system (e.g., chemical compositions for air quality modeling). In reality, the evolution of the state is continuous and described by a set of nonlinear differential equations in a set of $\Omega$ as follows:
$$ \frac{dx}{dt} = \mathcal{M}(x) \quad \text{in } \Omega \times [0, T] \qquad (\text{Eq.15}) \\ x(t=0) = x_0 $$

where the true initial condition is not known. Supposing that $y(t) \in \mathcal{O}$ are observations of the true state and $\mathcal{H}$ is observation operator, in the continuous context, **the objective J function** (also called **cost function**) defines the distance of the state of the atmosphere with respect to the background and the observations as follows:

$$ J(x) = \frac{1}{2} \int_{0}^{T} \|y(t) - \mathcal{H}(x(x_0, t))\|_{o}^{2} dt + \frac{1}{2} \|x - x^{b}\|_{x}^{2} \qquad (\text{Eq.16}) $$

While in the first term, integral calculates the distance between observations and modelled state under the observation operator in the infinite and continuous context, the second term comes from a regularization, as inverse problems are known to be ill posed and quantify the misfit to the background state. Adding the second term makes the problem well posed. The problem would then be to search for $x^a$ that minimizes the cost function $J$:

$$ x^a = \operatorname*{argmin}_{x} J(x) \qquad (\text{Eq.17}) $$

This minimization requires numerical optimization techniques. Quasi-Newtonian descending algorithms are used to solve this minimization problem iteratively with high dimensionality, starting with an initial state vector $x_0$. For the existence of a (local) minimum of Eq.17, the gradient of the cost function is necessarily $\nabla J(x^a) = 0$. However, the calculation of the derivatives of the cost function at each step with high dimensions can be computationally expensive. For this reason, these gradients are obtained by running an adjoint model simultaneously in parallel to the inverse model. This is why variational data assimilation is also known as ‘adjoint’ method.

### 2.1.4.2.2 4D-Var and its pioneering 3D-Var

This subsection introduces 4D-Var data assimilation formulation and its adjoint, as well as its pioneering 3D-Var method in a finite-dimensional and discrete context. In the finite-dimensional case, the two norms under the integral (Eq.16) are represented by the error covariance matrices $R$ and $B$, respectively. The observation operator $\mathcal{H}$ and the computational model $\mathcal{M}$ are assumed to be linear (thus the cost function is quadratic) according to the hypotheses made earlier and defined as their Jacobian matrix $H$ and $M$, respectively. Their
adjoint $\mathcal{H}^*, \mathcal{M}^*$ can be identified by the transpose for the Jacobian $\mathbf{H}^\mathrm{T}$ and $\mathbf{M}^\mathrm{T}$. Under these assumptions, the direct model $\mathbf{M}$ is defined by,

$$ \begin{aligned} \mathbf{x}_i &= \mathbf{M}_i(\mathbf{x}_{i-1}), \quad i = 1, n \\ \mathbf{x}_0 &= \mathbf{x}^b \end{aligned} \qquad (\text{Eq.18}) $$

where the state trajectory is $\mathbf{x} = (\mathbf{x}_0, \mathbf{x}_1, \dots, \mathbf{x}_{n-1}, \mathbf{x}_n)^\mathrm{T}$, $n$ is the number of time steps per assimilation window. Then, for a fixed time, the Eq.16 would define 3D-Var cost function as follows,

$$ J(\mathbf{x}) = \frac{1}{2}(\mathbf{x} - \mathbf{x}^b)^\mathrm{T}\mathbf{B}^{-1}(\mathbf{x} - \mathbf{x}^b) + \frac{1}{2}(\mathbf{y} - \mathbf{H}(\mathbf{x}))^\mathrm{T}\mathbf{R}^{-1}(\mathbf{y} - \mathbf{H}(\mathbf{x})) \qquad (\text{Eq.19}) $$

And the gradient of the $J$ function is given by:

$$ \nabla J(\mathbf{x}) = \mathbf{B}^{-1}(\mathbf{x} - \mathbf{x}^b) - \mathbf{H}^\mathrm{T}\mathbf{R}^{-1}(\mathbf{y} - \mathbf{H}(\mathbf{x})) \qquad (\text{Eq.20}) $$

where $\mathbf{B}$ and $\mathbf{R}$ are stated as the background and observation error covariance matrices respectively. Calculations using these matrices are not computationally storable and manipulable due to the high dimensionality of the problem. Therefore, the main idea of the variational method is to minimize the cost function $J$ in order not to deal with these difficulties. This minimization is achieved by an adjoint model involving the calculation of the gradient of the cost function and a suitable descent algorithm in the gradient direction. This minimization can be stopped by a stopping criterion that sets a decrease in the norm of the gradient of $J$ ($\|\nabla J(\mathbf{x})\|$) small enough or artificially limits the maximum number of iterations. When there is no time dependence, like in the case of 3D, the adjoint model is not necessary. Because the adjoint simply includes the transpose of $\mathbf{H}$ matrix (Eq.20). However, the 3D-Var method is a particular case of the 4D-Var method. While 3D-Var does not take into account temporal dimension in minimization process, 4D-Var method carries the minimization over a time interval. The process involves the observations, which are obtained at different times, and considers all the observations at once in the assimilation window (e.g., in the Figure 14 the assimilation window is from $\mathbf{t}_0$ to $\mathbf{t}_n$ (Lahoz and Schneider 2014)). Therefore, the adjoint model is required for the transpose of both the model and operator dynamics.
<figure><img src="image_18.png" /><figcaption>Figure 14 Illustration of 3D-Var & 4D-Var (Taken from Lahoz and Schneider, 2014)</figcaption></figure>

As a result, 4D-Var data assimilation scheme is defined by the minimization of a cost function $J(\mathbf{x}) = J^b(\mathbf{x}) + J^o(\mathbf{x})$ from time $i = 0$ to $i = n$ as follows:

$$J(\mathbf{x}) = \frac{1}{2}(\mathbf{x} - \mathbf{x}^b)^{\mathrm{T}}\mathbf{B}^{-1}(\mathbf{x} - \mathbf{x}^b) + \frac{1}{2}\sum_{i=0}^{n}(\mathbf{y}_i - \mathbf{H}_i(\mathbf{x}_i))^{\mathrm{T}}\mathbf{R}_i^{-1}(\mathbf{y}_i - \mathbf{H}_i(\mathbf{x}_i)) \quad (\text{Eq.21})$$

where the state $\mathbf{x}$ at time $i$ is obtained by an iterative composition of the model $\mathbf{M}$, e.g., $\mathbf{x}_i = \mathbf{M}_i \mathbf{M}_{i-1} ... \mathbf{M}_2 \mathbf{M}_1 \mathbf{x}$. Then the cost function $J$ is as follows,

$$J(\mathbf{x}) = \frac{1}{2}(\mathbf{x} - \mathbf{x}^b)^{\mathrm{T}}\mathbf{B}^{-1}(\mathbf{x} - \mathbf{x}^b) + \frac{1}{2}\sum_{i=0}^{n}\mathbf{d}_i^{\mathrm{T}}\mathbf{R}_i^{-1}\mathbf{d}_i \quad (\text{Eq.22})$$

where the innovation vector $\mathbf{d}_i = \mathbf{y}_i^0 - \mathbf{H}_i \mathbf{M}_i \mathbf{M}_{i-1} ... \mathbf{M}_2 \mathbf{M}_1 \mathbf{M}_i$. Then the gradient of the $\nabla J(\mathbf{x})$,

$$\nabla J(\mathbf{x}) = \mathbf{B}^{-1}(\mathbf{x} - \mathbf{x}^b) - \sum_{i=0}^{n} \mathbf{M}_i^{\mathrm{T}} \mathbf{H}_i^{\mathrm{T}} \mathbf{R}_i^{-1} \mathbf{d}_i \quad (\text{Eq.23})$$

By this way, the 4D-Var algorithm calculates first $J$ followed by $\nabla J$ one direct model integration and one adjoint model integration. One can find the evolution of Eq.22 and Eq. 23 from Eq.21 in Annex B: 4D-Var Cost and the Gradient Evolution.
2.1.4.2.3 4D-Var algorithm

* Initialization: $i = 0, \mathbf{x}^b = \mathbf{x}_0$

* While $\|\nabla J\| > \delta$ or $i \le i_{max} = n$:

* Compute $J$ using the direct model $\mathbf{M}$ and the observation operator $\mathbf{H}$ (forward mode)

* Compute $\nabla J$ using the backwards integration of the adjoint model $\mathbf{M}^T$ and the adjoint of the observation operator $\mathbf{H}^T$ (inverse mode)

* Gradient descent minimization and update of $\mathbf{x}$

* $i = i + 1$

At the initialization step, trajectory of the model is determined by the state of the starting point $\mathbf{x}_0$ and set it to the background $\mathbf{x}^b$. This is a first guess. At Step 1, Model runs forward in time to the end of the assimilation window to obtain the first-guess trajectory using Eq.18 and Eq.22. Using the first guess trajectory from Step 1, Step 2 integrates the adjoint model backwards in time from $t_i$ to $t_0$ using Eq.23. In Step 3, the optimal state of the cost function is calculated using the outputs from Step 1 and Step 2 via Gradient Descent Minimization (Figure 15, Bouttier and Courtier 1999) and set to the state variable. This procedure iteratively continues via Step 4 until it reaches the maximum number of steps ($i_{max}$) or the gradient of the cost function is sufficiently smaller than a delta ($\delta$) value determined by the user.

<figure><img src="image_19.png" /><figcaption>Figure 15 Visualization of gradient descent minimization for the state $\mathbf{x}$ in 2D. $J$ is a quadratic convex cost function. Its minimum is obtained $\mathbf{x}^a$. The image taken from, Bouttier and Courtier, ECMWF,1999</figcaption></figure>
### 2.1.4.2.4 Limitations of 4D-Var

Although the 4D-Var technique is powerful and robust, it is not easy to implement because it is highly dependent on the problem (linearity, dimensionality) and makes the calculations computationally expensive because the algorithm relies on the availability of derivative information (e.g., coding of the adjoint operators). The adjoint model is a mathematical trick used to calculate sensitivities between model parameters outside the forward model. However, developing an adjoint model is computationally costly and consists of more complex equations than forward model. At the same time, each adjoint model is specific to the forward model, so any changes in the forward model also require recoding in the adjoint model.

Also, the development of background **B** and the observation error covariance matrix **R** are challenging and require attention. Because they give perception about how much we trust the given information (e.g., how much we trust the observations), any misled specification of these matrices may strongly influence the results of the inversion, especially at high resolutions.

Noting that the model operator **M** and observation operator **H** are assumed to be linear here. It is known that, under these linearity assumptions, 4D-Var and Kalman filter methods provide the same results at the end of the assimilation period (Lorenc 1986). However, in the real world problems, they are non-linear, and non-linearity of the model **M** and observation operator **H** making the objective function non-convex, and its gradient complex. In this case, the BLUE formalism is not valid anymore. Thanks to the tangent linear hypothesis (TLH) and adjoint models (Lewis and Derber 1985), non-linear operators are replaced by their tangent linear operators (Annex C: Non-linear operators), resulting in the problem being an approximation of 4D-Var. Courtier, Thépaut, and Hollingsworth 1994 proposed the use of the TL model ($\delta$) using the increment vector in the classical formulation of 4D-Var (called incremental 4D-Var) and solving the problem in the neighborhood of the background state. In this case, the 4D-Var algorithm includes systematic and random errors raised by the linearization.

In the framework of this thesis, the variational approach of the CIF inverse modeling system is used. The next section is therefore devoted to an introduction about CIF and the development of CIF for our inversion problem.
## 2.2 CIF inverse modeling system

This thesis adapts the Community Inversion Framework (CIF), a state-of-the-art inverse modeling system. This system is developed to unify atmospheric inversion problems and bring together researchers interested in atmospheric inversions (Berchet et al. 2021). It is a user-friendly, easily accessible, transparent, and open-source Python-based tool on http://community-inversion.eu/ (Last access: 08/04/2023). It is used to estimate the atmospheric states, more specifically, fluxes of various greenhouse gases and reactive species by having the flexibility to run CTMs (Lagrangian or Eulerian) at different spatial (regional or global) and temporal scales using different types of observations (e.g., in-situ, satellite) in different data assimilation techniques (e.g., analytical, statistical, variational).

The next section explains how the CIF is structured in this thesis.

### 2.2.1 Development of CIF in this study

In this thesis, to inverse NOₓ emissions using satellite observations, we used variational approach of the CIF (more specifically 4D-Var). Although the variational approach offers a numerical approximation to the solution of the inverse problem, the formulation of the problem does not require complex matrix products. Therefore, the high dimensionality of the matrices is not an obstacle during the calculation process (Berchet et al. 2021). Furthermore, inversion within the assimilation window makes it possible to consider non-linear links between emissions and concentrations along trajectories on time scales equal to the chemical lifetime of the studied species. For this reason, due to the high diurnal characteristics of NO₂ concentrations and their non-linearity with NOₓ emissions, the 4D-Var data assimilation technique with 24h assimilation windows (that is daily temporal resolution) is used. Next, it is described the main elements of the 4D-Var used within CIF.

#### 2.2.1.1 The control vector ($\mathbf{x}$)

The control vector contains all the parameters related to the species to be optimized in the inversion problem (e.g., fluxes and/or initial and boundary concentration conditions). In this thesis as we aim to estimate NOₓ emissions, the control vector contains corrections coefficients (or multiplicative factors) to be applied to the relevant NOₓ emissions information (including anthropogenic and biogenic NOₓ emissions) at a temporal resolution of 1-day. The dimensions
of the control vector is as follows; hours in a day × level of altitude × number of latitudinal grid cells × number of longitudinal grid cellsThe daily corrections obtained from the inversions are applied to the background NOₓ emissions (xᵇ) as multiplicative factors. The relationship between the multiplicative factors of the control vector and background NOₓ emissions (prior) is then as follows:

Estimated NOₓ emissions (posteriori) = x. Prior NOₓ emissions

### 2.2.1.2 The background vector (xᵇ)

The background vector is an initial state (Eq.18) of the inversion problem. It contains a priori information to be considered by the inversion. It is linked to the background NOₓ daily emission fluxes (or prior NOₓ emissions). The background vector is defined with multiplication factors of 1 as these are reference emissions from bottom-up emission inventories, which will be described in detail in Section 2.3. Similar to the control vector, the background vector has the dimensions as follows; hours in a day × level of altitude × number of latitudinal grid cells × number of longitudinal grid cells.

### 2.2.1.3 The background error covariance matrix (B)

It is important to note that, the uncertainties in the prior NOₓ emissions are assumed to be unbiased, have a Gaussian distribution and are characterized by their B covariance matrix. Its diagonal terms are the weights to the components of the vector differences between the control vector and prior NOₓ emissions. On the other hand, non-diagonal terms of the covariance matrix give uncertainties about the spatial correlations of the NOₓ emission sources. In our configuration, they are calculated using exponentially decaying functions with the correlation length scale derived from the model resolution. As it is difficult to construct B matrix standard estimates of prior NOₓ emissions, uncertainties from the literature are used as a starting point in this thesis instead of complicated ones.

### 2.2.1.4 The observation operator (H(x))

The observation operator includes a computation of CTM, which is used to link NOₓ surface fluxes to NO₂ atmospheric concentration measurements and projection of simulated NO₂ concentrations to the NO₂ tropospheric column measurements by the satellite instrument. By
this way, the background emissions become comparable with the satellite instrument
measurements. The CTM model used in this thesis is explained in Section 2.3.

2.2.1.5 The observation vector (y)

The observation vector brings independent information to the inversion system using atmospheric observations. These observations can come from ground-based measurements, such as trace gas mixing ratios from fixed or mobile platforms, and remote sensing observations, such as satellite observations. In this thesis, satellite observations of the NO₂ tropospheric column are used because of their high spatial coverage in the area of interest. Section 2.4 explains the types of satellite instruments of which observations are used in this thesis.

2.2.1.6 The observation error covariance matrix (R)

As with the **B** matrix, the observation errors are assumed to be unbiased and Gaussian distributed. In this thesis, errors from the model are neglected and all the diagonal terms of the **R** matrix are taken as the measurement errors reported in the satellite data sets. Additionally, spatial error correlations between NO₂ satellite observations are neglected. Therefore, all non-diagonal terms of the **R** matrix are taken as zero.

2.2.1.7 The minimizer M1QN3

The minimization subroutine M1QN3, which uses the cost function and its adjoint as input, is
used to minimize the inversion problem. The minimizer written in Fortran was described by
(Gilbert and Lemaréchal 1989; Gilbert, Lemaréchal, and Simulator 2009). It is distributed under
the GNU General Public License. It is used to solve large-scale unconstrained minimization
problems using a quasi-Newton algorithm. Quasi-Newton methods are designed to use on linear
problems. It iteratively moves the cost function along the direction of its steepest descents,
which is, calculated an approximation of the Hessian matrices and by the adjoint. However,
non-linearity of the chemistry of NOₓ reactive species makes our inversion problem non-linear
and the cost function non-convex. For this reason, the quasi-Newton algorithm can reach the
local minima without reaching a global minimum. To solve this issue, this method needs an
initial regularization of the control vector that we want to minimize. In other words, instead of
minimizing $\boldsymbol{x}$ we minimize $\boldsymbol{\chi} = \mathbf{B}^{-\frac{1}{2}}(\boldsymbol{x} - \boldsymbol{x}^b)$. A simpler formulation of the adjoint of the cost function $J$ to be applied regulation is as follows:

$$ \nabla J(\boldsymbol{x}) = \mathbf{B}^{-1}(\boldsymbol{x} - \boldsymbol{x}^b) - \mathcal{H}^*(\mathbf{R}^{-1}(\mathcal{H}(\boldsymbol{x}) - \boldsymbol{y})) $$

The formulation of the adjoint of the cost function $J$ becomes:

$$ \nabla J(\boldsymbol{\chi}) = \boldsymbol{\chi} - \mathbf{B}^{-\frac{1}{2}} \mathcal{H}^* \left( \mathbf{R}^{-1} \left( \mathcal{H} \left( \mathbf{B}^{-\frac{1}{2}} \boldsymbol{\chi} + \boldsymbol{x}^b \right) - \boldsymbol{y} \right) \right) $$

The minimizer algorithm is required to meet the convergence criterion before reaching the stopping criterion. In cases where this criterion is not met, the algorithm is terminated by the stopping criterion and the optimization is considered non-convergent. In other words, this means that a suitable state vector that minimizes the cost function has not been found. The convergence criterion we adopt in this thesis is as follows; consider the cost function to have converged if, after minimization, the norm of its gradient has decreased by 95 % compared to the initial one. If such a vector that minimizes the cost function cannot be found, stop running the algorithm with at most 20 iterations.

### 2.2.1.8 Analysis ($\boldsymbol{x}^a$)

The optimal estimate of the NOx emissions is found by iteratively minimizing the cost function $J$ and is so-called the analysis, or posterior NOx emissions. Figure 16 presents the iterative structure of the CIF used for our inversion problem.

<figure><img src="image_20.png" /><figcaption>Figure 16 Scheme of the CIF</figcaption></figure>
## 2.3 Atmospheric Modeling

Atmospheric modeling consists of a set of mathematical equations representing the physical and chemical interactions of the atmosphere. It can simulate the evolution of atmospheric concentrations, taking into account all atmospheric physical and chemical phenomena within and at the boundaries of the study area (domain), but also taking into account emissions and meteorological parameters used as forcing data (e.g., precipitation, wind, temperature). The simulation of atmospheric concentrations is performed using CTMs based on two main approaches: (i) The Lagrangian approach and (ii) The Eulerian approach.

The Lagrangian approach considers a pollutant on the flow and simulates the trajectory of the air parcel with respect to its direction due to the wind. Therefore, the simulation is performed at different locations and different times. The Eulerian approach, on the other hand, considers the domain over which the pollutant moves as fixed and discretized (e.g., by boxes or cells), and simulations are carried individually within these discrete cells (so-called 'grid cells'). Chemical composition calculations are performed for each grid cell. The size of the grid cells determines the resolution of the CTM. A 3D CTM model is a model that performs simulations at the desired vertical and horizontal resolutions, taking into account the interactions between each grid cell (both vertically and horizontally), and determines the concentrations over the entire domain. It is necessary to take into account the transport between the grid cells and the physical/chemical phenomena in the grid cells in order to calculate the concentrations in each grid cell. For this reason, simulations at each grid cell are performed based on the following mass balance equation to determine the concentration c of the species s at time t.

$$ \frac{dc_s}{dt} = \left(\frac{dc_s}{dt}\right)_{\text{chemistry}} + \left(\frac{dc_s}{dt}\right)_{\text{advection}} + \left(\frac{dc_s}{dt}\right)_{\text{emission}} + \left(\frac{dc_s}{dt}\right)_{\text{dilution}} + \left(\frac{dc_s}{dt}\right)_{\text{deposition}} \quad (\text{Eq.24}) $$

Let us represent the chemistry balance term $(\frac{dc_s}{dt})_{\text{chemistry}}$ by $R_s$ which is the net chemical production rate of the species s at time t. The advection term can be expressed as follows,

$$ \left(\frac{dc_s}{dt}\right)_{\text{advection}} = \frac{u}{\Delta x} (c_s^0 - c_s) $$
where u is the wind speed, Δx is the length of the grid cell, c₀ is the background concentration of the species s. So the proportion of the wind speed to the length of the grid cell determines the residence time of the air parcel over the domain, which reduces some of the background concentration of the species s. Sometimes the source of the pollutant is inside of the grid cell so the emission term can be expressed as follows,

$$ \left(\frac{dc_s}{dt}\right)_{\text{emission}} = \frac{q_s}{h} $$

which represents the emission rate per unit area (kgm⁻²h⁻¹), where q_s is the quantity of the emission of the species s and h is the height of the grid cell. The dilution balance term (dc_s/dt) can be neglected when the grid cells have constant height because there will not be any transfer entraining or detraining of the concentrations into or out of the grid cell. The deposition term can be expressed using removal rate as follows,

$$ \left(\frac{dc_s}{dt}\right)_{\text{deposition}} = -\frac{v_{d_s} c_s}{h} $$

where v_d_s is dry or wet deposition velocity of the species s.

Overall, the balance of the species s can therefore be expressed with the simplified mass balance equation as follows:

$$ \frac{dc_s}{dt} = R_s + \frac{u}{\Delta x}(c_s^0 - c_s) + \frac{q_s}{h} - \frac{v_{d_s} c_s}{h} \quad (\text{Eq.25}) $$

This equation includes all source and sink terms of chemical constituents in the CTM; therefore, solving this equation gives the concentration of a given species s for a given time t (so-called output). Figure 17 visualizes the mass balance conception of the concentration of a species s at one grid cell and construction of the 3D CTM.
<figure><img src="image_21.png" /><figcaption>Figure 17 Illustration of simplified mass balance processes represented in each grid cell of a 3D CTM</figcaption></figure>

For each grid cell, solving the mass balance equation numerically consists of the calculation of the sources (e.g., emissions) and sinks (e.g., deposition) while considering the interaction between neighboring grid cells. To do this, some information (so-called input) must be provided to the CTMs before the simulations launch. For example, anthropogenic and biogenic emissions from emission inventories or other CTMs need to be provided as input. In order to determine atmospheric conditions and provide information about the atmospheric situation, meteorological conditions (e.g., temperature, wind, cloud, humidity) need to be defined through offline or online CTM modeling environments and meteorology models. While in an offline setting, a CTM is used to run after the meteorology model determines meteorological conditions in advance, in an online setting, there is back and forth relationship between the CTM and the meteorology model. Moreover, a chemical scheme has to be provided for chemical reaction types and rates of species including aerosol chemistry. Assuming that a CTM simulation covers only part of the globe, meaning that simulations run over a smaller domain, in this case, boundary conditions of the domain need to be presented to take into account transport from neighborhoods. Generally, these conditions are determined running another simulation over the larger domain (usually with lower resolutions compared to the smaller domain). These kinds of simulations are called nested-simulations. Moreover, initial conditions are presented to provide information about the existing atmospheric chemical conditions of the domain, which are fed to the model in the same manner as boundary conditions. Noting that boundary and initial conditions are necessary only for the first time-step. Finally, land use information needs to be presented.
CTMs have advantages that turn them into complementary tools for the measures observed. They provide continuous and potentially very complete information relative to physical sensing instruments. They make it possible to quantify interactions between different processes. On the other hand, like other man made tools, they face some limitations. They can only process the state of the information we have already defined. For example, they cannot take into account information from an as-yet-unknown chemical reaction. As such, they are only a representation of reality, and need to be validated to ensure that the representation is faithful enough for the use of simulations. Furthermore, the amount of information that needs to be used as input is detailed and these inputs also contain some limitations and uncertainties. For example, uncertainties in emissions used in CTMs. However, CTMs are still the only and necessary tool to reconstruct inaccessible information, to predict concentrations anywhere at any time in the past or the future.

In the framework of this thesis, in order to perform an inversion of NOₓ emissions on a regional scale, it is necessary to run a CTM to link NOₓ emission fluxes and their ambient concentrations. To do that, the regional CTM CHIMERE is used, since it is one of the CTMs whose adjoint has been coded (Laurent Menut 2003) and coupled with the CIF inversion system.

### 2.3.1 CHIMERE Regional CTM

CHIMERE regional chemistry and transport model has been co-developed since 1997 by Laboratoire de Météorologie Dynamique (LMD), Institut national de l'environnement industriel et des risques (INERIS) and Laboratoire Interuniversitaire des Systèmes Atmosphériques (LISA). It is a 3D offline (and online for the recent versions) Eulerian model designed to run long-term simulations for emission control scenarios and forecasting air pollution events within the atmospheric boundary layer. It uses external forcing (e.g., meteorological fields, pollutant emissions, boundary conditions) to simulate hourly concentrations of gaseous and aerosol species on a regular Cartesian grid from the surface up to the level of interest with resolutions from 1-2 km up to 100-200 km. CHIMERE simulated concentrations were compared and evaluated in different resolutions with various types of measurements from surface stations (Menut et al. 2000; Blanco-Ward et al. 2021; Menut et al. 2020), airborne or various satellite instruments for tropospheric columns (Huijnen et al. 2010; Konovalov et al. 2009).

This thesis uses the CHIMERE v2013b version, which is integrated into the CIF without activating the aerosol module (Menut et al. 2013). This is because the adjoint model of
CHIMERE (Laurent Menut 2003), which helps to calculate the sensitivity of simulated atmospheric concentrations to corrections in fluxes, was developed without including the aerosol module due to its high computational cost. Therefore, the CIF-CHIMERE only considers full gaseous chemistry, not aerosols. Considering aerosols in the model can lead to different simulated NO₂ concentrations as they significantly affect the photochemical process. As an example, we checked the impact of aerosols in the simulation only in the Paris region (within 15 km), where the consideration of aerosol leads to a reduction of the simulated NO₂ concentrations by 5 to 15%.

CHIMERE uses hybrid sigma-pressure ($\sigma - p$) coordinates which allows user to define any number of vertical layers. The formula to calculate the pressure (in hPa) at the top of each layer $k$, $p_k$ is as follows: $p_k = a_k p_0 + b_k p_{surf}$ where $p_{surf}$ is the surface pressure, $p_0$ is the reference state pressure, $a_k$ and $b_k$ are the hybrid coefficients determined by either the user or the model. More detailed information about CHIMERE v2013b and its adjoint can be found in (Menut et al. 2013; Laurent Menut 2003; Fortems-Cheiney et al. 2021).

### 2.3.1.1 Setup of the CHIMERE Simulations

In the framework of this thesis, the study domains include two principal urbanized regions of the Northern Hemisphere: Eastern China and Europe (studied in the Chapters 3 and 4 respectively); hence, CHIMERE's configurations vary depending on these regions. Therefore, the setup of the CHIMERE simulations and their inputs are explained separately according to these domains.
### 2.3.1.1.1 Domain : Eastern China

CHIMERE setup for the NOₓ emission inversion to simulate NO₂ concentrations over Eastern Asia covers the regular grid with horizontal resolution of 0.5° × 0.5° and 17 vertical levels from the surface up to 200 hectopascal (hPa) with 8 layers within the first three kilometers. It includes 65 longitude (lon) × 61 latitude (lat) grid cells that correspond to 18°N - 50°N (degrees north) and 102°E - 132°E (degrees east) and covers eastern and southern China, the South China Sea and Korean Peninsula domain as in the Figure 18.

<figure><img src="image_22.png" /><figcaption>Figure 18 Eastern China domain visualization</figcaption></figure>

To be able to perform simulations over Eastern China domain, the necessary ingredients used as inputs for CHIMERE CTM is summarized in the Figure 19.
<figure><img src="image_23.png" /><figcaption>Figure 19 Main inputs and their sources for CTM CHIMERE simulations over Eastern China</figcaption></figure>

Boundary and initial conditions of the domain are given by a climatology from the global Laboratoire de Météorologie Dynamique Zoom (LMDz) model that is coupled with the model Interaction with Chemistry and Aerosols (INCA) (Hauglustaine et al. 2004). Meteorological conditions are coming from the global meteorological forecasts of European Centre for Medium-Range Weather Forecasts (ECMWF) and hourly interpolated to fit the hourly simulations (Owens and Hewson 2018). MELCHIOR-2 gaseous chemical scheme which is the reduced MELCHIOR chemical mechanism (Lattuati 1997; Derognat 2003) comprising 44 species and 120 reactions is used to evaluate tropospheric chemical evolution of gaseous species in CHIMERE (Menut et al. 2013; Derognat 2003). Land use information is coming from the database Global Land Cover Facility (GLCF) of University of Maryland (CHIMERE User's Manual). Biogenic emissions are derived from the Model of Emissions of Gases and Aerosols from Nature (MEGAN) v2.1 model (Guenther et al. 2006). And finally, anthropogenic emissions are those from Hemispheric Transport of Air Pollution (HTAP) v2.2 bottom-up inventory (Janssens-Maenhout et al. 2015). It is important to note that, HTAP v2.2 anthropogenic emissions are taken as a reference (or prior) bottom-up emission inventory in the inversion to estimate NOₓ emissions for Eastern China; therefore, the following part is dedicated to HTAP v2.2 prior bottom-up inventory. It is worth noting that, in this study, NOₓ emissions are expressed as nitrogen oxides (NO₂ and NO) as budgets with the unit of kilotons (hereinafter kt NOₓ or Tg NOₓ referred to as kt or Tg).
### 2.3.1.1.1.1 Anthropogenic Bottom-up Emission Inventory: HTAP

HTAP are used as the prior inventory in our inversions (Henze, Hakami, and Seinfeld 2007). For the region of China, HTAP v2.2 is constructed using regional emission inventories in collaboration with the MICS Asia Scientific Community (MICS-Asia) and REAS experts with a horizontal resolution of 0.25° × 0.25° converted to 0.1° × 0.1°. Emissions are estimated from seven main sectors of anthropogenic activities: air, ships, energy, industry, transport, residential, and agriculture.

In this study, The EDGAR-HTAP V2.2 emissions inventory (Janssens-Maenhout et al. 2015) based on the year 2010 was used to generate anthropogenic emissions of gaseous species. The emissions are initially grouped by sector of activity (SNAP, Standardized Nomenclature for Air Pollutants); they are then broken down to correspond to the actual pollutants, and then aggregated to match the species represented in CHIMERE (a construction of emissions of the type (Middleton, Stockwell, and Carter 1990). Emissions are spatially distributed and daily hourly profiles emissions are allocated (in Coordinated Universal Time (UTC)) for each day of the week.

### 2.3.1.1.1.2 Anthropogenic Bottom-up Emission Inventory: MEIC

MEIC (Multi-resolution Emission Inventory for China) is a bottom-up emission inventory developed and operated at Tsinghua University since 1990 (http://meicmodel.org.cn/, Last access: 08/04/2023). It is designed to track anthropogenic emissions over continental China at different spatial resolutions. NOₓ emissions are mainly presented in five sectors: energy, transportation, industry, agriculture, and residential. For more details about MEIC inventory construction and sectors, see (Li, Liu, et al. 2017; Zheng et al. 2018).

In this thesis work, MEIC version 1.3 from 2015 (Zheng et al. 2018; 2021) is used as an independent inventory to evaluate the estimated NOₓ emissions for the same year. Note that this version of MEIC includes the impacts of the 12th Five-Year Plan mentioned in Section 1.6.
### 2.3.1.1.2 Domain : Europe

CHIMERE setup for the NOₓ emission inversion to simulate European NO₂ concentrations covers the regular grid with horizontal resolution of 0.5° × 0.5° and 20 vertical levels from the surface up to 200 hPa with 8 layers within the first three kilometers. It includes 101 lon × 85 lat grid cells that correspond to -15.25°N – 34.75°N and 31.75°E – 73.75°E and covers Europe (Figure 20).

<figure><img src="image_24.png" /><figcaption>Figure 20 European domain visualization</figcaption></figure>

To be able to perform the simulations, the inputs variables of CHIMERE are as follows (Figure 21): similar to the Eastern China domain (Figure 19), meteorological forecasts are obtained from ECMWF, the chemical scheme MELCHIOR-2 is used. Initial, lateral and boundary conditions are given by a climatology from LMZD-INCA global model. The biogenic emissions are those from MEGAN model as well. Moreover, anthropogenic emissions are provided by Copernicus Atmosphere Monitoring Service (CAMS). In this study, CAMS anthropogenic NOₓ emissions from 2016 are taken as a reference (or prior) bottom-up emission inventory in the inversion to estimate NOₓ emissions for Europe; therefore, the following part is dedicated to CAMS prior bottom-up inventory.
<figure><img src="image_25.png" /><figcaption>Figure 21 Main inputs and their sources for CTM CHIMERE simulations over Europe</figcaption></figure>

2.3.1.1.2.1 Anthropogenic Bottom-up Emission Inventory: CAMS

CAMS, conducted by ECMWF on behalf European Union, consists of global and regional
European anthropogenic emissions for various sectors. The regional CAMS (CAMS-REG)
inventory has a spatial resolution of 0.1° × 0.05°. It covers GHG emissions (referred to as
CAMS-REG-GHG) and air pollutants (referred to as CAMS-REG-AP). It consists of emissions
based on national reporting and a combination of other inventories (e.g., GAINS, EDGAR,
TNO inventory estimations) (J. J. P. Kuenen et al. 2014; J. Kuenen et al. 2022) and aggregated
across different sectors and fuel combinations (Table 2.2, Granier et al. 2019). Annual NOₓ
emissions are presented in the sectors such as; energy, residential, transportation, industry, and
agriculture. For the documentation describing CAMS inventories, construction and sectors, see
Granier et al. 2019.

In this study, the CAMS-REG-AP inventory v3.1 for 2016 published in 2019 was used as
anthropogenic emissions.

**2.4 Atmospheric observations**

In the framework of this thesis, NOₓ emissions are constrained using NO₂ tropospheric column
observations from the OMI and TROPOMI satellite instruments. The satellite observations
subsection is devoted to explaining the nature of these observations and the criteria by which
the observation data were selected for the inversions.
After performing the inversions, the obtained NOₓ emissions estimations were evaluated by considering independent surface observations. The surface observations subsection is devoted to the description of the ground-based networks at different domains (Eastern China and Europe) stated in this thesis.

## 2.4.1 Satellite observations

### 2.4.1.1 OMI

The Ozone Monitoring Instrument (OMI) is an ultraviolet-visible (UV-Vis) nadir-view solar backscatter spectrometer carried by the National Aeronautics and Space Administration's (NASA) Earth Observing System (EOS) Aura satellite launched in July 2004 (Figure 22). The instrument is a contribution of the Netherlands' Institute for Air and Space Development (NIVR), in collaboration with the Finnish Meteorological Institute (FMI).

<figure><img src="image_26.png" /><figcaption>Figure 22 Visualization and illustration of the Aura Spacecraft (images credit: Northrop Grumman Space Technology NGST)</figcaption></figure>

OMI has a wide field of view of 2600 km (Figure 23). This provides near global coverage with 14 orbits per day. The ground pixel size at nadir have a spatial resolution of (nadir) 13km × 24km (along × across-track) and a local overpass time of 13:40 local time (LT) (Levelt et al. 2006).
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /><figcaption>Figure 23 OMI measurement geometry (Taken from Levelt et al., 2006)</figcaption></figure>

It is a successor instrument to the European Global Ozone Monitoring Experiment (GOME) (with a spatial resolution of 320 km × 40 km) (Burrows et al. 1999), the Scanning Imaging Absorption Spectrometer for Atmospheric Cartography (SCIAMACHY) (with a spatial resolution of 30 km × 60 km) (Bovensmann et al. 1999) and NASA's Total Ozone Mapping Spectrometer (TOMS) (with a spatial resolution of 50 km × 50 km) (Wellemeyer et al. 1997) instruments. OMI is the most sophisticated instrument of its time in terms of high resolution and the ability to measure the entire UV-Vis-Near Infrared (NIR) spectrum. Figure 24 shows an example of the spectrum over the Netherlands in a clear sky.

Thanks to its two 2D charge coupled device (CDD) detectors, it covers UV1 (270–310 nm), UV2 (310–365 nm) and VIS-channel (365-500 nm) with a spectral sampling of 0.21 nm and a spectral resolution of 0.63 nm (Boersma et al. 2007). Signal-to-noise ratio (SNR) is estimated to be 400 for UV2 and 500 for VIS for 1 pixel (Zara et al. 2018, Table 1). OMI was initially designed to continue the TOMS total ozone record, so it includes TOMS-type of retrievals as well as other types of retrievals (e.g., Differential Optical Absorption Spectroscopy (DOAS) type, Raman scattering). It also provides mapping of air pollution products such as NO₂, SO₂ and aerosols from an urban to regional scale (Boersma et al. 2007; 2011). This study uses OMI
<figure><img src="image_29.png" /><figcaption>Figure 24 Solar irradiance and the Earth's radiance spectrum measured by OMI over the Netherlands. The image taken from Levelt et al., 2006</figcaption></figure>

NO₂ tropospheric columns v1.1 retrievals from the European Quality Assurance for Essential Climate Variables (QA4ECV) project (Boersma et al. 2018). This version of the OMI observations (so-called QA4ECV NO₂ data product) consists of an improved spectral fitting method for the retrieval of NO₂ slant columns (a photon path from the Sun to the atmosphere and satellite using DOAS spectral fitting procedure (Platt and Stutz 2008). Air mass factor (AMF) is another key element in the retrieval process. It is used to correct and convert the slant columns density (SCD) to vertical columns (VCD) density (AMF = SCD / VCD). However, it is calculated by using radiative transfer model and dependent on the model input parameters such as, albedo and clouds. It also requires a priori estimation of the vertical profiles of the species of interest. So, AMF calculations are subject to uncertainties. OMI QA4ECV retrieval includes data assimilation using the TM5-MP CTM model in 1° × 1° degree horizontal resolution to estimate stratospheric background NO₂ to strengthen the AMF and reduce the uncertainty associated with the AMF (Zara et al. 2018; Lorente et al. 2017) and thus improves the retrieval. Averaging kernels (AKs) correlate the derived quantities and the true profile during the retrieval (Rodgers 2000). They give information about the vertical sensitivity of the satellite retrievals, vary from one pixel to another and are used to minimize the discrepancies of model/satellite comparisons. AMFs used to calculate the tropospheric averaging kernels for below the tropopause as follows:

$$A^{\text{trop}} = \frac{M}{M^{\text{trop}}} A \qquad (\text{Eq.26})$$
where AK and tropospheric AK is denoted by the symbols *A* and *$A^{\text{trop}}$*; AMF and tropospheric AMF is denoted by *M*, and *$M^{\text{trop}}$*, respectively.

As in this thesis, we are interested in NO₂ tropospheric columns, it should be noted that OMI tropospheric NO₂ column data are defined as an integrated quantity from the surface to the tropopause for a specified surface and often expressed in molecule/cm². Validations studies show that OMI QA4ECV NO₂ column data is negatively biased or even gives good agreement (e.g., bias = -2 % over the polluted Tai’an (China) region in June 2006 (Boersma et al. 2018) compared to the surface observations from stations (Compernolle et al. 2020). For documentation of the dataset, see Boersma et al. 2017, user’s guide. In addition, the dataset is freely accessible via the Tropospheric Emission Monitoring Internet Service (TEMIS) (https://www.temis.nl/airpollution/no2.php, Last access: 08/04/2023).

### 2.4.1.1.1 Data selection and processing

In this thesis, we follow the OMI user's manual data quality recommendations to use the QA4ECV NO₂ data. The selection of OMI observation data is based on the following criteria: (i) pixels with a processing error flag equal to 0 are selected for successful processing; (ii) pixels with the solar zenith angle larger than 80 are selected to avoid observations from polar nocturnal darkness; (iii) pixels selected with the snow ice flag lower than 10 (ice coverage with max 10%) or equal to 255 (that is a flag for the ice-free ocean); (iv) pixels with the ratio of tropospheric AMF over geometric AMF is higher than 0.2 to avoid situations in which the retrieval is based on very low (relative) tropospheric AMFs; and (v) pixels with the cloud fraction is lower or equal than 0.5 meaning that most of the radiance originates from the clear-sky part of the pixel. These limit the contribution of a priori profile made for the calculation of AMF; (vi) In addition to these recommendations, we consider an additional criterion concerning the error. Data with an uncertainty of more than 100% are filtered out.

After the selection and filtering procedure, observation data needs to be processed to become comparable with the model (Jº observation cost function in Eq. 21). Since the satellite data are not on a regular grid (as in the model), and the spatial resolution of the OMI data is finer than the CHIMERE model grid cells, we sampled the observation data into super-observations. To do that, we chose the median of the OMI NO₂ tropospheric columns within the 0.5° × 0.5° grid cell and its corresponding AKs.
As we are dealing with the discrepancies between simulated and observed NO₂ concentrations, to compare with OMI super-observations, hourly-simulated NO₂ concentrations in the same grid cells are selected and are interpolated to OMI's pressure levels. The AKs are then applied to obtain a model estimate with the vertical sensitivity of the tropospheric NO₂ columns (Eq. 27).

$$y = \sum_{i=1}^{n} y_{i}^{s} A_{i} \qquad (\text{Eq. 27})$$

where y is the model equivalent of the satellite data, n is the number of levels of the satellite observations, $y_i^s$ is the simulated concentration interpolated at the i-th level and AK is denoted by the following symbol $A_i$ at the i-th level. Note that, these operations take place in the H observation operator, which includes the CHIMERE CTM to relate emissions to concentrations and transformations to calculate the simulated equivalent of observations as explained in the Section 2.2.

### 2.4.1.2 TROPOMI

TROPOMI is a nadir view solar backscatter instrument on board European Space Agency's (ESA) Copernicus Sentinel-5 Precursor (S5P) satellite, launched into a Sun-synchronous orbit in October 2017. Similar to OMI, TROPOMI has a wide field of view of 2600 km (Figure 25) and provides near global coverage with 14 orbits per day.

<figure><img src="image_30.png" /><figcaption>Figure 25 TROPOMI measurement geometry (Taken from TROPOMI user's manual, 2019)</figcaption></figure>
It is a nadir imaging spectrometer with a high spatial resolution of 3.5 km × 5.5 km (initially 3.5 km × 7 km changed on 6 August 2019) making it the most advanced instrument available today including, its predecessor, the OMI (Figures 26 and 27).

<figure><img src="image_31.png" /><figcaption>Figure 26 Spatial resolution comparison of satellite instruments. Note that TROPOMI pixel size further decreased now (Taken from BIRA-IASB, www.aeronomie.be/en/news/2021/three-years-tropomi-measurements)</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Figure 27 OMI (left) and TROPOMI (right) NO₂ tropospheric column comparison for 17/04/2018 (Taken from Lorente et al., 2019)</figcaption></figure>

The local overpass time is 13:30 LT, almost identical to OMI's. It covers the wavelength bands 270-495 nm and 710-775 nm with spectral resolution 0.55 nm and 2305-2385 nm with spectral resolution 0.25 (from ultraviolet to shortwave infrared) and it has a capacity to measure the atmospheric columns of trace gases such as NO₂, HCHO, SO₂, O₃, CO, CH₄ and aerosol and cloud properties (Figure 28). In this thesis, we are interested in NO₂ tropospheric columns.
<figure><img src="image_33.png" /><figcaption>Figure 28 Spectral Coverage of TROPOMI (Image credits: ESA)</figcaption></figure>

Compared to OMI, TROPOMI has 1-5 times higher SNR per-ground pixel (TROPOMI NO₂ data products user manual Van Geffen et al. 2019; Eskes et al. 2019). As with its predecessor OMI, the spectra are analyzed, using the DOAS technique to retrieve SCDs. Figure 29 gives an example of a NO₂ measurement spectral analysis (Borger et al. 2020).

<figure><img src="image_34.png" /><figcaption>Figure 29 Example of a spectral analysis of a TROPOMI NO₂ measurement spectrum (black line is the best fit while red line is noise) (taken from Borger et al., 2020)</figcaption></figure>

Similar to the OMI, TM5-MP CTM model is used to determine stratospheric and thus tropospheric SCDs. Then, the tropospheric SCDs are converted into a vertical column density using AMF, which is calculated using altitude-dependent AMF look-up tables, radiative transfer model, vertical distribution of NO₂ from an assimilation, cloud fraction, and pressure
information from the FRESCO-S algorithm (unlike to OMI) (Lorente et al. 2017; van Geffen et al. 2022). At the beginning of this thesis, the TROPOMI L2 NO2 OFFL (offline) v1.2 and later v1.3 data products were operationally available at the S5P Pre-Operations Data Hub. Since TROPOMI is a new satellite instrument, comparison of TROPOMI L2 NO2 OFFL v1.2 and v1.3 products with ground-based observations shows that underestimation of up to 51% and 37% in extremely and highly polluted areas, respectively (Verhoelst et al. 2021, Table 4). This underestimation may be due to the FRESCO-S cloud pressure algorithm or the use of low-resolution a priori NO₂ profiles from TM5-MP compared to the high resolution of TROPOMI. This problem necessitated an improvement of the data product and since the beginning of 2022; a new product has been made publicly available through the ESA's Sentinel-5P Product Algorithm Laboratory (S5P-PAL). This new PAL v2.2 product shows on average 10-40% higher tropospheric VCDs compared to v1.2 and v1.3 data, depending on pollution level and season (van Geffen et al. 2022). The latest version of the data was reprocessed between 1 May 2018 to 14 November 2021 and is labeled as PAL v2.3 which covers our study period in this thesis (Eskes et al. 2021). In this thesis, in subsection 4.1 and 4.2 we use v1.2 (at the time of the study PAL v2.3 products were not released) and in the rest of Chapter 4 we use PAL v2.3 products of TROPOMI.

### 2.4.1.2.1 Data selection and processing

The processing of TROPOMI observations is the same as for OMI observations (i) to (iv). To ensure the quality of the retrieval results, the data is distributed with different quality flags corresponding to each ground-pixel and ranging from 0 to 1. In this thesis, we have chosen to use a quality flag value (qa_value) of 0.75, which is a more restricted and also recommended filter. Thus we use pixels with a qa_value greater than 0.75. In addition to that, unlike OMI, data with more than 60% uncertainty are filtered out. This is due to the large amount of TROPOMI observations compared to OMI. We chose to be more restrictive when using TROPOMI data in the inversion to use less informative data. And last, the steps for building super observations are the same as in OMI.

## 2.4.2 Surface observations

Since true NOₓ emissions are not known (and cannot be) in the atmosphere, we used other observational tools currently available to assess our estimated NOₓ emissions, constrained by satellite NO₂ observations within the CIF-CHIMERE inversion system. To do that, in this
thesis, first we perform CHIMERE simulations using the prior and estimated NOₓ emissions (from the inversions) to obtain simulated NO₂ concentrations. Then, the simulated surface NO₂ concentrations are compared to surface observations to evaluate the inversion results in case these new estimate improve the comparisons. NO₂ surface observations are from different in-situ networks according to our study domain.

### 2.4.2.1 Chinese in-situ network

In this thesis, in-situ observations provided by the China National Environmental Monitoring Centre (CNEMC) network directly affiliated to the Chinese Ministry of Ecology and Environment (MEE) were used to evaluate our estimated NOₓ emissions for China (Chapter 3). The CNEMC network established in 1980, incorporates more than 2000 air quality-monitoring stations to evaluate the status of the national air quality by hourly measuring national NO₂, O₃, CO, SO₂, PM2.5 and PM10 since 2014 (Figure 30).

<figure><img src="image_35.png" /><figcaption>Figure 30 Spatial distribution of the CNEMC Network sites for Eastern China</figcaption></figure>

#### 2.4.2.1.1 Data selection and processing

Networks include different classifications of the stations, depending on their environment types. In this work, the types and environment of the stations follow the classification of Lachatre et al. 2020 based on the Flemming, Stern, and Yamartino 2005 method. Flemming et al. describes a new classification by using hourly concentrations of O₃ and their diurnal variability to
evaluate whether stations represent rural, suburban, urban, or traffic environments. In our analyses, we select in-situ data from stations in rural, suburban and urban environments and exclude stations in traffic environments. Because they are close to pollutant sources and cannot be compared to the model with its coarse model resolution. This study considers surface NO₂ concentration measurements from 648 stations in 2015 and 626 stations in 2019, located in Eastern China (also stated in the Figure 30).

These measurements are performed using a commercial chemiluminescence analyzer equipped with molybdenum converter. During the performance of this measurement analyzer catalyzing NO₂ to NO, other reactive oxidized nitrogen compounds (NOz) such as PAN and HNO₃ can also be converted to NO, resulting in an overestimation of ambient NO₂ concentrations (Lamsal et al. 2008). To compensate that, Lamsal et al. suggested using the following correction factor (CF) formula to alleviate and avoid this overestimation of NO₂ (Eq.28).

$$CF = \frac{NO_2}{NO_2 + \sum AN + 0.95PAN + 0.35HNO_3} \quad (Eq.28)$$

where the sum of all alkyl nitrates (AN), PAN, and HNO₃ at a resolution of 0.5° × 0.5° are calculated running CHIMEREv2013 simulations with the settings explained in Section 2.3.1.

The values of correction factors are between 0 and 1. In summer, the value is smaller, which means that larger corrections occur, while in winter, it is found to be closer to 1 (Lamsal et al. 2008) and Figure 31. This is because NOₓ has a longer lifetime in winter. These correction factors are then bilinearly interpolated by station locations and applied (by multiplication) to the hourly measured surface NO₂ concentrations at these stations.
<figure><img src="image_36.png" /><figcaption>Figure 31 Correction factor examples winter (left) vs summer (right) for 1 day</figcaption></figure>

### 2.4.2.2 European in-situ network

In-situ NO₂ concentrations from national monitoring networks were retrieved from the European Environmental Agency (EEA, https://www.eea.europa.eu, Last access: 30/03/2023) under the direction of EU. The network provides hourly monitoring data for various pollutants such as NO₂ and O₃ for a large number of stations (more than 3500 stations) since 2013, not only for EU member states but also for non-EU countries located on the European continent (Figure 32). The data provider has predetermined the classification of the stations depending on the predominant emission sources, and distribution of the building around the station. According to these, urban, suburban and rural classification represents the surrounding of the station according to the density of the building from higher to less dense. Traffic stations are located close to a main road (or roads), industrial stations are located close to industrial sources and background stations are located where the general population or vegetation is exposed to average pollution levels. The same correction factor as used for Chinese surface measurements is applied.
<figure><img src="image_37.png" /><figcaption>Figure 32 Spatial distribution of the EEA sites for Europe (the image is produced in https://www.eea.europa.eu/data-and-maps/explore-interactive-maps/up-to-date-air-quality-data, accessed on 30/03/2023)</figcaption></figure>

## 2.5 Sensitivity analysis of the CIF inverse modeling system

Since CIF is a new inverse modeling system, we were interested in assessing the sensitivity of some of its parameters, which are used to determine the optimization and convergence of the system and are critical for variational inversion. This allows us to identify which parameters have a higher influence on the outputs of the inversions and to guide our selection parameters more consistently when setting up the inversion system, as in the applications in Chapters 3 and 4. To this end, more than 40 inversion runs for the same day were done over Europe using OMI satellite observations. At the time of the study, the stabilized version of the TROPOMI product (PAL version v2.3.1) was not available. We chose a day in winter where NO₂ are usually higher. The date of 08/02/2020 was chosen because the number of OMI observations was reasonable (amongst the largest for this period) (Figure 33) and the observed NO₂ columns cover a representative range of values. In this section, an overview of these tests is presented. Either a set of parameters, related to the optimization or variational processes, which can influence the inversion results are selected to be tested. A reference test is performed with parameters fixed at values that we estimate to be suitable for a first try based on literature or expert knowledge. The other tests covering a range of values for the different parameters are compared with this
reference test by changing one parameter at a time. In the following, these parameters are described for inversion, firstly as optimization parameters (less effective on the solutions), then as **B** error parameters (more effective) and how they are selected is explained.

<figure><img src="image_38.png" /><figcaption>Figure 33 OMI tropospheric column for the reference day (08/02/2020)</figcaption></figure>

## 2.5.1 Parameters definition

### 2.5.1.1 Optimization parameters

CIF-CHIMERE uses the M1QN3 limited-memory quasi-Newton minimization algorithm to find the best estimate of NOₓ emissions that minimizes the cost function J (Gilbert and Lemaréchal 1989). The parameters related to the performance of the minimization algorithm can be classified under technical parameters. These parameters are as follows:

(i) stopping criteria (so-called ‘epsg’). It sets a threshold set on the ratio between the final and initial norms of the cost function gradient. For example, if less than 0.05 is taken as epsg, this equates to a requirement to reduce the final norm by at least 95%.

(ii) maximum number of iterations (so-called maxiter). It sets a limit to the number of iterations allowed in the algorithm. This allows the algorithm stop if epsg is not satisfied.

(iii) Step size (so-called ‘df1’). It is the expected decrease in the gradient norm of the cost function during the first iteration and determines the first step size.
Both (i) and (ii) allow the algorithm to stop. The inversion is considered convergent if it reaches
a certain amount of reduction in the norm of the gradient cost function before reaching the
maximum number of iterations. Otherwise, the system cannot find any minimum (case non-
convergent) in which case the criteria must be relaxed (should increase maxiter etc.).

It is important to note that the M1QN3 minimizer uses a line search before each iteration, so
that each iteration consists of at least one simulation. This number of simulations (so-called
nsim) can also be limited and taken as a criterion for performance. However, it is not considered
in our analysis. More information about the M1QN3 minimizer and its parameters can be found
in (Gilbert, Lemaréchal, and Simulator 2009).

2.5.2 Variational parameters

Besides the optimization parameters, the inversions are mainly driven by the error statistics considered in the variational processes: the **B** and **R** matrices and their respective weight (Anderson et al. 2000). The **R** matrix is usually well prescribed by the observation errors provided with the satellite products. Concerning the **B** matrix, attributing the errors and covariances to the prior emissions at the grid-cell size is more complicated. Indeed, this information is usually not provided with the emissions inventories. Uncertainties are mainly estimated for national and annual budgets, given as a range of values, and how they redistribute spatially and temporally is difficult to assess. In practice, the errors values considered for **B** are assigned from a range of values from the inventories, the literature, and the expert knowledge. They may also define the range of “freedom” one gives to the system to reach a reasonable solution. In this thesis work, we tested this range of possible values to guide our configuration choice accounting both the performances of the inversion and the computation time. This means to define range of parameter values one wants to test for the errors in NOₓ emissions and their spatial (or horizontal) correlations. These parameters are called variational parameters and are as follows:

(iv) Emission errors (**B** matrix variances).

(iv) Correlation of the emission errors (B matrix covariances). To consider these correlations, we did the assumption of Gaussian correlations of the errors and use the spatial correlation parameter defined in CIF as hcorr:sigma, which defines the horizontal correlation of grid cells by the sigma factor. It is formulated as follows;
$c(i, j) = \exp(-\text{dist}(i, j)/\sigma)$ where c, the correlation between i-th grid cell and j-th grid cell, is a decreasing exponential of the distances over sigma. Thus, an increasing sigma increases the correlation between grid cells. We consider the sigma values as the variable parameter in the tests. Note that, vertical correlation in the uncertainty is ignored.

For sensitivity analysis, the following configuration is then taken as a baseline in this thesis: 70 % anthropogenic emission error, biogenic emission errors are ignored, the spatial correlation is assumed 50 km, the stopping criterion is taken as 0.05 with a maximum of 10 iterations are allowed with the first step size of 0.01.

## 2.5.2 Sensitivity to optimization parameters

Table 2 summarizes the tests performed by changing the technical parameters on the day of 08/02/2020. Test A is taken as a reference with anthropogenic NOₓ emission error 70 %, no error in biogenic NOₓ emissions, spatial correlation as 50 km, and technical parameters; maxiter, epsg, df1 are as 10, 0.05 and 0.01, respectively. It can be deduced that the reference Test A converges in 4 iterations with the result that the norm of the gradient of the cost function J decreases by 99 %. According to the M1QN3 user manual, a too small value of df1 causes the algorithm to stop at the first step; however, a large value of df1 is reduced after a few iterations. Therefore, we tested for larger values at df1 (e.g., 0.1, 0.2 and 0.3), but did not obtain significant changes (from tests B to D) in terms of emission corrections. It almost only increases the computation duration. Next, we tested to go further in the minimization algorithm by stressing the stopping criterion. The epsg parameter is then decreased to 0.01 in Test E and to 0.001 in Test F. The maximum number of iterations is increased to avoid the system stops before reaching the epsg criterion. Reference Test A and Test E give the same results as the reference test stops with a reduction of the gradient norm close to 99 %. Test F goes further in the reduction of the gradient norm and reaches the convergence but without any significant impact on the differences of emissions between the prior and the posterior (Figure 34) and increasing the computation time (Table 2). Due to the limited impact of the optimization parameters on the inversion results but their implication in terms of computation time, these parameters are taken as in reference Test A for the next inversion tests.
<table><caption>Table 2 Inversion tests with changed optimization parameters for the day of 08-02-2020.</caption><thead><tr><th rowspan="2">Name of the test</th><th rowspan="2">Anthropogenic emission error (%)</th><th rowspan="2">Spatial correlation (km)</th><th rowspan="2">Biogenic emission error (%)</th><th colspan="3">Optimization Parameters</th><th rowspan="2">Duration (in hour)</th><th rowspan="2">Convergent? (Yes / No)</th></tr><tr><th>Max iteration</th><th>Stop: epsg</th><th>Step size: df1</th></tr></thead><tbody><tr><td>A</td><td>70</td><td>50</td><td>-</td><td>10</td><td>0.05</td><td>0.01</td><td>3</td><td>Y</td></tr><tr><td>B</td><td>70</td><td>50</td><td>-</td><td>10</td><td>0.05</td><td>0.1</td><td>2.5</td><td>Y</td></tr><tr><td>C</td><td>70</td><td>50</td><td>-</td><td>10</td><td>0.05</td><td>0.2</td><td>4</td><td>Y</td></tr><tr><td>D</td><td>70</td><td>50</td><td>-</td><td>10</td><td>0.05</td><td>0.3</td><td>10</td><td>Y</td></tr><tr><td>E</td><td>70</td><td>50</td><td>-</td><td>20</td><td>0.01</td><td>0.01</td><td>3</td><td>Y</td></tr><tr><td>F</td><td>70</td><td>50</td><td>-</td><td>20</td><td>0.001</td><td>0.01</td><td>5</td><td>Y</td></tr></tbody></table>

<figure><img src="image_39.png" /><figcaption>Figure 34 NOₓ emission corrections obtained by Test A (left), Test E (middle), Test F (right). Note that, the corrections are the same.</figcaption></figure>

### 2.5.3 Sensitivity to variational parameters

In this section, we tested different parameters related to the **B** error covariance matrix, which is a key parameter in variational inversion (Sections 2.1 & 2.2), with the aim of building **B** matrix and analyzing how the CIF-CHIMERE inversion system reacts to different estimations of these values. To do so, we need to know the uncertainties and correlations on emissions, but we have low knowledge about these uncertainties as discussed previously. For this reason, we perform different tests over a range of uncertainties.

#### 2.5.3.1 **B** variances

Variances of **B** are the diagonal terms of the matrix representing the different anthropogenic emission errors. In Europe, the uncertainty of national anthropogenic NOₓ emission estimates stays at 20-30 % (Citepa 2021), but uncertainty values are more likely to be larger on a smaller scale (grid-cell scale). In literature, Müller and Stavrakou 2005 reported 50 % of uncertainties in their inversions for Europe. Souri et al. 2020 assumed these uncertainties to be 50 % for their
inversions in Asia. Therefore, to evaluate the impact of the choice of the values in correcting NOₓ emissions without pushing the limits of computational resources, we considered B variances for our tests between 30 % and 100 %. Table 3 summarizes the performed tests. It is important to note that these errors indicate how much confidence we have in the prior emissions and which freedom we allow the system to have. A higher error value gives the inversion system more freedom to correct emissions. Hereunder, test I, which has the largest error value we tested, is not convergent with the given parameters.

<table><caption>Table 3 Inversions tests with different errors in anthropogenic NOₓ emissions for the day of 08-02-2020.</caption><thead><tr><th rowspan="2">Name of the test</th><th colspan="3">Variational Parameters</th><th rowspan="2">Max iteration</th><th rowspan="2">Stop: epsg</th><th rowspan="2">Step size: df1</th><th rowspan="2">Duration (in hour)</th><th rowspan="2">Convergent? (Yes / No)</th></tr><tr><th>Anthropogenic emission error (%)</th><th>Spatial correlation (km)</th><th>Biogenic emission error (%)</th></tr></thead><tbody><tr><td>A</td><td>70</td><td>50</td><td>-</td><td>10</td><td>0.05</td><td>0.01</td><td>3</td><td>Y</td></tr><tr><td>I</td><td>100</td><td>50</td><td>-</td><td>10</td><td>0.05</td><td>0.01</td><td>12</td><td>N</td></tr><tr><td>J</td><td>60</td><td>50</td><td>-</td><td>10</td><td>0.05</td><td>0.01</td><td>2.5</td><td>Y</td></tr><tr><td>K</td><td>50</td><td>50</td><td>-</td><td>10</td><td>0.05</td><td>0.01</td><td>2</td><td>Y</td></tr><tr><td>L</td><td>40</td><td>50</td><td>-</td><td>10</td><td>0.05</td><td>0.01</td><td>2</td><td>Y</td></tr><tr><td>M</td><td>30</td><td>50</td><td>-</td><td>10</td><td>0.05</td><td>0.01</td><td>2</td><td>Y</td></tr></tbody></table>

Tests J to M have an emission error that decreases by 10%, respectively. They are all convergent. Tests with smaller errors require shorter computation time to converge compared to the reference test A. Figure 35 shows the spatial distribution of the emission corrections obtained from these tests. It is seen that the range of corrections decreases with respect to the error value. More precisely, these corrections range from ±90 %, ±70 %, ±50 %, ±30 % and ±20 %, corresponding to the errors of 70 %, 60 %, 50 %, 40 % and 30 %, respectively. While tests L and M are too restrictive, test K provides sufficient freedom when the error parameter B is 50 % and the corrections are within the 50 % range. This means that, the inversion system is well responding to the error changes in the anthropogenic NOₓ emissions.
<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>

Figure 35 Increments according to the different errors in prior NOₓ emissions (From left top to bottom left, test A and tests J to M, respectively) (differences<1% is neglected).

### 2.5.3.2 B covariances

Defining the covariances or non-diagonal terms of **B** matrix is even more difficult to evaluate.
First NOₓ inversion studies neglected to spatial correlations of emissions errors making the
assumption that at medium/coarse resolution the transport of NO₂ between neighbouring grid cells can be neglected (e.g., Martin 2003). Konovalov et al. 2006 discussed this assumption with a CHIMERE at 0.5° resolution and concluded that considering two layers of neighbouring grid cell were sufficient to represent the NO₂ transport. Based on this, we define the range of variation of the sigma parameter used to prescribe the horizontal correlation in our case as an exponentially decaying function as mentioned previously. We tested different correlation lengths based on the CHIMERE resolution (50 km) from 0 to 100 km. Table 4 summarizes these tests for different spatial error correlations of emissions.

<table><caption>Table 4. Inversions tests with different spatial correlations for the day of 08-02-2020</caption><thead><tr><th rowspan="2">Name of the test</th><th colspan="3">Variational Parameters</th><th rowspan="2">Max iteration</th><th rowspan="2">Stop: epsg</th><th rowspan="2">Step size: df1</th><th rowspan="2">Duration (in hour)</th><th rowspan="2">Convergent? (Yes/No)</th></tr><tr><th>Anthropogenic emission error (%)</th><th>Spatial correlation (km)</th><th>Biogenic emission error (%)</th></tr></thead><tbody><tr><td>A</td><td>70</td><td>50</td><td>-</td><td>10</td><td>0.05</td><td>0.01</td><td>3</td><td>Y</td></tr><tr><td>G</td><td>70</td><td>100</td><td>-</td><td>10</td><td>0.05</td><td>0.01</td><td>2.5</td><td>Y</td></tr><tr><td>H</td><td>70</td><td>-</td><td>-</td><td>10</td><td>0.05</td><td>0.01</td><td>2.5</td><td>Y</td></tr></tbody></table>

Figure 36 shows the difference between estimated anthropogenic NOₓ emissions and prior emissions in %, from left to right, without considering any correlation, a correlation taken as 50 km (same as the CHIMERE model resolution) and 100km (correspond to the tests in the Table 4).

<figure><img src="image_45.png" /><figcaption>Figure 36 Anthropogenic NOₓ emissions increments (in %) provided by inversions test from left to right; with no correlation (Test H), sigma=50 (Test A, reference test) and sigma=100 (Test G), respectively (differences<1% is neglected).</figcaption></figure>

According to these tests, neglecting correlation in correction of the emissions (leftmost) leads to corrections only localized compared to the tests where correlations is considered (middle and rightmost). According to the study of Konovalov et al. 2006, this assumption is likely too strict for the considered CHIMERE resolution. Accounting for a spatial correlation with a Gauss
function with the half width (sigma) equal to the CHIMERE resolution fills the recommendation of Konovalov et al. 2006 study of considering two layers of neighboring grid cells. Test G, with sigma equal to 100 km gives too diffusive corrections. Another “advantage” of considering horizontal correlation is somehow to counterbalance the lack of observations such it can be the case with OMI. Therefore, it seems reasonable to take spatial correlations with the correlation length derived from the model resolution as in the reference test A (middle).

### 2.5.3.3 Biogenic emission error

So far, we have performed inversion tests without taking into account the error in biogenic emissions. In winter months NOₓ emission is dominated by anthropogenic activities but in summer months biogenic emissions (biomass burning, lightning and soil) are also a substantial source (Vinken et al. 2014). Like anthropogenic emissions, estimates of biogenic emissions are subject to uncertainties; e.g., lack of good quality datasets due to highly variable agricultural activities such as fertilization of soils, vegetation types and soil management practices (Skiba et al. 2021). These estimates are often larger than these for anthropogenic emissions are. Therefore, we consider to errors between 100 % and 200 % in biogenic emissions in our tests and we include a summer day where the impact could be more important (Table 5).

<table><caption>Table 5 Inversions tests with different errors in biogenic NOₓ emissions for two days 08-02-2020 and 28-08-2019</caption><thead><tr><th rowspan="2">Name of the test</th><th rowspan="2">Day</th><th colspan="3">Variational Parameters</th><th rowspan="2">Duration (in hour)</th><th rowspan="2">Convergent ? (Yes/No)</th></tr><tr><th>Anthropogenic emission error (%)</th><th>Spatial correlation (km)</th><th>Biogenic emission error (%)</th></tr></thead><tbody><tr><td>N</td><td>08/02/2020</td><td>50</td><td>50</td><td>100</td><td>2</td><td>Y</td></tr><tr><td>O</td><td>08/02/2020</td><td>50</td><td>50</td><td>200</td><td>2.5</td><td>Y</td></tr><tr><td>P</td><td>28/08/2019</td><td>50</td><td>50</td><td>100</td><td>2</td><td>Y</td></tr></tbody></table>

Note that the anthropogenic emission error is assigned to 50 % as discussed in 2.5.3.1, 2005 for Europe. Figure 37 shows the biogenic NOₓ emission differences from tests N and O (respectively). According to the figure, we can conclude that, doubling the biogenic emissions errors results in only minor changes in the correction of the emissions – the corrections are smaller than 1 % in winter. On the other hand, on a summer day (Figure 38), when taking the same error (100 %) in biogenic emissions, the correction of the emissions is ten times more significant compared to a winter day but remains small with ±10 %. Therefore taking 100 % on error in biogenic emissions would be sufficient for our next inversion configuration.
<figure><img src="image_46.png" /></figure>

Figure 37 Increments according to the 100% (Test N, left) and 200% (Test O, right) errors in biogenic NOₓ emissions (differences<1e-5% is neglected).

<figure><img src="image_47.png" /></figure>

As a conclusion, these are preliminary tests to evaluate how the state-of-the-art CIF-CHIMERE inversion system responds to a number of changing parameters and to prepare for the inversions for Europe and China in the following chapters. These tests show that our choices based on expert knowledge are reasonable and consistent with the literature. For a comprehensive sensitivity study, more tests should be conducted, as one day is not representative. e.g., for seasonality.

The final configuration we choose is 50 % anthropogenic emission errors with 50 km spatial correlation and 100 % biogenic emission errors for Europe. For China, NOₓ emissions uncertainties reported in literature are slightly larger than in Europe, we then decide to increase the B variances to 70%.

<figure><img src="image_48.png" /></figure>

Figure 38 Increments using 100% biogenic error emissions in a winter day (Test N, left) and in a summer day (Test P, right) (differences<1e-5% is neglected).
# Third part : Scientific
# Applications

## Chapter 3. Inverse modeling of NOₓ
emissions over China

This chapter is an adapted and expanded version of the scientific paper with its supplementary
elements published in MDPI's journal Atmosphere under the title "Anthropogenic NOₓ
Emission Estimates over East China for 2015 and 2019 Using OMI Satellite Observations and
the New Inverse Modeling System CIF-CHIMERE" and can be found as Savas et al., 2023 on
https://doi.org/10.3390/atmos14010154 (Last Access: 08/04/2023). We acknowledge the open
access license policy of the publisher MDPI.

## Summary

China is the world's largest emitter of air pollutants (Lin et al. 2014). Over the last few decades, NOₓ emission levels in China have further increased due to growth in industrial activities, urbanization and weak pollution regulations (Jin, Andersson, and Zhang 2016; Shaw and Van Heyst 2022). However, in March 2011, the Chinese government for the first time introduced regulations to reduce the level of NOₓ pollutants and control these emissions, called the 12th Five-Year Plan (2011-2015). The target of this regulation was to reduce the annual Chinese NOₓ emissions nationwide by 10% by 2015 compared to the levels of 2010 (de Foy, Lu, and Streets 2016). Studies have shown that these regulations can reduce NOₓ emissions and the Chinese government has achieved its emission reduction target by reducing total NOₓ emissions by 18.6% (Liu et al. 2016, 13th Five-Year Plan). Since then, the changes in NOₓ emissions have been assessed using various approaches to evaluate the impact of the regulations (e.g., (Zheng et al. 2018; Xi Lu et al. 2020). Complementary to the previous studies, this study estimates anthropogenic NOₓ emissions in 2015 and 2019 over Eastern China using as a reference the Hemispheric Transport of Air Pollution (HTAP) v2.2 emission inventory for 2010 and the new variational inversion system the Community Inversion Framework (CIF) interfaced with the CHIMERE regional chemistry transport model and OMI satellite observations. We also compared the estimated NOₓ emissions with the independent Multi-Resolution Emission
Inventory for China (MEIC) v1.3, from 2015. The inversions show a slight global decrease in NOₓ emissions (in 2015 and 2019 compared to 2010), mainly limited to the most urbanized and industrialized locations. In the locations such as Baotou, Pearl River Delta, and Wuhan, the estimations in 2015 compared to 2010 are consistent with the target reduction (10%) of the 12th Five-Year Plan. Comparisons between our emission estimates and MEIC emissions in 2015 suggest that our estimates likely underestimate the emission reductions between 2010 and 2015 in the most polluted locations of Eastern China. However, our estimates suggest that the MEIC inventory overestimates emissions in regions where MEIC indicates an increase of the emissions compared to 2010.

## 3.1 Introduction

Nitrogen oxides (NOₓ=NO₂+NO) are a family of highly reactive trace gases that cause air quality deterioration in the troposphere. They play a key role in the formation of ozone (O₃) through photochemistry and the formation of secondary organic aerosols (Seinfeld and Pandis 2016; Ng et al. 2007). Most tropospheric NOₓ are emitted as NO and rapidly converted to NO₂ by various reactions (Atkinson 2000). Anthropogenic activities such as fossil fuel combustion are dominant NOₓ emission sources, especially in industrialized regions (Houghton et al. 2001). On the other hand, natural and biogenic sources can be a major source of NOₓ, such as soil NOₓ emissions in agricultural areas (Almaraz et al. 2018).

Over the past few decades, NOₓ emission levels in China increased due to growth in industrial activities and weak pollution regulations (Claire Granier et al. 2011; Shaw and Van Heyst 2022; Jin, Andersson, and Zhang 2016). The control of NOₓ emissions was added to the regulations in March 2011 with the 12th Five-Year Plan (2011–2015) to reduce the annual Chinese NOₓ emissions nationwide by 10% compared to 2010 (de Foy, Lu, and Streets 2016). Studies have shown that the Chinese Government met its emission reduction target due to these regulations in this period and even exceeded the plan's target e.g., (CAA (Clean Air Asia) 2016; Anger et al. 2016). This plan is followed by a 13th Five-Year Plan (2016–2020) with a target of reducing Chinese NOₓ emissions nationwide by 15% in 2020 compared to 2015 (“Notice of the State Council on Printing and Distributing the ‘13th Five-Year’ Ecological and Environmental Protection Plan, Government Information Disclosure Column” 2016).

Chemistry transport models (CTMs) have been widely used to study the efficiency of environmental regulations, assess their impacts, and make future predictions through different
scenarios (Anger et al. 2016; Colette et al. 2012; Zhao, Zhang, and Nielsen 2014; Permadi, Kim Oanh, and Vautard 2018; Nellie Elguindi et al. 2020). Emission inventories are one of the main inputs of CTMs, and their spatial and temporal accuracies are essential for simulations and forecasts. One common way to build an emission inventory is the bottom-up approach. This approach quantifies emissions based on the local distribution of activities, emission factors, and annual self-declarations by emitters. Numerous studies estimated China's NOₓ emissions using these activity rates and species' emission factors (Xia, Zhao, and Nielsen 2016; Shi et al. 2014). For example, Li, Zhang, et al. 2017 calculated the annual anthropogenic Chinese NOₓ emissions for 2008 as 26.55 Tg using the Multiresolution Emission Inventory for China (MEIC) v1.0 and 20.66 Tg using Hemispheric Transport of Air Pollution (HTAP) v2.2 bottom-up emission inventories. For 2010, Li, Zhang, et al. 2017 and Zheng et al. 2018 estimated annual NOₓ emissions as 26.5 Tg and 29.1 Tg, respectively, using two different versions of MEIC bottom-up inventory (v1.3 and v1.0, respectively). However, even though the information comes firsthand, significant uncertainties remain. Using Monte Carlo simulations, Zhao, Zhang, and Nielsen 2013 predicted Chinese annual NOₓ emissions' uncertainties to be 35 % of estimated emissions in 2010. Crippa et al. 2018 estimated NOₓ emission uncertainties in 2012 as 56.2 % for top emitting regions in China. It is worth noting that the uncertainties are given on an annual and national basis. They might be larger when considering grid points or smaller timescales. Moreover, the collection of data to build bottom-up inventories is time-consuming; then, updating these inventories may take several years (Ding et al. 2017).

On the other hand, using inverse modeling techniques, which constrain emissions by atmospheric observations, brings complementary information and, in some regions such as China, captures current emission trends better than traditional bottom-up inventories (Nellie Elguindi et al. 2020). In addition, unlike bottom-up inventories, these top-down inventories provide broad spatial coverage and have a fast update ability, as emissions can be estimated when observations are ready to interpret. Nevertheless, top-down inventories have some limitations and uncertainties related to the performances of CTMs on transport and chemical processes (de Foy et al. 2014), and the observations themselves (e.g. limited sensitivity to the surface and cloud coverage for satellite observations) (Boersma, Vinken, and Eskes 2016). Thus, different inversion techniques have been developed to estimate NOₓ emissions. For example, earlier studies used a mass balance approach to constrain NOₓ emissions on various spatial scales (Leue et al. 2001; Cooper et al. 2017; Han, Kim, and Song 2020; Itahashi et al. 2019; Martin 2003). Most of the studies assume a linear relationship between NOₓ prior
emissions and NO₂ columns, and do not consider horizontal emission transport. With the increase in computational capacity, advanced inversion techniques such as the Kalman filter (Napelenok et al. 2008) and variational data assimilations (Kurokawa et al. 2009) have been developed to take into account the transport of emissions and nonlinearities in chemical reactions. Ronald J. van der A et al. 2017 derived anthropogenic NOₓ emissions on a 25 × 25km² provincial level using the DECSO v4 algorithm (Mijling and van der A 2012; Ding et al. 2015; Ding et al. 2017), which estimates emissions using Kalman filter, the regional CHIMERE v2013 CTM (Menut et al. 2013), MEIC v1.0 inventory for 2010 (Zheng et al. 2018; Li, Liu, et al. 2017) and Ozone Monitoring Instrument (OMI) NO₂ observations using DOMINO v.2 algorithm (Boersma et al. 2011) for the period 2005–2015 in Eastern China. They found an increasing trend in total NOₓ emissions until 2012, except for 2009. Even though the total NOₓ emission peaks occur in 2012, different peaks occur between 2011 and 2014 on a provincial level. Qu et al. 2022 co-assimilated anthropogenic and biomass burning NOₓ, SO₂, and CO emissions from bottom-up inventory HTAP v2.2 2010 Janssens-Maenhout et al. 2015 in China using a sector-based variational algorithm and GEOS-Chem Adjoint model (Henze, Hakami, and Seinfeld 2007) for January 2005–2012. Unlike Ronald J. van der A et al. 2017, their inversion captures the NOₓ emissions peak for 2011. Moreover, they found that all estimated NOₓ emission top-down inventories (sector-based & species based) are smaller than HTAP v2.2 bottom-up inventory by 21–26 % in January 2010. Zheng et al. 2018 quantified anthropogenic NOₓ emissions using a combination of MEIC v1.3 bottom-up inventory and index decomposition analysis. They compared their bottom-up emission estimates with other studies, which include OMI observations and top-down emission estimates between 2010 and 2015 for China and Eastern China. According to their study, all top-down NOₓ emission inventories calculated by inverse modeling using OMI NO₂ columns show smaller relative differences between these two years than bottom-up estimations and OMI observations (Zheng et al. 2018)(p.14106, Table 2).

The previous studies show a significant range of uncertainties regarding the NOₓ emissions from bottom-up and top-down approaches and a significant spread in the emissions derived from different approaches (Ding et al. 2017). The changes in the emissions such as the emission reduction targets of governmental mitigation strategies, which the inventories need to assess, are usually smaller than the uncertainties on the emissions and then an important challenge for the inventories. A variety of tools to estimate emissions is essential to assess the emissions and their related uncertainties better. Our study proposes the use of a new inverse modeling system,
the Community Inversion Framework (CIF) (Berchet et al. 2021), interfaced with the CHIMERE chemistry-transport model (Menut et al. 2013) and its adjoint (Fortems-Cheiney et al. 2021) and OMI satellite observations to estimate Chinese NOₓ emissions in 2015 and 2019. The objectives of this study are twofold: (i) testing the new CIF-CHIMERE over China and evaluating if the derived emissions and concentrations simulated using them are closer to independent inventory and measured concentrations, respectively; (ii) evaluating if the NOₓ emission regulations planned in the two last Five-Year Plans are visible in the emissions derived with the new CIF-CHIMERE system. We choose to focus on the last year (2015) of the 12th Five-Year Plan. We do not consider the last year of the 13th Five-Year Plan, 2020 as it was perturbed by the COVID-19 pandemic. We analyze the NOₓ emission estimates at the scale of Eastern China, at the scale of the most urbanized, and industrialized Chinese locations. We evaluate the estimated NOₓ emissions by comparing with the independent inventory MEIC v1.3 and surface NO₂ observations from the Chinese air quality network. Section 3.2 describes the methodology and tools used in this study. Results and discussions are presented in Section 3.3.

## 3.2 Materials and methods

### 3.2.1 Inversion-related tools and data

#### 3.2.1.1 Inverse modeling and the CIF system

The atmospheric inverse modeling approach applied to emission estimates aims at constraining a prior knowledge of emissions with atmospheric observations (concentrations or columns). A chemical transport model combined with an observation operator is used to map the emissions into the observation space. The emissions are usually constrained at the model grid point and at a temporal resolution better than annual. Due to the high dimension of the problem, sophisticated approaches (Bocquet et al. 2015) are necessary. The evaluation of the estimated or posterior emissions can only be performed indirectly by comparing concentrations simulated using these posterior emissions (named Conc_Post in the following) to independent measured concentrations. This step is often named “experiment run”. To evaluate if the inversion succeeds to better assess the emissions compared to our prior knowledge used, it is also necessary to perform a “control run” to simulate the concentrations using the prior emissions (named Conc_Prior in the following) and quantify if the Conc_Post are closer to the measured
concentrations than the Conc_Prior. Figure 39 summarizes the different steps needed to perform atmospheric inversions and to evaluate them.

<figure><img src="image_49.png" /><figcaption>Figure 39 Schematic view of the atmospheric inversion and evaluation.</figcaption></figure>

As the dimension of the problem to solve is high, we use a variational approach based on a Bayesian framework and normal distributions of the errors and uncertainties. The system aims to find an optimal solution of the control vector x by minimizing the cost function J (Equation 29).

$$J(x) = \frac{1}{2}(x - x_b)^T B^{-1}(x - x_b) + \frac{1}{2}(H(x) - y_o)^T R^{-1}(H(x) - y_o) \quad (\text{Eq.29})$$

where $x_b$ is a prior estimate of the control vector, B is the prior error covariance matrix, $H(x)$ is the observation operator matrix which links control variables to atmospheric measurements at specified locations and time using a CTM (CHIMERE in this study), $y_o$ is the observation vector (OMI NO₂ columns in this study), and R is the covariance matrix of the observation errors. In our case, the control variables are not emissions, but rather daily multiplicative factors applied to NOₓ emissions. During a given day, the corresponding multiplicative factor is used to scale homogeneously the hourly emission values, retaining the relative variations of the diurnal profile. This relationship between control variables and emissions can be stated as follows: