Characterizing lightning NO_x production in the Ebro Valley with meteosat 12 LI and LMA observations

Lightning is one of the main sources of NO_x in the Earth's atmosphere, yet there is considerable variability in NO_x production during thunderstorms. In this study, we combined optical lightning data from the Meteosat Third Generation (MTG) Lightning Imaging (LI) instrument, which provides information on the occurrence rate and physical characteristics of lightning flashes, measurements from an extended Lightning Mapping Array (LMA) consisting of 24 sensors, lightning measurements from the Lightning Location System (LLS) of the Meteorological Service of Catalonia (MSC), and TROPOspheric Monitoring Instrument (TROPOMI) cloud and NO₂ research products to investigate lightning NO_x emissions in the Ebro Valley (Spain). WRF-Chem simulations were used to calculate the necessary air mass factors and to subtract the NO_x not produced by lightning from the TROPOMI measurements, allowing a more accurate estimation of lightning-induced NO_x. We found positive relationships between lightning-produced NO_x per flash and several lightning properties, including flash channel length, radiance, footprint, duration and the length of the optical continuous signal detected from space. These results highlight the importance of combining detailed lightning observations with atmospheric modeling to quantify NO_x emissions from thunderstorms.

1 Introduction

Lightning produces nitrogen oxides (NO_x= NO + NO₂) by following the Zeldovich mechanism (Zel'dovich, 1946). Lightning accounts for about 10 % of total emissions of NO_x (Schumann and Huntrieser, 2007) in the atmosphere, producing between 3 and 7 Tg N in NO_x per year (Pickering et al., 2026), or between 150 and 350 mol NO_x per flash. These emissions are even more important in the tropics, accounting for up to 20 % of the total NO_x production (Schumann and Huntrieser, 2007; Pickering et al., 2026, and references therein). Lightning-produced NO_x (LNO_x) is about 6 times more efficient in driving ozone production than anthropogenic NO_x emissions, producing about 100 molecules of ozone per molecule of lightning-emitted NO_x molecule (Schumann and Huntrieser, 2007; Dahlmann et al., 2011). In turn, LNO_x affects the oxidizing capacity of the atmosphere and the tropospheric budget of carbon monoxide and methane (Wu et al., 2007; Murray et al., 2012; Gordillo-Vázquez et al., 2019; Pérez-Invernón et al., 2025a).

The LNO_x Production Efficiency (PE), defined as the total production of NO_x mol or molecules per lightning flash, exhibits substantial variability across different measurements, likely reflecting both measurement uncertainties and influences of lightning flash characteristics. Precise determination of the variability of LNO_x PE is essential to determine the role of lightning in the chemical composition of the atmosphere (Price et al., 1997; Allen and Pickering, 2002; Tost et al., 2007; Murray et al., 2012; Jöckel et al., 2016; Gordillo-Vázquez et al., 2019; Luhar et al., 2021; Pérez-Invernón et al., 2023a, 2025a). Airborne campaigns have significantly contributed to set the current estimations of LNO_x PE, such as The Deep Convective Clouds and Chemistry (DC3) field campaign (Barth et al., 2015) in the USA, the TROCCINOX campaign in the tropics (Huntrieser et al., 2008) and the EULINOX campaign in Europe (Huntrieser et al., 2002; Fehr et al., 2004). Combination of airborne campaign with ground-based instruments suggest that there is a proportional relationship between LNO_x PE and the flash channel lengths (Stith et al., 1999; Schumann and Huntrieser, 2007; Huntrieser et al., 2008; Pickering et al., 2024, 2026). In turn, laboratory and modeling studies have reported a positive relationship between lightning flash energy and flash channel length with LNO_x PE (Stark et al., 1996; Wang et al., 1998; Ripoll et al., 2014; Jenkins et al., 2021; Pérez-Invernón et al., 2025b). According to estimates reported in the literature, the NO_x production per unit discharge of energy ranges from 1.4×10¹⁶ to $\mathrm{30}\times {\mathrm{10}}^{\mathrm{16}}\mathrm{molec}.{\mathrm{J}}^{-\mathrm{1}}$, while the NO_x production per unit flash channel length spans from 0.5×10²¹ to $\mathrm{10}\times {\mathrm{10}}^{\mathrm{21}}\mathrm{molec}.{\mathrm{m}}^{-\mathrm{1}}$ (Schumann and Huntrieser, 2007). Estimations of LNO_x PE from triggered lightning reported by Rahman et al. (2007) found a production of LNO_x per flash channel length in three measured flashes, with values ranging between 2.0×10²² and $\mathrm{2.4}\times {\mathrm{10}}^{\mathrm{22}}\mathrm{molec}.{\mathrm{m}}^{-\mathrm{1}}$. According to Rahman et al. (2007), slow discharge processes occurring on time scales ranging from milliseconds to hundreds of milliseconds, such as continuing currents and other steady currents lasting for tens to hundred milliseconds, account for approximately 80 % of the total NO_x production per lightning flash. However, the possible relationship between continuing current and LNO_x production has never been confirmed using systematic measurements in storms, mostly because of difficulties in monitoring the occurrence of lightning with continuing currents over large areas. Recent efforts have been made to detect the presence of continuing current in lightning from space (Bitzer, 2017; Fairman and Bitzer, 2022). However, space-based optical sensors have a low detection efficiency for continuing currents occurring below clouds (Ding et al., 2024; Roncancio et al., 2025; Wemhoner et al., 2025). The combination of space-based measurements with ground radio detections could serve to monitor and characterize the continuing currents over large areas (Bozóki et al., 2025; Camino-Faillace et al., 2026). Recent results have already shown that global Extremely Low Frequency (ELF) magnetic field signals, which are strongly influenced by the presence of continuing current in lightning, can serve as a proxy for estimating inter-annual variations in LNO_x (Wang et al., 2025).

Nadir-viewing satellites have provided significant advances in the estimations of LNO_x PE (e.g., Pickering et al., 2016; Laughner and Cohen, 2017; Bucsela et al., 2019; Allen et al., 2019; Lapierre et al., 2020; Allen et al., 2021a; Zhang et al., 2022; Pérez-Invernón et al., 2022). In particular, systematic retrievals of LNO_x from the Ozone Monitoring Instrument (OMI) by Bucsela et al. (2019); Allen et al. (2019) and from the Sentinel-5P TROPOspheric Monitoring Instrument (TROPOMI) by Allen et al. (2021a) have reported an inverse relationship between lightning flash rates and the LNO_x PE. Bucsela et al. (2021) presented a revised evaluation of TROPOMI-based LNO_x PE estimates by incorporating lightning observations from the Geostationary Lightning Mapper (GLM) and a new set of atmospheric chemistry simulations to account for the influence of background-NO_x not produced by fresh lightning on the retrievals. This study reported a weaker relationship than that found by Bucsela et al. (2019), particularly for weak to moderately active thunderstorms (defined as fewer than 3000 flashes per hour and degree).

Using a combination of GLM lightning observations and TROPOMI measurements, Allen et al. (2021b) and Allen et al. (2021a) reported a weak positive relationship between LNO_x PE and lightning flash optical energy. In addition, Pérez-Invernón et al. (2023b) and Pickering et al. (2024) estimated LNO_x PE per flash by combining Lightning Mapping Array (LMA) observations with satellite-based and aircraft-based NO_x measurements, respectively. These studies found that thunderstorms with higher lightning flash rates tend to produce shorter lightning channel lengths and lower LNO_x PE per flash, consistent with earlier findings. This behavior can be reconciled by the inverse relationship between lightning channel length and flash occurrence frequency within storms (Bruning and MacGorman, 2013; Bruning and Thomas, 2015).

More recently, Seiler et al. (2025) assessed the uncertainty of satellite-based LNO_x PE estimates by comparison with aircraft measurements, reporting differences between satellite-retrieved and aircraft-retrieved lightning-NO_x columns ranging from approximately 20 % to 50 %.

The Ebro Valley, located in the North of Spain, is an area of high lightning activity. Airflows over the Ebro Valley are shaped by the Mediterranean Sea and Atlantic Ocean, the Pyrenean mountain range, cold fronts traversing Europe, and a thermal low over the Iberian Peninsula (Pineda et al., 2010). Pérez-Invernón et al. (2022) provided the first TROPOMI-based LNO_x PE estimations over the Ebro Valley using different lightning location systems and TROPOMI products. Pérez-Invernón et al. (2023b) combined lightning measurements from the Ebro LMA (eLMA) with TROPOMI cloud and NO₂ products, as well as Earth Networks Lightning Total Network (ENTLN) lightning measurements, to report a positive relationship between lightning flash channel lengths and LNO_x PE. Pérez-Invernón et al. (2023b) estimated the background-NO_x not produced by lightning by using TROPOMI-derived tropospheric NO_x over TROPOMI pixels without lightning in the previous 5 h. In this study, to reduce uncertainties in estimates of LNO_x PE, we investigate relationships between lightning flashes and LNO_x PE in the Ebro Valley to contribute to reduce current uncertainties in the measured variability of LNO_x PE. We use optical measurements of lightning provided by the Meteosat Third Generation (MTG) Lightning Imager (LI) onboard Meteosat-12 satellite (Dobber and Grandell, 2014; Holmlund et al., 2021) and Very High Frequency (VHF) Lightning Mapping Arrays (LMAs) (Rison et al., 1999) measurements of lightning from the recently upgraded extended Catalonia LMA system (XCALMA) (van der Velde and Montanyà, 2013; Rodriguez et al., 2025) together with the Deutsches Zentrum für Luft- und Raumfahrt (DLR) TROPOMI operational cloud and TROPOMI research NO₂ products, also known as TROP-DLR (Liu et al., 2021; Loyola et al., 2018). Meteosat 12 LI lightning space-based measurements include geometric and temporal metrics of lightning that are influenced by cloud cover, such as the area illuminated by the discharge, named flash footprint, the detected total radiance of the flash, the optical duration and the temporal evolution of the optical signal detected from space. In turn, XCALMA provides estimations of the total flash channel length, volume and duration of flashes from Very High Frequency (VHF) sources without the effect of cloud distortion. Additionally, we employ the Weather Research and Forecasting (WRF) coupled with chemistry (WRF-Chem) (Grell et al., 2005), version 4.3., to estimate the case-based background-NO_x not produced by lightning and to extract the average LNO₂ and LNO_x vertical profiles that are needed to calculate the air mass factor with simulations nudged towards lightning data from Meteosat 12 LI. The WRF-Chem model has demonstrated effectiveness in calculating UV-visible air mass factors, which are crucial for deriving vertical column densities of NO_x from nadir-viewing satellite measurements of NO₂ slant column densities (Laughner and Cohen, 2017).

This paper is organized as follows: Sect. 2 describes the data and methodology; Sect. 3 presents the main results, discusses their implications and limitations; and Sect. 4 concludes with final remarks and future research directions.

2 Data sets and methods

This Section describes the data and methodology. Section 2.1 describes the TROPOMI operational and research cloud and NO₂ products. The lightning measurements are described in Sect. 2.2. Section 2.3 presents the case study selection. Section 2.4 describes the air mass factor and bakcground-NO_x calculations based on WRF-Chem model simulations. Finally, Sect. 2.5 presents the method to calculate the LNO_x PE. The complete methodology is summarized in Fig. 1.

Figure 1 Overview of the methodological framework used to investigate lightning-produced NO_x using satellite observations, ground-based lightning measurements, and numerical modeling.

2.1 TROP-DLR Operational and Research Products

TROPOMI, operating since 13 October 2017 from a low Earth polar orbit onboard the European Space Agency Sentinel-5 Precursor satellite, is a passive imaging spectrometer with eight spectral bands from ultraviolet (UV) to shortwave infrared (SWIR). It provides daily global measurements of cloud properties and several trace gases (Veefkind et al., 2012), with overpasses at approximately 13:30 local solar time.

In this study, we used the TROP-DLR operational cloud and NO₂ research products (Liu et al., 2021; Loyola et al., 2018), which were recently applied by Pérez-Invernón et al. (2022, 2023b) for LNO_x PE estimates. The TROP-DLR operational cloud and NO₂ research products provide the slant column density (SCD) of NO₂, the associated SCD error, the quality assurance (QA) value, the stratospheric vertical column density (VCD) of NO₂, the stratospheric air mass factor (AMF) of NO₂, the cloud fraction (CF), the cloud optical thickness (COT), and the optical centroid pressure (OCP), with a horizontal resolution at nadir of 3.5 km×5.5 km. The TROP-DLR NO₂ research product is better suited for estimating LNO_x than the operational TROPOMI NO₂ product, as it provides more reliable cloud properties and NO₂ retrievals over bright surfaces, such as thunderclouds. Cloud properties are derived using the OCRA/ROCINN algorithms within the Clouds-As-Layers (CAL) model (Loyola et al., 2018), which treats clouds as optically uniform layers and employs a more realistic cloud scattering approach than the Cloud as Reflecting Boundaries (CRB) model implemented in the operational TROPOMI NO₂ product (Lindfors et al., 2018; Pérez-Invernón et al., 2022; Van Geffen et al., 2022).

The TROP-DLR NO₂ research product further applies the Directionally dependent STRatospheric Estimation Algorithm from Mainz (DSTREAM) to separate the contributions of the troposphere and stratosphere to the NO₂ column density (Liu et al., 2021), without relying on atmospheric model inputs. Following Pérez-Invernón et al. (2022, 2023b), we selected TROPOMI pixels with SCD NO₂ errors below $\mathrm{2}\times {\mathrm{10}}^{\mathrm{15}}\mathrm{molec}.{\mathrm{cm}}^{-\mathrm{2}}$. Pixels with deep convection were defined as those with an effective cloud fraction greater than 0.95 and an OCP below 534 hPa (Pérez-Invernón et al., 2022, 2023b).

2.2 Lightning measurements

The MTG-LI onboard Meteosat-12 satellite (Dobber and Grandell, 2014; Holmlund et al., 2021) was launched on 13 December 2022. Pre-operational Meteosat 12 LI (previously Meteosat Third Generation, MTG) L2 data were publicly released 4 July 2024, and are declared operational on 31 October 2024 at 10:00 UTC (Enno et al., 2025). Meteosat 12 LI detects continuously optical signal emitted by lightning from geostationary orbit focused around 0° longitude, providing near-continuous views of the Earth's full disk, including Europe, Africa, the Atlantic Ocean and parts of South America (Gedik et al., 2025; Bližňák and Sokol, 2026). The CMOS detector of Meteosat 12 LI detects light at 1 kHz sampling frequency within a 1.9 nm wide band centered on the atomic oxygen triplet at 777 nm, with a 4.5 km resolution at sub-satellite point (0° latitude, 0° longitude) (Enno et al., 2025). Meteosat 12 LI L2 data include lightning detections at three hierarchical levels: events, groups, and flashes. Events are pixels illumined by lightning during one frame. Events are clustered into a group if they are within one LI frame and spatially contiguous. Finally, flashes are formed by groups separated by less than 330 ms and within about 16 km of each other. Relevant physical information about flashes and groups provided by Meteosat 12 LI L2 data is used in this study, such as flash and groups id, flash and group time, flash duration, flash and group coordinates, total flash radiance, number of groups in a flash, and flash footprint. It is important to note that flash and group properties, such as radiance, duration, and footprint, are influenced not only by the intrinsic characteristics of the lightning but also by the cloud structure (Luque et al., 2020). Apart from the flash and group explicitly included in the Meteosat 12 LI L2 data, we have produced two new variables to account for the time continuity of the optical signal of flashes detected by Meteosat 12 LI. These variables are the total number of time-contiguous groups per flash and the ratio of flashes containing at least 5 ms of time-contiguous optical signal to the total number of flashes reported by Meteosat 12 LI in each flash. The reason for defining these new variables is that the duration of flashes reported by Meteosat 12 LI does not contain information about the total time in which the discharge is optically active. It is important to note that time contiguous optical signal of flashes reported by Meteosat 12 LI could be a proxy for continuing current in CG, but also for in-cloud lightning activity. Bitzer (2017) classified lightning flashes reported by the Lightning Imaging Sensor (LIS) as CG flashes with continuing current using a criterion of at least 10 ms of time-continuous signal. However, the GLM, in geostarionary orbit as Meteosat 12 LI, has a low detection efficiency for the optical signal emitted by continuing currents (Montanyà et al., 2021; Ding et al., 2024; Roncancio et al., 2025; Camino-Faillace et al., 2026) and is more sensitive to in-cloud activity associated with CG lightning (Fairman and Bitzer, 2022; Wemhoner et al., 2025). Therefore, we cannot use Meteosat 12 LI solely to distinguish the LNO_x produced by continuing currents or by time contiguous in-cloud activity to determine if continuing currents lasting more than 5 ms produce more LNO_x in CG strokes than shorter CG strokes, as proposed by Rahman et al. (2007). However, we can analyze possible relationships between LNO_x and the time continuity of the optical signal of lightning detected by Meteosat 12 LI to better understand the role of continuing currents and in-cloud continuous activity in the emissions of LNO_x. The total number of time-contiguous groups per flash and the ratio of flashes containing at least 5 ms of time-contiguous optical signal to the total number of flashes are determined by using Meteosat 12 LI L2 group data. A time-contiguous group is defined as a set of groups within the same flash where at least one group occurs with a time difference of less than 1.1 ms.

The Ebro Lightning Mapping Array (eLMA) (van der Velde and Montanyà, 2013; Montanyà et al., 2014) is a VHF time-of-arrival network (Rison et al., 1999) that was initially established in the Ebro Valley in 2011 and has been operational since then. Between 2023 and 2025, the eLMA network expanded up to 30 stations. The now so-called XCALMA (eXtended CAtalonia Lightning Mapping Array) extends the eLMA coverage mostly to the northern and eastern regions of Catalonia, with 15 new rev. 5 stations deployed by the Meteorological Service of Catalonia (Rodriguez et al., 2025). Furthermore, the initial 15 eLMA stations were enhanced to the LMA rev. 5 version (Van Der Velde et al., 2023; Montanya et al., 2024; van der Velde et al., 2024). In this study, we use lightning data reported by XCALMA between 4 July 2024 and 4 May 2025, based on measurements from the 24 active stations shown in Fig. 2. Lightning sources are located using a minimum of 5 or 6 stations. The 3-dimensional locations of VHF sources in time reported by XCALMA are clustered into flashes in a different way than the previous study Pérez-Invernón et al. (2023b), which used eLMA data. The data are initially processed using a simple filter to remove scattered sources (fewer than 3 sources in a $\mathrm{7}\times \mathrm{7}\times \mathrm{2}$ km grid cell within a 100 ms interval). A DBSCAN-based flash clustering method (Fuchs et al., 2016) was refined (Van Der Velde et al., 2023; Montanya et al., 2024; van der Velde et al., 2024) using two passes in order to perform better for both convective and stratiform flashes. The first pass groups sources at 300 ms and 4 km separation criteria. The second pass identifies convective and stratiform flash regions, and proceeds to recluster the associated sources with longer time and distance criteria to prevent the large stratiform flashes from being broken into multiple small flashes. At greater distances from the LMA network continuous leader activity within a flash (Van Der Velde et al., 2023) is no longer detectable, and the vertical position accuracy is reduced. For this reason, maps of quality metrics were produced to select storms used for this study. Pineda et al. (2025) compared lightning detection between XCALMA and Meteosat 12 LI, revealing that the DE of Meteosat 12 LI and XCALMA are comparable for thunderstorms with low lightning activity (below 15 flashes per minute). However, the DE of Meteosat 12 LI decreases significantly for more active thunderstorms. Conversely, the study also reported an increase in the DE of Meteosat 12 LI for higher intracloud (IC) flash rates.

Figure 2 Left panel: 98th percentile vertical density score (vdensity98) averaged over all the flashes analyzed in this study and contained in each 0.1°×0.1°. Right panel: Masked data showing only pixels with vdensity98 ≥15, showing the region with the best quality XCALMA data. Locations of the XCALMA sensors are included (red triangles).

The variables of each flash extracted from XCALMA for this study include the median flash latitude and longitude, flash start time, flash channel length, flash volume, flash duration, and the 98th percentile vertical density score (vdensity98). The vdensity is calculated using the standard deviation of altitude in 10 ms time steps. It clusters horizontally within every 10 ms, and for each cluster, the metric is defined as the square root of the number of sources in the cluster divided by the variance of their altitudes. This metric correlates inversely with the horizontal and, in particular, the vertical dispersion of points around the true location. Subsequently, vdensity98 is calculated per flash, indicating the quality of the best part of the flash. The vdensity98 ranges approximately up to 100, with higher values indicating better quality. In Fig. 2, we show the 98th percentile vertical density score averaged over all the flashes analyzed in this study and contained in each 0.1°×0.1° latitude and longitude grid boxes.

We incorporate CG lightning data provided by the Lightning Location System (LLS) operated by the Meteorological Service of Catalonia (SMC) (San Segundo et al., 2020). In particular, we calculate the ratio of CG flashes detected by the SMS-LLS to the total number of flashes reported by XCALMA. The SMC-LLS is composed of four LS8000 and three LS-7002 stations (Vaisala Inc.) covering the region of Catalonia, including the Ebro Valley. CG strokes are detected by a Low Frequency (LF) sensor and located using a combination of the Time-of-Arrival/Magnetic Direction Finding techniques (Cummins et al., 1998; Cummins and Murphy, 2009). Like other Vaisala's LLSs based on this technology, the flash algorithm implemented in the SMC-LLS relies on the criteria of the U. S. National Lightning Detection Network (NLDN). Throughout the years of operation, the SMC-LLS performance has been experimentally evaluated in successive campaigns by means of electromagnetic field measurements and video recordings of natural lightning (Montanyà et al., 2006; Pineda and Montanyà, 2009; Montanyà et al., 2012), with a CG flash detection efficiency around 90 % and a location accuracy below 1 km. The method selected for grouping CG strokes into flashes significantly impacts the annual count of lightning flashes. Grouping algorithms primarily depend on two key characteristics of a lightning flash: its spatial extent and total duration. Since the mid-1990s upgrade of the U.S. NLDN, other LLSs based on the same technology have adopted the “NLDN” grouping criteria (Cummins et al., 1998). Under these criteria, a subsequent stroke is grouped with the first return stroke to form a flash if the stroke occurs within 1 s of the first return stroke, its location is within 10 km of the first return stroke, and the time interval between successive strokes is less than or equal to 500 ms. These criteria have remained largely unchanged through subsequent processor upgrades and changes in NLDN ownership (Cummins and Murphy, 2009; Schulz et al., 2005; Murphy and Nag, 2015; Schulz et al., 2016). More recently, the International Standard IEC 62858 (International Electrotechnical Commission, 2019) has also adopted this flash algorithm.

2.3 Case study selection

The case study selection is done by the following method. First, we identify candidate dates by filtering days between 4 July 2024, and 4 May 2025, with at least 650 flashes reported by Meteosat 12 LI within the latitude range of 39.5 to 42.9° N and the longitude range of 1.5° W to 4.5° E, and between 07:00 and 14:00 UTC. This spatial region is chosen as it encompasses most of the area where the LMA vdensity98 exceeds 15 (see Fig. 2). The 07:00–14:00 UTC time window is selected because the TROPOMI overpass occurs approximately between 12:00 and 14:00 UTC in all cases. This yields 35 candidate dates, as shown in Table S1 in the Supplement.

After this preselection of dates, we analyze the spatial distribution of lightning activity and discard dates where no electrical activity cores are centered in the area where the LMA vdensity98 exceeds 15. Twelve dates are excluded based on this criterion, as detailed in Table S1. We further analyze the remaining dates and exclude those with fewer than 650 flashes reported by Meteosat 12 LI in the 5 h preceding the TROPOMI overpass, or where the average flash age exceeds 3.5 h, thereby excluding an additional 4 dates (see Table S1). The age of each flash defined as the hours between the flash occurrence and the overpass of TROPOMI (Allen et al., 2021a). Finally, we exclude the date 17 October 2024, as TROPOMI NO₂ observations indicated that the LNO_x was advected outside the LMA-covered area, where most lightning occurred. This results in the selection of the 18 dates shown in Table S1. Finally, we define the region of interest (ROI) for each date based on the more active thunderstorm within the latitude range of 39.5 to 42.9° N and the longitude range of 1.5° W to 4.5° E from the spatial distribution of Meteosat 12 LI lightning flash activity, CF and OCP.

2.4 Air mass factor and background-NO_x calculations

2.4.1 Simulated LNO_x emissions

The model WRF‐Chem (Grell et al., 2005), version 4.3., is employed to estimate the case-based background-NO_x not produced by lightning and to extract the average LNO₂ and LNO_x vertical profiles needed for calculating the air mass factor of LNO_x in each TROPOMI pixel (${\mathbf{AMF}}_{{\mathrm{LNO}}_x}$).

Simulations were performed using two two-way nested domains. The outer domain (d01) had a horizontal grid spacing of 10 km (300 × 200 grid points) between 32 and 47° N latitude, and between 14.5° W and 16.5° E degrees longitude. The inner domain (d02) had a resolution of 2 km (406 × 366 grid points) between 37 and 42.9° N degrees latitude, and between 3.2° W and 5° E degrees longitude. To improve the vertical resolution of the simulation in both domains, the WRF preprocessing system was used to interpolate the meteorological input data from 34 to 37 vertical levels. This adjustment enhances the ability of the model to resolve atmospheric processes, particularly in the boundary layer and upper atmosphere, leading to more accurate and physically consistent simulations. The parent-child grid ratio was 5:1 in both space and time, and feedback from the inner to the outer domain was enabled. The model time step was 50 s in the outer domain. The base simulations with LNO_x spanned three days, with the day of the TROPOMI overpass included as the final day of the integration, in order to ensure at least two days of model spin-up prior to the analysis period. Initial and boundary meteorological conditions were obtained from analyses of the NCEP GDAS/FNL 0.25 Degree Global Tropospheric Analyses and Forecast Grids (National Centers for Environmental Prediction, 2015) by using the WRF Preprocessing System (Skamarock et al., 2019). Boundary conditions were updated every 6 h. Model outputs were saved every 60 min (d01) and 10 min (d02).

Microphysics was represented using the WRF Single-Moment 6-class scheme (WSM6) scheme (Hong and Lim, 2006). Longwave and shortwave radiation were computed with the RRTM (Mlawer et al., 1997) and Goddard scheme (Chou and Suarez, 2026), respectively, with a radiation time step of 2 min. Surface layer fluxes were calculated using Monin-Obukhov similarity theory (Janjić, 1994), land surface processes were represented with the Noah Land Surface Model (Ek et al., 2003), and the planetary boundary layer was parameterized using the Yonsei University (YSU) scheme (Noh et al., 2003; Hu et al., 2013). Deep convection was parameterized only in the 10 km domain using the Grell 3D ensemble cumulus scheme (Grell and Dévényi, 2002), while convection was explicitly resolved in the 2 km domain. Upper-level damping was applied above 5 km. Grid nudging (FDDA) was applied during the first 2 d of simulation with relaxation coefficients of 3 × 10⁻⁴ s⁻¹ for wind, temperature, and water vapor, excluding the planetary boundary layer.

Gas-phase and aerosol chemistry were simulated using MOZART gas-phase chemistry (Brasseur et al., 1998) coupled with GOCART aerosols (Chin et al., 2000). Photolysis rates and biogenic emissions were updated every 15 min, and gas-phase chemistry time steps were 5 s (d01) and 1 s (d02). Dry deposition of gases and aerosols, wet scavenging, vertical mixing, and convective transport (outer domain only) were enabled. The “wesely” utility from the WRF-Chem preprocessor tool (National Center for Atmospheric Research (NCAR), 2021) was used for calculating dry deposition rates. In turn, the “exo-coldens” utility was employed for calculating O₂ and O₃ column densities. Aerosol-radiation feedback was disabled. Anthropogenic emissions were prescribed hourly from external files interpolated by the “anthro-emiss” utility from EDGAR-HTAP anthropogenic emission files for the year 2010 (Janssens-Maenhout et al., 2015). Biogenic emissions were prepared by the “bio-emiss” usility and were computed online (Guenther et al., 2012). Biomass burning emissions were created by the “fire-emiss” utility from Fire INventory from NCAR (FINN) input (Wiedinmyer et al., 2011) and included with plume-rise parameterization (30 min update frequency). Chemical initial and lateral boundary conditions were produced by the “mozbc” utility from the Whole Atmosphere Community Climate Model (WACCM) (National Center for Atmospheric Research (NCAR), 2024; Buchholz et al., 2019).

Lightning nudging was applied in the nested domain with a 600 s interval (Fierro et al., 2012) by using Meteosat 12 LI lightning flash measurements to ensure that deep convection occurred at the location of observed flashes and to place lightning-NO_x emissions. Lightning NO emissions are parameterized using the lightning nudging scheme in the nested domain, with an emission rate of 311 moles of NO per flash. We have defined a production of 311 moles of NO per flash, based on the widely used values of 1112 moles per flash for CG flashes (Price et al., 1997), 111.2 moles per flash for IC flashes (Price et al., 1997), and a CG-to-total flash ratio of 0.2 (Boccippio et al., 2001). Previous studies have reported substantially different estimates of the amount of NO produced per lightning flash (Schumann and Huntrieser, 2007). In the present study, the specific choice of the absolute LNO_x production per flash is not expected to significantly affect our main conclusions, as our analysis focuses primarily on the relative variability of LNO_x PE among individual thunderstorms, rather than on its absolute magnitude. Therefore, while the assumed emission factor influences the absolute LNO_x concentrations simulated by WRF-Chem, it should have a limited impact on the relative differences in LNO_x PE between the selected storm cases. In another set of simulations, the quantity of moles of NO emitted per flash is scaled using the ratio of the total number of flashes detected by the LMA to the total number of flashes detected by the Meteosat 12 LI within the region of interest. Given that the DE of XCALMA varies across the simulation domain, which extends beyond the region where vdensity98 exceeds 15, we adopted a simplified approach: rather than performing separate simulations with identical PE but differing total flash counts, we maintained a consistent total number of flashes in both simulations and adjusted the PE values accordingly. This method ensures a direct comparison of the impact of detection efficiency on the simulated outcomes. This approach enables the derivation of the background-NO_x and the ${\mathbf{AMF}}_{{\mathrm{LNO}}_x}$, which are subsequently used in the retrieval of LNO_x PE using both the Meteosat 12 LI and LMA flash products. The vertical distribution of NO emissions follows the profiles described in Ott et al. (2010). To isolate the impact of fresh LNO_x, we restarted the base simulations at 07:00 UTC (approximately 5 h before the TROPOMI overpass) and disabled LNO_x emissions while preserving lightning nudging. This approach ensures that meteorological fields remain effectively unchanged.

The choice of a fixed amount of NO_x per flash in the simulations may influence the mean LNO_x per flash (LNO_x PE) results obtained in this study, as the background NO_x estimation will be strongly affected by this quantity. However, the primary goal of this study is to analyze the variability in LNO_x production per flash among different storms. Therefore, any potential bias introduced by assuming a fixed amount of NO_x per flash in the simulations will likely affect all storms in a similar manner, allowing us to explore correlations between this variability and the observed characteristics of the flashes.

Figure 3 WRF-Chem simulation results for 6 July 2024 at 13:50 UTC, focused on the region of interest. (a) Assimilated lightning flashes between 08:50 and 13:50 UTC reported by Meteosat 12 LI. (b) Cloud Top Pressure (CTP). (c) Vertical Column Density of NO_x. (d) Vertical column density of NO_x in the simulation without fresh LNO_x, i.e., with LNO_x deactivated after 07:00 UTC. Meteosat 12 LI lightning data © EUMETSAT.

The “wrf-python” tool (Visualization & Analysis Systems Technologies, 2017) is used to obtain the cloud fraction, the cloud top temperature and the cloud top pressure at each grid point. In this study, the tropopause is defined according to the World Meteorological Organization (WMO) criteria as the lowest atmospheric level at which the vertical temperature lapse rate decreases to 2 °C km⁻¹ or less, provided that this lapse rate is maintained over a minimum vertical extent of 2 km above this level (World Meteorological Organization, 1957). The tropopause height was calculated for each grid point in the model domain, integrating vertical profiles of temperature from the surface to the upper troposphere.

Figure 3 shows the simulation results for 6 July 2024 at 13:50 UTC, focused on the region of interest. The majority of the pixels with lightning activity in the previous 5 h (Fig. 3a) are located in regions with CTP below 534 hPa (Fig. 3b), indicating high convection. However, some lightning activity occurred in regions without significant convective activity at 13:50 UTC. The signal of LNO_x can be clearly seen by comparison of Fig. 3c and d. The simulation results for other cases analyzed in this study can be seen in the Supplement (Figs. S1–S51).

2.4.2 Background-NO_x and advection effects

The pairs of simulations with and without the effect of fresh LNO_x emissions described in Sect. 2.4.1 were used to estimate the background-NO_x not produced by lightning and the effect of advection in LNO_x in each studied case, both for the case of Meteosat 12 LI and LMA flashes. We extracted simulation outputs from the nested domain for the 10 min window closest to the TROPOMI overpass and 5 h prior. We then calculated the VCD of LNO_x for all pixels where the cloud fraction was ≥0.95 and the cloud top pressure was ≤534 hPa, integrating from the surface to the tropopause level. Finally, we compared simulations with and without fresh lightning NO_x (LNO_x) pixel by pixel to determine the ratio of the VCD attributable to fresh LNO_x in pixels with lightning activity, named as R_f, and in pixels without lightning activity, named as R_nf, during the preceding 5 h. The ratios R_f and R_nf averaged over all the pixels with and without lightning will subsequently be used to estimate LNO_x from TROPOMI measurements.

Figure 4 Vertical mixing ratio profiles of NO (a), NO₂ (b), NO_x (c), LNO_x and LNO₂ (d) extracted from WRF-Chem simulations with (solid lines) and without (dashed lines) fresh LNO_x on 6 July 2024 at 13:50 UTC (close to the TROPOMI overpass (background: bck)).

Figure 4 shows the vertical mixing ratio profiles of NO_x and LNO_x from WRF-Chem simulations with and without fresh LNO_x on 6 July 2024 at 13:50 (close to the TROPOMI overpass). These vertical profiles have been obtained by averaging over simulation grid boxes where the cloud fraction was ≥0.95 and the cloud top pressure was ≤534 hPa, and only over the ROI, that is, between 42 and 42.9° N degrees latitude, and between 0° W and 3.5° E degrees longitude. The tropopause height in this case is 16.7 km, and it is located at 99.8 hPa pressure level. The vertical mixing ratio profiles of NO_x and LNO_x from WRF-Chem simulations with and without fresh LNO_x for other dates analyzed in this study can be seen in the Supplement (Figs. S1–S51).

Table 1 Tropopause levels from the WRF-Chem simulation, ratio of the VCD attributable to fresh LNO_x in pixels with lightning activity (R_f), and in pixels without lightning activity (R_nf) for each analyzed date.

The obtained ratios R_f and R_nf, as well as the tropopause level, are shown in Table 1. The pressure at the tropopause is lower during summer months, as expected. The ratio of the tropospheric NO_x produced by lightning is higher in simulation grid boxes with lightning than in boxes without lightning. In turn, the ratio is larger in simulations with LNO_x emissions escalated with the total number of flashes reported by the LMA, which is larger than the total number of flashes detected by Meteosat 12 LI. Additionally, the ratio of tropospheric NO_x produced by lightning is elevated during the summer months, potentially due to an increased planetary boundary layer height, which facilitates greater vertical mixing and accumulation of NO_x in the troposphere.

2.4.3 Air mass factor of LNO_x

The ${\mathbf{AMF}}_{{\mathrm{LNO}}_x}$ is required to calculate the VCD of LNO_x from the SCD of NO₂ (Beirle et al., 2009) provided by the TROP-DLR NO₂ research product. The ${\mathbf{AMF}}_{{\mathrm{LNO}}_x}$ is calculated by following the method proposed by Bucsela et al. (2013) and later used by, for example, Pickering et al. (2016); Laughner and Cohen (2017); Allen et al. (2021a); Pérez-Invernón et al. (2022). We used the algorithm of Bucsela et al. (2013, 2019) to calculate the scattering weights for each of the studied cases using the viewing geometry and the cloud properties for each TROPOMI pixel, based on a look up table generated by Bucsela et al. (2013) with TOMRAD forward vector radiative transfer model. The calculated scattering weights at each TROPOMI pixel are combined with the LNO₂ and LNO_x vertical profiles from the WRF-Chem simulations (see Sect. 2.4.1) to calculate the ${\mathbf{AMF}}_{{\mathrm{LNO}}_x}$ at each TROPOMI pixel. The ${\mathbf{AMF}}_{{\mathrm{LNO}}_x}$ is calculated by integrating Eqs. (2) and (3) of Bucsela et al. (2013) between the surface and the tropopause level over all simulation grid points where the cloud fraction was ≥ 0.95 and the cloud top pressure was ≤534 hPa.

2.5 LNO_x PE estimations

The calculation of the LNO_x PE is based on a modification of the algorithm proposed by Allen et al. (2021a) and later used by Pérez-Invernón et al. (2022) and Pérez-Invernón et al. (2023b). The LNO_x PE in units of mol of NO_x produced per flash is calculated as

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathrm{LNO}}_x\mathrm{PE}& \left({\mathrm{mol}\mathrm{m}}^{-\mathrm{1}}\right)=\left[V_{{\mathrm{tropLNO}}_x}\times A\right]/\\ & \left[{\mathrm{DE}}^{-\mathrm{1}}\times N_{\mathrm{A}}\times \sum _i^F\left(\exp (t_i/\mathit{\tau })\right)\right],\end{array}\end{array}$$

