Deriving regional and point source nitrogen oxides emissions in China from TROPOMI using the directional derivative approach with nonlinear chemical lifetime fitting

An appropriate representation of the NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio and NO_x lifetime is essential for estimating NO_x emissions from satellite NO₂ observations. We introduce a satellite-based, data-driven approach that applies variable NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio and derives the nonlinear chemical lifetime using a piecewise fitting method based on the directional derivative approach (DDA). This method enables the estimation of both regional and point-source NO_x emissions across China, representing the first application of a lightweight, satellite-driven method to directly capture nonlinear NO_x lifetime for emission estimation over large, topographically complex region. The incorporation of a variable NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio enhances the accuracy of source divergence and emission estimates and the improved fitting scheme captures the nonlinear behavior of NO_x chemistry. Anthropogenic contributions are isolated by subtracting natural sources from satellite-derived total emissions, with natural NO_x identified using a seasonal criterion and further constrained by Nighttime Light (NTL) data. Estimated anthropogenic NO_x emissions in China from 2019 to 2024 are 20.2, 18.5, 19.4, 18.9, 20.7 and 18.8 Tg, respectively, with annual uncertainties of 27 %–30 %. These values show good agreement with both bottom-up inventories and top-down inversions, with national scale discrepancies ranging from −11.8 % to 0.8 %. The DDA captures key spatial and temporal emission patterns, including consistent decline in NO_x emissions in megacities and provincial disparities linked to urbanization and economic development. The DDA estimates are consistent with previous studies on power plant emissions, and emissions from 124 plants vary between 0.02–2.13 kg s⁻¹ for 2019–2024, with uncertainties spanning 4 %–78 %, averaging 16 %. This satellite-based, lightweight method enables low-latency, long-term monitoring of NO_x emissions and offers a promising alternative to bottom-up inventories and resource-intensive top-down models. The data are publicly available at https://doi.org/10.5281/zenodo.18923337 (Chen et al., 2026).

1 Introduction

Nitrogen oxides (NO_x= NO + NO₂) are key reactive trace gases and major air pollutants in the troposphere, influencing ozone photochemistry, particulate matter and acid rain formation (Galloway et al., 2004). Accurate and near-real-time estimation of NO_x emissions is crucial for air pollution control and provides a practical approach for estimating co-emitted CO₂ (Reuter et al., 2014, 2019; Miyazaki and Bowman, 2023; Li and Zheng, 2024).

NO₂ tropospheric vertical column densities (TVCDs) retrieved from new generation satellite instruments at unprecedented high spatiotemporal resolution and accuracy (Verhoelst et al., 2021; van Geffen et al., 2022) have accelerated the development of lightweight, mass conservation-based inversion approaches for top-down NO_x emission estimation. Methods such as the cross-sectional flux (CSF) (Reuter et al., 2019; Santaren et al., 2025), integrated mass enhancement (IME) (Santaren et al., 2025) and flux divergence approach (FDA) (Beirle et al., 2019, 2021, 2023) have advanced rapidly due to their low computational cost and minimal latency. These methods rely on NO₂ data and approximate meteorological transport velocities to estimate emissions. The CSF and IME methods first identify sources and then integrate emissions within plume contours based on enhancements above background. They are highly sensitive to wind speed and best suited for strong, distinct point sources with synchronous wind observations (Frankenberg et al., 2016; Koene et al., 2024; Krings et al., 2011; Varon et al., 2018). In contrast, the FDA generates a regional emission map and identifies sources by calculating the divergence of the horizontal flux. It is especially effective for detecting point sources without additional a priori knowledge (Ayazpour et al., 2025; Beirle et al., 2019, 2021, 2023; Koene et al., 2024).

Since its proposal by Beirle et al. (2019), the FDA has been widely applied and refined for NO_x emissions estimation, primarily driven by TROPOspheric Monitoring Instrument (TROPOMI) NO₂ data (Beirle et al., 2021, 2023; De Foy and Schauer, 2022; Koene et al., 2024). It has also been extended to long-lived gases like CH₄ (Veefkind et al., 2023) and CO₂ (Santaren et al., 2025). The FDA theory recommends explicit background subtraction (Koene et al., 2024), even for short-lived species like NO₂ (Cifuentes et al., 2025). However, a well-defined NO₂ background field may not exist (Koene et al., 2024). In practice, the NO₂ background is either empirically removed (Beirle et al., 2019; De Foy and Schauer, 2022; Rey-Pommier et al., 2022, 2023) or the entire column is treated as the enhanced field (Beirle et al., 2021, 2023). An alternative to the FDA is the directional derivative approach (DDA) as proposed by Sun (2022), which replaces the flux divergence term by a directional derivative, i.e., the inner product between horizontal wind and gradient of column amounts. Both the FDA and DDA establish the link between satellite-observed column amounts and emissions by vertically integrating the three-dimensional continuity equation, incorporating reasonable assumptions and approximations (Ayazpour et al., 2025; Koene et al., 2024). The primary difference between the FDA and DDA is the upper limit of vertical column integration (Ayazpour et al., 2025; Lonsdale and Sun, 2023). Integration in DDA is performed from surface to an intermediate altitude that does not need to be explicitly defined. As demonstrated in Ayazpour et al. (2025), the different upper limit used in the DDA ultimately leads to the directional derivative term and a topography term, which in combination replace the flux divergence term in the FDA. By spatially differentiating the column amounts, the DDA implicitly removes the background field that is invariant at the spatial scale of single grid cells. Another resultant benefit of the DDA is that the adopted horizontal wind does not have to approximate the full wind profile but only the lower levels, which leads to higher accuracy and error tolerance.

Two key challenges remain for NO_x emission estimation using the FDA and DDA, including the NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio and the nonlinear NO_x lifetime. When emitted, NO_x is dominated by NO, which is rapidly oxidized by ozone (O₃) to form NO₂. As the plume mixes with ambient air, the balance between NO and NO₂ is determined by the availability of oxidants that oxidizes NO to NO₂ and radiations that photolyze NO₂ to NO. The dominant sink of NO_x is the reaction between NO₂ and OH, which is modulated by the complex interplay between NO_x, O₃, and volatile organic compounds (VOCs) (Laughner and Cohen, 2019). This nonlinear photochemical evolution causes significant variations in the NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio and NO_x lifetime, both over time and distance from the emission source (Krol et al., 2024; Meier et al., 2024). Constrained by the availability of observational data and the rationale for maintaining algorithmic efficiency, both existing FDA and DDA applications have necessarily involved substantial simplifications. Previous studies have accounted for the variability of NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratios across different pixels using auxiliary data (Beirle et al., 2021, 2023; Ayazpour et al., 2025; Cifuentes et al., 2025; Meier et al., 2024), rather than assuming a constant value (Beirle et al., 2019; De Foy and Schauer, 2022; Lonsdale and Sun, 2023; Sun, 2022). Efforts have also been made to consider factors such as NO_x concentration, latitude, season (Beirle et al., 2023; Lange et al., 2022) and distance from the emission source (Krol et al., 2024; Meier et al., 2024) to account for the variability of NO_x lifetime. However, capturing these nonlinear variations (Laughner and Cohen, 2019) remains difficult. Meanwhile, due to the challenges above, most existing studies focus on point sources rather than emissions from large regions.

