Influence of convection on the upper-tropospheric O₃ and NO_x budget in southeastern China

Thunderstorms can significantly influence the air composition via strong updraft and lightning nitrogen oxides (LNO_x). In this study, the ozonesondes and TROPOMI (TROPOspheric Monitoring Instrument) nitrogen dioxide (NO₂) observations for two cases are combined with a model to investigate the effects of typical strong convection on vertical redistribution of air pollutants in Nanjing, southeastern China. The ozonesonde observations show higher O₃ and water vapor mixing ratios in the upper troposphere (UT) after convection, indicating the strong updraft transporting lower-level air masses into the UT and the possible downward O₃-rich air near the top of UT over the convective period. During the whole convection life cycle, the UT O₃ production is driven by the chemistry (5–10 times the magnitude of dynamic contribution) and reduced by the LNO_x (−40 %). Sensitivity tests demonstrate that neglecting LNO_x in standard TROPOMI NO₂ products causes overestimated air mass factors over fresh lightning regions and the opposite for outflow and aged lightning areas. Therefore, a new high-resolution retrieval algorithm is applied to estimate the LNO_x production efficiency. Our work shows the demand for high-resolution modeling and satellite observations on LNO_x emissions of both active and dissipated convection, especially small-scale storms.

1 Introduction

Convection can transport the surface pollutants and moisture from the planetary boundary layer to the upper troposphere (UT) in a short time, where the gaseous pollutants have a longer lifetime due to the slower reaction rates, except for photolysis, in the colder environment (Dickerson et al., 1987). As trace gases remain for more than 1 week in the UT, they are distributed by the upper-level winds around the globe (Ridley et al., 2004). Meanwhile, the vertical profiles of trace gases are reshaped by the updraft and downdraft on timescales of hours (Huntrieser et al., 2016; Barth et al., 2019). On a global scale, the chemical reactions of transported ozone (O₃) and its precursors can increase the amount of UT O₃ (Lawrence et al., 2003; Murray, 2016).

Nitrogen oxides (NO_x), volatile organic compounds (VOCs), and hydrogen oxide radicals (HO_x) are the three main O₃ precursors which can be uplifted or produced by thunderstorms (Bozem et al., 2017). Lightning-produced nitrogen oxides (LNO_x) are the dominant natural source of UT NO_x, contributing as much as 35 %–45 % of global free-tropospheric ozone (Allen et al., 2010; Liaskos et al., 2015). Several lightning parameterizations have been developed for chemistry transport and climate models to evaluate the relationship between LNO_x and other trace gases (Murray et al., 2012; Gordillo-Vázquez et al., 2019; Luhar et al., 2021). Therefore, the precise estimation of LNO_x is crucial for the global O₃ trend and feedback between lightning and climate change (Finney et al., 2016; Chen et al., 2021). Globally, the LNO_x is estimated as 2–8 Tg N yr⁻¹, which is substantially less well quantified than the flash rate (Schumann and Huntrieser, 2007).

The large uncertainty in LNO_x estimation might be reduced by using cloud-resolving chemistry models in combination with satellite and aircraft observations. It is beneficial to take advantage of sonde and satellite observations for exploring the convection effects, especially for the estimation of the LNO_x production efficiency (PE). Recently, many studies have used different satellites to quantify the LNO_x PEs (Beirle et al., 2009; Pickering et al., 2016; Zhang et al., 2020). Two main methods have been proposed to distinguish LNO_x from the NO₂ background pollution: (1) subtracting the weighted temporal average NO₂ of areas with few flashes before the satellite passing time (Pickering et al., 2016; Bucsela et al., 2019; Allen et al., 2019; Lapierre et al., 2020) and (2) directly using customized lightning air mass factors (AMFs) for each convection event (Beirle et al., 2009; Zhang et al., 2020). Recently, Allen et al. (2021) proved the potential of deriving LNO_x PEs by the geostationary lightning instruments (e.g., Lightning Mapping Imager, LMI; Yang et al., 2017; Geostationary Lightning Mapper, GLM; Rudlosky et al., 2019) and NO₂ observations such as Tropospheric Emissions: Monitoring of Pollution (TEMPO; Chance et al., 2019).

Furthermore, aircraft observations and chemical models indicate that the transport from the stratosphere to the troposphere can also increase the UT O₃ besides the chemical production from LNO_x (Pan et al., 2014). As revealed in the Deep Convective Clouds and Chemistry (DC3) 2012 studies and mesoscale convective system simulations, the compensation of subsidence and differential advection beneath the convective core can lead to the anvil wrapping effects (Huntrieser et al., 2016; Phoenix et al., 2020). The different mechanisms of stratosphere–troposphere exchange and the effects on the tropospheric chemistry have been discussed in Holton et al. (1995) and Stohl (2003).

At present, most aircraft observations and model simulations of convection effects focus on the tropics or the United States (Vaughan et al., 2008; Barth et al., 2019). Little is known about the role of convection in southeastern China (Murray, 2016; Guo et al., 2017), where thunderstorms and lightning have increased significantly because of urbanization during recent decades (Yang and Li, 2014; Pérez-Invernón et al., 2021). This is likely due to the increasing aerosol concentration, which can invigorate storms in a moist and convectively unstable environment (Koren et al., 2008; Rosenfeld et al., 2008; Tao et al., 2012). In this study, we combine ground observations and model simulations to investigate the origin of higher UT O₃ and water vapor mixing ratio (Qv) after convection, and we try to distinguish the contributions of physical processes, chemical reactions, and LNO_x. For the first time, the TROPOspheric Monitoring Instrument (TROPOMI; Veefkind et al., 2012) NO₂ observations are used to identify LNO_x PEs in southeastern China. Section 2 describes the datasets used with a brief introduction of the cloud-resolved chemistry model and the LNO_x retrieval method. Section 3 evaluates the model simulations, and the physical and chemical effects of convection are analyzed in Sects. 4 and 5. We apply new a priori NO₂ profiles into the retrieval algorithm to explore the sensitivity of AMFs to LNO_x in Sect. 6. Conclusions are summarized in Sect. 7.