where F is the total number of flashes reported by Meteosat 12 LI or the LMA in the analyzed region and up to 5 h before the overpass of TROPOMI, t_i is the age of each flash defined as the hours between the flash occurrence and the overpass of TROPOMI, A is the area (cm⁻²) of the thunderstorm, defined as the area of TROPOMI pixels satisfying the Deep Convective Constraint (DCC), which are the pixels with a cloud fraction equal or greater than 0.95 (Allen et al., 2021a; Pérez-Invernón et al., 2022, 2023b), and OCP equal or lower than 534 hPa (Pérez-Invernón et al., 2022, 2023b) or with undefined OCP. N_A is the Avogadro's number (molec. mol⁻¹). τ is the lifetime of NO_x in the near-field of convection, in hours. The value of τ is uncertain and can vary between 2 h and 2 d (Pickering et al., 1998; Schumann and Huntrieser, 2007; Beirle et al., 2010; Nault et al., 2017; Allen et al., 2021a) depending on the altitude of LNO_x emissions within the particular thunderstorm and subsequent transport. In this study, we primarily assume τ to be 5 h, as it is a commonly employed value in the literature (Allen et al., 2021a; Pérez-Invernón et al., 2022, 2023b), which also consider τ=3 h, and aligns with the temporal window of flashes included in Eq. (1). To account for the variability in results due to this choice, we also performed calculations using τ=3 h. DE stands for the flash detection efficiency, ranging between 0 and 1. We assume a DE of 1 for the LMA and 0.75 for Meteosat 12 LI (Enno et al., 2025; Gedik et al., 2025; Bližňák and Sokol, 2026). Finally, $V_{{\mathrm{tropLNO}}_x}$ is the VCD NO_x produced by recent lightning (molec. cm⁻²). Allen et al. (2021a); Pérez-Invernón et al. (2022) and Pérez-Invernón et al. (2023b) calculated $V_{{\mathrm{tropLNO}}_x}$ as the median VCD NO_x over TROPOMI pixels satisfying the DCC criteria minus a given value of background derived from measurements or TROPOMI pixels without flashes in the same day or before, or from previous aircraft measurements close to the investigated region. In this study, we calculate $V_{{\mathrm{tropLNO}}_x}$ by using the WRF-Chem simulations as

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{V}_{{\mathrm{tropLNO}}_x}& =\left(N_f\times V_{{\mathrm{tropNO}}_x}^{\mathrm{f}}\times R_f+N_{nf}\times V_{{\mathrm{tropNO}}_x}^{\mathrm{nf}}\times R_{\mathrm{nf}}\right)/\\ & \left(N_f+N_{nf}\right),\end{array}\end{array}$$

