Introducing new lightning schemes into the CHASER (MIROC) chemistry–climate model

. The formation of nitrogen oxides (NO x ) associated with lightning activities (hereinafter designated as LNO x ) is a major source of NO x . In fact, it is regarded as the dominant NO x source in the middle to upper troposphere. Therefore, improving the prediction accuracy of lightning and LNO x in chemical climate models is crucially important. This study implemented three new lightning schemes with the CHASER (MIROC) global chemical transport and climate model. The ﬁrst lightning scheme is based on upward cloud ice ﬂux (ICEFLUX scheme). The second one (the original ECMWF scheme), also adopted in the European Cen-tre for Medium-Range Weather Forecasts (ECMWF) forecasting system, calculates lightning ﬂash rates as a function of Q R (a quantity intended to represent the charging rate of collisions between graupel and other types of hydrometeors inside the charge separation region), convective available potential energy (CAPE), and convective cloud-base height. For

Abstract. The formation of nitrogen oxides (NO x ) associated with lightning activities (hereinafter designated as LNO x ) is a major source of NO x . In fact, it is regarded as the dominant NO x source in the middle to upper troposphere. Therefore, improving the prediction accuracy of lightning and LNO x in chemical climate models is crucially important. This study implemented three new lightning schemes with the CHASER (MIROC) global chemical transport and climate model. The first lightning scheme is based on upward cloud ice flux (ICEFLUX scheme). The second one (the original ECMWF scheme), also adopted in the European Centre for Medium-Range Weather Forecasts (ECMWF) forecasting system, calculates lightning flash rates as a function of Q R (a quantity intended to represent the charging rate of collisions between graupel and other types of hydrometeors inside the charge separation region), convective available potential energy (CAPE), and convective cloudbase height. For the original ECMWF scheme, by tuning the equations and adjustment factors for land and ocean, a new lightning scheme called the ECMWF-McCAUL scheme was also tested in CHASER. The ECMWF-McCAUL scheme calculates lightning flash rates as a function of CAPE and column precipitating ice. In the original version of CHASER (MIROC), lightning is initially parameterized with the widely used cloud-top height scheme (CTH scheme). Model evaluations with lightning observations conducted using the Lightning Imaging Sensor (LIS) and Optical Transient Detector (OTD) indicate that both the ICEFLUX and ECMWF schemes simulate the spatial distribution of lightning more accurately on a global scale than the CTH scheme does. The ECMWF-McCAUL scheme showed the highest prediction accuracy for the global distribution of light-ning. Evaluation by atmospheric tomography (ATom) aircraft observations (NO) and tropospheric monitoring instrument (TROPOMI) satellite observations (NO 2 ) shows that the newly implemented lightning schemes partially facilitated the reduction of model biases (NO and NO 2 ), typically within the regions where LNO x is the major source of NO x , when compared to using the CTH scheme. Although the newly implemented lightning schemes have a minor effect on the tropospheric mean oxidation capacity compared to the CTH scheme, they led to marked changes in oxidation capacity in different regions of the troposphere. Historical trend analyses of flash and surface temperatures predicted using CHASER (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020) show that lightning schemes predicted increasing trends of lightning or no significant trends, except for one case of the ICEFLUX scheme, which predicted a decreasing trend of lightning. The global lightning rates of increase during 2001-2020 predicted by the CTH scheme were 17.69 % • C −1 and 2.50 % • C −1 , respectively, with and without meteorological nudging. The un-nudged runs also included the short-term surface warming but without the application of meteorological nudging. Furthermore, the ECMWF schemes predicted a larger increasing trend of lightning flash rates under the short-term surface warming by a factor of 4 (ECMWF-McCAUL scheme) and 5 (original ECMWF scheme) compared to the CTH scheme without nudging. In conclusion, the three new lightning schemes improved global lightning prediction in the CHASER model. However, further research is needed to assess the reproducibility of trends of lightning over longer periods.

