Impact of Lightning NO x Emissions on Atmospheric Composition and Meteorology in Africa and Europe

: NO x emissions from lightning have been added to the CHIMERE v2020r1 model using a parameterization based on convective clouds. In order to estimate the impact of these emissions on pollutant concentrations, two simulations, using the online coupled WRF-CHIMERE models with and without NO x emissions from lightning, have been carried out over the months of July and August 2013 and over a large area covering Europe and the northern part of Africa. The results show that these emissions modify the pollutant concentrations as well as the meteorology. The changes are most signiﬁcant where the strongest emissions are located. Adding these emissions improves Aerosol Optical Depth in Africa but has a limited impact on the surface concentrations of pollutants in Europe. For the two-month average we ﬁnd that the maximum changes are localized and may reach ± 0.5 K for 2 m temperature, ± 0.5 m s − 1 for 10 m wind speed, 10 W m − 2 for short wave radiation surface ﬂux, and 50 and 2 µ g m − 3 for dust and sea salt surface concentrations, respectively. This leads to maximum changes of 1 µ g m − 3 for surface concentrations of PM 2.5 .


Introduction
Lightning is an important source of NO x in the troposphere. This process has been studied for many years, mainly at the global scale, and the basics of the physics and chemistry are presented in review articles such as those in [1][2][3]. This source of NO x is estimated to represent 10% of all NO x emitted over the world, with hot spots at equatorial latitudes [4]. The relative part of NO in NO x is about 75-95% [5]. The main principle of NO x emissions is detailed in [6]: NO is produced when the heating due to the flash induces a thermal dissociation of O 2 then a recombination with dinitrogen N 2 , similar to what occurs during combustion processes. NO is quickly oxidized to form NO 2 and an equilibrium is reached between NO and NO 2 . In models, parameterizations consider an instantaneous equilibrium between the two chemical species. The impact of this production on atmospheric chemistry is detailed in [3].
To quantify these emissions, numerous processes must be understood and then parametrized. (i) First, it is necessary to estimate the number of flashes. For that, a link between the meteorology, the studied area, and the flashes must be established. (ii) It is important to know the type of flash: emissions will depend on the flash nature: Cloud-to-ground (CG) or Intra-cloud (IC). (iii) Then, it is necessary to estimate the amount of NO x emitted by one flash. (iv) Finally, it is also necessary to estimate the vertical distribution of the emissions along the flash (from the surface to the top of the clouds). All these points are very uncertain and very few observations are available to constrain the calculations [7].
For the calculation of flash number, several types of parameterizations exist, mainly based on three approaches. (i) The "Cloud Top Height" scheme (CTH) by [8], is based on a linear relationship between the cloud height and the frequency of flashes over the continent. (ii) The "convective precipitation" schemes by the authors of [9][10][11][12], for example, with the Community Multiscale Air Quality (CMAQ) model and (iii) the "ice flux" scheme by the authors of [13]. These schemes are used in many models, with some adaptations and scaling factors depending on model resolutions, the studied regions, and periods of the year. The uncertainty is so large that many studies are devoted to implement several schemes and compare them to each other, as in [14,15]. These comparisons rarely conclude on one better scheme, the number of flashes being always tuned after comparisons to observations such as the satellite observations of TRMM-LIS and OTD, presented in [16,17]. It is noticed that these observations themselves remain very uncertain [18].
The abundance of NO x nonlinearly changes O 3 and HO x concentrations depending on the chemical regimes [1]: close to the emitting sources (such as urban centers), freshly emitted NO x (mainly because of NO) can reduce O 3 by titration, while downwind emission sources, O 3 can be produced when NO x and VOCs are present with sunlight. Reactions with OH induce a NO x increase and cause a positive radiative forcing leading to a warming via O 3 . The NO x concentrations are also strongly dependent on their conversion into nitric acid (HNO 3 ) and particulate nitrate (NO − 3 ): it corresponds to the main sink for NO x through dry and wet deposition. The impact of NO x on ozone will also depend on the altitude of the emissions via the chemical regimes. As pointed out in [3], the NO x production in the free troposphere will increase concentrations of OH, NO − 3 , H 2 O 2 , and ozone. More atmospheric NO − 3 also increases its deposition. Lightning production has also an impact on aerosol, NO x being oxidized with OH to nitric acid (HNO 3 ), then condense into ammonium nitrate (NH 4 NO 3 ) in the presence of excess ammonia (NH 3 ), after neutralization of sulfate aerosol. The net impact is an increase downwind in aerosol due to faster oxidation rates of SO 2 to sulfate and NO x to nitrate. It also increases the conversion of DMS to SO 2 , then sulfate aerosol in marine environments. Regarding nitrate, complex interactions between chemistry, climate, and lightning are described in [19]. They first show that lightning NO x is an important source of nitrate in the upper troposphere. It has also an impact on gas-phase sulfate formation, then the aerosol size distribution, then the extinction properties of aerosol and clouds, then a decrease of shortwave radiation. On the other hand, some studies are describing the impact of aerosols on NO x production by lightning through direct and indirect areosol effects [20].
In this study, a simple scheme of NO x lightning emissions is implemented in the CHIMERE chemistry-transport model [21], coupled to the WRF regional meteorological model. This scheme, based on the CTH approach, is presented in Section 2, along with the rest of the modeling framework. To analyze the impact of these emissions in a regional model, several simulations are performed including those with or without lightning NO x emissions. Comparisons between the simulations are presented in Section 3.1 for ozone concentrations, Sections 3.2 and 3.3 for other meteorological and chemical variables with time-averaged maps, and Section 3.4 as time-series. Section 3.5 presents results as statistical scores. Discussions and conclusions are finally presented in Section 4.