where N_f and N_nf are the total number of TROPOMI pixels satisfying the DCC criteria that contain and do not contain Meteosat 12 LI flash groups, respectively. R_f and R_nf are derived from WRF-Chem simulations (Table 1). The advantage of using high-resolution WRF-Chem simulations over other methods is that, by emitting LNO_x at the exact locations where lightning has been reported, the model faithfully reproduces the advection of LNO_x and its interaction with other chemical species. This allows for the determination of the LNO_x ratio relative to the total tropospheric NO_x at each grid point. Although the horizontal resolution of the simulations and TROPOMI pixels is not exactly the same as the horizontal resolution of WRF-Chem simulations, we utilized the coefficients from the table because TROPOMI pixels do not exhibit significant variations from case to case. Additionally, their resolutions are of the same order of magnitude (2 km×2 km vs 3.5 km×5.5 km). $V_{{\mathrm{tropNO}}_x}^{\mathrm{f}}$ and $V_{{\mathrm{tropNO}}_x}^{\mathrm{nf}}$ are the median values of VCD NO_x over pixels satisfying the DCC criteria that contain and do not contain Meteosat 12 LI flash groups, respectively. Given that ${\mathbf{V}}_{{\mathrm{tropNO}}_x}$ is the union of $V_{{\mathrm{tropNO}}_x}^{\mathrm{f}}$ and $V_{{\mathrm{tropNO}}_x}^{\mathrm{nf}}$, it is calculated as