2 Datasets

2.1 Ozonesonde data

Five ozonesondes were launched from the Nanjing National Reference Climatological Station (31.93^∘ N, 118.90^∘ E) on 25 July 2019 and 1 September 2020; both days experienced strong convection. Both pre-convection and during-convection/post-convection campaigns were designed to investigate the convection effects. The convection and ozonesonde trajectories are illustrated in Fig. 1a and b.

Three Institute of Atmospheric Physics (IAP) ozonesondes were launched near the air mass convection that developed on 25 July 2019. The IAP ozonesonde uses an electrochemical concentration cell (ECC). The complete parameters and performance are described in Zhang et al. (2014). The average bias is less than 0.3 mPa from the surface up to 2.5 km, close to zero below 9 km, and less than 0.5 mPa between 9 and 18 km. The first IAP ozonesonde was launched at 05:35 UTC on a sunny day (23 July), and the other two were at 05:10 UTC (pre-convection) and 06:35 UTC (post-convection) on 25 July. Because of water leakage, the pre-convection one lost signal just a few seconds after the release, and instead the ozonesonde launched on 23 July is chosen. Although the time interval is 2 d, the largest relative difference of forecast O₃ profiles above 10 km is usually smaller than 25 % (Fig. S2); therefore the daily variation cannot explain the observed difference of more than 65 %.

The two Vaisala ECC ozonesondes were launched successfully at 23:45 UTC on 31 August (pre-convection) and 06:10 UTC on 1 September (during-convection), following the standard manual to ensure the precision is better than 5 % and the accuracy is within ± (5–10) % below 30 km (Smit et al., 2007). The captured squall line was developing from the convergence of cold air and typhoon Maysak's outer region circulation. Note that the during-convection ozonesonde entered directly into the cloud, providing a unique opportunity of exploring the ozone affected by the convective clouds.

Figure 1 (a, b) The convection detected by the FY-4A Advanced Geostationary Radiation Imager (AGRI) visible channel (0.65 µm) field at the time when the post-convection/during-convection ozonesondes reached around 10 km. The pre-convection ozonesonde trajectories are colored blue, while the others are in red. The white star symbol stands for the observation station, and the thin yellow lines are the TROPOMI swath pixels. Panels (c) and (d) are the observed O₃ (red) and Qv (black) profiles in the pre-convection (dashed) and post-convection/during-convection (solid) periods. The initial (dashed) and simulated post-convection or during-convection (solid) O₃ profiles are in blue. The dark gray shading is the 50 % confidence interval, while the light one is the 90 % confidence interval. The gray lines are the lapse rate tropopauses.

2.2 Lightning data

Three lightning datasets were used in this study: the China National Lightning Detection Network (CNLDN; Yang et al., 2015), the Earth Networks Total Lightning Network (ENTLN; Marchand et al., 2019), and the World Wide Lightning Location Network (WWLLN; Rodger et al., 2006). The detection efficiency (DE) of cloud-to-ground (CG) flashes is about 90 % for the CNLDN data in Jiangsu Province (Li et al., 2017a), while ENTLN and WWLLN detect both intra-cloud (IC) and CG flashes with a specific detection frequency (1 Hz–12 MHz for ENTLN and 3–30 kHz for WWLLN). In the ENTLN data, groups of pulses are classified as a flash if they are within 700 ms and 10 km. Both strokes and lightning flashes composed of one or more strokes are included in the preprocessed data obtained from the ENTLN. The detailed processing algorithm of the WWLLN is given by Rodger et al. (2004). The WWLLN strokes and pulses are combined with ENTLN into one dataset (ENGLN) within 10 km and 0.7 s as mentioned in Virts and Goodman (2020). To increase the lightning data coverage in our study, the CG flashes of ENGLN and CNLDN datasets are combined using spatial and temporal clustering criteria of 10 km and 0.5 s (Zhao et al., 2020). Merging these three datasets should provide a sufficiently high CG flash detection efficiency for this analysis. Because the IC DE of all these lightning data is low in China, we conservatively use the merged CG data with a constant IC $/$ CG ratio of 3:1 based on Wu et al. (2016) and Bandholnopparat et al. (2020). IC data will become more accurate if more Chinese total lightning networks, such as Beijing Lightning Network (BLNET; Srivastava et al., 2017), are available to be compared with lightning imaging sensors (Rudlosky and Shea, 2013; Poelman and Schulz, 2020).

2.3 TROPOMI data

On 13 October 2017, TROPOMI on the Sentinel-5 Precursor satellite was launched successfully (Veefkind et al., 2012). For the current study, we used the Royal Netherlands Meteorological Institute (KNMI) standard product v2.1 test as input to our LNO_x retrieval algorithm. A spike removal is included in the product to better deal with the detector saturation and blooming effects, which enables more valid data over the bright clouds generated by convection (Ludewig et al., 2020; van Geffen et al., 2022). Each pixel includes a slant column density quality flag (no2_scd_flag = 0, see Appendix A) that can be used to get the data without known retrieval errors (Allen et al., 2021; van Geffen et al., 2022).

The official NO₂ columns are retrieved using the near-ultraviolet and visible (UV-VIS, 405–465 nm) spectrometer backscattered solar radiation measurements on TROPOMI (van Geffen et al., 2015). The retrieval consists of three main procedures for each measured Level-1b spectrum.

1. Total NO₂ slant column density (SCD) is determined by the DOAS (differential optical absorption spectroscopy) method.

2. The stratospheric and tropospheric SCDs are separated by data assimilation of slant columns in the tracer model, version 5, tailored for the application of satellite retrievals (TM5-MP; Williams et al., 2017).

3. The stratospheric and tropospheric NO₂ vertical column densities (VCDs) are obtained via the air mass factor (AMF) lookup tables (Lorente et al., 2017).