Materials and Methods
In this section, the whole modeling system is presented. This includes the WRF and CHIMERE regional models and the way they are coupled. This also includes how the calculation of NO x lightning emissions was implemented and the design of the test case defined for this study.

The WRF and CHIMERE Models
The regional models WRF 3.7.1 [22] and CHIMERE 2020r1 [23] are used. This version of CHIMERE is designed to run in offline or online modes, corresponding to the simulation of the direct and indirect effects of aerosols on radiation and cloud properties. How these effects are taken into account is described in [24] (for the direct effects) and in [25] (for the indirect effects). The model configuration is exactly the same as in [26]: it includes emissions from anthropogenic, biogenic, sea salt, biomass burning, and mineral dust sources. The chemical evolution of gaseous species is calculated using the MELCHIOR2 scheme [27]. The photolysis rates and Aerosol Optical Depth (AOD) are calculated using the FastJX radiation module (version 7.0b) [28]. The aerosols are modeled using the scheme developed by [29]. This module takes into account sulfate, nitrate, ammonium, primary organic matter and elemental carbon (EC), secondary organic aerosols, sea salt, dust, and water. The aerosol size is represented using ten bins, from 40 nm to 40 µm, in mean mass median diameter. Their life cycle is completely represented with nucleation of sulfuric acid, coagulation, absorption, wet and dry deposition and scavenging (in-cloud and sub-cloud scavenging). The inorganic part constitutes the major part of the particulate matter in the fine mode (for D p < 2.5 µm). To determine the gas particle partitioning of these semi-volatile species, the ISORROPIA model is used [30]. At the boundaries of the regional domain, climatologies from global model simulations with LMDz-INCA [31] for all gaseous and aerosols species and GOCART for mineral dust [32]. The only difference with the previous version is the addition of lightning NO x emissions in CHIMERE. Note that a lightning parameterization already exists in WRF since version 3.5 [7,33]. However, only the calculation of the flash rate is done, not the NO x emission. To be consistent with the deep convection scheme implemented in CHIMERE, the choice was made for this study to implement the whole lightning NO x scheme in CHIMERE, from the flash rate estimation to the NO x emissions.