$$\begin{array}{}\text{(3)}& {\mathbf{V}}_{{\mathrm{tropNO}}_x}=\left[{\mathbf{S}}_{{\mathrm{NO}}_{\mathrm{2}}}-\mathrm{avg}\left({\mathbf{V}}_{{\mathrm{stratNO}}_{\mathrm{2}}}\times {\mathbf{AMF}}_{\mathrm{strat}}\right)\right]/{\mathbf{AMF}}_{{\mathrm{LNO}}_x},\end{array}$$

where ${\mathbf{S}}_{{\mathrm{NO}}_{\mathrm{2}}}$ and AMF_strat are the SCD of NO₂ and the stratospheric AMF from the TROP-DLR NO₂ research product, respectively. Only TROPOMI pixels with SCD NO₂ errors below $\mathrm{2}\times {\mathrm{10}}^{\mathrm{19}}\mathrm{molec}.{\mathrm{m}}^{-\mathrm{2}}$ are used (Pérez-Invernón et al., 2022, 2023b). Finally, ${\mathbf{AMF}}_{{\mathrm{LNO}}_x}$ is the AMF of LNO_x derived from WRF-Chem simulations and using the algorithm of Bucsela et al. (2013, 2019) (see Sect. 2.4.3). $\mathrm{avg}\left({\mathbf{V}}_{{\mathrm{stratNO}}_{\mathrm{2}}}\times {\mathbf{AMF}}_{\mathrm{strat}}\right)$ stands for the average of the product ${\mathbf{V}}_{{\mathrm{stratNO}}_{\mathrm{2}}}\times {\mathbf{AMF}}_{\mathrm{strat}}$.