We replaced the tropospheric AMF (AMF_trop) with a new AMF called AMF for LNO_x (AMF${}_{{\mathrm{LNO}}_x}$) to derive the tropospheric LNO_x vertical column density (VCD${}_{{\mathrm{LNO}}_x}$). The concept of AMF${}_{{\mathrm{LNO}}_x}$ is inherited from the AMF_trop derived by a function of several parameters (solar zenith angle, viewing zenith angle, relative azimuth angle, surface albedo, surface pressure, cloud fraction, cloud height, and a priori trace gas profile). Briefly, the numerator is the modeled tropospheric NO₂ slant column density (SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$) and the denominator is the modeled VCD (VCD${}_{{\mathrm{NO}}_{\mathrm{2}}}$ or VCD${}_{{\mathrm{LNO}}_x}$). In detail, these two AMFs can be calculated as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathrm{AMF}}_{\mathrm{Trop}}={\displaystyle \frac{\begin{array}{c}\left(\mathrm{1}-f_{\mathrm{r}}\right){\int }_{p_{\mathrm{surf}}}^{p_{\mathrm{tp}}}{w}_{\mathrm{clear}}\left(p\right){\mathrm{NO}}_{\mathrm{2}}\left(p\right)\mathrm{d}p\\ +f_{\mathrm{r}}{\int }_{p_{\mathrm{cloud}}}^{p_{\mathrm{tp}}}{w}_{\mathrm{cloudy}}\left(p\right){\mathrm{NO}}_{\mathrm{2}}\left(p\right)\mathrm{d}p\end{array}}{{\int }_{p_{\mathrm{surf}}}^{p_{\mathrm{tp}}}{\mathrm{NO}}_{\mathrm{2}}\left(p\right)\mathrm{d}p}},\text{(2)}& {\displaystyle }{\mathrm{AMF}}_{{\mathrm{LNO}}_x}={\displaystyle \frac{\begin{array}{c}\left(\mathrm{1}-f_{\mathrm{r}}\right){\int }_{p_{\mathrm{surf}}}^{p_{\mathrm{tp}}}{w}_{\mathrm{clear}}\left(p\right){\mathrm{NO}}_{\mathrm{2}}\left(p\right)\mathrm{d}p\\ +f_{\mathrm{r}}{\int }_{p_{\mathrm{cloud}}}^{p_{\mathrm{tp}}}{w}_{\mathrm{cloudy}}\left(p\right){\mathrm{NO}}_{\mathrm{2}}\left(p\right)\mathrm{d}p\end{array}}{{\int }_{p_{\mathrm{surf}}}^{p_{\mathrm{tp}}}{\mathrm{LNO}}_x\left(p\right)\mathrm{d}p}},\end{array}$$

where p_surf is the surface pressure, p_tp is the tropopause pressure, p_cloud is the cloud optical pressure, f_r is the cloud radiance fraction in the NO₂ window, w_clear and w_cloudy are respectively the pressure-dependent scattering weights from the lookup Table (Lorente et al., 2017) for clear and cloudy parts, and NO₂(p) is the NO₂ vertical profile simulated by WRF-Chem (Weather Research and Forecasting model coupled with chemistry). Moreover, LNO_x(p) is the LNO_x vertical profile calculated by the differences in the vertical profile between WRF-Chem simulations with and without lightning. All other parameters in the KNMI v2.1 test product, including the total SCD, stratospheric SCD, total VCD, stratospheric VCD, surface albedo, and scattering weights, remain unchanged.

In comparison with this study, Pickering et al. (2016), Allen et al. (2019), Bucsela et al. (2019), and Allen et al. (2021) defined another AMF${}_{{\mathrm{LNO}}_x}$ to convert SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ to the tropospheric NO_x vertical column density (VCD${}_{{\mathrm{NO}}_x}$). Then, their VCD${}_{{\mathrm{LNO}}_x}$ can be calculated as the VCD${}_{{\mathrm{NO}}_x}$ subtracted from a tropospheric NO_x background. Because our AMF${}_{{\mathrm{LNO}}_x}$ converts the SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ to VCD${}_{{\mathrm{LNO}}_x}$ directly, the additional estimation of background NO₂ is not needed for calculating LNO_x PE in Sect. 6.2.

2.4 Model simulations

This study uses WRF-Chem version 4.1.4. The initial and boundary conditions of meteorological parameters are provided by the hourly ECMWF atmospheric reanalysis (ERA5) data (Hersbach et al., 2020). Simulations are performed with one‐way nesting with 75 vertical levels and a 50 hPa model top. The domain settings are illustrated in Fig. S1. The microphysical processes are computed with the WRF Single-Moment 6-class scheme (WSM6; Hong and Lim, 2006), while the shortwave and longwave radiation is calculated by the Rapid Radiative Transfer Model for GCMs (general circulation models) scheme (RRTMG; Iacono et al., 2008). The land surface processes are simulated by the Noah scheme (Koren et al., 1999). However, we use different planetary boundary layer (PBL) parameterizations to simulate the convection. Specifically, the 2019 case uses the Yonsei University scheme (YSU; Hong et al., 2006), while the quasi-normal scale elimination (QNSE; Sukoriansky et al., 2005) is applied to the 2020 case.