In this work, we augment the DDA by applying a variable NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio from a global high-resolution chemical transport model and deriving a more realistic nonlinear NO_x lifetime using an improved satellite data-driven piecewise fitting approach. Based on the augmented DDA, we estimate NO_x emissions in China from 2019 to 2024. We evaluate the estimated emissions by comparing regional and point-source emissions with bottom-up inventories and top-down datasets. The paper is structured as follows: Sect. 2 describes the input datasets. Section 3 outlines the methods with a focus on NO_x lifetime fitting, the isolating anthropogenic NO_x from total emissions, and uncertainties. Section 4 presents the validation, along with the distribution and variations of emissions. Section 5 discusses limitations, followed by the conclusions (Sect. 7).

2 Data

2.1 TROPOMI NO₂

The operational offline TROPOMI NO₂ TVCD product (S5P_L2__NO2___HiR2) (van Geffen et al., 2024) from NASA GES DISC (https://daac.gsfc.nasa.gov/datasets/, last access: 17 April 2026) for 2019–2024 is used in this study. TROPOMI is on board of ESA's Sentinel-5 Precursor (S5P) early-afternoon LEO satellite with a high signal-to-noise ratio (Veefkind et al., 2012). It provides global daily coverage with a spatial resolution of 7.0 km × 3.5 km before 6 August 2019 and updated to 5.5 km × 3.5 km thereafter. The data are filtered according to the following criteria: qa > 0.75 and cloud fraction (CF) < 0.3 to remove very cloudy scenes, ice-snow cover scenes and erroneous retrievals (van Geffen et al., 2022; Verhoelst et al., 2021); solar zenith angle (SZA) < 65° to exclude observations with low solar elevation, and viewing zenith angles (VZA) < 56° to minimize unfavorable viewing conditions at the edges of the swath, following Beirle et al. (2023).

2.2 Wind fields data

Wind fields are provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 reanalysis data (https://cds.climate.copernicus.eu/datasets/, last access: 17 April 2026). The 500 m wind fields are interpolated from hourly data on pressure levels, while the 10 m wind fields are directly obtained from hourly data on single levels. Both datasets have a spatial resolution of 0.25° × 0.25°. Details on the selection of wind heights and data processing are provided in Sect. 3.1.

2.3 Chemical data

The chemical data for NO and NO₂ are sourced from the Goddard Earth Observing System composition forecast (GEOS-CF) system (https://gmao.gsfc.nasa.gov/gmao-products/geos-cf/, last access: 16 April 2026) to derive the NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio. GEOS-CF integrates the offline GEOS-Chem chemistry module into the GEOS weather and aerosol modeling system, enabling global near real-time estimates (hindcasts) and 5 d forecasts of atmospheric constituents at a high spatial resolution of 25 km × 25 km (Keller et al., 2021; Knowland et al., 2022). The NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratios are highest near emission sources because freshly emitted NO rapidly consumes local ozone and is subsequently oxidized to NO₂ as the plume mixes with background air (Meier et al., 2024). Following Beirle et al. (2021, 2023), NO₂ columns are converted to NO_x using the near-surface photostationary state, with the NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio from the near-surface model layer with a 1-hour temporal resolution (chm_tavg_1hr_g1440x721_v1) used in this study to better represent near-source conditions and the resulting NO_x gradients (Ayazpour et al., 2025; Sun, 2022).

2.4 Inventories

Two inventories for natural and anthropogenic sources are used. Soil NO_x emissions are obtained from the CAMS global emission inventory (https://ads.atmosphere.copernicus.eu/datasets/cams-global-emission-inventories/, last access: 17 April 2026), which provides monthly estimates from fertilizer/manure application and atmospheric deposition at a 0.5° resolution (Hoesly et al., 2018; Simpson et al., 2014; Yienger and Levy II, 1995). Since the inventory is updated only until 2018 and soil NO_x variability is relatively small compared to total emissions, the 2018 data serve as a proxy for 2019–2024.

Biomass burning and vegetation fires NO_x emissions are derived from the CAMS global biomass burning emissions (GFAS) (https://ads.atmosphere.copernicus.eu/datasets/cams-global-fire-emissions-gfas/, last access: 17 April 2026). GFAS v1.2 provides near-real-time daily averaged fire NO_x fluxes using satellite observations of fire radiative power (FRP) on a global 0.1° × 0.1° grid (Kaiser et al., 2012). The data covers 2019–2024.

Anthropogenic NO_x emissions for validation are sourced from the Multi-resolution Emission Inventory model for Climate and air pollution research (MEIC, v1.4, 2019–2020, http://meicmodel.org.cn, last access: 17 April 2026) and the Emissions Database for Global Atmospheric Research (EDGAR, v8.1, 2019–2022, https://edgar.jrc.ec.europa.eu, last access: 17 April 2026). MEIC provides high-resolution, multi-scale databases of anthropogenic emissions for China at a 0.25° × 0.25° resolution (Geng et al., 2024; Li et al., 2017). EDGAR delivers global coverage at 0.1° × 0.1° resolution, with emissions derived through statistical downscaling of national inventories using high-resolution spatial proxies (Crippa et al., 2024; Solazzo et al., 2021).

2.5 Point source data

The point source emissions derived from DDA are evaluated against the point source catalog provided by Beirle et al. (2023). For point source detection, their catalog uses coal, gas, and oil power plants with capacities ≥ 100 MW from the Global Power Plant Database (Byers et al., 2019), identifying over 1100 NO_x point sources worldwide. Their validation shows good agreement in Germany and the United States, demonstrating the catalog's reliability.

2.6 Nighttime light data

NASA's Black Marble nighttime lights (NTL) product suite, derived from the Visible Infrared Imaging Radiometer Suite (VIIRS) Day/Night Band (DNB) onboard the Suomi National Polar-orbiting Partnership (SNPP), serves as a constraint on anthropogenic emissions (Román et al., 2018; Wang et al., 2021). The Lunar BRDF-Adjusted Nighttime Lights Yearly L3 Global 15 arcsec Linear Lat/Lon Grid product (VNP46A4) provides high spatial resolution at 500 m and is available from January 2012 onward (https://viirsland.gsfc.nasa.gov/Products/NASA/BlackMarble.html, last access: 17 April 2026). This study uses data from 2019 to 2023. Due to the current unavailability of VNP46A4 data for 2024, the 2023 dataset is applied as a substitute.

2.7 National accounts data

National accounts data of annual Gross Domestic Product (GDP) and provincial Gross Regional Product (GRP) are from National Bureau of Statistics of China (https://data.stats.gov.cn, last access: 17 April 2026).

3 Methods

3.1 Framework of the directional derivative approach

The mathematical framework of the DDA is detailed in Sun (2022) and further developed by Ayazpour et al. (2025). Based on satellite-observed NO₂ TVCDs, NO_x emissions E can be estimated as:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}E& ={\displaystyle \frac{\partial \left(f\mathrm{\Omega }\right)}{\partial t}}+f\mathit{u}\cdot \left(\mathrm{\nabla }\mathrm{\Omega }\right)+\mathrm{\Omega }\mathit{u}\cdot \left(\mathrm{\nabla }f\right)+Xf\mathrm{\Omega }{\mathit{u}}_{\mathrm{0}}\\ & \cdot \left(\mathrm{\nabla }{z}_{\mathrm{0}}\right)+{\displaystyle \frac{f\mathrm{\Omega }}{\mathit{\tau }}},\end{array}\end{array}$$

where ∇ denotes the two-dimensional horizontal vector differential operator in Cartesian coordinates, i.e., ∇= ($\partial /\partial x$, $\partial /\partial y)$. No terrain-following coordinate system is used in this formulation. Ω is NO₂ TVCD, f is the NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio, u is the profile-weighted horizontal wind from surface to a conceptual altitude (z₁) which is not explicitly needed and defined, u₀ is surface wind, and z₀ is surface altitude. X and τ are fitting parameters representing the inverse of scale height and NO_x lifetime, respectively. The tendency term ($\partial \left(f\mathrm{\Omega }\right)/\partial t$) becomes negligible when averaging over a month or longer or under the steady-state approximation. It should be noted that this approximation may not hold for a single satellite overpass, potentially leading to significant errors (Koene et al., 2024). Assuming X and τ stay constant over the averaging period, the spatiotemporal averaged emissions 〈E〉 for a given time period and horizontal resolution can be expressed as:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}〈E〉& =〈f\mathit{u}\cdot \left(\mathrm{\nabla }\mathrm{\Omega }\right)〉+〈\mathrm{\Omega }\mathit{u}\cdot \left(\mathrm{\nabla }f\right)〉+X〈f\mathrm{\Omega }{\mathit{u}}_{\mathrm{0}}\cdot \left(\mathrm{\nabla }{z}_{\mathrm{0}}\right)〉\\ & +{\displaystyle \frac{〈f\mathrm{\Omega }〉}{\mathit{\tau }}},\end{array}\end{array}$$

Here, the operator 〈〉 denotes spatiotemporal averaging. The first and second terms of Eq. (2) are referred as the DDf estimator (Ayazpour et al., 2025), reflecting the contributions of advection transport to local NO_x emissions. The third component is the topography correction as emphasized by Ayazpour et al. (2025), it is directly derived from the continuity equation rather than serving merely as an empirical adjustment (Koene et al., 2024). The sum of the first three terms is referred to as DDf_topo estimator. The last term describes the NO_x chemistry approximated by first-order loss.

Integration in DDA is performed from z₀ to an intermediate altitude z₁. While the ideal approach would involve winds at all vertical levels, the full wind profile is typically unavailable. Instead, a single-layer wind, known as the effective wind field, is used to approximate the average state of the full profile. Ayazpour et al. (2025) utilized the effective wind field height using WRF-CMAQ simulations and found that the 500 m wind is most suitable, which aligns with Beirle et al. (2023). Consequently, the horizontal wind u is interpolated to 500 m above the surface, while u₀ directly obtained as the 10 m wind from ERA5.

To apply Eq. (2), the physics-based oversampling approach (Sun et al., 2018) is used to resample the level 2 pixels into 0.025° × 0.025° grid cells with appropriate weighting, then coarsened to 0.05° × 0.05°. Winds in ERA5, the NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio (f) in GEOS-CF data are resampled to match this spatial resolution and temporally aligned with satellite overpasses, allowing for the calculation of all bracketed terms in Eq. (2). Yearly NTL data from SNPP/VIIRS are also resampled at the same spatial resolution. The spatial differentiation in gradient calculation is conducted on 0.05° × 0.05° grid cells with second order central difference. All bracketed terms are averaged at a monthly scale before X and τ are fitted (see in Sect. 3.2) to derive emissions 〈E〉 in mol m⁻² s⁻¹, and the conversion to mass assumes NO_x as NO₂.

3.2 Nonlinear NO_x lifetime fitting

To avoid additional assumptions and external computations (Beirle et al., 2021, 2023), DDA performs a data-driven fitting approach based on monthly fluxes to determine X and τ across grids with negligible emissions (〈E〉≈ 0) (Ayazpour et al., 2025; Lonsdale and Sun, 2023; Sun, 2022), where Eq. (2) can be rewritten as:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}〈f\mathit{u}\cdot \left(\mathrm{\nabla }\mathrm{\Omega }\right)〉& +〈\mathrm{\Omega }\mathit{u}\cdot \left(\mathrm{\nabla }f\right)〉={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}〈f\mathrm{\Omega }{\mathit{u}}_{\mathrm{0}}\cdot \left(\mathrm{\nabla }{z}_{\mathrm{0}}\right)〉\\ & +{\mathit{\beta }}_{\mathrm{2}}〈f\mathrm{\Omega }〉+\mathit{\epsilon },\end{array}\end{array}$$

Here, −β₁ is an estimate of X, −β₂ is an estimate of the inverse of τ, β₀ and ε represent the offset and random error in the predicted variable, which cannot be accounted for by any linear combinations of predictors.

Both the FDA and DDA methods in prior studies have exhibited substantial negative emissions and/or underestimation of total emissions (Ayazpour et al., 2025; Beirle et al., 2019, 2021, 2023; De Foy and Schauer, 2022; Lonsdale and Sun, 2023; Sun, 2022). A single τ fitted in relatively clean region tends to overestimate the chemistry lifetime, then the NO_x flux is underestimated. Therefore, a more realistic representation of NO_x lifetime is critical for accurate emission estimation (Beirle et al., 2023; Krol et al., 2024; Laughner and Cohen, 2019; Meier et al., 2024). The photochemical reactions of NO_x–O₃–VOCs depends on the concentrations of NO_x, O₃, and VOCs, as well as on photolysis and meteorological conditions (Pusede et al., 2015; Sillman et al., 1990; Souri et al., 2023). Even under identical NO_x concentrations and at the same latitude (Beirle et al., 2023; Lange et al., 2022), the lifetime still vary considerably due to the complex interplay of NO_x and VOC chemistry (Laughner and Cohen, 2019).

To derive more realistic lifetime, we fit β₂ across stratified column amount intervals 〈Ω〉 within subregions instead of fitting a single lifetime across the domain to account for the NO_x reactivity nonlinearity. We partition the study domain into four subregions by employing NO_x concentration as the primary criterion due to its dominant control over chemical lifetime (Pusede et al., 2015) with additional consideration of latitude and terrain effects. Figure 1a shows the four subregions. The east (E) and west (W) sections are divided by the Hu line (HL) that extends from Heihe (127.54° E, 50.25° N) to Tengchong (98.50° E, 25.03° N). The E and W sections are further subdivided into north (N) and south (S) subregions by latitudes 35° N for the west and 30° N for the east. The resultant four subregions are denoted as WN, WS, EN, ES thereafter. The HL is a well-known natural geographical boundary that divides China into two contrasting regions in terms of terrain, climate, population density, and economic activity. The eastern part, encompassing most of China's plains and rivers (Fig. 1b), is dominated by the East Asian monsoon, supports 96 % of the population, contributes the majority of national productivity, and faces severe air pollution. In contrast, the western region is characterized by elevated terrain, a westerly climate, and sparse population (Hu, 1935; Zhang et al., 2022). So, the above factors implicitly account for variations in O₃ and its precursors influenced by human activities (Jin et al., 2020; Martin et al., 2004; Souri et al., 2023), differences in natural VOCs emissions from vegetation under varying climatic and geographical conditions (Guenther et al., 1995; Sprengnether et al., 2002; Palmer et al., 2006), and meteorological influences on transport and photochemistry (Duncan et al., 2010; Li et al., 2020; Pusede et al., 2015).

Figure 1 Four subregions divided by the Hu line and (a) average NO₂ TVCDs of 2019 as an example; (b) altitude.

Given the small interannual variability and limited impact on reactive species emission estimation (Ayazpour et al., 2025), β₁ is fitted using the same climatological months, whereas β₂ is fitted for each individual month using NO₂ TVCD percentile bins constructed independently for that month and subsequently averaged over the same months for each subregion during 2019–2024. First, β₁ is fitted in grid cells with rough terrain (0.001 m s${}^{-\mathrm{1}}<〈{\mathit{u}}_{\mathrm{0}}\cdot \left(\mathrm{\nabla }{z}_{\mathrm{0}}\right)〉<$ 0.1 m s⁻¹), where emissions are negligible (DDf < 5 × 10⁻⁹ mol m⁻² s⁻¹). Then, β₂ is fitted in flat terrain (where $〈{\mathit{u}}_{\mathrm{0}}\cdot \left(\mathrm{\nabla }{z}_{\mathrm{0}}\right)〉<$ 0.001 m s⁻¹) without strong NO_x emission sources (DDf_topo < 10⁻¹⁰ mol m⁻² s⁻¹). The fitting is performed across piecewise bins defined by NO₂ TVCD percentiles. TVCD percentiles are binned at 10 % intervals for each month in WN and EN, while in WS and ES, due to narrower lower percentile ranges, bin widths are expanded to 20 % below the 80th percentile and remain 10 % above it to ensure robust fitting performance. For month-bins where fitting fails (p>0.01 or β₂>0), appropriate adjustments are made by merging neighboring bins or revising percentile thresholds. A total of 1393 month-bins are successfully fitted across the four subregions. The number of grid cells per bin ranges from 3 to 9388 (10th percentile: 17; median: 91; mean: 199), indicating that most bins contain sufficient data for robust lifetime fitting, while bins with very small sample sizes are rare and have negligible influence on the overall regression. For each subregion, average TVCD and β₂ per month-bin are calculated from successfully fitted bins. TVCDs can be grouped into intervals as described above, the mean TVCD of each interval is matched to the nearest month-bin and the corresponding β₂ is used to calculate flux.

3.3 Anthropogenic NO_x emissions separation

The DDA quantifies total NO_x emissions, from which anthropogenic contributions to the total emission rates over a certain area need to be isolated by subtracting natural sources. Globally, anthropogenic emissions dominant, while natural sources, such as soil emissions, biomass burning, and lightning, account for approximately 30 %–40 % (Jaeglé et al., 2005; Müller and Stavrakou, 2005). Regional variation of natural sources contributions is substantial, with around 14 % in East Asia (Zhao and Wang, 2009), and 8 % in eastern China (Lin, 2012). Nevertheless, these estimates involve significant uncertainties (Ding et al., 2017; Rey-Pommier et al., 2022), including potential underestimation (Song et al., 2021). Previous studies typically designated regions dominated by anthropogenic emissions as mean NO₂ TVCDs higher than 1.0 × 10¹⁵ molec. cm⁻² (Li and Zheng, 2024; Liu et al., 2016a). However, due to the presence of unexpectedly high or poorly understood natural NO_x emissions over this threshold (Kong et al., 2023; Song et al., 2021), it is not applied to filter natural sources in this study. Since natural NO_x is primarily emitted during the summer and remains low in winter (Liu et al., 2016a; van der A et al., 2006), while anthropogenic NO_x emissions typically peak in winter (Lonsdale and Sun, 2023), we identify natural-source grid cells based on seasonal criterion and further constrain them using NTL data. Grid cells with either the highest averaged NO₂ TVCDs in summer or the lowest values in winter comparing to other seasons, and with NTL < 0.01 nW cm⁻² sr⁻¹, are classified as natural NO_x sources and excluded from the areal integration to get anthropogenic NO_x emission rates. For the remaining grid cells, we subtract soil and biomass burning emissions (from CAMS data) from the satellite-derived emissions to isolate anthropogenic contributions.

3.4 Uncertainty analysis

To quantify uncertainties, DDf is calculated as the mean directional derivative along the zonal/meridional ($\mathit{x}/\mathit{y}$) and diagonal directions ($\mathit{r}/\mathit{s}$) (Li et al., 2026) on a 0.05° × 0.05° grid:

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{rl}{\mathrm{DDf}}_{x/y}& =f\mathit{u}\cdot {\left(\mathrm{\nabla }\mathrm{\Omega }\right)}_{x/y}+\mathrm{\Omega }\mathit{u}\cdot {\left(\mathrm{\nabla }f\right)}_{x/y}=f{u}_x\cdot {\displaystyle \frac{\partial \mathrm{\Omega }}{\partial x}}\\ & +f{u}_y\cdot {\displaystyle \frac{\partial \mathrm{\Omega }}{\partial y}}+\mathrm{\Omega }{u}_x\cdot {\displaystyle \frac{\partial f}{\partial x}}+\mathrm{\Omega }{u}_y\cdot {\displaystyle \frac{\partial f}{\partial y}},\end{array}\text{(5)}& {\displaystyle }\begin{array}{rl}{\mathrm{DDf}}_{r/s}& =f\mathit{u}\cdot {\left(\mathrm{\nabla }\mathrm{\Omega }\right)}_{r/s}+\mathrm{\Omega }\mathit{u}\cdot {\left(\mathrm{\nabla }f\right)}_{r/s}=f{u}_r\cdot {\displaystyle \frac{\partial \mathrm{\Omega }}{\partial r}}\\ & +f{u}_s\cdot {\displaystyle \frac{\partial \mathrm{\Omega }}{\partial s}}+\mathrm{\Omega }{u}_r\cdot {\displaystyle \frac{\partial f}{\partial r}}+\mathrm{\Omega }{u}_s\cdot {\displaystyle \frac{\partial f}{\partial s}},\end{array}\end{array}$$

Standard deviation of the difference between DDf${}_{x/y}$ and DDf${}_{r/s}$ is used to estimate random error σ of DDf:

$$\begin{array}{}\text{(6)}& \mathit{\sigma }=\mathrm{0.5}\cdot \mathrm{SD}({\mathrm{DDf}}_{x/y}-{\mathrm{DDf}}_{r/s}),\end{array}$$

The random errors at daily, monthly, and annual scales for 2019–2024 consistently decrease as the mean satellite data coverage (N) increases (Fig. 2). This observed scaling follows the theoretical relationship $\mathit{\sigma }={\mathit{\sigma }}_{\mathrm{0}}/\sqrt{N}$ (black line in Fig. 2), in agreement with the central limit theorem for independent random errors. Here, σ₀ reflects the precision of a single satellite overpass, calculated at the monthly scale as:

$$\begin{array}{}\text{(7)}& {\mathit{\sigma }}_{\mathrm{0}}=\exp (\langle \log {\mathit{\sigma }}_i\rangle +\mathrm{0.5}\cdot \langle \log N_i\rangle ),\end{array}$$

where σ_i and N_i denote the random error and mean satellite coverage for month i, respectively. The results further demonstrate that temporal averaging effectively reduces random errors in emission quantification.

Figure 2 The relationship between the random error (σ) and mean satellite data coverage (N) at different time scales.

Since the topography correction and chemical loss terms in Eq. (2) are determined through fitting, only the overall Root Mean Square Error (RMSE) can be evaluated, rather than grid-level precision. Therefore, the random error of DDf is employed to characterize the uncertainties in DDA emission estimation, encompassing both satellite retrieval noise and the random error arising from the DDf estimator.

4 Results

4.1 Effects of NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio correction and lifetime fitting

To illustrate the developments in this study, Fig. 3 compares three DDA-based results before and after applying the NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio correction and improved fitting scheme: (1) a constant NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio of 1.32 and monthly single-lifetime fitting (fixed_f and single_τ), corresponding to the original DDA framework by Sun (2022); (2) the variable NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio and monthly single-lifetime fitting (variable_f and single_τ), based on the modification by Ayazpour et al. (2025) and marking the first application of GEOS-CF chemical data in satellite-based emission estimation; (3) a combination of the variable NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio and piecewise fitting with nonlinear NO_x lifetime for each month (variable_f and nonlinear_τ), while the latter represents the major improvement in this study. Using a variable NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio better captures strong NO_x gradients near point sources, improving the accuracy of point source emission estimates. However, this approach still leads to notable underestimation of regional emissions. The nonlinear lifetime fitting more effectively accounts for the balance among local emissions, horizontal transport, and chemical loss, reducing negative emission grids and increasing regional emission estimates. Consequently, the improved fitting scheme minimizes artifacts in mountainous and remote regions compared to earlier results (Ayazpour et al., 2025; Beirle et al., 2023; Lonsdale and Sun, 2023; Sun, 2022).

4.2 NO_x lifetime

The lifetime of NO_x is known as a complicated function of NO_x chemistry regimes (Laughner and Cohen, 2019). At very low NO_x concentration, NO_x lifetime increases with NO_x concentration. As NO_x concentration rises, the lifetime decreases because increased NO enhances the chain reactions involving organic compounds (RH) and HO_x (HO_x= OH + HO₂+ RO₂), accelerating RH oxidation to produce O₃ (“NO_x-limited” regime) and further oxidizing NO. At high NO_x concentration, as NO_x reaches saturation, the reaction between OH and NO₂ becomes much faster than the reaction between OH and RH, dominating the fate of HO_x and slowing O₃ production (“NO_x-suppressed” regime), resulting in the opposite trend in NO_x lifetime (Laughner and Cohen, 2019; Pusede et al., 2015).

Figure 3 Comparison of three DDA-based results, including (a) anthropogenic NO_x emissions estimates in China for 2019 and (b) spatial distribution in Shanxi province for 2019, shown as an example of complex topography.

Figure 4 shows the monthly climatology of NO_x lifetimes from the improved fitting scheme described in Sect. 3.2. For each subregion and each climatological month, the fitted lifetime is shown as bubbles corresponding to bins of NO₂ TVCD. The size of bubbles scales with the mean NO₂ TVCD for the bin. The results clearly demonstrate the nonlinear variability of NO_x lifetime as a function of TVCDs and show significant discrepancies between subregions. The results closely match the theoretical calculated NO_x lifetime versus NO_x concentrations under different VOC reactivities by Laughner and Cohen (2019), capturing the turning points marked by an increase in lifetime at low NO_x concentrations (region I), the subsequent decrease with rising NO_x (region II), and the eventual increase under NO_x-saturated conditions (region III). The range of lifetimes vary from 0.71–26.47 h, and the average values across the subregions range from 3.17–7.85 h. Due to fitting failures in July–August, the number of bins in WN is substantially reduced, requiring broad bin merging and resulting in lifetimes that are likely overestimated and unrepresentative. The results are consistent with the 2–8 h range reported by Lange et al. (2022) as well as the 2 h (low NO₂) to over 7 h (high NO₂) given by Laughner and Cohen (2019).

Figure 4 Monthly NO_x lifetimes for (a) WN, (b) WS, (c) EN, and (d) ES. The bubble size indicates monthly mean NO₂ TVCD (µmol m⁻²) per bin. The τ (h), R² (coefficient of determination), and RMSE (root mean square error, nmol m⁻² s⁻¹) for each subregion represent the month-bin averages. N denotes the count of fitting grids across all months.

The consistent lifetime patterns highlight the dominant role of NO_x concentration in determining τ. However, even at comparable NO_x levels and over the same periods, τ exhibits subregional variations driven by distinct ambient conditions (e.g., O₃ and VOCs concentrations, meteorological parameters) and differences in NO_x emission sources. WN is located in the northern China, sees increased NO_x emissions during the colder half of the year due to heating demand. WS is situated on the sparsely populated Tibetan Plateau, has primarily natural NO_x sources, including unexpectedly high NO emissions from lakes (Kong et al., 2023). EN exhibits distinct anthropogenic emissions, with higher NO_x and longer τ in winter due to heating, and lower NO_x in summer, where intense photochemical reactions result in a shorter τ. ES is located in southern China with smaller annual temperature variations, shows less pronounced seasonal discrepancies in τ. Additionally, high natural VOCs emissions during the growing season of vegetables (Cao et al., 2022) in ES contribute to a longer τ (Laughner and Cohen, 2019).

We provide detailed comparisons of the lifetime fitting parameters for the three DDA-based approaches in Table S1 in the Supplement. The variable NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio correction improves the accuracy of source divergence and emission estimates, while the piecewise fitting approach captures nonlinear NO_x chemistry and yields a shorter overall lifetime with lower fitting RMSE. Across the four subregions, NO_x lifetimes without ratio correction and fitting scheme improvement are approximately 2 to 3 times those derived in this study, specifically 1.7 times in ES, 2.8–2.9 times in EN and WN, and 3.3 times in WS (figure omitted), although the values primarily represent the mean state associated with the nonlinear characterization of lifetime. These results highlight the great importance of accounting for the variable NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio and nonlinear NO_x lifetime, particularly in clean and heavily polluted regions, while the influence is comparatively less pronounced in moderately polluted areas.

In addition, the scale height ($\mathrm{1}/X$) fitting results for each subregion are presented in Fig. S1. Scale height shows clear seasonal and regional variability, generally higher in the cleaner and more remote regions (WN and WS) than in the more polluted regions (EN and ES). Its variability is primarily controlled by boundary layer mixing, as X links surface concentration to column density, reflecting the effective vertical mixing depth (Lonsdale and Sun, 2023; Sun, 2022).

4.3 Regional NO_x emissions

4.3.1 Comparison with inventories

The spatial distribution of national anthropogenic NO_x emissions based on DDA in China for 2019 is shown in Fig. 5a, with the corresponding total emissions provided in Fig. S2. The identified hotspots align closely with areas of intensive human activity, covering both urban and rural regions as well as transportation routes. Traces of faint anthropogenic emissions remain clearly visible even in the remote, sparsely populated Northwest China (Fig. S3).

The derived anthropogenic emissions of China from 2019 to 2024 are compared with: (1) two bottom-up emission inventories (MEIC and EDGAR), and (2) two top-down datasets, the Inversed Emission Inventory for Chinese Air Quality (CAQIEI) (Kong et al., 2024) and results from Li and Zheng (2024). CAQIEI assimilates surface observations using an ensemble Kalman filter (EnKF) and the Nested Air Quality Prediction Modeling System, while the results provided by Li and Zheng (2024) are based on TROPOMI NO₂ combined with GEOS-Chem. Note that these datasets are only available in limited years, while DDA covers all years. National scale comparisons show that DDA agrees well with other inventories and produces slightly lower estimates (Fig. 5b). Based on DDA calculations, the anthropogenic NO_x emissions of China are estimated to be 20.2, 18.5, 19.4, 18.9, 20.7 and 18.8 Tg from 2019 to 2024, with uncertainties of 27 %–30 %. The corresponding total emissions are 29.8, 28.8, 28.7, 28.2, 28.7 and 27.4 Tg, respectively, with uncertainties of 29 %–32 %. These results indicate fluctuating trends in anthropogenic NO_x emissions: a sharp decline in 2020 due to the COVID-19 lockdowns (Lonsdale and Sun, 2023; Miyazaki et al., 2021; Cooper et al., 2022), a modest rebound in 2021–2022, a peak in 2023 (surpassing 2019 levels despite total NO_x remaining lower), and a subsequent drop in 2024 to 7.1 % below 2019 levels. The differences between DDA and EDGAR, MEIC, CAQIEI and Li and Zheng are −8.7 % to 0.8 %, −6.3 % to −3.4 %, −7.2 % to −4.5 %, and −11.8 % to −3.8 %, respectively. Notably, the DDA approach provides long-term emission estimates with low latency, demonstrating a key advantage of this satellite-based lightweight estimation method over conventional bottom-up inventories and computationally intensive top-down data assimilation systems.

Figure 5 National anthropogenic NO_x emissions of China for (a) spatial distribution at 0.05° × 0.05°, with 2019 as an example, and (b) comparisons with previous inventories from 2019–2024.

The primary cause of the lower estimates is the systematic low bias in TROPOMI NO₂ column amounts (Beirle et al., 2023; Liu et al., 2024), most likely due to a priori assumptions and a tropospheric negative bias of approximately 30 % from ground-based validation (van Geffen et al., 2022). Meanwhile, we attribute a larger proportion of total emissions to natural sources compared to the latter two top-down methods, varying between 28.1 %–35.6 %, with an average of 32.1 % (Fig. S4). CAQIEI assigned soil and biomass burning NO_x from CAMS and GFAS as natural sources (approximately 2.0 Tg) and subtracts them from the total NO_x estimates (23.8 Tg in 2019 and 21.4 Tg in 2020, relative to 29.8 and 28.8 Tg in DDA). While Li and Zheng (2024) designated regions dominated by natural emissions as grids with mean NO₂ TVCDs less than 1.0 × 10¹⁵ molec. cm⁻², this threshold is not applied in this study, as described in Sect. 3.3. Based on satellite observations, Kong et al. (2023) found unexpectedly high NO emissions from remote lakes on the Tibetan Plateau, with per-unit-area emissions exceeding those from crop fields in summer. Lin et al. (2024) highlighted a severe underestimation of soil NO_x in the current CAMS inventory, while Opacka et al. (2025) also identified underestimated soil and lightning NO_x. Using natural nitrogen isotopes in precipitation to trace atmospheric NO_x sources, Song et al. (2021) showed that the relative contributions of natural NO_x average 57 ± 13 % in East Asia, a value that had long been underestimated. In this work, we identify natural sources based on seasonal emission patterns and constrain them with NTL data. From 2019 to 2023, the number of grid cells indicative of human activity increased by 48.6 % (Fig. S5). As a result, anthropogenic NO_x emissions in 2023 are 0.5 Tg higher than in 2019, despite total NO_x in 2023 being 1.1 Tg lower.

Table S2 shows NO_x emissions from sensitivity tests of the anthropogenic NO_x filter threshold, and Fig. S6 presents the probability density function (PDF) of NTL to support the selection of the NTL threshold. The results indicate that without applying the NTL constraint, anthropogenic NO_x emissions would be underestimated. The differences among results obtained with different NTL thresholds are relatively small (−4.9 % to 3.6 %). Under the selected threshold of 0.01 nW cm⁻² sr⁻¹, more than 90 % of NTL grids are identified as anthropogenic sources, and this threshold also helps minimize the resampling effect from 500 m to 0.05° in dark regions. Overall, the small differences among different thresholds and the good agreement with other datasets demonstrate the robustness of the anthropogenic NO_x source filtering approach used in this study. Although separating anthropogenic and natural sources in dimly lit grids remains challenging, it is clear that natural emissions in China are indeed underestimated, as noted by Song et al. (2021) and isolating these contributions from total emissions continues to require further work.

For subnational scale (Fig. 6a, b), DDA shows excellent agreement with inventories in the Central, North, and Southwest regions; slightly higher values in the Northeast and Northwest; somewhat lower values in the South; and significantly lower values in East China during 2019–2020. The results also demonstrate a stronger consistency between DDA and MEIC (The Pearson correlation coefficient (r) ranges from 0.86–0.87) compared to EDGAR (r 0.76–0.79). MEIC uses high-resolution localized data, while EDGAR relies on coarser global datasets, resulting in discrepancies in spatial allocation accuracy for China (Crippa et al., 2024; Liu et al., 2016b). Given the good agreement of DDA with inventories in 2019–2021 (Fig. 5b), the discrepancy with EDGAR in 2022 does not undermine the reliability of the estimates in this work, despite the current lack of inventories for comparison.

Figure 6 Comparisons of anthropogenic NO_x emissions between DDA, MEIC and EDGAR for (a) national and subnational scales in 2019, and (b) same as (a) but in 2020, (c) provincial scale in 2019.

At the provincial scale (Fig. 6c), NO_x emissions derived from DDA align well with bottom-up inventories, with about 80 % of provinces falling within ±60 % of MEIC and EDGAR. DDA shows consistently lower emissions in economically developed provinces (mostly in East and South China) and higher in less developed western provinces. The significant difference in Xizang's emissions arises from DDA identifying substantially more human activity than reported in inventories.

4.3.2 Differentiated patterns in provincial NO_x emissions changes

The anthropogenic NO_x emissions map reveals a significant decline in urban areas, particularly in heavily pollutant regions, from 2019 to 2024 (Fig. 7a). However, certain localized areas or point sources, particularly those in the Northeast, Northwest, and South, exhibit noticeable increases. At provincial scale (Fig. 7b), NO_x emissions increase by less than 10 % in Gansu, Hunan, and Inner Mongolia; by 10 %–30 % in Chongqing, Guangxi, Qinghai, Yunnan, and Xinjiang; and by more than 50 % in Xizang. At the subnational scale (Fig. 7c), which aggregates total emissions from provincial groups, NO_x emissions decline by more than 10 % in the North and East, and by less than 10 % in the Northeast, Central, and South, while slight increases under 10 % are observed in the Southwest and Northwest. Year-on-year changes in anthropogenic emissions from 2020 to 2024 are shown in Fig. S7. In other words, NO_x emissions at the provincial scale display a differentiated temporal pattern, with declines in the central and eastern regions but increases in the west.

Figure 7 Changes in anthropogenic NO_x emissions from 2019 to 2024, (a) grid scale at 0.05° × 0.05°, (b) provincial scale, and (c) subnational scale.

Air pollutant emissions evolve in response to urbanization and macroeconomic development, which influence policy regulations and industrial restructuring, forming a feedback loop that further shapes economic growth at both national and city scales (Miyazaki and Bowman, 2023; Wang et al., 2019). Anthropogenic NO_x emissions constrained by nighttime lights reflect the extent of human activity, and their share of total NO_x emissions can serve as a useful indicator of regional urbanization (Wang et al., 2021). Figure 8 presents the relationships among provincial NO_x emission changes, their contributions to total NO_x and GRP. Provincial positive changes in anthropogenic NO_x emissions often coincide with increases in their share of total NO_x (Fig. 8a), particularly in economically less developed regions with relatively low GRP, where anthropogenic sources constitute a small to moderate portion of total NO_x (Fig. 8b) and NO₂ TVCDs are generally low (Fig. S8). In the early stages of economic development, a resource-intensive model driven by fuel consumption dominates, with GRP growth accompanied by rising NO_x emissions at the expense of environmental quality. As the economy progresses, shifts in industry sectors and air quality mitigation measures contribute to a decline in NO_x emissions, with regional economic levels potentially playing a key role in driving these changes (Miyazaki and Bowman, 2023; Wang et al., 2019). The disparities in industrial structure and economic development levels across regions may account for the current differentiated patterns in provincial NO_x emissions.

Figure 8 (a) Provincial changes in anthropogenic NO_x emissions from 2019 to 2024 (ΔNO_x) and their contributions to total NO_x in 2024 ($A/T$), with bubble color indicating changes in the $A/T$ proportion ($\mathrm{\Delta }A/T$), provinces with positive ΔNO_x are labeled; (b) Same as (a), but showing the relationship between $A/T$ and GRP in 2024, with bubble color indicating ΔNO_x.

4.3.3 Consistent downward trend in NO_x emissions across megacities

Six megacities of China with GDP rankings in the top ten and populations exceeding ten million are selected to present the time series of monthly NO_x emissions and changes (Fig. 9). Values in the bottom right corner represent the percentage of months with negative relative changes. The results are in good agreement with the decreased NO_x during the strict COVID-19 lockdowns (e.g., February–May 2020 in Wuhan, February–May 2022 in Shanghai) (Miyazaki et al., 2021; Cooper et al., 2022), and the increase during major holidays (e.g., Spring Festival in January–February, summer holiday in July–August, National Day holiday in October). Seasonal variability of NO_x emissions in megacities is influenced by meteorological conditions. Heating demand leads to markedly higher emissions in autumn and winter than in spring and summer (Miyazaki et al., 2021), with more pronounced variability in northern cities such as Beijing. In contrast, emission fluctuations associated with rising electricity demand during hot seasons (Lange et al., 2022) are more pronounced in the five southern cities. Meanwhile, as China's clean air actions progress, megacities exhibit a consistent downward trend in NO_x emissions, with the proportion of negative changes higher than positive. Shanghai and Beijing have relatively high emissions, but they show a greater share of months with declines from 2020–2024 compared to 2019, reflecting stronger efforts in emission reduction.

Figure 9 Monthly NO_x emissions for 2020–2024 and relative changes from 2019 in China's megacities, with Beijing, Shanghai, Wuhan, Hangzhou, Shenzhen and Chengdu as examples.

4.4 Point source emissions

Point source emissions in DDA (from power plants in this study) are quantified by integrating over a 15 km radius, following the approach and locations of 124 power plants reported by Beirle et al. (2023). Figure 10a presents the estimated NO_x emissions from these plants in China for 2019. From 2019 to 2024, emissions from the 124 plants range from 0.02–2.13 kg s⁻¹, with uncertainties between 4 %–78 %, averaging 16 %. Overall, the plants show an average emission decline of 23 % over this period.

Figure 10 Comparisons of NO_x point source emissions between the DDA and Beirle et al. (2023) data.

The resulting estimates are validated against Beirle et al. (2023) for 2019–2021, as shown in Fig. 10b. The two datasets show good agreement, with an R² of 0.81, both indicating a decline in NO_x emissions from 2019 to 2021. Generally, NO_x emissions from the DDA are slightly lower than those from Beirle et al. (2023), with a slope of 0.83. For comparison, results from DDA test without ratio correction and fitting scheme improvement (fixed_f and single_τ) are also shown in Fig. S9, further illustrating the improvement in point source quantification achieved in this study. It should be noted that point source emissions include all fossil fuel emission sources within the defined radius, leading to a positive bias (Beirle et al., 2023).

5 Discussion

We estimate NO_x emissions using the DDA by applying a variable NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio and deriving a more realistic nonlinear NO_x lifetime through a piecewise fitting approach in each subregion. This scheme combines multiple linear fits based on NO₂ TVCD levels, providing a simple and robust way to represent nonlinear NO_x chemistry in the DDA framework.

Region partitioning mainly depends on NO_x concentration, latitude, and terrain. By using HL as the east-west dividing line, the results implicitly consider variations in O₃ and its precursors driven by human activities, differences in natural VOCs emissions from vegetation across different climates and geographies, and meteorological effects on transport and photochemistry. It is worth noting that HL is used as a boundary for partitioning because it provides a suitable division, where regional divisions improve fitting performance but have a limited impact compared to grading NO_x concentration levels.

The random error of DDf in Eq. (2) described in Sect. 3.4 is used to characterize the uncertainties in DDA emission estimation. Annual anthropogenic emissions uncertainties range from 27 %–30 % at the national scale and 15 %–40 % at the subnational scale. At provincial scale, uncertainties vary more widely, from 8 %–59 %. The average uncertainty for point sources is 16 %. Note that the error does not include the fitting uncertainties related to the topographic correction and chemical loss terms in Eq. (2), indicating that the reported uncertainty is likely underestimated.

Additionally, the anthropogenic NO_x emissions are derived by subtracting natural sources from total emissions, a step that may introduce unknown uncertainties. For point sources, spatial integration that includes all fossil fuel sources within the defined radius may lead to a positive bias.

The DDA adopts a data-driven fitting approach to derive parameters independently, eliminating the need for additional assumptions or external calculations. However, this method requires regions with sufficient terrain diversity (including both rough and flat areas) and an adequate number of observations that satisfy the fitting criteria. The fitted scale heights and chemical lifetimes represent effective subregional parameters derived from grid cells in clean regions, which may lead to overestimation of scale heights in heavily polluted areas and do not fully reflect the sensitivity of NO_x lifetimes to environmental conditions, although the regional averages still capture nonlinear NO_x chemistry. As these are fitting parameters, their physical interpretation should be treated cautiously to avoid over-interpretation (Lonsdale and Sun, 2023). Future observations with higher spatial resolution could improve the representation of scale height and chemical lifetime across heterogeneous regions, while future studies could focus on developing appropriate thresholds for automated global partitioning.

Since satellite-based inversion quantifies total NO_x emissions, isolating natural sources is crucial for accurately estimating anthropogenic emissions. However, the current understanding of natural sources remains limited, with substantial discrepancies and uncertainties across previous studies (Jaeglé et al., 2005; Kong et al., 2023; Li and Zheng, 2024; Lin, 2012; Lin et al., 2024; Müller and Stavrakou, 2005; Song et al., 2021; Zhao and Wang, 2009). We attribute 28.1 %–35.6 % of China's total NO_x emissions to natural sources, indicating that natural emissions in China have been underestimated, as noted by Song et al. (2021).

This work primarily aims to propose a practical and insightful perspective for addressing nonlinear NO_x chemistry in satellite-based emission estimation, rather than focusing on improvements in data quality itself (e.g., air mass factor correction in satellite NO₂ retrievals).

6 Code and data availability

Code relevant to this paper can be found in Sun (2022) at https://doi.org/10.5281/zenodo.7987812 (Sun, 2023). The data can be found at https://doi.org/10.5281/zenodo.18923337 (Chen et al., 2026), Table S3 provides a description of the dataset.

7 Conclusions

We present an improved satellite-based framework for estimating NO_x emissions across China, leveraging the directional derivative approach (DDA) alongside TROPOMI NO₂ observations, ERA5 wind fields, and variable NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratios derived from GEOS-CF. The DDA addresses several challenges in satellite-based emission quantification. It corrects divergence artifacts induced by terrain-driven flows, reduces biases from using single-level wind fields, and suppresses background signal contamination implicitly.

In this work, by incorporating a spatially variable NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio and implementing a data-driven, piecewise fitting strategy, we account for nonlinear NO_x chemistry and improve the estimation of NO_x lifetimes across diverse emission regimes. The NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratio correction improves the accuracy of source divergence and emission estimates, with the major advancement being the piecewise fitting approach, which captures the nonlinear NO_x chemistry. The fitting scheme clusters NO₂ concentration levels within localized regions, reduces fitting errors, mitigates artifacts in mountainous and remote areas, and improves the overall reliability of the estimates. It enables the estimation of both point-source and regional NO_x emissions across China, representing the first application of a lightweight, satellite-driven NO_x emissions estimator in such a large and topographically complex region.

NO_x lifetimes vary from 0.71–26.47 h with NO_x concentrations, and average values across the four subregions range from 3.17–7.85 h, reflecting the nonlinear variability of lifetime as a function of NO₂ TVCDs. Significant discrepancies in regional ambient conditions, such as O₃ and VOCs concentrations, meteorological parameters, and differences in NO_x emission sources, also contribute to this variability. Without ratio correction and fitting improvements, lifetimes across the four subregions are about 2 to 3 times those found in this study, reaching 1.7 times in ES, 2.8–2.9 times in EN and WN, and 3.3 times in WS, highlighting the critical need to consider variable NO${}_x/{\mathrm{NO}}_{\mathrm{2}}$ ratios and nonlinear lifetimes, particularly in clean and heavily polluted regions.

Anthropogenic contributions are isolated by subtracting natural sources from satellite-derived total emissions, with natural NO_x identified using a seasonal criterion and further constrained by NTL data. Estimated anthropogenic NO_x emissions in China are 20.2, 18.5, 19.4, 18.9, 20.7 and 18.8 Tg from 2019 to 2024, with annual uncertainties of 27 %–30 %. The corresponding total emissions are 29.8, 28.8, 28.7, 28.2, 28.7 and 27.4 Tg, respectively, and natural sources account for 28.1 %–35.6 % of the totals in this study.

Spatial and temporal trends show consistent NO_x reductions in megacities, while provincial-level patterns reflect regional differences in urbanization and economic development.

From 2019 to 2024, emissions from the 124 power plants range from 0.02–2.13 kg s⁻¹, with uncertainties between 4 %–78 %, averaging 16 %. Overall, the plants show an average emission decline of 23 % over this period.

Validation against established inventories, including MEIC, EDGAR, and other top-down inversions, demonstrates strong agreement with national level discrepancies ranging from −11.8 % to 0.8 %. At the provincial scale, the DDA shows consistently lower emissions in economically advanced regions and higher in less developed Northeast and western areas. For point sources, our estimates closely match those from previous study, with an R² of 0.81 and a slope of 0.83.

Looking forward, this framework holds promise for global-scale application and for separating natural and anthropogenic NO_x sources. Its low-latency, data-driven nature offers critical value for air quality management, CO₂ co-emission estimation, and international efforts such as the Global Stocktake.