Additionally, we have estimated the averaged LNO_x PE in terms of molec NO_x per meter, per second or per unit of radiance (mW m⁻² sr⁻¹) by modifying Eq. (1) to better account for the metrics of individual flashes and the influence of their ages as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\mathrm{LNO}}_x\mathrm{PE}& \left({\mathrm{molec}.\mathrm{m}}^{-\mathrm{1}}\right)=\left[V_{{\mathrm{tropLNO}}_x}\times A\right]/\\ & \left[{\mathrm{DE}}^{-\mathrm{1}}\times \sum _i^F{M}_i\left(\exp (t_i/\mathit{\tau })\right)\right],\end{array}\end{array}$$

where M_i is the metric representing the flash channel length, the radiance, or the duration of the flash i reported by Meteosat 12 LI or the LMA.

As an alternative approach, we modified Eq. (2) to calculate $V_{{\mathrm{tropLNO}}_x}$ as the difference between the median VCD NO_x over TROPOMI pixels satisfying the DCC criteria and the 10th and the 30th percentiles of the tropospheric VCD NO_x over TROPOMI non-flashing pixels (i.e., pixels without detected lightning activity), as previously proposed by Allen et al. (2021a). This alternative formulation facilitates comparison with the results from Pérez-Invernón et al. (2022) and Pérez-Invernón et al. (2023b) for the same region.

3 Results and discussion

This section presents the key results of the study, along with a detailed discussion of their implications. Section 3.1 presents the vertical column densities of LNO_x derived from TROP-DLR product, alongside the corresponding lightning observations. Section 3.2 presents the LNO_x estimates and analyzes their relationships with lightning characteristics. Finally, Sect. 3.3 discusses the uncertainties and limitations of this study, with a focus on how key input parameters in the LNO_x retrieval process influence the observed relationships between LNO_x and lightning characteristics.

Figure 5 TROP-DLR product of orbit 34875, Meteosat 12 LI and LMA lightning data for the case 6 July 2024. (a) Slant Column Density (SCD) NO₂ provided by the TROP-DLR NO₂ research product. (b) SCD NO₂ minus the stratospheric vertical column density multiplied by the stratospheric air mass factor of NO₂ (AMF_strat) obtained from the TROP-DLR NO₂ research product. (c) and (d) Cloud Fraction (CF) and Optical Centroid Pressure (OCP) from the TROP-DLR NO₂ research product, respectively. (e) VCD NO_x ($V_{{\mathrm{tropNO}}_x}$ from Eq. 3) over pixels satisfying the Deep Convective Constraint (DCC). (f) and (g) Lightning flashes during the 5 h period per TROPOMI pixel before the TROPOMI overpass reported by Meteosat 12 LI and the LMA, respectively. (h) TROPOMI pixels containing at least one lightning group detected by Meteosat 12 LI. Sentinel-5P TROPOMI data © ESA, processed by DLR. Meteosat 12 LI lightning data © EUMETSAT.