The chemical initial and boundary conditions are defined using the output from the Whole Atmosphere Community Climate Model (WACCM; https://www.acom.ucar.edu/waccm/, last access: 22 April 2022). The initial O₃ profile of the 2020 case is replaced by the O₃ profile from the ozonesonde. Anthropogenic emissions are driven by the 2016 Multi-resolution Emission Inventory for China (MEIC) version 1.3 (http://www.meicmodel.org/, last access: 22 April 2022). The Model of Emissions of Gases and Aerosol from Nature (MEGAN; Guenther et al., 2006) is used for biogenic emissions. The chemical mechanism is the Model for Ozone and Related chemical Tracers (MOZART) for gas-phase chemistry and Goddard Chemistry Aerosol Radiation and Transport (GOCART) model for aerosols (Pfister et al., 2011). The photolysis rates are adjusted by the presence of aerosols and clouds using the new TUV (tropospheric ultraviolet and visible) photolysis option with the scaled cloud optical depth (cloud_fraction^1.5). Note that the bimodal profile modified from the standard Ott et al. (2010) profile is employed as the vertical distribution of LNO in WRF-Chem (Laughner and Cohen, 2017), while the LNO_x parameterization is activated as 500 mol NO per flash (Zhu et al., 2019). The resulting lightning nitrogen monoxide (LNO) and lightning nitrogen dioxide (LNO₂) profiles are defined as the difference of vertical profiles between simulations with and without lightning.

To simulate the convection and LNO_x realistically, a lightning data assimilation (LDA) technique is applied to WRF-Chem. The details of the LDA technique can be found in Fierro et al. (2012) and Li et al. (2017b) and are briefly illustrated here. The water vapor mass mixing ratio is increased at constant temperature layers in columns where flashes occur:

$$\begin{array}{}\text{(3)}& Q_v=A{Q}_{\mathrm{sat}}+B{Q}_{\mathrm{sat}}\tanh \left(CX\right)\left[\mathrm{1}-\tanh \left(D{Q}_g^{\mathit{\alpha }}\right)\right],\end{array}$$

where Q_sat is the water vapor saturation mixing ratio (g kg⁻¹), Q_g is the graupel mixing ratio (g kg⁻¹), and X is the flash rate. In our simulations, the layer between 263.15 and 290.15 K is chosen to let the convection root itself in the PBL quickly as the Q_v in the lower troposphere is the deeper layer (Marchand and Fuelberg, 2014; Finney et al., 2016; Li et al., 2017b). Parameter settings follow Li et al. (2017b): A=0.94, B=0.2, C=0.001, and D=0.25 and α=2.2. The resampled total lightning flash data are read through the auxiliary input stream of WRF every 10 min. For example, if the beginning time of LDA is 05:00 UTC with a time step of 10 min, all flashes in a specific grid between 05:00 and 05:10 UTC are summed as the contribution during this period. At the next time step, the flashes are classified as the next new group. Therefore, the flash count is the flash rate density (units: flashes per 10 min per dx km per dy km), where dx and dy are the resolutions of model grids in the x and y direction, respectively. This leads to flash rate densities that are the same through all nested domains, as done in Fierro et al. (2012) and Li et al. (2017b), and these are used as the flash counts in WRF-Chem directly following the method employed in CMAQ (Community Multiscale Air Quality; Kang et al., 2019a, b, 2020).

3 Model evaluation

Compared with the radar observations, the simulated initiations of convection are ahead by 60 and 30 min for the 25 July 2019 and 1 September 2020 cases, respectively (Figs. 2 and 3). The lightning assimilations are applied with the same time step forwards. For comparisons below, we choose the matched stages instead of the same time.

The air mass convection on 25 July 2019 was initialized as isolated cells. WRF-Chem reproduces the position and intensity of isolated convections at the initial stage (Fig. 2a and d). At 05:40 UTC, the cells were presented with the northeast–southwest orientation, and the column maximum radar reflectivity (CRF) reached 60 dBZ (Fig. 2b), which is stronger than the simulated convection (CRF = 55 dBZ, Fig. 2e). The vertical cross section of simulated radar reflectivity across the core of cells is compared with the observations (Fig. S3). Although there were much missing data caused by the long distance between convection and radar, the horizontal and vertical structures of isolated cells are roughly shown without artificial interpolation. While the simulated 45 dBZ contour extends to 12 km, the observed one only reaches 10 km because of reduced data quality above 10 km.

The squall line on 1 September 2020 was formed in the north, strengthened, and moved towards the observation site (Fig. 3). The strongest convective stage with a CRF of 60 dBZ was at 05:50 UTC, which is the TROPOMI overpass time (Fig. 3b and e). Although the highest level reached was lower than that of the 2019 case, the reflectivity in the lower troposphere (2–8 km) was larger and broader (Fig. S4). Note that the simulated dissipated cells deviate from the radar observation, and this leads to the region for the ozonesonde comparison moving to the west of the station (Fig. 3c and f).

Figure 2 Observed radar composite reflectivity at (a) 04:20, (b) 05:40, and (c) 06:50 UTC. (d–e) WRF-Chem simulated composite reflectivity 1 h before the radar observation times. The A–B solid lines in (b) and (e) are cross-section lines for Fig. S3. The C–D solid line in (e) is the cross-section line for Fig. 4b. The black rectangle is the region for the comparison with ozonesonde.

The measured O₃ and Qv profiles at different convection stages are shown in Fig. 1c and d. Generally, the observed UT O₃ and Q_v are higher with convection, while the largest enhancements are between 10 and 16 km. However, the 2020 case showed a larger increase in the lower troposphere (2–8 km, LT). Additionally, a two-valley shape of O₃ profile exists in both cases but at different levels: 2 and 8 km for the 2019 air mass and 4 and 10 km for the 2020 squall line. Although the WRF-Chem model tends to underestimate the O₃ concentration in the LT and UT for the 2019 and 2020 cases, respectively, it can reproduce the detailed O₃ structures and provide the opportunity to analyze the mechanisms of convection.

Three possible sources can explain the enhancements of O₃ in the UT: convective transport, chemical production, and O₃ directly produced by lightning. Only the first two factors are discussed in detail in Sect. 4 as lightning O₃ is beyond the scope of this study and still uncertain, as shown by limited observations and model simulations (Morris et al., 2010; Ripoll et al., 2014).

Figure 3 Same as Fig. 2 but for the case on 1 September 2020. The simulation time is 30 min ahead of each radar observation.

4 Convection impacts

To evaluate the higher UT O₃ concentration after the convection, the mean vertical profiles of O₃ in the atmospheric regions passed by the ozonesondes are illustrated at three stages of convection: initiation, development, and dissipation (Fig. 4a and d). The UT O₃ increases continuously during the whole cycle of the 2020 case, but for the 2019 case, it declines at the developing stage because of the uplift of O₃-poor air, while later it starts rising again. This phenomenon can be explained by the O₃ vertical cross sections at the developing stages (Fig. 4b and e). For the 2019 case, the cells of low O₃ concentration reach 16 km by the updraft, and then the high O₃ air wrapped behind the convection moves into the region. However, the observed increasing O₃ of the 2020 case is mostly from the vertically transported background O₃.

Figure 4 (a, d) The mean O₃ profiles in the regions passed by the ozonesondes at three stages: initiation (black), development (orange), and dissipation (blue). (b, e) Vertical O₃ distribution within the convective periods along the line crossing the convective core (Figs. 2e and 3f). The dashed blue lines stand for the boundaries of regions passed by the ozonesondes, and the lapse rate tropopause is shown as the white dots. The cloud boundaries (the sum of cloud liquid water mixing ratio [q_cloud] and ice mixing ratio [q_ice] ≥ 0.01 g kg⁻¹) are shown with white lines. (c, f) The vertical distributions of the O₃ net production rate and tendency due to horizontal advection (advh), vertical advection (advz), and chemistry (chem) during the convective periods.

To determine the processes causing the differences between the two cases, we analyzed the outputs of mean integrated physical rates (IPRs) from 10 to 14 km during the convective period (Fig. 4c and f). Generally, the opposite trend of the horizontal advection (advh) and vertical advection (advz) governs the net decrease in the UT O₃ production rate of the 2019 case. Note that the advz contribution is negative between 10 and 11.5 km and positive between 11.5 and 13.8 km. This is due to the uplifted O₃-poor air and downward O₃-rich air caused by the stronger updraft compared with the 2020 case. As indicated by the higher Q_v after convection (Fig. 1c) and tropopause height (Fig. 4b), the updraft of the 2020 case is not strong enough to wrap the stratospheric O₃ like the mesoscale convective system (Phoenix et al., 2020). While the dynamic processes play an important role in the O₃ production, the positive chemistry contribution cannot be neglected in both cases and leads to the net increase in UT O₃ during the convective period of the 2020 case. Specifically, the chemistry increases O₃ in both cases, and the magnitude of the effect is 5–10 times that of dynamic effects. This demonstrates the dominant chemistry role in the overall effects of convection.

5 Impact of lightning NO_x on the O₃ profile

Furthermore, the IPR outputs including LNO emissions are compared with these excluding LNO to explore the effects of LNO_x on O₃ (Table 1). The LNO_x reduces the net O₃ production by 25 % during the convective period and life cycle of the 2019 case. The decreased chemistry contribution is less significant (≤ 1 %) for the 2020 case which has a smaller lightning density near the station. Note that the LNO_x can certainly enhance the downwind ozone production on the scale of days (Pickering et al., 1996; DeCaria et al., 2005). Therefore, it is necessary to estimate the LNO_x PE accurately (discussed later in Sect. 6.2).

Additionally, the convection was divided into three regions by TROPOMI data: fresh lightning, downwind of fresh lightning, and aged lightning (See Sect. 6.1 and Fig. 8 for details). Firstly, the difference of O₃ (ΔO₃) profiles is obtained with different LNO PE assumptions (Fig. 5a–c). In contrast with the net loss of ozone (< 4 ppbv) over all height levels in Ott et al. (2007), the ΔO₃ is mostly positive (< 1 ppbv) between 2 and 5 km and negative ($>-\mathrm{3}$ ppbv) between 5 and 12 km in our cases. The higher PE (700 mol per flash) slightly reduces the O₃ concentration by less than 1 ppbv at all levels compared with the default PE (500 mol per flash), and it even leads to negative ΔO₃ between 2 and 5 km downwind of fresh lightning (Fig. 5b). The maximum O₃ loss is between 8 and 10 km due to the peak of LNO_x (up to 2.6 ppbv) introduced in the model.

Then, the integrated reaction rate (IRR) is applied to determine the chemistry mechanism and the effect of LNO on the O₃ variation for two layers in which the ΔO₃ is opposite: 800–500 hPa (ΔO₃>0) and 500–200 hPa (ΔO₃<0). The tropospheric O₃ is mainly controlled by five reaction rate terms (Pickering et al., 1990; Bozem et al., 2014):

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left[{\mathrm{O}}_{\mathrm{3}}\right]& =k_{\mathrm{1}}\left[\mathrm{NO}\right]\left[{\mathrm{HO}}_{\mathrm{2}}\right]+\sum _i{k}_i\left[\mathrm{NO}\right]\left[{\mathrm{R}}_{\mathrm{i}}{\mathrm{O}}_{\mathrm{2}}\right]\\ & -k_{\mathrm{3}}\left[{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\right]\left[\mathrm{O}{(}^{\mathrm{1}}\mathrm{D})\right]-k_{\mathrm{4}}\left[{\mathrm{HO}}_{\mathrm{2}}\right]\left[{\mathrm{O}}_{\mathrm{3}}\right]\\ & -k_{\mathrm{5}}\left[\mathrm{OH}\right]\left[{\mathrm{O}}_{\mathrm{3}}\right],\end{array}\end{array}$$

where k_i is the rate coefficient of the reaction between peroxy radicals (R_iO₂) and NO. The time series of each contribution to O₃ production are illustrated in Fig. 5d–i. Overall, the time series of IRR are more variable for the 2019 case due to the stronger activity as clarified in Sect. 4. The total net chemistry production of O₃ remains positive for both layers. In detail, the reaction between NO and HO₂ always dominates the production, while the oxidation of NO by RO₂ is about 40 %–60 % of that production. The dominant loss of O₃ is the photolysis (described by the reaction of O(¹D) and H₂O), while the reaction between O₃ and OH is comparable during the convective period. The lowest contribution to O₃ loss, O₃ + HO₂ → OH + 2O₂, is reduced during the convection because of the production of LNO which captures the HO₂ reacting with O₃. Note that although the increase in total IRR induced by LNO_x can reach 1.36 × 10⁷ and 2.60 × 10⁶ molec. cm⁻³ s⁻¹ in the low layer and high layer over these three regions, respectively, the net O₃ production actually decreases in the high layers (Fig. 5a–c) due to the combination of dynamic transport and chemical production related to LNO_x.

Figure 5 (a–c) Changes in O₃ profiles due to LNO_x at TROPOMI overpass time in three regions (fresh lightning region, downwind of fresh lightning, and aged lightning area) as defined in Fig. 7. (d–f) Time series of the mean integrated reaction rate (IRR) between 800 and 500 hPa. The legend shows detailed species and reactions. The “total IRR” is the O₃ loss IRR subtracted from the O₃ production IRR (red and purple lines). The solid line shows the IRR with LNO (500 mol per flash), while the dashed line is without LNO. (g–i) Same as (d)–(f) but between 500 and 200 hPa.

6 TROPOMI products over the convection

6.1 Relation between lightning and TROPOMI products

As the LNO_x production estimation from TROPOMI depends upon the SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$, we compare the SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ distributions with the observed lightning flashes (Fig. 6a–d). Although the SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ over the most active pixels is not valid due to the detector saturation and blooming effect, the nearby or outflow regions still have useful data. While the flashes occurred less than 30 min before the TROPOMI overpass time in the 2019 case, the 2020 case had both fresh and aged LNO₂ (Fig. 6d).

Figure 6 Events on 25 July 2019 (a, c) and 7 September 2020 (b, d). (a, b) The tropospheric NO₂ slant column density (SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$, filled color) and cloud optical pressure (line color). These white grid cells stand for missing TROPOMI data (no2_scd_flag ≠ 0), as defined in Appendix A. The solid black border is Jiangsu Province. (c, d) The cloud radiance fraction in the NO₂ window and flashes whose color depends on the occurrence time relative to the TROPOMI overpass time.

Specifically, the SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ of the convective pixels (f_r ≥ 0.7) is smaller than that in other regions. This is opposite to previous studies of large-scale convective systems with high flash density (Beirle et al., 2009). Four factors could lead to this unexpected result: the cloud top heights, flash counts, flash occurrence time, and background NO₂. Either the inadequate flash or weak convection could lead to a smaller SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ over pixels with f_r ≈ 1 because TROPOMI can only see the LNO₂ above the clouds. In other words, the polluted NO₂ below the broken or thinner clouds is partially exposed if f_r < 1. The sensitivity tests of the WRF-Chem a priori SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ can help explain this phenomenon clearly (Fig. S5). The pixels of high SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ with low cloud fraction belong to background NO₂ pollution (Fig. S5a and e), but the SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ increased by UT LNO₂ is still visible compared with the lower SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ without LNO₂ (Fig. S5b–d and f–h).

An LNO PE upper limit of 700 mol NO per flash (Ott et al., 2010) is applied to WRF-Chem for investigating the importance of LNO_x for the AMF_trop and AMF${}_{{\mathrm{LNO}}_x}$ calculations. The changes in the retrieved AMFs are examined by replacing profiles in three tropospheric layers independently: middle troposphere (MT, 800 to 400 hPa), upper troposphere (UT, 400 to 150 hPa), and troposphere (surface to tropopause). Unless otherwise specified, the changes in AMFs are obtained by increasing LNO_x in this section. Figure 7 shows that the AMF changes are mostly controlled by the LNO_x in the UT layer where the detection sensitivity is high (Beirle et al., 2009; Laughner and Cohen, 2017) and the LNO_x production reaches the peak (Fig. S6). While the AMF${}_{{\mathrm{LNO}}_x}$ decreases by 5 %–40 % for both cases, the changes in AMF_trop (ΔAMF_trop) are regionally specific and can be classified by the lightning activity: fresh lightning (MT ΔAMF${}_{\mathrm{trop}}<-\mathrm{20}$ %), downwind of fresh lightning (MT ΔAMF_trop>20 %), and aged lightning (UT ΔAMF_trop>20 %). Figure 8a illustrates the relationship between p_cloud and f_r over these three regions. The clouds are higher than 400 hPa (p_cloud<400 hPa), and f_r is larger than 0.8 over fresh lightning pixels, but both aged lightning and downwind of fresh lightning areas have clouds lower than 400 hPa. This coincides with the mean cloud pressures in Fig. S6 and explains why UT ΔAMF_trop>20 % exits in Fig. 7bi and biii, indicating the possibility of LNO_x estimations over the aged lightning regions (Sect. 6.2).

Table 1 Process analysis table for the mean O₃ integrated tendencies (10–14 km).

^* The unit is 10¹⁰ molec. cm⁻³. The percentage is the proportion of each part in the net O₃ change.

As defined in Appendix B, SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}^{{\mathrm{LNO}}_x}$ $/$ SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}^{{\mathrm{noLNO}}_x}$ and VCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}^{{\mathrm{LNO}}_x}$ $/$ VCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}^{{\mathrm{noLNO}}_x}$ can be used to determine which parameter controls ΔAMF_trop: enhanced a priori SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ or a priori VCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ (Fig. 8b–d). Briefly, the dominant one belongs to the larger ratio. First, if the LNO₂ is included in the tropospheric layer of a priori NO₂ profiles (Fig. 7ci), the AMF_trop decreases over most fresh lightning pixels because of the increased a priori VCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ (Fig. 8b). The situation is opposite for the downwind of the fresh lightning region. There the AMF_trop is larger no matter which layer is chosen (Fig. 7ai–ci) because the LNO₂ was convected, reached above the cloud top, and led to the larger a priori SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ (Fig. 8c). Interestingly, for the 2020 case, the AMF_trop of aged lightning pixels increases more than 50 % for the UT layer (Fig. 7biii). It demonstrates the important role of advected UT LNO₂ and that the cloud exists as a barrier, causing the difference between a priori SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ and a priori VCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ (Fig. 8d). Although the difference is smaller than that of the other two regions due to the LNO₂ lifetime, it is still useful for retrieving the LNO_x. Moreover, considering the region-specific LNO_x effects on AMFs, we need to include the representation of LNO₂ in the TROPOMI NO₂ retrievals better, especially outflow regions. The comparisons (Fig. S6) between the TROPOMI standard NO₂ profiles from TM5-MP and WRF-Chem also illustrate the importance of LNO_x, resolved convection transport, and emissions. Aircraft observations of NO and NO₂ will be useful to determine the exact roles (Laughner and Cohen, 2017).