Introduction
Nitric oxide (NO) can be formed during lightning activities. Also, NO can be oxidized quickly to nitrogen dioxide (NO 2 ). An equilibrium between NO and NO 2 can be reached during daytime. Those gases are known collectively as NO x (Finney et al., 2014). Actually, LNO x is estimated as contributing approximately 10 % of the global NO x source. Regarded as the dominant NO x source in the middle to upper troposphere (Schumann and Huntrieser, 2007;Finney et al., 2016b), NO x is associated with many chemical reactions in the atmosphere. Most importantly, NO reacts with the peroxy radical to reproduce the OH radical. Photochemical dissociation of NO 2 engenders the production of ozone (Isaksen and Hov, 1987;Grewe, 2007). The primary oxidants in the atmosphere, which are the OH radical and ozone, control the oxidation capacity of the atmosphere. Results of several studies have indicated that global-scale LNO x emissions are an important contributor to ozone and other trace gases, especially in the upper troposphere (Labrador et al., 2005;Wild, 2007;Liaskos et al., 2015). Consequently, LNO x influences atmospheric chemistry and global climate to a considerable degree (Schumann and Huntrieser, 2007;Murray, 2016;Finney et al., 2016a;Tost, 2017). However, large uncertainties remain in predicting lightning and LNO x in chemical climate models (Tost et al., 2007). Therefore, improving lightning prediction accuracy and quantifying LNO x in chemical climate models is crucially important for future atmospheric research.
Global chemical climate models (CCMs) such as CHASER (MIROC) (Sudo et al., 2002;Sudo and Akimoto, 2007;Watanabe et al., 2011) most often use the convective cloud-top height to parameterize the lightning flash rate (Price and Rind, 1992;Lamarque et al., 2013). The Earth system models (ESMs) recently used in the sixth Coupled Model Intercomparison Project (CMIP6) all used the convective cloud-top height to calculate the lightning flash rates (Thornhill et al., 2021). Not only with global CCMs but also studies of LNO x with regional-scale models have made significant progress in recent years (Heath et al., 2016;Kang et al., 2019aKang et al., , b, 2020. The spaceborne Lightning Imaging Sensor (LIS) and Optical Transient Detector (OTD) lightning observation data (Cecil et al., 2014) are often utilized to evaluate the performance of different lightning schemes. A new lightning scheme proposed by Finney et al. (2014), which is based on upward cloud ice flux, has shown better spatial and temporal correlation coefficients as well as root mean square errors (RMSEs) than the cloud-top height scheme compared against the LIS/OTD lightning observations. Another lightning scheme also showed more accurate lightning prediction than the cloud-top height scheme, which is also adopted in the ECMWF forecasting system (Lopez, 2016). This lightning scheme uses Q R (a quantity intended to represent the charging rate of collisions between graupel and other types of hydrometeors inside the charge separation region), con-vective available potential energy (CAPE), and convective cloud-base height to compute the lightning flash rate (Lopez, 2016). The two new lightning schemes (Finney et al., 2014;Lopez, 2016) mentioned above have only been evaluated in a few chemical transport and climate models. The new lightning schemes are expected to be evaluated and compared in more chemical transport and climate models, such as CHASER. To achieve better prediction accuracy for lightning and better quantification of LNO x in chemical climate models, comparing and optimizing the existing lightning schemes and evaluating them with various observation data are also important.
Lightning simulations are also fundamentally important in chemical climate model studies for predictions of atmospheric chemical fields and climate. Nevertheless, different lightning schemes respond very differently on decadal to multi-decadal timescales under global warming. Some lightning schemes such as those using cloud-top height or CAPE × precipitation rate as a proxy for calculating lightning indicate that lightning increases concomitantly with increasing global average temperature. By contrast, other lightning schemes, such as those using convective mass flux or upward cloud ice flux as a proxy for lightning, indicate that lightning will decrease as the global average temperature increases (Clark et al., 2017;Finney et al., 2018). Several studies (Price and Rind, 1994;Zeng et al., 2008;Jiang and Liao, 2013;Banerjee et al., 2014;Krause et al., 2014;Clark et al., 2017) have found 5 %-16 % increases in lightning flashes per degree of increase in global mean surface temperatures with the lightning scheme based on cloud-top height. Over the contiguous United States (CONUS), the CAPE × precipitation rate proxy predicted a 12 ± 5 % increase in the CONUS lightning flash rate per degree of global mean temperature increase (Romps et al., 2014). Compared to the findings reported by Romps et al. (2014), Finney et al. (2020 found a relatively small response of lightning to climate change (2 % K −1 ) over Africa using a cloud-ice-based parameterization for lightning. By contrast, Finney et al. (2018) found a 15 % global mean lightning flash rate decrease with the lightning scheme based on upward cloud ice flux in 2100 under a strong global warming scenario. Furthermore, a 2.0 % decrease in global mean lightning flashes per degree of increase in the global mean surface temperature with the lightning scheme based on convective mass flux has been reported by Clark et al. (2017). Although it remains unclear which lightning scheme is best suited to predicting future lightning (Romps, 2019), comparing different lightning schemes in different chemical climate models is valuable for consideration of the sensitivity of lightning to global warming.
This study introduced three new lightning schemes into CHASER (MIROC). The first lightning scheme (Finney et al., 2014) is based on upward cloud ice flux. The second one (Lopez, 2016), also adopted in the ECMWF forecasting system, calculates lightning flash rates as a function of Q R (defined in Sect. 2.2), CAPE, and convective cloud-base height. In the case of the second lightning scheme, by tuning the equations and adjustment factors based on a study reported by McCaul et al. (2009), a new lightning scheme called the ECMWF-McCAUL scheme was also tested for CHASER (MIROC). The ECMWF-McCAUL scheme calculates lightning flash rates as a function of CAPE and column precipitating ice. Our study conducted a detailed evaluation of lightning and LNO x by LIS/OTD lightning observations, NASA/ATom aircraft observations, and TROPOMI satellite observations. The effects of different lightning schemes on the atmospheric chemical fields were evaluated. Also, 20year (2001-2020) historical trend analyses of lightning densities and LNO x emissions simulated by different lightning schemes were conducted. Based on the results, the effects of LNO x emissions during 2001-2020 on tropospheric NO x and O 3 column trends were estimated and discussed.
Research methods, including the model description and experiment setup, are described in Sect. 2. In Sect. 3.1, the evaluation of lightning schemes using LIS/OTD lightning observations is explained. In Sect. 3.2, LNO x emission simulation by different lightning schemes is evaluated with aircraft and satellite measurements. Section 3.3 presents a discussion of the effects of different lightning schemes on the atmospheric chemical fields. Historical trends of lightning simulated by different lightning schemes are analyzed and discussed in Sect. 3.4. Section 3.5 discusses how LNO x emission (2001-2020) trends influence the tropospheric NO x and O 3 column trends. Section 4 presents the discussion and conclusions of this study.

Chemistry-climate model
The model used for this study is the CHASER (MIROC) global chemical transport and climate model (Sudo et al., 2002;Sudo and Akimoto, 2007;Watanabe et al., 2011;Ha et al., 2021), which incorporates consideration of detailed chemical and transport processes in the troposphere and stratosphere. CHASER calculates the distributions of 94 chemical species and reflects the effects of 269 chemical reactions (58 photolytic, 190 kinetic, 21 heterogeneous). Its tropospheric chemistry incorporates consideration of non-methane hydrocarbon (NMHC) oxidation and the fundamental chemical cycle of O x -NO x -HO x -CH 4 -CO. Its stratospheric chemistry simulates chlorinecontaining and bromine-containing compounds, chlorofluorocarbons (CFCs), hydrofluorocarbons (HFCs), carbonyl sulfide (OCS), and N 2 O. Furthermore, it simulates the formation of polar stratospheric clouds (PSCs) and heterogeneous reactions on their surfaces. CHASER is coupled to the MIROC AGCM version 5.0 (Watanabe et al., 2011). Grid-scale large-scale condensation and cumulus convection (Arakawa-Schubert scheme) are used to simulate cloud and precipitation processes. Aerosol chemistry is coupled with the SPRINTARS aerosol model (Takemura et al., 2009), which is also based on MIROC.

Lightning NO x emission scheme
The lightning flash rate in CHASER is originally parameterized by cloud-top height Rind, 1992, 1993), with a "C-shaped" NO x vertical profile adopted (Pickering et al., 1998). The equations used to calculate the lightning flash rate by cloud-top height are where F represents the total flash frequency (fl. min −1 ), H stands for the cloud-top height (km), and subscripts "l" and "o" respectively denote the land and ocean (Price and Rind, 1992). For this study, three new lightning schemes are implemented into CHASER (MIROC). One is based on upward cloud ice flux. It calculates the lightning flash rate by the following equations, as described by Finney et al. (2014).
Therein, f l and f o respectively represent the flash density (fl. m −2 s −1 ) over land and ocean. Also, φ ice is the upward cloud ice flux at 440 hPa (kg ice m −2 cloud s −1 ) as calculated using where q denotes the specific cloud ice water content at 440 hPa (kg ice kg −1 air ), mass stands for the updraft mass flux at 440 hPa (kg air m −2 cell s −1 ), and c represents the fractional cloud cover at 440 hPa (m 2 cloud m −2 cell ). The 440 hPa pressure level is chosen because it is a representative pressure level of fluxes in deep convective clouds (Finney et al., 2014). Moreover, Romps (2019) has proposed an alternative approach to 5630 Y. He et al.: Introducing new lightning schemes into the CHASER (MIROC) chemistry-climate model applying the ICEFLUX scheme by using the upward cloud ice flux at 260 K isotherms instead of at 440 hPa isobars. As suggested by Romps (2019), the isotherm alternative is more appropriate for climate change simulations because the charge separation zone will follow the isotherms instead of the isobars with climate change. The 260 K isotherm is chosen because it is close to the 440 hPa isobar based on a present-day tropical sounding and it lies within the mixedphase regions of clouds (Romps, 2019). To distinguish the two different approaches to applying the ICEFLUX scheme, the isobar approach is abbreviated as ICEFLUX_P and the isotherm alternative is abbreviated as ICEFLUX_T.
Another new lightning scheme, also adopted in the ECMWF forecasting system, calculates lightning flash rates as a function of the Q R (defined in Eq. 8), CAPE, and convective cloud-base height (Lopez, 2016) as where f T is the total lightning flash density (fl. m −2 s −1 ), z base is the convective cloud-base height in meters, and α (fl. kg −1 m −3 ) is a constant obtained after calibration against the LIS/OTD climatology, which is set to 1.11×10 −15 in this study. As explained by Lopez (2016), the number 1800 used in Eq. (6) is a constraint to let the term z base remain constant after it exceeds 1800 m. Note that Eq. (6) is standardized on base SI units. CAPE (m 2 s −2 ) is diagnosed from the following equation.
In that equation, g is the constant of gravity. Also, T u v and T v respectively denote the virtual temperatures in the updraft and the environment. The integral in Eq. (7) starts at the level of free convection z LFC and stops at the level at which negative buoyancy is found (w = 0). Quantity Q R (kg m −2 ) is intended to represent the charging rate of collisions between graupel and other types of hydrometeors inside the charge separation region. It is empirically calculated as where z 0 and z −25 signify the heights (m) of the 0 and −25 • C isotherms, and q cond denotes the mass mixing ratio of cumulus cloud liquid water (kg kg −1 ). The respective amounts of graupel (q graup ; kg kg −1 ) and snow (q snow ; kg kg −1 ) are computed from the following equations for each vertical level of the model.
In those equations, P f denotes the vertical profile of the frozen precipitation convective flux (kg m −2 s −1 ), ρ stands for the environmental air density (kg m −3 ), and V graup and V snow respectively express the typical fall speeds for graupel and snow set to 3.0 and 0.5 m s −1 . The dimensionless coefficient β is set as 0.7 for land and 0.45 for ocean to account for the observed lower graupel contents over oceans. For the original ECMWF scheme, by tuning the calculation equations based on findings reported by McCaul et al. (2009) and the adjustment factors for land and ocean, the lightning prediction accuracy is improved further, as explained in Sect. 3.1. We named the new lightning scheme the ECMWF-McCAUL scheme, and it simulates the lightning flash rate by the following equations.
Therein, f l and f o respectively denote the total flash density (fl. m −2 s −1 ) over land and ocean. Also, α l and α o are constants (fl. s 1.6 kg −1 m −2.6 ) obtained after calibration against LIS/OTD climatology, respectively, for land and ocean. For this study, α l and α o are respectively set to 2.67 × 10 −16 and 1.68 × 10 −17 . Then CAPE is computed in the same way as the original ECMWF scheme. In addition, Q Ra (kg m −2 ) is a proxy for the charging rate resulting from the collisions between graupel and hydrometeors of other types inside the charge separation region (from 0 to −25 • C isotherm), as reported by McCaul et al. (2009). Also, Q Ra represents the total volumetric amount of precipitating ice in the charge separation region, calculated as where q graup , q snow , and q ice respectively represent the mass mixing ratios (kg kg −1 ) of graupel, snow, and cloud ice. In this study, q graup and q snow were respectively computed by Eqs. (9) and (10). For the ECMWF-McCAUL scheme, V graup and V snow are respectively set to 3.1 and 0.5 m s −1 . Then q ice was diagnosed using the Arakawa-Schubert cumulus parameterization. Table 1 presents all the lightning schemes examined for this study. As described in this paper, the original ECMWF scheme and the ECMWF-McCAUL scheme are collectively designated as ECMWF schemes. Based on recent studies, the intra-cloud (IC) lightning flashes are as efficient as the cloudto-ground (CG) lightning flashes in NO x generation, and the A fourth-order polynomial is used to calculate the proportion of total flashes that are cloud-to-ground (p) based on the cold cloud depth, as described in an earlier report (Price and Rind, 1993).
In that equation, D represents the depth of cloud above the 0 • C isotherms in kilometers.

Lightning observations
The LIS/OTD gridded climatology datasets are used for this study, consisting of climatologies of total lightning flash rates observed using the Lightning Imaging Sensor (LIS) and spaceborne Optical Transient Detector (OTD): OTD aboard the MicroLab-1 satellite and LIS aboard the Tropical Rainfall Measuring Mission (TRMM) satellite (Cecil et al., 2014). Both sensors detect lightning by monitoring pulses of illumination produced by lightning in the 777.4 nm atomic oxygen multiplet above background levels. Both sensors, in low Earth orbit, view an Earth location for about 3 min as OTD passes overhead or for 1.5 min as LIS passes overhead. Actually, OTD and LIS circle the globe 14 times a day and 16 times a day, respectively. OTD collected data between +75 and −75 • latitude from May 1995 through March 2000, whereas LIS observed between +38 and −38 • latitude from January 1998 through April 2015. The product used throughout this paper is the LIS/OTD 2.5 • Low Resolution Time Series (LRTS). The LRTS includes the daily lightning flash rate on a 2.5 • regular latitude-longitude grid from May 1995 through April 2015.

Atmospheric tomography (ATom) aircraft observations
To evaluate the LNO x emissions calculated by different lightning schemes, we used NO observations by the atmospheric tomography (ATom) aircraft missions (Wofsy et al., 2018). By deploying an extensive gas and aerosol payload on the NASA DC-8 aircraft, ATom is designed to sample the atmosphere systematically on a global scale, performing continuous profiling from 0.2 to 12 km of altitude. Flights took place in each of the four seasons of 2016 through 2018. Since most of the lightning occurs over land regions during summer, ATom1 (July-August 2016) and ATom2 (January-February 2017) were used to evaluate LNO x emissions (corresponding to summer in the Northern and Southern Hemisphere, respectively). Both ATom1 and ATom2 originate from the Armstrong Flight Research Center in Palmdale, California, USA, fly north to the western Arctic, south to the South Pacific, east to the Atlantic, north to Greenland, and return to California across central North America. Figure 1 exhibits the respective flight tracks of ATom1 and ATom2. To evaluate the model-simulated NO against the ATom observations, we have sampled the specific flight track and timings from the modeled data.

TROPOMI satellite observations
The Tropospheric Monitoring Instrument (TROPOMI) is the payload aboard the Sentinel-5 Precursor (S5P) satellite of the European Space Agency (ESA), which was launched in October 2017. TROPOMI has been providing observations of important atmospheric pollutants (NO 2 , O 3 , CO, CH 4 , SO 2 , CH 2 O) with an unprecedented horizontal resolution of approximately 7 × 3.5 km 2 since August 2017 (changed to 5.5 × 3.5 km 2 after August 2019) (Verhoelst et al., 2021). For the direct comparisons between TROPOMI level-2 products and CHASER results, the following procedures were conducted to pre-process the TROPOMI data and CHASER modeled fields.
1. The TROPOMI retrievals with quality assurance (QA) values larger than or equal to 0.75 were selected.
3. The modeled results were sampled based on the TROPOMI overpass time. The CHASER tropospheric NO 2 columns were calculated by using the sampled modeled results, the averaging kernels retrieved from the TROPOMI retrievals, and the temperature and pressure profiles provided by TROPOMI retrievals. The averaging kernels are applied to each layer of the CHASER outputs following the Eq. (16).
4. The pre-processed data described above were used to produce the monthly averaged data.

OMI satellite observations
The Ozone Monitoring Instrument (OMI) is a key instrument aboard NASA's Aura satellite for measuring criteria pollutants such as O 3 , NO 2 , SO 2 , and aerosols. OMI has been providing observations with spatial resolution varying from 13 km × 25 km to 26 km × 128 km since October 2004 (Goldberg et al., 2019). The NO 2 product used in this study is the level-3 daily global gridded (0.25 • × 0.25 • ) nitrogen dioxide product (OMNO2d) . The O 3 product used in this study is the monthly mean tropospheric column O 3 product developed from OMI in combination with the Aura Microwave Limb Sounder (MLS) with the detailed method described by Ziemke et al. (2006).

Experiment setup
For this study, all the introduced lightning schemes were implemented into CHASER (MIROC). Six sets of experiments were conducted for this study, and the detailed settings of all experiments are presented in Table 2. For each set of experiments, the same initial conditions and chemical emissions were used except for LNO x emissions. The set of experiments that applied meteorological nudging also has the same meteorological conditions. The monthly varying soil NO x emissions used are constant each year for all experiments derived from Yienger and Levy (1995). All experiments used the "backward C-shaped" LNO x vertical profile (Ott et al., 2010). The LNO x PE values of IC and CG used in all experiments are set to the same value (250 mol per flash), which is based on the recent literature (Ridley et al., 2005;Cooray et al., 2009;Ott et al., 2010;Allen et al., 2019). It is noteworthy that there are still large uncertainties in determining the LNO x PE values Bucsela et al., 2019), and the choice of different LNO x PE values may influence the simulated LNO x emissions and chemical fields. A more sophisticated parameterization of LNO x PE values needs to be implemented and verified in chemistry-climate models in future research.
The first set of experiments was conducted for the years 2001-2020. It was used to evaluate the distribution of the lightning flash rate against LIS/OTD lightning observations and to derive the historical lightning trend. The second set of experiments is the same as the first set of experiments but uses daily mean LNO x emission rates of 2001 calculated using lightning schemes for each year. This set of experiments is used to produce results for comparison with those of the first set of experiments to estimate the effects of LNO x emission trends on tropospheric NO x and O 3 column trends. The third set of experiments gives results for 2011-2020. These experiments are used to estimate the effects of different lightning schemes on atmospheric chemical fields. To normalize the different annual LNO x emission amounts by different lightning schemes, temporally and spatially uniform adjustment factors were applied to adjust the mean LNO x production (2011-2020) to 5.0 Tg N yr −1 . Note that the 10-year (2011-2020) mean LNO x production was adjusted to 5.0 Tg N yr −1 , but the LNO x production in each year is not exactly 5.0 Tg N yr −1 . This adjustment was achieved by first conducting the simulations without any adjustment; the 2011-2020 mean LNO x production (P LNO x ) was calculated, and then the corresponding adjustment factor (adj_factor) can be calculated by using the following equation.
Similarly, we also adjusted the LNO  (1996)(1997)(1998)(1999)(2000) within a broader range of ±75 • latitude. We have evaluated the potential uncertainties associated with the inconsistency of the time period of simulated lightning and observed lightning (2007-2011 and 1996-2000). The statistical analysis between LIS (2007-2011) and LIS/OTD (1996-2000 within ±38 • latitude exhibits an extremely high spatial correlation coefficient (R = 0.99) and relatively small relative bias (0.65 %), which supports the reasonability of comparing model results with the observation data within different time ranges.
The distribution of lightning observed by LIS/OTD and simulated by CHASER (MIROC) with different lightning schemes is depicted in Fig. 2. Figure 2 shows that lightning over the ocean is not well reproduced by the original CTH scheme. Actually, it is improved considerably by the new lightning schemes. Compared with the CTH scheme, the original ECMWF scheme better represents the lightning distribution in South Asia including the Indian region. The ECMWF schemes and the ICEFLUX_P scheme reduced negative biases in North America compared to the CTH scheme. In Australia, the ECMWF schemes better simulate the horizontal distribution of lightning. All lightning schemes failed to capture the worldwide maximum value found over the Congo Basin, although all lightning schemes captured the active region in central Africa.
To directly estimate the prediction accuracy of all lightning schemes, the Taylor diagrams are displayed in Fig. 3. In Fig. 3a, the overall prediction accuracy of the ICE-FLUX_P and original ECMWF schemes evaluated against the LIS 2007-2011 lightning climatology is slightly improved compared to the CTH scheme. This improvement is more obvious when considering land and ocean separately (Fig. 3b and c). In the case of Fig. 3a-c, the ECMWF-McCAUL scheme has shown the best prediction accuracy among all lightning schemes. In Fig. 3d, comparison of the annual mean lightning flash rate of LIS/OTD 1996-2000 and the CHASER calculation for 2007-2011 yields spatial correlation coefficients of 0.80 and 0.79 for the ICE-FLUX_P and original ECMWF schemes, respectively, which are slightly higher than that found for the CTH scheme (0.78). The overall RMSE of the ICEFLUX_P scheme is 3.32 fl. km −2 yr −1 , which is slightly less than that of the CTH scheme of 3.44 fl. km −2 yr −1 . Among all lightning schemes, the ECMWF-McCAUL scheme exhibits the highest spatial correlation coefficient (0.83) and the lowest RMSE (3.20 fl. km −2 yr −1 ) as depicted in Fig. 3d. As displayed in Fig. 2, the prediction accuracy of lightning over the ocean is significantly improved, which can also be verified in Fig. 3f.
To estimate whether the improvement of prediction accuracy discussed in Fig. 3 is significant, a significance test is conducted by considering the uncertainties in the LIS/OTD observations. Based on the uncertainties in the LIS/OTD observations, the probability density distributions (PDDs) of spatial correlation coefficients (R) and RMSE between the model and observations are derived by using a Monte Carlo method and displayed in Fig. 4. The uncertainties in the LIS/OTD observations are determined based on the uncertainties of the instrument bulk flash detection efficiency of LIS (88 ± 9 %) and OTD (54 ± 8 %) (Boccippio et al., 2002). The R and RMSE shown in Fig. 4 are all normally distributed, which is determined by the Kolmogorov-Smirnov test. Based on the probability density functions of R and RMSE derived from Fig. 4, the order of R between the model and observations is estimated to be ECMWF-McCAUL > ICEFLUX_P > ECMWF-original > CTH with a confidence limit larger than 99.9 %. Moreover, the order of RMSE between the model and observations is estimated to be ECMWF-McCAUL < ICEFLUX_P < ECMWF-original and CTH with a confidence limit larger than 95 %. According to the significance test described above, we can conclude that the newly implemented lightning schemes have improved the lightning prediction accuracy compared to the original CTH scheme. Figure 5 displays a comparison of seasonal and annual meridional average lightning flash densities from simulations ( -2011( ) and LIS satellite observations (2007( -2011. As Fig. 5 shows, the pairs of curves are usually in good agreement, even though the annual plot (Fig. 5e) a Nudging off means the meteorological fields (u, v, T ) are free-running instead of nudging towards the NCEP FNL data. b LNOx is interactively calculated by using different lightning schemes. c The climate change is simulated by prescribed sea surface temperature (SST) and sea ice fields as well as prescribed varying concentrations of GHGs (CO 2 , N 2 O, methane, chlorofluorocarbons -CFCs, and hydrochlorofluorocarbons -HCFCs) utilized only in the radiation scheme. The SST and sea ice fields are obtained from the HadISST dataset (Rayner et al., 2003).  have made improvements within Africa. Also, the ICE-FLUX_P and the original ECMWF schemes have slightly reduced the biases over North America. A noticeable underestimation over the Americas in JJA and overestimation in MAM can be respectively observed in Fig. 5c and b. Lightning densities over Africa are generally underestimated to varying degrees in different seasons, with the greatest underestimation occurring in JJA (Fig. 5c). Lightning densities over Asia (from 60 to 120 • E) are slightly underestimated in MAM (Fig. 5b). The ICEFLUX_P scheme has reduced the biases. Figure 6 is the same as Fig. 5, but for the zonal mean distributions. The curves of the model results and the observation results in Fig. 6 generally show good agreement. Figure 6f shows that, overall, the ICEFLUX_P and the ECMWF-McCAUL schemes slightly overestimated the lightning densities near the Equator (10 • S-10 • N). All lightning schemes underestimated the lightning densities within 15-38 • N and 20-38 • S. Figure 6f also shows that the ICEFLUX_P scheme has reduced the biases within 10-38 • N and 15-38 • S compared to the CTH scheme. In DJF (Fig. 6a), all lightning schemes overestimated the flash densities over the lowlatitude regions but slightly underestimated the flash densi-ties over the middle-latitude regions in the Southern Hemisphere. In MAM (Fig. 6b), lightning densities are overestimated near the Equator and underestimated over 15-38 • N and 15-38 • S by all lightning schemes to varying degrees. In JJA (Fig. 6c), noticeable overestimation around 10 • N by the original ECMWF scheme is apparent. Moreover, the CTH and the original ECMWF schemes respectively facilitated reduction of model biases over 15-38 • S and 15-38 • N. As Fig. 6d shows, the model-predicted lightning maximum value is shifted approximately 15 • to the north in SON compared to the lightning observations. Figure 6d also shows that all lightning schemes underestimated the lightning densities over 15-38 • N and 0-38 • S. The ICEFLUX_P scheme has shown improvement over these regions.

Evaluation of LNO x emissions by ATom1 and ATom2 observations
To evaluate the LNO x emissions of different lightning schemes, we used ATom1 and ATom2 aircraft measurements (NO) for comparison against model results. All lightning schemes, when implemented in CHASER, produce flash  Table 3 presents model biases of different lightning schemes against the ATom1 and ATom2 observations. Figure 7 displays the vertical profile of biases divided by the ATom observations in percentage terms. In Table 3 and Fig. 7, case ZERO is the case with the lightning flash, with LNO x emissions completely switched off. Comparisons between model results and ATom observations were conducted within two specific regions (South America region and western Pacific region) in which LNO x is the major source of NO x (Fig. 8). As Table 3 and Fig. 7 show, the model generally tends to underestimate the NO concentrations. The model biases are reduced considerably by including lightning NO x sources. For ATom1, overall, the ICEFLUX_P scheme has the smallest model bias. The original ECMWF scheme also reduced the model biases compared to the CTH scheme ( Table 3). The ICEFLUX_P and the ECMWF-McCAUL schemes reduced the model biases substantially within 0-30 • N latitude where the lightning activities are most dominant during the ATom1 observation period (29 July to 23 Au-gust 2016). In the range of 30 • S to 80 • N in ATom1, overall the ICEFLUX_P scheme reduced the model biases considerably and the ECMWF schemes slightly reduced or extended the model biases compared to the CTH scheme (Table 3, Fig. 7b-e). However, in the range of 30-80 • S, the model biases were extended by the ICEFLUX_P and the ECMWF schemes compared to the CTH scheme (Table 3, Fig. 7f and g).
For ATom2, overall, the ECMWF schemes slightly reduced the model biases over the upper troposphere compared to the CTH scheme (Fig. 7j). During the ATom2 observation period (26 January to 21 February 2017), the lightning activities are most dominant within the range of 0-30 • S, where the model biases were reduced significantly by newly implemented lightning schemes. A hotspot of lightning activities during the ATom2 observation period is the South America region, where the model biases were reduced dramatically by the ECMWF schemes. The model biases were mostly reduced by the newly implemented lightning schemes within the low-latitude and middle-latitude regions but slightly extended within the high-latitude regions. The model biases were mostly reduced or extended over the middle to upper troposphere (Fig. 7). This is true because most LNO x was distributed over the middle to upper troposphere. Also, NO x has a longer lifetime over the middle to upper troposphere. In the western Pacific region, results obtained from comparisons with ATom1 and ATom2 indicate that all lightning schemes overestimated LNO x emissions in the upper troposphere; also, both the ICEFLUX_P scheme and ECMWF schemes reduced the total model biases considerably more than the CTH scheme did.

Evaluation of LNO x emissions by TROPOMI satellite observations
TROPOMI satellite observations of tropospheric NO 2 columns were used to evaluate LNO x emission results ob-tained using the CHASER model. To eliminate differences in annual total LNO x emissions attributable to the different lightning schemes, we adjusted the annual LNO x emissions of all lightning schemes to 5.0 Tg N yr −1 using different adjustment factors. For direct comparison between CHASER and TROPOMI tropospheric NO 2 columns, the averaging kernel information from TROPOMI observations was used.
In that equation, X chaser represents the CHASER tropospheric NO 2 column after averaging kernels applied, A tropomi denotes the TROPOMI averaging kernels, x chaser denotes the CHASER NO 2 partial column at layer i, and N denotes the number of tropospheric layers. Comparison between TROPOMI observations and CHASER outputs indicates that the CHASER model tends to underestimate tropospheric NO 2 columns. Overall, the newly implemented lightning schemes have not shown improvements of model biases of tropospheric NO 2 columns at an annual global scale. To minimize the uncertainties of model biases of tropospheric NO 2 columns caused by other   (Fig. 9), where LNO x is the major source of NO x (as shown in Fig. 8). Figure 10 presents a comparison of smoothed CHASER and TROPOMI tropospheric NO 2 columns over four target regions in 2019. The spatial average values of each month in 2019 are shown in Fig. 10. That figure shows that the model generally captured the temporal variation of tropospheric NO 2 columns in the four regions. Actually, the temporal variations of modeled tropospheric NO 2 columns are close to each other. For the Amazon region, lightning activities are most dominant during MAM and SON, when the ECMWF-McCAUL scheme has shown noticeable improvements in model biases (Fig. 10a). Figure 10b reveals that all the newly implemented schemes slightly reduced the model biases with the original ECMWF scheme, showing the smallest model biases during the most prevailing season of lighting (DJF). Figure 10c is for the South Atlantic region where the most prevailing season of lighting is also DJF. Figure 10c shows that the ECMWF schemes slightly reduced the model biases compared to the CTH scheme. Referring to Fig. 10d, the dominant season of lightning is JJA, when the ECMWF-original scheme considerably reduced the model biases and the ECMWF-McCAUL scheme also slightly reduced the model biases.

Effects of different lightning schemes on tropospheric chemical fields
In the tropospheric chemical field, LNO x has an important role. The LNO x effects on the tropospheric chemical fields vary along with differences in the horizontal distribution of LNO  were conducted with the 10-year mean LNO x production of all lightning schemes adjusted to 5.0 Tg N yr −1 (Sect. 2.4). The CTH scheme with a "backward C-shaped" profile is regarded as the base scheme. The effects of different lightning schemes on the atmospheric chemistry are calculated as shown in Eq. (17).
Therein, Impact ij represents the effects of the ith lightning scheme on the concentrations of target atmospheric component j . Also, LS ij denotes the concentrations of target atmospheric component j simulated by the ith lightning scheme. Base j stands for the concentrations of target atmospheric component j as simulated using the base scheme. Figure 11 presents the respective effects of the ECMWF-McCAUL, original ECMWF, and ICEFLUX_P schemes on the atmospheric chemical fields (NO x , O 3 , OH, CO) relative to the base scheme CTH. The ECMWF-McCAUL scheme led to an increase (approximate maximum is 12 %) in NO x concentration at low-latitude regions and a decrease (approximate maximum is 15 %) at middle-to high-latitude regions. In the case of the ECMWF-McCAUL scheme, the concentration of ozone and the OH radical mostly increased at lowlatitude regions and decreased at middle-to high-latitude regions in the Southern Hemisphere, which corresponds to the changing pattern of NO x . The effects of the original ECMWF scheme on the atmospheric chemical fields are similar to that of the ECMWF-McCAUL scheme. However, the original ECMWF scheme led to a higher total tropospheric CO burden compared to the ECMWF-McCAUL scheme. As Fig. 11 shows, the three lightning schemes led to a marked decrease in NO x , O 3 , and OH radical concentrations over the South Pole region. This decrease occurred because the lightning densities and the LNO x emissions simulated by the CTH scheme are markedly higher than those simulated using other lightning schemes at this latitude band (Fig. 2). Moreover, NO x can engender the formation of ozone and the OH radical. In the case of the ICEFLUX_P scheme, the concentrations of NO x , ozone, and the OH radical mostly increased in the tropics and decreased at middle-to high-latitude regions in the Southern Hemisphere.
Methane lifetime is an indicator reflecting the tropospheric oxidation capacity. The global mean tropospheric lifetime of methane against the tropospheric OH radical spanning 2011-2020 with the CTH, original ECMWF, ECMWF-McCAUL, and ICEFLUX_P schemes is respectively estimated as 9. 226, 9.299, 9.256, and 9.229 years. Compared to the CTH scheme, the ECMWF schemes led to a slight increase in methane's global mean tropospheric lifetime. In contrast, methane's global mean tropospheric lifetime simulated by the ICE-FLUX_P scheme is almost the same as that simulated by the CTH scheme. Although little difference exists in the total tropospheric oxidation capacity simulated by different lightning schemes, the ECMWF schemes and ICEFLUX_P scheme led to marked changes in oxidation capacity in different regions of the troposphere.

Historical trend analysis of lightning density
The accuracy of predicting the simulated lightning distribution under the current climate is only one aspect of lightning scheme evaluation. The ability of the lightning scheme to reproduce the trend of lightning under climate change is also an important factor. For this study, 20 years of (2001-2020) experiments were conducted to analyze the historical trends of lightning flash rates simulated using different lightning schemes (Sect. 2.4). Figure 12 shows the global anomaly of lightning flash rates calculated from the simulation results. Because nudging to meteorological reanalysis data cannot be used when predicting lightning trends under future climate changes, we also showed the results without nudging. The un-nudged runs also represented the short-term surface warming like the experiments with nudging. The only differences between the unnudged and nudged experiments are whether the meteorological fields are nudged towards the 6-hourly NCEP FNL data. We used the Mann-Kendall rank statistic to ascertain whether the lightning trends in Fig. 12 are significant (Hussain and Mahmud, 2019). From results of the Mann-Kendall rank statistic test (significance set as 5 %), all the trends in Fig. 12 were inferred to be significant except for the trends shown in Fig. 12a, e, and i. As Fig. 12 shows, all lightning schemes predicted increasing trends or no significant trends of lightning except the ICEFLUX_P scheme without nudging, which predicted a decreasing lightning trend. The isotherm alternative application of ICEFLUX (ICEFLUX_T) led to slightly enhanced lightning trends toward positive lighting trends compared to the ICEFLUX_P scheme. As explained by Romps (2019), the ICEFLUX_P approach is based on a fixed isobar, which makes it less convenient for climate change studies. Therefore, at least the lightning trends simulated by the ICEFLUX_T approach are expected to be closer to the real situation than the ICEFLUX_P approach.
As displayed in Fig. 12, the positive lightning trends are generally enhanced by application of meteorological nudging. A further investigation of the trends of CAPE during 2001-2020 reveals that the trends of global averaged CAPE are also enhanced by application of meteorological nudging. Higher CAPE means higher buoyancy in the updrafts, which led to the higher lightning densities calculated by the lightning schemes used in this study. It is worth noting that even though the meteorological fields (u, v, T ) of nudged simulations are expected to be closer to the real situations, we can-  not analogously deduce that the lightning trends predicted by the nudged runs are also closer to the real situations. This is because the predicted lightning trends are not only controlled by the meteorological fields but also controlled by many other physical processes (e.g., cumulus parameterization). Few studies have specifically examined the lightning trends predicted by the ECMWF schemes under shortterm surface warming. When nudging was not applied, the ECMWF schemes predicted the increasing trends of lightning flash rates under short-term surface warming by factors of 4 (ECMWF-McCAUL scheme) and 5 (original ECMWF scheme) compared to the CTH scheme (Table 4). Figure 13 shows a global map of changes in the lightning flash rate (% yr −1 ) during 2001-2020. In Fig. 13, the area in which the trend was found to be significant by the Mann-Kendall rank statistic test (significance level inferred for 5 %) is marked with hatched lines. As Fig. 13 shows, the distribution of trends simulated by the same lightning scheme is similar whether meteorological nudging was applied or not. As displayed in Fig. 13, in the Arctic region of the Eastern Hemisphere, both the CTH scheme and the ECMWF schemes showed an increasing trend of lightning. Earlier studies based on the World Wide Lightning Location Network (WWLLN) observations have indicated that lightning densities in the Arctic increase concomitantly with increasing global mean surface temperature (Holzworth et al., 2021). Earlier studies indicate that the results of the CTH scheme and the ECMWF schemes are reasonable for the Arctic region of the Eastern Hemisphere. In the highlatitude region of the Southern Hemisphere (60-70 • S), both the CTH scheme and the ECMWF schemes showed decreasing lightning trends. Lightning is rarely observed south of 60 • S (Kelley et al., 2018). Moreover, the trends of lightning in this region expected to occur with the short-term surface warming remain highly uncertain. In some parts of the northern Pacific Ocean, the ECMWF schemes and ICEFLUX scheme results showed increasing trends of lightning, which is consistent with results obtained from an earlier study (Petersen and Buechler, 2008). All schemes show decreasing trends for lightning flash rates in Indonesia, although only the ICEFLUX scheme explicitly passed the significance test. In the North Atlantic, all schemes showed increasing lightning trends. Only the CTH scheme failed the significance test.  Figure 14 presents trends of annual global LNO x emissions calculated from the simulation results (2001-2020) obtained using different lightning schemes. As Fig. 14 shows, the annual global LNO x emission trends correspond to the trends of lightning presented in Fig. 12. Similar to the trends found for lightning, the trends of annual global LNO x emissions are also increased by application of meteorological nudging. Only the ICEFLUX scheme simulated decreasing trends of annual global LNO x emissions. Figure 15 portrays trends of global mean tropospheric NO x columns calculated from the first and second set of experiments ( Table 2). As Figs. 14 and 15 depict, when the trends of annual global LNO x emissions are not strong (e.g., Fig. 14a), their effects on the trends of global mean tropospheric NO x columns are negligible. The marked increasing trends of annual global LNO x emissions (Fig. 14f-h) led to great increases (12.12 %-20.59 %) of the increasing trends of tropospheric NO x columns (Fig. 15f-h). In the case of the ICEFLUX_P scheme without nudging, because of the decreasing trends of LNO x emissions, the increasing trends of the tropospheric NO x columns decreased by around 10 %. Figure 16 is similar to the results shown in Fig. 15, but for tropospheric O 3 columns. Because NO x causes the formation of O 3 by the fundamental chemical cycle of O x -NO x -HO x , the trends of the global mean tropospheric O 3 columns are strongly affected by trends of the global mean tropospheric NO x columns. In some cases, the simulated trends of tropospheric O 3 columns are almost identical, as portrayed in Fig. 16a, b, e, i, and j, because the trends of tropospheric NO x columns simulated by the two sets of experiments are very similar (Fig. 15a, b, e, i, j). As Figs. 14 and 16 show, the marked increasing trends of annual global LNO x emissions led to increases in the increasing trends of tropospheric O 3 columns by around 15 % (Fig. 16f-h). In the case of ICE-FLUX_P without nudging, because of the decreasing trend of LNO x emissions, the increasing trend of the tropospheric O 3 columns decreased by around 10 % (Fig. 16d). Note that the tropospheric NO x or O 3 columns in 2001 simulated by the first set of experiments and the second set of experiments are not exactly the same. This is because the blue lines show results with interactively calculated LNO x emission rates (the time resolution is 10-30 min). But the orange lines show results calculated by reading daily mean input data for LNO x  emission rates, which inhibits interaction of LNO x with meteorology in the model.

Effects of LNO
In conclusion, because the ICEFLUX scheme predicts opposite trends of LNO x emissions from the other lightning schemes, they simulate opposite effects on the historical trends of global mean tropospheric NO x and O 3 columns. Furthermore, an evident trend of annual global LNO x emissions has a strong effect on the trend of global mean tropospheric NO x and O 3 columns.

Discussion and conclusions
Three new lightning schemes, the ICEFLUX, the original ECMWF, and the ECMWF-McCAUL schemes, were implemented into CHASER (MIROC), a global chemical climate model. Using LIS/OTD lightning observations as validation data, both the ICEFLUX_P and ECMWF schemes simulated the spatial distribution of lightning more accurately on a global scale than the CTH scheme did, and the lightning distribution in the ocean region was especially improved. The ECMWF-McCAUL scheme showed the highest prediction accuracy for the spatial distribution of lightning on a global scale. It is noteworthy that whilst the ice-based parameterizations showed superb prediction accuracy of lightning distribution under today's climate, they have greater uncertainties associated with inputs, especially regarding the microphysics scheme used (Charn and Parishani, 2021).
To verify the LNO x emissions of different lightning schemes, we used NO observations from ATom1 and ATom2. Overall, both the ICEFLUX_P scheme and the ECMWF schemes partially reduced model biases, typically over the dominant regions of lightning activities, compared to the CTH scheme. We also used TROPOMI tropospheric NO 2 columns to verify the LNO x emissions of different lightning schemes. Although the ICEFLUX_P and the ECMWF schemes have not shown improvements of model biases of tropospheric NO 2 columns at an annual global scale, they generally led to an obvious reduction of model biases in the prevailing seasons of lightning within the regions where LNO x is a dominant source of NO x . Several studies have pointed out that the TROPOMI data used in this study are biased negatively compared to airborne or ground-based observation data (Tack et al., 2021;Verhoelst et al., 2021;van Geffen et al., 2022). The TROPOMI data used are generally negatively biased and the simulated tropospheric NO 2 columns are underestimated compared to the TROPOMI observations. Therefore, the uncertainties that existed in the TROPOMI data have negligible impacts on the conclusions of our study.
Effects of the newly implemented lightning schemes on the tropospheric chemical fields are evaluated compared to the CTH scheme. Compared with the CTH scheme, the ECMWF schemes mainly led to a slight increase in NO x , ozone, and OH radical concentrations at low-latitude regions and a decrease at middle-latitude to high-latitude regions. Effects of the ICEFLUX_P scheme on the tropospheric chemical fields slightly differ from those of the ECMWF schemes. The ICEFLUX_P scheme mainly causes a slight increase in NO x , ozone, and OH radical concentrations from the tropics to the Northern Hemisphere and a decrease in the concentrations of the three chemical species in the Southern Hemisphere except the tropics. The commonality between the ECMWF schemes and the ICEFLUX_P scheme is that they both result in decreasing concentrations of NO x , ozone, and the OH radical at the middle-to high-latitude regions of the Southern Hemisphere. Although the newly implemented lightning schemes have little effect on the total oxidation capacity of the troposphere compared to the CTH scheme, they led to marked changes in oxidation capacity in different regions of the atmosphere.
This study also analyzed the historical trends of lightning simulated by different lightning schemes under short-term surface warming during 2001-2020. The Mann-Kendall rank statistic was used to ascertain whether the lightning trends were significant. Use of Mann-Kendall rank statistic tests revealed that all the simulated historical lightning trends are significant, except the CTH and ICEFLUX_T schemes without nudging and the ICEFLUX_P scheme with nudging, for significance at 5 %. All the lightning schemes predicted increasing lightning trends or no significant trends except the ICEFLUX_P scheme without nudging, which predicted a decreasing lightning trend. The ICEFLUX_T scheme predicted a decreasing trend without nudging even though the trend failed the significance test. If it is accepted that the noninductive charging mechanism is an appropriate basis for a lightning parameterization, then the implication is that in the future if cloud ice (and cloud ice fluxes) is reduced then electrical charging will be reduced too. This provides a line of scientific reasoning to explain why lightning may be reduced in the future. Moreover, findings showed that when nudging was not applied, the ECMWF schemes predicted an increasing trend of lightning flash rate under short-term surface warming by factors of 4 (ECMWF-McCAUL scheme) and 5 (original ECMWF scheme) compared to the CTH scheme. Although a considerable degree of uncertainty remains in determining the sensitivity of lightning activity to changes in surface temperature on the decadal timescale (Williams, 2005), the majority of past estimates show that the sensitivity tends average close to 10 % K −1 (Betz et al., 2009, p. 521). This value is most consistent with the lightning increase rate predicted by the ECMWF-McCAUL scheme without nudging in this study. Future research should be undertaken for specific examination of development of lightning schemes that both accurately predict the global distribution of LNO x and predict the changes in lightning that are expected to occur concomitantly with global climate change. Finally, we quantitatively estimated the LNO x emission effects on tropospheric NO x and O 3 column trends during 2001-2020. Results showed that a marked trend of annual global LNO x emissions significantly affects the trend of global mean tropospheric NO x and O 3 columns.
Code availability. The source code for CHASER to reproduce results in this work is obtainable from the repository at https://doi.org/10.5281/zenodo.5835796 (He et al., 2022).
Author contributions. YH introduced new lightning schemes into CHASER (MIROC) by adding new codes to CHASER (MIROC), conducted all simulations, interpreted the results, and wrote the paper. KS developed the model code, conceived of the presented idea, and supervised the findings of this work and the paper preparation. HMSH provided the TROPOMI data and the relevant codes to preprocess the TROPOMI data.