3.1 Estimated Vertical column density of LNO_x from TROP-DLR and lightning observations

Figure 5 shows the TROP-DLR products, lightning measurements and calculated $V_{{\mathrm{tropNO}}_x}$ for the case 6 July 2024, while the figures showing the rest of the analyzed dates are provided in Figs. S1–S51. In the majority of the analyzed cases, a distinct enhancement in the SCD NO₂ (Fig. 5a, b) is observed in or near TROPOMI pixels characterized by high lightning activity in the last 5 h, as in the TROPOMI orbits analyzed in other studies (Allen et al., 2021a; Pérez-Invernón et al., 2022, 2023b). However, this enhancement is not evident in all cases, particularly on 6 July, 23 August, 25 August, 2 September, and 31 October 2024. Nevertheless, when the VCD NO_x (Fig. 5e) is calculated using the ${\mathbf{AMF}}_{{\mathrm{LNO}}_x}$, which incorporates cloud properties (Fig. 5c, d), observation geometry, and simulated LNO_x vertical profiles, the enhancement in TROPOMI pixels with or near lightning activity becomes apparent. In turn, a clear similarity between the spatial distribution of the OCP (Fig. 5d) and the lightning activity (Fig. 5f, g) can be seen. The spatial distribution of flashes detected by Meteosat 12 LI and the LMA do not show large discrepancies (Fig. 5f, g), with the small differences observed likely attributable to the parallax correction applied in Meteosat 12 LI (Pineda et al., 2025). In fact, larger differences emerge when comparing the spatial distribution of individual flashes (Fig. 5f, g) and groups detected by Meteosat 12 LI (Fig. 5h). These differences suggest that relying solely on the flash product may introduce inaccuracies in distinguishing between TROPOMI pixels with and without lightning activity, as lightning flashes coordinates are estimated as the radiance-weighted geometrical centroid of the groups detected by Meteosat 12 LI (Dobber and Grandell, 2014; Holmlund et al., 2021; Enno et al., 2025) or by the median coordinates of the sources detected by the LMA. However, each lightning flash may not be fully contained within a single TROPOMI pixel. By using the flash groups reported by Meteosat 12 LI, we can access information about the TROPOMI pixels that actually contain at least a portion of a flash. Note that the VCD NO_x (Fig. 5e) can have negative values, indicating that the average stratospheric column exceeds the local vertical column (Pérez-Invernón et al., 2022). The TROP-DLR products, lightning measurements and calculated $V_{{\mathrm{tropNO}}_x}$ for other cases analyzed in this study can be seen in the Supplement (Figs. S1–S51).

Figure 6 Total number flashes detected by Meteosat 12 LI, LMA and CG flashes detected by LLS-SMC up to 5 h before the overpass of TROPOMI (a), mean age of the flashes detected by Meteosat 12 LI and the LMA (b), total number of valid TROPOMI pixels satisfying the Deep Convective Constraint (DCC), total number of TROPOMI pixels with Meteosat 12 LI or LMA flashes satisfying the DCC (c), and ratio of the total number of CG flashes to all flashes detected by the LMA.

The comparison of WRF-Chem simulation and TROPOMI measurements indicates that the CTP (Fig. 3b) and the OCP (Fig. 5d) are nearly in agreement, with a slight displace of the simulated high clouds westward. The simulated tropospheric VCD NO_x in areas with high convection (Fig. 3c) is consistently displaced westward, and lower in about 1 order of magnitude than the tropospheric VCD NO_x calculated from measurements (Fig. 5e). The obtained disagreement in the maximum simulated and measurements-based tropospheric VCD NO_x can be due to different horizontal resolutions in simulations and measurement, as well as the inherent uncertainties in the WRF-Chem model, such as emission inventories, chemical reaction rates, meteorological inputs, or parameterizations of physical processes.

Figure 6 and Table S2 shows the basic measurements for each analyzed date. The LMA detects 1.3 to 6.5 times more flashes than the Meteosat 12 LI. The ratio of CG flashes reported by the LLS-SMC to total flashes reported by the LMA ranges between 0.01 and 0.39. There are no significant differences in the average age of flashes reported by Meteosat 12 LI and the LMA. The average OCP ranges between 323.07 and 499.34 hPa, which are typical values in thunderstorms (Allen et al., 2021a). The percentage of valid TROPOMI pixels satisfying the DCC that contain Meteosat 12 LI lightning flashes with respect to pixels without flashes ranges between 9 % and 46 %. In turn, the ratio of the total number of TROPOMI pixels satisfying the DCC and containing flashes according to the Meteosat 12 LI to the pixels containing flashes according to the LMA ranges between 0.7 and 1.4, indicating that there are no significant differences in the spatial distribution of flashes detected by both instruments.

Table 2 Analyzed dates, average footprint, average radiance, average duration from Meteosat 12 LI, average total number of time contiguous groups per flash, average flash channel length, average flash volume, average flash duration from LMA, LNO_x PE based on Meteosat 12 LI and LMA lightning flash measurements with τ=5 and τ=3 h.

3.2 Relationships between LNO_x and lightning characteristics

The characteristics of lightning flashes analyzed in this study and the LNO_x PE estimates are shown in Table 2. There are notable variations in the characteristics of lightning flashes detected by Meteosat 12 LI. The average footprint ranges from 4.88 to 14.59 Meteosat 12 LI pixels, with a mean value of approximately 7.71 pixels, indicating spatial variability in the optical extent of detected events. The average radiance spans a broad range from 152.69 to 910.57 mW m⁻² sr⁻¹, reflecting substantial differences in the optical intensity of the flashes. This variability may be attributed to differences in the energy released during lightning events and to differences in the cloud opacity. Regarding the average duration of Meteosat 12 LI detected flashes, values range from 159 to 374 ms, with a mean duration of around 230 ms. The average flash channel length reported by the LMA ranges from 12.33 to 75.48 km, and the average flash volume spans from 8.60 to 64.93 km³ (number of $\mathrm{1}\times \mathrm{1}\times \mathrm{0.25}$ km boxes), indicating substantial variability in the spatial extent of lightning channels. Additionally, the average duration of LMA-detected flashes ranges from 170 to 420 ms, with a mean duration of approximately 295 ms, slightly longer than the average duration provided by Meteosat 12 LI. The obtained average LNO_x PE ± the standard deviation with Meteosat 12 LI flash data and τ=5 h is 37 ± 25 mol NO_x per flash, while it is 31 ± 24 mol NO_x per flash when using the LMA lightning flash data. When setting τ=3 h, we obtain 41 ± 29 mol NO_x and 38 ± 31 per flash when using Meteosat 12 LI and LMA lightning data, respectively. The approximately similar LNO_x PE obtained with both lightning data sources suggests that the WRF-Chem simulations used in the estimation of LNO_x PE help mitigate the impact of the difference in the total number of flashes reported by Meteosat 12 LI and the LMA (see Fig. 6). The averaged LNO_x PE obtained in this study falls at the lower end of the range of previous estimations in the same region, ranging between 42 and 108 mol NO_x per flash (Pérez-Invernón et al., 2022, 2023b). In turn, the obtained LNO_x PE falls at the lower end of the widely consensus value, ranging between 100 and 400 mol NO_x per flash (Schumann and Huntrieser, 2007). However, it is important to mention that the obtained averaged LNO_x PE value is significantly influenced by the production of NO_x per flash assumed in WRF-Chem simulations, as shown in the last four columns of Table 2. Alternatively, we approximated the background-NO_x by considering the 10th and the 30th percentiles of the tropospheric VCD NO_x over TROPOMI non-flashing pixels (i.e., pixels without detected lightning activity) (Allen et al., 2021a). Based on LMA lightning data, we obtained 116 ± 89 and 49 ± 47 mol NO_x per flash by considering the 10th and the 30th percentiles, respectively. These values are in agreement with Pérez-Invernón et al. (2022, 2023b), who employed similar methods to estimate the background-NO_x.

Table 3 Spearman and Pearson correlation coefficients (ρ and r, respectively) and corresponding p-values between $V_{{\mathrm{tropLNO}}_x}$ (as defined in Eq. 2) and various lightning flash characteristics, computed across TROPOMI pixels included in the 18 analyzed dates. These results are based on calculations with τ=5 h.

Table 3 shows the Spearman and Pearson correlation coefficients and corresponding p-values between $V_{{\mathrm{tropLNO}}_x}$ (as defined in Eq. 2) and various lightning flash characteristics, computed across all TROPOMI pixels included in the 18 analyzed dates. Only pixels with values of $V_{{\mathrm{tropLNO}}_x}$ between the 5th and the 95th percentiles are included in the calculation of the correlation coefficients in order to exclude outliers. The lightning characteristics for each TROPOMI pixel were derived by averaging the respective values over all lightning flashes occurring within that pixel. The Spearman correlation coefficients presented in Table 3 reveal measurable, albeit weak, associations between lightning flash characteristics and $V_{{\mathrm{tropLNO}}_x}$. Among the analyzed variables, footprint exhibits the weakest correlation with $V_{{\mathrm{tropLNO}}_x}$, suggesting a less pronounced influence compared to other flash characteristics. Nevertheless, the results indicate that the properties of lightning flashes, such as length, volume, duration, and radiance, play an appreciable role in the production of LNO_x. The data employed to calculate the correlations presented in Table 3 are shown in Fig. S52.