Figure 7 The percent differences of AMFs by replacing the a priori NO₂ profiles at three layers: middle troposphere (left), upper troposphere (middle), and troposphere (right). ΔAMF_trop is the comparison of the AMF_trop with 500 mol NO per flash relative to 0 mol NO per flash. ΔAMF${}_{{\mathrm{LNO}}_x}$ is the comparison of the AMF${}_{{\mathrm{LNO}}_x}$ with 700 mol NO per flash relative to 500 mol NO per flash. Three regions are annotated: fresh lightning (blue), downwind of fresh lightning (red), and aged lightning (black). Because of the quite large AMF${}_{{\mathrm{LNO}}_x}$ values in pixels with little lightning, ΔAMF${}_{{\mathrm{LNO}}_x}$ is shown over pixels where 0 < AMF${}_{{\mathrm{LNO}}_x}$ < 10.

Figure 8 (a) The relationship between cloud pressure and cloud radiance fraction for three regions defined in Fig. 7: fresh lightning region, downwind of fresh lightning, and aged lightning area. (b–d) The a priori SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}^{{\mathrm{LNO}}_x}$/SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}^{{\mathrm{noLNO}}_x}$ and a priori VCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}^{{\mathrm{LNO}}_x}$/VCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}^{{\mathrm{noLNO}}_x}$ of pixels in these three regions. The LNO_x superscript indicates that the a priori variable is calculated with LNO_x (500 mol NO per flash), and the noLNO_x superscript is without LNO_x.