The Lighning No x Emissions
The lightning NO x emission calculation follows the scheme of [34]. It is added in the model following [13,35]. The prerequisite is to use a deep convection diagnostic scheme to evaluate the occurrence of convective clouds. If convective clouds are diagnosed, their top altitude is evaluated and directly used in the parameterization to estimate flash localization and frequency. This is the so-called Cloud Top Height (CTH) approach: with F c and F m the flash frequencies (# mn −1 25 km −2 ) [36], with c for continental and m for maritime.
H is the top of the convective cloud (in km). A scaling factor [37], designed to take into account the model grid cell size is applied as where ∆λ and ∆Φ are the cell dimensions in latitude and longitude respectively, expressed in degrees. In practice, except for very coarse resolutions, this factor evolves around 0.97 and has in practice a low impact on the results, compared to the uncertainties of all others parameters. These uncertainties are generally compensated by adjustment factors to bring the model closer to the observations, as described in [15]. This kind of CTH scheme was extensively tested (e.g., in [35]). For example, the authors of [13] used an additional scaling factor of α=0.05 to fit the number of flash after comparison to global yearly satellite data. For each model cell, the relative percentage of sea x sea is estimated using the land-sea mask from the land use database. Then, the flash frequency is estimated as with F the final number of flash, divided by the surface of the cell to have number per km 2 , then in (# mn −1 km 2 ). For each flash, it is needed to estimate the relative part "in cloud" (IC) and "cloud to ground" (CG). This estimate is first based on the cloud ice depth estimation, H f r , proposed by [34] and calculated as H f r (km) = −6.64 × 10 −5 |L| 2 − 4.73 × 10 −3 |L| + 7.34 (4) with L the latitude in degrees. Note that their fit was only provided for northern hemisphere. For this study, and in order to have a symmetrical relationship for the two hemispheres, the absolute value of latitude is used. H f r is used to calculate the β parameter: p is defined as the relative part of CG in the total (IC+CG) and is defined as Note that β is used with boundary values and must be 1 < β < 50. This corresponds to 5.5 km < H f r < 14 km and it enables to avoid unrealistic values. Finally, knowing the flash frequency and the vertical distribution of IC and CG for each one, the NO x concentration is estimated following the work in [5] as P(CG, NO) = 6.7 × 10 26 molecules/flash P(IC, NO) = 6.7 × 10 25 molecules/flash (7) This scheme is simple but also rather uncertain. In addition, it has been designed after simulations of several years at global scale: all constants and parameters have been tuned for other modeling conditions than those used in this study. For the vertical shape of the emission profiles, the simplest approach is implemented by considering the emission as constant in altitude: the CG flux from surface to H f r and the IC value from H f r to the cloud top. Some limitations are taken into account: the cloud depth must be at least 5.5 km [34], and it is mandatory to have H f reeze < H cloud_top < H domain_top to consider a possible emission. This simple approach could be improved later as, for example, following the recommendations in [38] for the vertical profile shape. Finally, and following schemes as implemented in CMAQ, for example, in [7], flashes are considered only in model grid cell where a convective precipitation was diagnosed during the current hour. This mix between CTH method and convective precipitation was already used by the authors of [9,39], for example.

The Studied Region and Period
To evaluate the impact of these emissions, a simulation ranging from 15 June to 31 August 2013 is performed. The first fifteen days are considered as a spin-up period and the results are only analyzed for the period from 1 July to 31 August 2013. The modeled domain has a horizontal resolution of 60 × 60 km and covers Europe and the northern part of Africa. Domain and period correspond to studies already performed with the WRF-CHIMERE model, where comparisons with measurements were performed [40,41]. With this large horizontal domain, it is expected to catch numerous convective clouds with different sizes and lifetimes, depending if they are in Europe or along the coast of Guinea during the monsoon.
Two simulations are performed with the same domain and period, and using the online model capability for the calculation of direct and indirect effects of aerosols on radiation and microphysics. The first simulation has all emissions except NO x by lightning and is called noLiNOx. The second one is identical to the first one but with NO x lightning emissions enabled and is called LiNOx. By comparing these two simulations, it is possible to discuss the following questions. (i) What is the impact of lightning NO x on pollutants concentrations? (ii) Is there a long range transport of modified pollutants concentrations? (iii) Is there a substantial impact on meteorological variables? Figure 1 presents the time-averaged number of flashes during the period from 1 July to 31 August 2013. For the model, the results are presented in number of flashes per minute and per km 2 . The simulation is compared with observational data from satellite detecting the number of flash events. The data are available over all the tropical and subtropical band between 38 • S and 38 • N at 3 to 6 km horizontal resolution. They are also averaged over the two months of July and August 2013. The TRMM Lightning Imaging Sensor (LIS) dataset was collected by the LIS instrument on the Tropical Rainfall Measuring Mission (TRMM) satellite used to detect the distribution and variability of total lightning occurring in the Earth tropical and subtropical regions. These data can be used for severe storm detection and analysis, as well as for lightning-atmosphere interaction studies. The LIS instrument makes measurements during both day and night with high detection efficiency [42][43][44].
Satellite WRF-CHIMERE model Model outputs and satellite data are not directly comparable as the model displays all calculated flashes while the satellite displays only events observed by the satellite when passing over a specific region. Having no information to extract the model data corresponding to the satellite leg, the comparison may be done only to see that modeled flashes are well distributed. If the distribution is correct for the main emitting area (Gulf of Guinea coast and the Sahelian area, some flashes are not properly reproduced by the model (northern Africa and Atlantic ocean). There is several reasons for this discrepancy: (i) as previously said, the data and model does not represent exactly the same thing, as the satellite has partial spatial measurements. To have a more significant comparison, it would be needed to make the comparison during several years to accumulate satellite data. However, in our case, we are making a two-month duration study with hourly data, and we are thus not in the same time scale. (ii) the diagnostic of flashes is done based on the deep convection scheme, itself uncertain. Differences may thus be due to errors in the deep convection diagnostic.

Impact on Ozone Concentrations
Results are presented in this section as maps of impact on ozone. This impact is estimated mainly by presenting the results as differences between the two simulations as (LiNOx-noLiNOx).
The first expected impact on pollutant concentrations in the troposphere is linked to the ozone-NO x chemistry. By adding NO x directly in the troposphere where no other sources of emissions exist (except aircraft, which are closer to tropopause), this additional source should be clearly visible. Maps are presented in Figure 2 Differences (LiNOx-noLiNOx) are also presented. At both altitudes (10 and 4000 m a.g.l.) the additional NO x emissions from lightning cause a net increase in ozone concentrations. The maximum impact is located where more convective clouds were diagnosed, and therefore where more NO x emissions from lightning were parameterized: a latitudinal band ranging between 0 and 15 • N. The range of the differences is 10 ppb. All differences are positive and have the same order of magnitude as in other studies such as in [10] with CMAQ, for example. The simulated impact is stronger and more spread in the free troposphere than close to the ground, due to stronger wind and more efficient advection aloft. The Wilcoxon signed-rank test, explained in the Appendix B, is applied to these maps of differences. The "+' symbols is added on the map where the differences are found to be statistically significant. Except for some areas in northern Africa and above the Atlnatic sea, the increase of ozone due to lightning NO x is significant.  For 2 m temperature (K) and 10 m wind speed (m s −1 ), the differences are alternatively negative and positive over the whole domain. No clear spatial relation between wind speed differences and temperature differences appear. Wind and temperature differences exist over the whole domain, not only where the lightning flashes were modeled. Warming up to 0.5 K is simulated over Africa and Saudi Arabia, and cooling down to −0.5 K appears over Africa and Europe. For 10 m wind speed, differences range from −0.5 to +0.5 m s −1 with maximum over seas. As for 2 m temperature, there is no spatial regularity and differences show patchy structures. The Wilcoxon signed-rank test is also applied for the meteorological parameters. For temperature, the two largest spots of differences (around latitude +20 • N) are found to be significant. Moreover, in Europe, many non-negligible impacts are found to be significant. On the other hand, the increase of temperature at latitude +30 • N is not significant and is certainly the results of a nonlinear perturbation on the meteorological model, a direct impact of the online c coupling of these simulations. For the wind speed, the largest changes due to the emissions are found to be all significant.

Impact on Other Variables at the Surface
The impact on AOD is rather weak: differences never exceed ± 0.05. The most important differences are inland and in the center of Africa, far from the lightnings and probably due to stronger wind over this area (Figure 3, upper-right panel) leading to more emissions of mineral dust over these desert areas. Only a few areas are found to be significant with the Wilcoxon test. Precipitation rates increase around 10 • N, with +0.05 kg m −2 h −1 . Over western Europe, more precipitation is simulated when lightning emissions are included, but on the contrary precipitation is reduced over eastern Europe. All in all, precipitation sensitivity over Europe is lower (≈±0.01 kg m −2 h −1 ) than over Africa. Spatially, the precipitation rate is much more impacted in the Sahelian region that the cloud fraction: it means that the surface coverage of clouds is not really modified and therefore the droplet size is probably impacted with more associated precipitation. Note that the cloud fraction and cloud mixing ratio were also examined (but not shown here). The cloud fraction was unchanged and the cloud mixing ratio shows the same differences than the precipitation rate, mainly linked to the change in wind speed, then latitudinal position of the clouds. As for AOD, only a few areas are found to be significant with the Wilcoxon test: but the area of the precipitation due to the monsoon is found to be significant for many locations. Figure 4 presents the impact of NO x lightning emissions on pollutant concentrations. For these results, time-averaged maps of concentrations differences are presented at the first model level (representing the near-surface). The Wilcoxon test is applied to the maps of differences as for ozone and the meteorological parameters. In this case, all changes are found to be statistically significant.
The impact on mineral dust is contrasted with negative and positive differences, but positive differences dominate. In Saudi Arabia, the differences are only positive and over a large area. In Africa, the differences are spatially correlated to the wind speed differences, showing the direct impact on the emission process. Values are important with changes about ± 30 µg m −3 . Sea salt concentrations are increased along the north-western African coast. Other hot spots of increase are also modeled in the Mediterranean. Compared to dust, the impact on concentrations is less important with extrema of differences around ± 1 µg m −3 . Changes on sea salt are spatially correlated to changes on wind speed, which reflects the fact that wind speed is the man driver of sea salt emissions.  For NO x , the most important changes are located where more flashes were diagnosed (Figure 1), which reflects the fact that flashes are directly emitting NO x following Equation (7). The largest impact is in the Sahelian belt, for latitudes between 5 and 15 • N. Differences are also diagnosed in Western Europe, where fewer flashes were simulated. However, over Europe, the differences are negative and this is probably due to the transport of ozone from Africa to Europe then local reaction with NO.
For PM 2.5 , the spatial pattern is close to the NO x one. However, the most important increase is located close to the megacities of the Gulf of Guinea (Abidjan, Lomé, Lagos). In this area, there is no change of temperature and wind speed, but the model diagnoses less precipitation. This change may explain this PM increase. In average the changes in PM are not very important with maximum around ± 2 µg m −3 .
Results are also presented for nitric acid HNO 3 and nitrate NO − 3 . Impacts are low for the two species, with maximum changes around ± 0.2 µg m −3 . For HNO 3 , the most important differences are calculated over the whole Africa and are partially co-located with the increase in ozone and NO x . For NO − 3 , changes are about ±0.5 µg m −3 with opposite variability than sea salt.

Vertical Cross Section of Impact on Ozone Concentrations
To better understand the impact of NO x emissions on pollutants and meteorology, latitude-altitude plots of wind speed, mineral dust, NO x , O 3 , nitric acid HNO 3 , and nitrate NO − 3 concentrations are presented in Figure 5. As for maps, modeled values are time-averaged during the two months of simulation. The values are also spatially averaged here for longitudes between 10 • W and 10 • E.
For wind speed, one can note two different vertical structures of net increase, located around latitudes +10 and +30 • N. Positive or negative, the impact is not very important with maximal changes around ±0.5 m s −1 . The impact on mineral dust concentrations may also be negative or positive, with ± 30 µg m −3 .
For NO x , and as already analyzed with the surface maps, the net change is an increase, directly due to NO x production by the lightning. Values are 0.05 ppb, maximal where most flashes are located, between 5 and 20 • N. A very substantial increase in ozone concentrations, maximal between 5 and 20 • N. This increase extends down to the surface and up to the top of model, with maximal values 10 ppb between 3000 and 7000 m a.s.l. A similar structure is simulated for HNO 3 , with an increase of 0.25 ppb around 4000 m AGL and for latitudes between +5 and +20 • N. For NO − 3 , the impact is mainly an increase of concentrations, stronger in the lower troposphere, but with very low values.

Comparison of Time Series
The previous sections showed impact of lightning NO x emissions with time-averaged values providing a synthetic view of the results. As time-averaging hides variability, we compare in this section time series of the two simulations for selected locations where surface measurements are available. Coordinates are presented in Table 1. Units for ozone and NO 2 is now in¯g m −3 to be directly compared with surface measurements. Measurement data and model results are daily averaged.  Figure 6 presents comparisons of ozone and NO 2 surface concentrations for the sites of Campisabolos and Barcarrola (Spain). For ozone, surface concentrations vary between 60 and 120 ppb during the two months and there is no daily averaged major peak or event. The time variability is close between the two sites, and one can distinguish the two local minima around 1 and 10 August. The comparison between the two simulations NoLiNOx and LiNOx shows that the differences are low (as already shown on the maps) and all along the period, even if they are larger in August than in July. The LiNOx simulation shows larger ozone concentration values but lower NO 2 due to increased formation of ozone in Africa followed by ozone transport towards Europe. The reduced NO 2 concentrations in Europe are also due to increased oxidation of NO 2 into nitrates nitrate NO − 3 .

AOD and PM 10
For AOD, Figure 7, the two AERONET stations of Banizoumbou and Dakar are used. Results show that the model slightly underestimates AOD, but, globally, is able to reproduce the observed background values and the time variability. As for surface concentrations of the gaseous species, there is no systematic difference between the two time-series. The changes on AOD are never larger than ≈ ± 0.1. For PM 10 surface concentrations also, small differences exist between both simulations but no systematic pattern or bias is visible. With LiNOx emissions, the most important effect is modeled at the end of the period, in August, with differences up to ± 2 µg m −3 . For the selected locations, the effect of NO x emissions by lightning is small compared to the model negative bias.  Figure). Crosses indicate the model grid cell where the result is found significant, using the Wilcoxon signed-rank test.

Statistical Scores
All previous results are summarized as statistical scores in Table 2. These scores are defined in Appendix A. For AOD, we used data from 32 AERONET stations located in Africa and Western Europe. For pollutants surface concentrations, we used data from 35 stations located in Western Europe. The complete list of the stations used is presented in [23]. R s is the spatial correlation, R t is the temporal correlation, RMSE is the Root Mean Square Error, and the bias is the difference (LiNOx-NoLiNOx). Details on these calculations are presented in Appendix A. Results are presented for ozone, NO 2 and PM 10 surface concentrations, and AOD. Results show good spatial correlation for the four variables. For the temporal correlations, based on daily mean values, the correlation is low for NO 2 (0.28 and 0.31), PM 10 (0.27 and 0.26) and AOD (0.30 and 0.31), but with a better correlation for ozone (0.64). The difference between the two simulations is not important and these scores show there is no significant improvement of surface concentrations correlation by adding NO x lightning emissions. A slight improvement is shown for ozone and spatial correlation, meaning that the emissions are well located and the transport done in a right way. For the RMSE and bias, there is also no significant improvement with these additional emissions.

Conclusions
This study has presented the implementation of NO x emissions from lightning in the CHIMERE regional chemistry transport model. The parameterization, based on the work in [34], is classical and used in many models. To evaluate this parameterization, a large simulation domain covering Africa and Europe is designed. Two simulations are done on July and August 2013: one with lightning emissions, the other one without. This corresponds to a period where the monsoon is active in the Sahelian zone and large convective clouds create NO x emissions from lightning. Emissions are also present in Europe due to summertime thunderstorms.
The impact of these emissions is analyzed by calculating the difference between the two simulations. The WRF and CHIMERE models are online coupled, including radiative and microphysical effects of the aerosols. The addition of NO x emissions in the troposphere can therefore modify gas and aerosol concentrations but can also impact meteorology.
The results show that ozone has differences of ± 10 ppb throughout the troposphere and averaged over the two months. These differences are mainly in the tropical band, where emissions are strongest (latitude 5-15 • N). Concentrations of NO x and nitrates also increase following the same pattern, but the simulated increase is small (0.05 ppb for NO x , 0.3 ppb for HNO 3 ). Maps of the other variables show that there is an impact on all gas and aerosol species and on meteorology but with changing signs and less well-defined patterns. 2 m temperature is modified by about ± 0.5 K, 10 m wind speed by about ± 0.5 m s −1 , precipitation rate by 0.05 kg m −2 h −1 (especially where is the monsoon), and AOD by 0.05 (on the continent and the ocean). Changes in temperature and wind induce a feedback on emissions with changes in dust and sea salt emissions. This corresponds to a change of 50 µg m −3 for dust and 1 µg m −3 for sea salt.
Surface concentrations of NO x are changed but mostly where emissions occur. The changes in PM 2.5 correspond mainly to changes in sea-salt and the NO x lightning emissions area close to the Gulf of Guinea. They are about ±1¯g m −3 As a long-lived species, the additional ozone produced in the tropics can be transported towards Northern Africa and Europe, where a positive impact on ozone concentrations and negative impact on NO 2 concentrations is simulated. From a statistical point of view, even if we see changes in meteorology and concentration results, this does not particularly improve the scores when comparing simulations to surface measurements. The impact of NO x lightning on surface pollution and on meteorology is not very pronounced but is not negligible. The present study confirms that NO x emissions for lightning bring large changes in the tropospheric ozone cycle in the tropics, small changes for NO x and HNO 3 concentrations, and additional variability for other variables. It should therefore be included if possible in simulations of atmospheric composition, particularly in the tropics and/or in zones with active atmospheric convection.
As perspective, some improvements could be done for this parameterization. An important point is that the model resolution is too large to well take into account the fast NO x /O 3 chemistry. This could induce an underestimation of ozone production by ≈3 ppb in the emission region, as found in [45]. A possible model improvement could be to add to our scheme a plume-in-grid approach such as the one proposed by these authors and successfully used in the Meso-NH and GeosChem models.

Appendix A. Definition of Statistical Scores
Three statistical indicators are used: the spatial Pearsons' correlation, the temporal Pearsons' correlation, and the normalized RMSE. The temporal correlation, R t , is computed for each station and is directly related to the hourly variability. O t,i and M t,i represent the observed and modeled values, respectively, at time t and for the station i, for a total of T days, and a total of I stations. The mean time averaged value X i is The temporal correlation R t,i for each station i is calculated as The mean temporal correlation, R t , used in this study is thus The spatial correlation, noted R s , uses the same formula type except it is calculated from the temporal mean averaged values of observations and model for each location where observations are available.
The spatio-temporal mean averaged value is calculated as and the spatial correlation, R s : The normalized Root Mean Square Error (nRMSE) is expressed as for all stations i and all times t.

Appendix B. The Mann-Whitney-Wilcoxon Test
To quantify the statistical significance of differences, the Mann-Whitney test is applied (also called the Wilcoxon signed-rank test) [46]. This test is nonparametric: there is less restrictive assumptions as, for example, the fact that the distribution has not to be normal. This test examines the two sets of data (in our case, two different simulations) by combining all data and by sorting them in ascending order.
The first and second simulations, noted x 1,i and x 2,i , have here the same dimension, N. We first calculate the difference between the two datasets: In the original scheme, a reduced dataset with dimension N r is built by removing data where d i = 0. In our case, to avoid noise, an value is defined as =0.01 ×max(d i ). Values are then removed from the dataset if − < d i < . The remaining data d i are sorted, in ascending order and by their rank, R i , and stored. The statistic test W is calculated as A z score is the calculated as z = W σ w (A9) with, when N r ≥ 20, σ w = N r (N r + 1)(2N r + 1) 6 (A10) The null hypothesis H 0 is rejected (i.e the differences are significant and not due to hazard) if |z| > z critical . For a level of significance of 0.05, we have the value z critical ≈ 1.645. In the results section, Figures present crosses at the points where the difference was found to be significant.