Figure 7 Spearman correlation heatmap of the most important lightning characteristics and LNO_x PE estimates. These results are based on caculations with τ=5 h.

Figure 7 shows a Spearman monotonic correlation heatmap of the most important lightning characteristics and LNO_x PE estimates with τ=5 h. There are strong direct relationships between the average flash channel length (LMA – av. Length), median flash channel length (LMA – med. length) flash volume (LMA – volume) and flash duration (LMA – duration) estimates by the LMA and the average footprint (MTG – footprint), radiance (MTG – rad) and duration (MTG – duration) reported by Meteosat 12 LI. These correlations indicates consistency in the methods used to cluster LMA sources into flashes and Meteosat 12 LI groups into flashes. In turn, the high correlation between durations and geometric variables, such as flash channel length, flash volume and duration, suggests that larger flashes then to have longer durations. There is a significant direct correlation between the flash length and the age of flashes (Flash age), showing that flashes that occur near the sunrise are longer than flashes that occur near the mid-day, which is consistent with previous findings of shorter flashes in more active thunderstorms (Bruning and Thomas, 2015; Chronis et al., 2015). In the same manner, there is an even stronger correlation between the age of the flashes and the footprint, radiance and duration reported by Meteosat 12 LI, indicating that Meteosat 12 LI is capable of reporting a longer and brighter signal of flashes during early morning than during the mid-day, due to the lower DE when the sun is close to the zenith. Conversely, a higher DE during the early morning could allow Meteosat 12 LI to capture more weak flashes that would otherwise be missed when the DE is lower. However, the effect of observing a larger portion of bigger flashes when the DE is higher ultimately dominates. As shown by the LMA measurements reported in Chronis et al. (2015), small flashes are typically less numerous than large flashes early in the day. The total number of time contiguous groups (MTG – Cont. Groups) and the ratio of flashes with at least 5 ms of time contiguous optical signal to the total number of Meteosat 12 LI (MTG – CC (5 ms)/f) flashes is moderately correlated with the flash channel length, the flash volume and the flash duration provided by the LMA, while they are strongly correlated with the footprint, radiance and duration measured by Meteosat 12 LI.

The ratio of CG flashes to total flashes (CG/f) is strongly correlated with the flash channel length, flash volume, duration, footprint and radiance (Fig. 7), indicating that CG flashes are longer in duration, larger and brighter from space than IC flashes. Similar results were obtained by Ringhausen et al. (2021) for GLM lightning measurements. The ratio of CG flashes to total flashes reported is moderately correlated with the total number of time contiguous groups and with the ratio of flashes with at least 5 ms of time contiguous optical signal, indicating that CG flashes reported by Meteosat 12 LI tend to have longer time contiguous optical signals detected by Meteosat 12 LI. This may be due to the presence of continuing currents in CG, or to longer in-cloud stroke activity in CG flashes than in IC flashes.

The ratio of the total number of flashes reported by the LMA (N LMA flashes) to the total flashes reported by Meteosat 12 LI (N MTG flashes) is inversely correlated with the OCP (TROP-DLR OCP) (Fig. 7), indicating that the detection efficiency of Meteosat 12 LI with respect to the LMA is higher when the cloud tops are lower, possibly due to lower cloud opacity in the travel of the light between the emitting source and Meteosat 12 LI. In turn, the ratio of the total number of flashes reported by the LMA to the total flashes reported by Meteosat 12 LI is inversely correlated with the average footprint and with the ratio of CG flashes to total flashes, indicating that the detection efficiency of Meteosat 12 LI with respect to the LMA is higher for brighter flashes, particularly for CG flashes. The positive correlation observed between the total number of flashes detected by the LMA and the ratio of this quantity to the total number of flashes reported by Meteosat 12 LI suggests that the DE of Meteosat 12 LI decreases for more active thunderstorms, in agreement with Pineda et al. (2025).

The analysis of the flash initiation altitude (LMA – start altitude) reveals that lightning initiated at lower altitudes exhibits greater spatial extent, longer duration, and higher radiance as observed by MTG, likely due to the predominance of CG flashes, which tend to be more energetic and prolonged (Fig. 7). Additionally, lower initiation altitudes correlate with a higher number of flashes lasting at least 5 ms and elevated OCP values (indicating lower cloud tops). Conversely, lightning initiated at higher altitudes, though typically smaller in size (Pérez-Invernón et al., 2023b), is more effectively detected by LMA than MTG, suggesting higher DE of LMA than MTG for smaller flashes.

The LNO_x PE is positively correlated with the flash channel length, the flash volume and the flash duration estimates by the LMA and the footprint, radiance, duration and total number of contiguous groups per flash reported by Meteosat 12 LI (Fig. 7), indicating that larger, brighter and longer duration flashes produce more LNO_x. The stronger correlations are obtained when comparing the LNO_x PE based on LMA flash measurements with the average flash channel length, the flash volume and the flash duration, in agreement with previous studies (Pérez-Invernón et al., 2023b). There is a weak positive correlation between the LNO_x PE and the total number of Meteosat 12 LI time contiguous groups, as well as with the ratio of flashes with at least 5 ms of time contiguous optical signal, indicating that time contiguous emissions could contribute to a larger LNO_x PE.

The lack of significant correlations between the LNO_x PE and the average flash age indicates that the assumption of a temporal window of 5 h and a τ=5 h do not introduce a significant bias in the LNO_x PE estimates (Fig. 7). In turn, the lack of significant correlations between the LNO_x PE and the average OCP suggests that the calculation of the ${\mathbf{AMF}}_{{\mathrm{LNO}}_x}$ does not introduce a significant bias (Figure 7). The Spearman monotonic correlation heatmap with τ=3 h is shown in Fig. S53, indicating minor differences with the results obtained by using τ=5 h.

Figure 8 Pearson correlation between the LNO_x PE and the most important lightning characteristics included in this study. The average LNO_x PE (blue points) is obtained as the average between the LNO_x PE estimates by using both Meteosat 12 LI and the LMA for each analyzed case, while each pair of grey bars are the values of the LNO_x PE estimates by using Meteosat 12 LI and the LMA separately. The flash channel length, the flash volume and the flash durations in panels (a), (b) and (c) are the average values calculated from LMA measurements. The footprint, the radiance, the flash durations and the total number of contiguous groups per flash in panels (d), (e), (f) and (g) are the average values calculated from Meteosat 12 LI measurements. The average ratio of CG to total flashes in panel (h) is obtained from LLS-SMC and LMA lightning measurements. Each subplot includes the Pearson correlation coefficient r and its corresponding p-value, as well as a linear fitting of the blue points (red line). These results are based on caculations with τ=5 h.

Figure 8 shows the Pearson linear correlations between the average LNO_x PE and the most important lightning characteristics with τ=5 h, while nearly similar results are shown in Fig. S54 with τ=3 h. The average LNO_x PE is obtained as the average between the LNO_x PE estimates by using both Meteosat 12 LI and the LMA. Figure 8a, b and c shows a positive linear correlation between the LNO_x PE and the average flash channel length, volume and duration estimated from LMA measurements, respectively. In turn, Figure 8d, e and f show a positive but weaker linear correlation between the LNO_x PE and the average footprint, radiance and duration based on the Meteosat 12 LI flash product, respectively. Allen et al. (2021b, a) reported that the correlation between LNO_x PE and the flash energy reported by GLM ranged between 0.18 and 0.57, which is consistent with the 0.3578 correlation reported in this study (Fig. 8e).

There is a positively and moderate correlation between the total number of contiguous Meteosat 12 LI groups per flash and the LNO_x PE Fig. 8d. The correlation between the LNO_x PE and the total number of contiguous groups per flash is lower than the correlation with the total duration measured by the LMA (Fig. 8c, g). Geostationary lightning detection optical systems detect more efficiently in-cloud activity than cloud-to-ground processes, especially during the occurrence of continuing currents in CG strokes (Wemhoner et al., 2025). The lower correlation of the LNO_x PE with the total number of contiguous groups per flash than with the total duration reported by the LMA (0.36 vs 0.68) indicates that there may be lightning activity below clouds not detected by Meteosat 12 LI that contributes to the production of LNO_x. Finally, Fig. 8h indicates a positive linear correlation between LNO_x PE and the ratio of CG to total flashes.

The linear positive correlations between the LNO_x PE and the characteristics of lightning measured by the LMA (Fig. 8a–c) are stronger than the correlations with the lightning metrics reported by Meteosat 12 LI (Fig. 8d–g) especially due to differences in two particular cases, that are the 14 August 2024 and the 4 May 2025. These cases are characterized by medium flashes in terms of flash channel length, volume, duration from the LMA and LNO_x PE, but small and short in terms of the optical signal detected by Meteosat 12 LI. This suggests that the lightning characteristics measured by the LMA are better proxies for the LNO_x PE than Meteosat 12 LI flash properties, potentially due to the effect of cloud cover in the optical signals detected from space.