6.2 Estimations of LNO_x

Satellite observations with fresh convection are usually used to estimate the LNO_x PE. However, it is difficult to apply the same method to regions with small convection like the 2019 case because of the pixel saturation of TROPOMI and limited coverage of convective area (Fig. 6c). Instead, we focus on dissipated convection, which is the southern part of the 2020 case, while the northern part has multiple missing values regardless of the adequate flashes and high cloud fractions. As shown in Fig. 9b–c, the time differences between flashes and TROPOMI overpass time are longer than 30 min but shorter than 3 h. Since the lifetime of NO₂ is ∼ 3 h in or near the field of convection (Nault et al., 2016) and ranges from 2 to 12 h depending on the convective location, these pixels can still be used for the LNO_x estimation. Equation (5) is applied to determine the mean LNO_x PE (mol per flash):

$$\begin{array}{}\text{(5)}& {\mathrm{PE}}_{{\mathrm{LNO}}_x}=\sum _p{V}_i{A}_i/\sum _N{F}_j{e}^{-\left(t_{\mathrm{0}}-t_j\right)/\mathit{\tau }},\end{array}$$

where p stands for pixels affected by LNO_x, V_i (mol m⁻²) is the LNO_x vertical column densities (VCD${}_{{\mathrm{LNO}}_x}=$ SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$ $/$ AMF${}_{{\mathrm{LNO}}_x}$) over pixel i with an area called A_i (m²), N is the total number of flashes contributing to VCD${}_{{\mathrm{LNO}}_x}$, and the exponential component considers the lifetime of NO_x for each flash (F_j). Specifically, t₀ is the time of TROPOMI overpass, t_j is the time of the lightning flash, and τ represents the NO_x 3 h lifetime near convection.

As the dissipated convection produced enough lightning and the UT winds within the storm were blowing from west-northwest to east-southeast (Fig. 9a), the pattern of VCD${}_{{\mathrm{LNO}}_x}$ can still be clearly identified (dashed rectangle in Fig. 9). Fortunately, there is a low VCD${}_{{\mathrm{LNO}}_x}$ strip separating the northern and southern convection. With the careful selection, the LNO_x PE is estimated as 60 mol NO_x per flash. Although there are a few lightning flashes related to the VCD${}_{{\mathrm{LNO}}_x}$ outside of the region selection, it only affects the LNO_x PE by ≈ 2 mol, which is within the uncertainty discussed below.

Figure 9 The background is the distribution of LNO_x vertical column densities (VCD${}_{{\mathrm{LNO}}_x}$). The white rectangles are manually selected regions for the LNO_x PE estimation. The overlaid wind arrows in (a) are the 500 hPa horizontal wind simulated by WRF-Chem. The lightning dots in (b) and (c) are the flashes whose color depends on the occurrence time relative to the TROPOMI overpass time.

Following Allen et al. (2019) and Zhang et al. (2020), the uncertainty of LNO_x is determined by the LNO_x lifetime, lightning DE, NO / NO₂ ratio, LNO profile, and other sources (summarized in Table 2). The lifetime (τ) of NO₂ is replaced by 2 and 6 h to evaluate the uncertainty as 27 %, while another uncertainty is also 27 %, which is related to lightning DE by changing the ratio of IC to CG to 2:1 and 4:1. The uncertainty caused by modeled NO / NO₂ ratios is assumed to be 30 % based on Allen et al. (2019), and the uncertainty related to the LNO profile is 26 % by using the a priori NO₂ profile with 330 and 700 mol NO per flash. The uncertainty associated with the stratospheric vertical column is considered as 7 % by applying a bias of ± 10¹⁴ molec. cm⁻² (van Geffen et al., 2022). The uncertainty caused by other possible error sources including systematic errors in slant columns, optical cloud pressure, and NO₂ redistributed by convection is difficult to quantify and is assumed to be 10 % following Allen et al. (2021). Assuming no correlation between errors, the total uncertainty (56 %) is estimated as the square root of the sum of the squares of all individual uncertainties. As a result, the LNO_x PE is 60 ± 33 mol NO_x per flash. It is less than that of our previous work (90 ± 50 mol NO_x per flash) over the continental United States (Zhang et al., 2020) and at the lower end of 120 ± 65 mol NO_x per flash obtained in Allen et al. (2021). Thus, more studies over China are necessary for the estimation of region-dependent LNO_x PEs.