The obtained LNO_x PE in molecules of NO_x per meter obtained as the LNO_x PE averaged over the 18 analyzed cases divided by the averaged flash channel length with τ=5 h (Table 2) is 0.56 × 10²¹ molec. NO_x m⁻¹, which falls in the lower range of previous estimates based on laboratory and aircraft measurements (Wang et al., 1998; Stith et al., 1999; Huntrieser et al., 2002; Skamarock et al., 2003), with theoretical results based on electrodynamical models ranging between 0.045 and 0.106 × 10²¹ molec. NO_x m⁻¹ (Pérez-Invernón et al., 2025b), and with TROPOMI based estimates in the same region, ranging between 0.8 × 10²¹ and 1.6 × 10²¹ molec. NO_x m⁻¹ (Pérez-Invernón et al., 2023b). In turn, we have estimated the averaged LNO_x PE in terms of molec NO_x per meter from Eq. (4), obtaining 0.48 × 10²¹ molec. NO_x m⁻¹. Additionally, we have estimated the averaged LNO_x PE in terms of molec NO_x per second of discharge by using Eq. (4), yielding 6.53 × 10²⁵ molec. NO_x s⁻¹ of duration reported by Meteosat 12 LI, and 6.02 ×  10²⁵ molec. NO_x s⁻¹ of duration reported by the LMA. Finally, we have calculated the averaged LNO_x PE in terms of molec NO_x per flash radiance seen by Meteosat 12 LI from Eq. (4), obtaining a production of 9.79 × 10²² molec. NO_x per mW m⁻² sr⁻¹.

3.3 Limitations and uncertainties

3.3.1 Uncertainty in the obtained LNO_x PE and lightning observations

The main limitations and sources of uncertainties related with nadir-based estimations of LNO_x discussed by Allen et al. (2021a); Pérez-Invernón et al. (2022) and Seiler et al. (2025), ranging from approximately 20 % to 50 %, applies to this study. However, in our present study some of them have been partially mitigated by introducing modifications to the methods developed to estimate the background-NO_x, the ${\mathbf{AMF}}_{{\mathrm{LNO}}_x}$, and the total number of flashes.

In this study, the uncertainty in the LNO_x PE estimate was evaluated by addressing key sources of error identified by Pérez-Invernón et al. (2022). Specifically, we targeted the 29 % uncertainty attributed to the background-NO_x estimate, the 30 % uncertainty attributed to the parameterization of lightning activity in atmospheric models, which affected the a priori profiles of LNO_x and LNO₂, scattering weights, and DCC definitions, and the 7 % uncertainty due to the lightning dataset DE. By running WRF-CHEM simulations nudged towards observed lightning data instead of relying on parameterizations, we expect the 30 % source of uncertainty to be significantly reduced. Additionally, we mitigated the 29 % influence of the background-NO_x estimate by using Meteosat 12 LI reported groups combined with high horizontal spatial resolution simulations. Furthermore, the use of the upgraded Lightning Mapping Array (LMA) over the Ebro Valley, with its high detection efficiency (DE), allowed us to reduce the 7 % uncertainty associated with the lightning dataset. Assuming independent error sources, the total uncertainty was recalculated as the square root of the sum of the remaining squared contributions provided by Pérez-Invernón et al. (2022). If we only account for the maximum reduction in background-NO_x and lightning dataset uncertainties, the overall uncertainty is reduced by up to 20 %. However, if we also consider the elimination of the 30 % uncertainty from lightning parameterization and related factors, the total uncertainty reduction could reach about 41 %. This means that the total uncertainty, ranging from approximately 20 % to 50 % according to Seiler et al. (2025), would be reduced to a range of 12 % to 40 % .This improvement underscores the critical role of accurate lightning and background estimations in constraining LNO_x PE estimates.

3.3.2 Uncertainty in the obtained relationships between LNO_x PE and lightning characteristics

More significant than the absolute LNO_x PE values are the relationships between the case-to-case variability of LNO_x PE and lightning characteristics. Figures S53 and S54 indicate that variation in the assumed lifetime of NO_x in the near-field of convection (τ) does not significantly influence the obtained relationships between LNO_x PE and lightning characteristics.

We analyze the influence of the background-NO_x estimate in the obtained relationships between LNO_x PE and lightning characteristics by applying different methods to estimate the background-NO_x. In addition to estimating background NO_x from WRF-Chem simulations (Table 1), we also approximated the background NO_x by considering the 10th and the 30th percentiles of the tropospheric VCD NO_x over TROPOMI non-flashing pixels (i.e., pixels without detected lightning activity) (Allen et al., 2021a), as well as without substracting the background-NO_x. The results, presented in Fig. S55, indicate that the choice of method for calculating background NO_x has no significant impact on the observed relationships between lightning characteristics and LNO_x PE, as all obtained correlations are positive regardless of the method used to estimatethe background-NO_x.

Regarding the uncertainty in the characteristics of lightning flashes, the upgraded LMA network XCALMA (van der Velde and Montanyà, 2013; Rodriguez et al., 2025) provides better estimations than the results reported by Pérez-Invernón et al. (2023b). In particular, Pérez-Invernón et al. (2023b) reported the average flash channel lengths for only a portion of the studied thunderstorms due to the limited and inhomogeneous coverage of the previous LMA network. However, in this study we have been able to explore the characteristics of lightning in the entire analyzed thunderstorm cells. In turn, the strong correlations obtained between the geometric and temporal lightning characteristics derived from both LMA and Meteosat 12 LI observations further reinforce the robustness of our results.

4 Conclusions

This study examined the correlations between the geometric and temporal characteristics of lightning flashes in the Ebro Valley in Spain, as measured by LSS, LMA and Meteosat 12 LI systems, and their impact on LNO_x production efficiency (LNO_x PE). The results reveal positive correlations between lightning properties, such as flash channel length, volume, duration, and radiance, the ratio of CG to total flashes, and LNO_x production. While these correlations are moderate, their statistical significance underscores the relevance of lightning characteristics in tropospheric LNO_x emissions.

Estimates of LNO_x production per flash showed slight differences depending on the measurement system and the assumption of the lifetime of NO_x in the near-field of convection (τ). For τ=5 h, the mean value is 37 ± 25 mol NO_x per flash for Meteosat 12 LI flashes, while it is 31 ± 24 mol NO_x per flash when using the LMA lightning flash data. When setting τ=3 h, we obtain 41 ± 29 and 38 ± 31 mol NO_x per flash when using Meteosat 12 LI and LMA lightning data, respectively. The averaged LNO_x PE derived in this work lies toward the lower bound of previously reported estimates for the same region, which span from 42 to 108 moles of NO_x per flash (Pérez-Invernón et al., 2022, 2023b). However, it is important to note that this study primarily focuses on determining the storm-to-storm variability of LNO_x. It does not aim to retrieve an absolute value for NO_x produced per flash. This is because the background NO_x was calculated from WRF-Chem simulations, which already incorporate an assumed LNO_x PE as an input parameter.

The strong correlation between LMA and Meteosat 12 LI geometry and duration measurements reinforces the consistency and robustness of the results, validating the combined use of both systems for studying atmospheric electrical activity and its chemical impact. According to the results of this study, the lightning characteristics measured by the LMA are better proxies for the LNO_x PE than Meteosat 12 LI flash properties, potentially due to the effect of cloud cover in the optical signals detected from space. Positive but weak correlations between contiguous optical signals detected from space and LNO_x PE are lower than the correlation between the total duration of lightning flashes and LNO_x PE. This suggests that there may be in- or below-cloud lightning processes with highly attenuated optical signals that are contributing to LNO_x emissions. Overall, these findings contribute to a deeper understanding of LNO_x emissions and their dependence on lightning properties, providing a solid foundation for improving atmospheric models and environmental impact assessments related to thunderstorms.

In this study, we have refined the method to estimate the LNO_x PE and used an upgraded LMA (XCALMA) with respect to Pérez-Invernón et al. (2023b), as well as including new Meteosat 12 LI lightning data. We analyzed a new set of thunderstorms, obtaining results that are consistent with those reported by Pérez-Invernón et al. (2023b). However, significant sources of uncertainty remain that may influence the obtained results. The relationships between LNO_x PE and lightning characteristics, such as flash channel length, volume, radiance, footprint, and duration, are consistently positive regardless of the approach used to estimate LNO_x PE. Nevertheless, these uncertainties complicate the quantification of these relationships. Future research could address these challenges by exploring additional factors, such as thunderstorm temporal evolution, the vertical structure of LNO_x emissions, or vertical profiles of relative humidity, to further refine these estimates. Combination of LNO_x estimates with systematic monitoring of lightning with continuing current could provide valuable insights into the observed variability of LNO_x PE (Wang et al., 2025; Camino-Faillace et al., 2026). The launch of the geostationary Tropospheric Emissions: Monitoring of Pollution (TEMPO) (Zoogman et al., 2017; Marchenko et al., 2025; Dang et al., 2026) provides air composition measurements over the Americas. Additionally, the Meteosat Third Generation (MTG) Sounder (Holmlund et al., 2021), launched in 2025, is equipped with the Copernicus Sentinel-4 Ultraviolet-Visible-Near-Infrared (UVN) imaging spectrometer, which will report air composition, including tropospheric NO₂ estimates. Meanwhile, the MTG Infrared Sounder (IRS) will deliver vertical profiles of humidity and temperature over Europe. Together, these instruments will enable a deeper understanding of the relationships between thunderstorm evolution, lightning activity, and LNO_x emissions. In particular, the combination of MTG-IRS and MTG-UVN (Holmlund et al., 2021) with the upgraded XCALMA network across the Ebro Valley could significantly advance our knowledge of the key factors contributing to LNO_x emissions, thereby improving their representation in atmospheric models.