Table 2 Uncertainties for the estimation of LNO_x production efficiency.

7 Conclusions

Both the 2019 and 2020 cases saw enhanced UT O₃ concentrations according to the ozonesonde observations. As revealed by the WRF-Chem model, the dynamic contribution of O₃ variation in the 2019 case was generated by mixing with the UT O₃-rich air due to strong updraft, while it was caused by vertical advection of high background O₃ in the 2020 case. The detailed analysis of integrated physical rates shows that the dynamic processes dominate the UT O₃ decrease during the convective stage of both cases. However, in the convection life cycle, the contribution of chemistry reactions to the UT O₃ production is 5–10 times that of dynamics. Moreover, the UT O₃ enhancement of the 2019 case decreases by 40 % if the lightning nitrogen oxides (LNO_x) are included in the model, indicating the importance of the LNO_x. Utilizing the outputs of integrated reaction rates improves our understanding of the chemical reactions and effects of lightning nitrogen oxides (LNO_x) on O₃. While the reaction of NO and HO₂ dominates the production of O₃, the O₃ loss processes are suppressed because the LNO_x reacts with the HO₂ which is the reactant in the production.

The WRF-Chem results are incorporated into the retrieval algorithm to explore how LNO_x affects the official TROPOMI products and prove that TROPOMI is useful for LNO_x studies over small-scale convection regions. The sensitivity tests imply that air mass factors of tropospheric NO₂ are smaller for fresh lightning regions and the opposite for the aged lightning pixels or downwind of fresh lightning. The air mass factors of tropospheric LNO_x always decrease with increased LNO_x. Since the LNO_x affects the variation in tropospheric NO₂ in the outflow regions, better consideration of LNO_x is essential for studies of the tropical and mid-latitude regions in summer.

Because the saturation of TROPOMI pixels and blooming effects lead to the failure of detecting LNO_x over active convection, we focus on the dissipation regions where aged LNO_x still exists. The production efficiency is estimated to be 60 ± 33 mol LNO_x per flash, which is less than but still within the range of previous studies. Although the current results are limited, the technique can be implemented worldwide to region-specific LNO_x retrievals. To quantify and refine the LNO_x production estimates, TROPOMI data over both active and dissipated convection could provide valuable information.

Appendix A: Flag definitions used in this study

The no2_scd_flag is introduced to make the usage of the NO₂ slant column density (SCD) data easier by gathering information from a few variables into one flag (Table A1). This flag can thus be used for filtering, though with care as it probably does not cover all possible situation. Here “δ” refers to the SCD error (in mol m⁻²), and “pqf” stands for processing quality flag.

Table A1 Definition of no2_scd_flag.

Appendix B: Contributions to ΔAMF_trop

In Fig. 8b–d, the contribution of LNO_x to ΔAMF_trop is divided into two parts: SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}^{{\mathrm{LNO}}_x}$ $/$ SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}^{{\mathrm{noLNO}}_x}$ and VCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}^{{\mathrm{LNO}}_x}$ $/$ VCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}^{{\mathrm{noLNO}}_x}$, where the LNO_x superscript indicates that the a priori variable is calculated with LNO_x (500 mol NO per flash), and the noLNO_x superscript is without LNO_x. The two contributions are derived by taking the logarithms of Eqs. (B1) and (B2) and then subtracting them in Eq. (B3). Here, several abbreviations are defined to simplify the symbols: S is SCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$, V is VCD${}_{{\mathrm{tropNO}}_{\mathrm{2}}}$, subscript 1 is with LNO_x, and 0 is without LNO_x.

$$\begin{array}{}\text{(B1)}& {\displaystyle }{\mathrm{AMF}}_{\mathrm{1}}={\displaystyle \frac{S_{\mathrm{1}}}{V_{\mathrm{1}}}}\text{(B2)}& {\displaystyle }{\mathrm{AMF}}_{\mathrm{0}}={\displaystyle \frac{S_{\mathrm{0}}}{V_{\mathrm{0}}}}\text{(B3)}& {\displaystyle }\begin{array}{rl}\log \left({\mathrm{AMF}}_{\mathrm{1}}\right)-\log \left({\mathrm{AMF}}_{\mathrm{0}}\right)& =\log \left({\displaystyle \frac{S_{\mathrm{1}}}{V_{\mathrm{1}}}}\right)-\log \left({\displaystyle \frac{S_{\mathrm{0}}}{V_{\mathrm{0}}}}\right)\\ & =\log \left({\displaystyle \frac{S_{\mathrm{1}}}{S_{\mathrm{0}}}}\right)-\log \left({\displaystyle \frac{V_{\mathrm{1}}}{V_{\mathrm{0}}}}\right)\end{array}\end{array}$$

Therefore, if $\frac{S_{\mathrm{1}}}{S_{\mathrm{0}}}$ is larger than $\frac{V_{\mathrm{1}}}{V_{\mathrm{0}}}$ (the solid line is higher than the dashed line in Fig. 8b–d), then AMF₁ is larger than AMF₀. In other words, these two variables determine how a priori LNO_x affects the retrieval of NO₂.