Evaluation of a new lightning-produced NOx parameterization for cloud resolving models and its associated uncertainties

A new parameterization of the lightning-produced NO x has been developed for cloud-resolving models. This parameterization is based on three unique characteristics. First, the cells that can produce lightning are identiﬁed using a vertical velocity threshold. Second, the ﬂash rate in each cell is estimated from the non-precipitation and pre- 5 cipitation ice mass ﬂux product. Third, the source location is ﬁlamentary instead of volumetric as in previous parameterizations. This parameterization has been tested on the 10 July 1996 Stratospheric-Tropospheric Experiment: Radiation, Aerosols and Ozone (STERAO) storm. Comparisons of the simulated ﬂash rate and NO mixing ratio with observations at di ﬀ er- 10 ent locations and stages of the storm show a good agreement. An individual ﬂash produces on average 121 ± 41 moles of NO (7.3 ± 2.5 × 10 25 molecules NO) for the simulated high cloud base, high shear storm. Sensitivity tests have been performed to study the impact of the ﬂash rate, the cloud-to-ground ﬂash ratio, the ﬂash length, the spatial distribution of the NO molecules, and the production rate per ﬂash on the NO 15 concentration and distribution. Results show a strong impact from the ﬂash rate, the spatial placement of the lightning-NO x source and the number of moles produced per ﬂash.

cloud flashes. The way the LNO x is horizontally distributed also differs from one parameterization to another. The simplest way consists of distributing the NO molecules horizontally in the whole cloud or within the 20 dBZ contour (Pickering et al., 1998;De-Caria et al., 2000. Fehr et al. (2004) chose randomly one column in the cloud with mixing ratio higher than a threshold where NO would be deposited. The most 15 sophisticated horizontal distribution is the one of Ott et al. (2007) who attempted to mimic the tortuous aspect of the lightning channel. For each vertical level, the number of points where the NO is deposited follows the bimodal distribution, and the points where the NO molecules are distributed is chosen randomly among possible points located in a predetermined area downwind of the convective core. With this approach, 20 the NO is no longer instantly diluted in a large volume of the cloud, but it is distributed in a region where lightning flashes are expected to propagate.
The last step of a LNO x parameterization is prescribing the magnitude of NO produced by lightning flashes. Pickering et al. (1998) adopted the Price et al. (1997) values (6.7×10 26 molecules NO per CG=1113 moles NO per CG and 6.7×10 25 molecules per 25 IC=111 moles NO per CG) assuming that a cloud-to-ground flash produces ten times more NO molecules than an intra-cloud flash. Wang and Prinn (2000) tested both the Price et al. (1997) values and the Franzblau and Popp (1989) DeCaria et al. (2005) used the method of Price et al. (1997) to deduce the production rate per CG in the 12 July 1996 STERAO storm. To deduce the production rate per IC, they performed several simulations with different ratios between the production rate per IC and the production rate per CG. They concluded that taking a production rate of NO per IC that is 75 to 100% of the production rate 5 of NO per CG best matches with observations. Using the same method, Ott et al. (2007) deduced that a production rate of 360 moles of NO per IC and CG compares favorably with observations for the 21 July 1998 European Lightning Nitrogen Oxides Project (EULINOX) storm. For this same storm, Fehr et al. (2004) concluded that the flash production rates suggested by Price et al. (1997) are not supported by the anal- 10 ysis. They estimated the CG and IC production rates to be 2.1×10 26 molecules NO (349 moles NO) and 2.9×10 26 molecules NO (482 moles NO), respectively. The recent results of DeCaria et al. (2000DeCaria et al. ( , 2005, Fehr et al. (2004) and Ott et al. (2007) are in agreement with Ridley et al. (2005)  Interactive Discussion duction rates are the ones using explicit electrical schemes, which placed the lightningproduced NO source in a small volume along the simulated lightning flash and reproduced the observed peaks and flux quite well for this storm. Since the other models distribute the NO in a large volume, these results tend to suggest that the position of the flash and its spatial distribution (vertical and horizontal) is fundamental in modeling 5 LNO x at the cloud scale. Despite some studies that investigated the impact of the relative contribution of intracloud and cloud-to-ground discharges to the LNO x production (DeCaria et al., 2000(DeCaria et al., , 2005Ott et al., 2007), few have been made to understand the relative impact of different choices for the different steps of the LNO x parameterizations. The objectives of 10 this paper are to introduce a new LNO x parameterization and to evaluate the sensitivity of the LNO x parameterization to lightning parameters (flash rate, flash length, spatial distribution of the NO molecules, presence of short duration flashes) through a series of sensitivity studies. A new parameterization for the LNO x production is presented in Sect. 2. The model framework and configuration and the studied storm are described 15 in Sect. 3. In Sect. 4, the reference simulation is studied, and the sensitivity tests are presented in Sect. 5.

Description of the lightning-produced NO x parameterization
The goal of this new parameterization is to predict the temporal and spatial distribution of individual lightning flashes without using a computationally explicit electrical scheme.
Except for the explicit electrical schemes (Zhang et al., 2003;Barthe et al., 2007b), LNO x parameterizations distribute the NO molecules either where the radar reflectivity exceeds 20 dBZ or within cloudy regions where the temperature is colder than −15 • C (Pickering et al., 1998;DeCaria et al., 2000). A first approach has been made by Ott et al. (2007) to simulate the filamentary aspect of a lightning flash. The parameteriza-25 tion presented here builds on this approach by estimating the lightning activity in each individual convective cell rather than the entire storm. This LNO x parameterization is Introduction When several convective cells are present, they may not all produce lightning flashes despite having radar reflectivity greater than 20 dBZ or cloud top colder than −15 • C. An algorithm has been developed to detect potentially electrified cells. In order for 5 lightning flashes to occur, the updraft speed must exceed 15 m s −1 which is similar to the 10-12 m s −1 threshold estimated by Zipser and Lutz (1994) for three different environmental regimes in Australia.
Unique to this parameterization is the prediction of lightning flash rate based on the fluxes of non-precipitating and precipitating ice. Through theoretical and obser-10 vational investigations, Blyth et al. (2001), Deierling (2006) and Latham et al. (2007) have shown a strong correlation between the total flash rate and the precipitation and non-precipitation ice mass flux product which is called the flux hypothesis. Barthe et al. (2007a) (hereinafter referred to as BDB07) show that simulated ice mass flux product for the entire storm is quite similar to the ice mass flux product derived from radar 15 observations. In this study, the parameterization is improved by calculating the nonprecipitation and precipitation ice mass flux product for each individual convective cell and is associated to a total flash rate per cell. Thus, the initial equation of BDB07 has been slightly modified to take into account the computation of the total flash rate per individual cell.
where F MF is the total flash rate (fl. min −1 ) computed from the precipitation and non-precipitation ice mass fluxes, f p (kg m s −1 ) and f np (kg s −1 ), respectively. The Interactive Discussion be triggered in the region downwind of the maximum vertical velocity (Proctor, 1981;Christian et al., 1999;Ushio et al., 2003;Ott et al., 2007). This region is defined by the convective core and by the region extending 10 km downwind of the maximum vertical velocity (C trig ) (Fig. 1). The downwind direction is assumed to be the mean wind direction at the altitude where both non-precipitation and precipitation ice particles 5 are encountered. The center of the region where an individual flash can propagate is chosen randomly among all the points of the cylinder that are in the glaciated part of the cloud. Another cylinder C prop where a lightning flash can propagate (radius 4 km but could depend on the flash length) is centered on the randomly chosen point. The points of the cylinder where the discharge can propagate are restricted to the region of 10 the cloud where ice particles can be found since the hydrometeors that carry most of the electric charges are the ice particles (Barthe and Pinty, 2007a). The vertical distribution of the flash channel follows a bimodal distribution (DeCaria et al., 2000(DeCaria et al., , 2005Ott et al., 2007). This kind of structure is the most commonly observed (Shao and Krehbiel, 1996;Krehbiel et al., 2000;Rison et al., 1999;Thomas 15 et al., 2001;Wiens et al., 2005;Bruning et al., 2007) and simulated by explicit electrical schemes (Mansell et al., 2002;Barthe and Pinty, 2007b). For each altitude level, the grid points reached by the lightning channel are chosen randomly among the possible points in the cylinder C prop to mimic the filamentary and tortuous aspect of a lightning flash (Ott et al., 2007). The flash length of the lightning flash is prescribed either to 20 be constant or to have a lognormal distribution Pinty and Barthe, 2008).
The amount of NO produced per flash is assumed to depend on the flash length and on the altitude : with n NO the number of NO molecules produced per flash length (molecules m −1 ), P the pressure (Pa). Wang et al. (1998) (Höller et al., 1999;Stith et al., 1999;Huntrieser et al., 2002;Skamarock et al., 2003;Ott et al., 2007). The aim of this new parameterization is to reproduce the global morphology of a lightning flash in terms of spatial distribution and length in order to avoid the instan-5 taneous dilution of the NO in the storm. This is important for the redistribution of the chemical species and for the comparison between model results and observations.

Experimental and model design
The lightning-produced NO x parameterization described above has been placed in the Weather Research and Forecasting (WRF) model. Simulations of the 10 July 1996 10 STERAO storm have been conducted to evaluate and assess its sensitivities.

The WRF model
The WRF model solves the conservative (flux-form), non-hydrostatic compressible equations using a split-explicit time-integration method based on a 3rd order Runge-Kutta scheme (Skamarock et al., 2005;Wicker and Skamarock, 2002). Scalar transport 15 is integrated with the Runge-Kutta scheme using 5th order (horizontal) and 3rd order (vertical) upwind-biased advection operators. Transported scalars include water vapor, the different hydrometeor categories and the chemical species. The cloud microphysics is described by the single moment (bulk water) approach (Lin et al., 1983). Mass mixing ratios of cloud water, rain, ice, snow, and graupel/hail are predicted. Cloud water 20 and ice are monodispersed and rain, snow, and hail have prescribed inverse exponential size distributions. For the graupel/hail category, the intercept parameter of the exponential distribution is 4×10 4 m −4 and the density is 917 kg m −3 which corresponds to characteristics of hail particles.
The model predicts the mixing ratios of methane (CH 4 ), carbon monoxide (CO), Introduction , hydroxyl radical (OH), hydroperoxy radical (HO 2 ), methyl hydroperoxy radical, nitrogen dioxide (NO 2 ), nitric oxide (NO), nitric acid, hydrogen peroxide (H 2 O 2 ), methyl hydrogen peroxide, formaldehyde (CH 2 O), formic acid, sulfur dioxide, ammonia, and aerosol sulfate . The gas-phase chemistry represents O 3 -NO x -CH 4 chemistry, but not non-methane hydrocarbon chemistry which likely has 5 a negligible effect on NO x mixing ratios in the anvil. Species are partitioned between the gas and aqueous phases via Henry's law equilibrium for low solubility species (e.g. CO) or via diffusion-limited mass transfer for high solubility species (Barth et al., 2001). The aqueous chemistry occurring in the cloud water and rain represents OH and HO 2 oxidation of dissolved species plus the oxidation of S(IV) by O 3 and H 2 O 2 . The pH of the 10 drops is iteratively calculated via a charge balance assuming CO 2 is 360 µmol mol −1 . The chemical mechanism is solved with an Euler backward iterative approximation using a Gauss-Seidel method with variable iterations. A convergence criterion of 0.01% is used for all the species. The dissolved species are transferred to other cloud hydrometeors via microphysics processes (Barth et al., 2001). When cloud or rain drops 15 freeze, it is assumed that all of the dissolved species is retained in the frozen hydrometeor. In addition, adsorption of gas-phase nitric acid, sulfur dioxide, formaldehyde, and hydrogen peroxide is represented using the Langmuir equilibrium model approach (Tabazadeh et al., 1999;Popp et al., 2004). In addition to the chemically active species, a tracer of NO x from lightning (LNO x ) 20 is included in WRF. The LNO x tracer corresponds to the NO mixing ratio produced by lightning flashes. LNO x is transported only and does not undergo any chemical reactions.

The 10 July 1996 STERAO storm
The LNO x production in the 10 July 1996 STERAO storm has been widely studied 25 (Stith et al., 1999;Skamarock et al., 2003;Barthe et al., 2007b;Barth et al., 2007a,b (Barth et al., 2007b). Thus, it is interesting to investigate the origin of such discrepancies and to evaluate the uncertainties associated with each step of the LNO x parameterization. The 10 July 1996 storm has also been chosen because there is a unique set of data for this storm: storm structure and kinematics from radar data, lightning 5 flash characteristics and in-situ chemical species measurements from two aircraft . During the STERAO-A experiment the ONERA VHF interferometeric mapper (ITF) measured the total (IC+CG) lightning activity while the National Lightning Detection Network (NLDN) documented the CG activity . NLDN provided lo-10 cations of the ground connections from the measurements of the electric and magnetic field due to the high current of return strokes (Cummins et al., 1998). ITF was designed to detect and locate VHF radiation emitted during both IC and CG flashes (Defer et al., 2001). VHF radiation is recorded during stepped leaders, dart leaders, recoil streamers and return strokes of negative CGs (Defer et al., 2001). Comparison between optical radiation detected by NASA Optical Transcient Detector (OTD) and VHF radiaton recorded by ITF for one passage over the STERAO-A domain (9 July 1996) showed consistent observations between OTD and ITF (Defer et al., 2006), suggesting that a flash sensed by OTD corresponds to a flash sensed by ITF. However, Boccippio et al. (2002) estimated the flash detection efficiencies to be 93% (nighttime) 20 and 73% (local noon) for LIS and 56% (nighttime) and 44% (local noon) for OTD. Thus, extrapolations of the results reported here must be done cautiously as the detection of flashes by satellite instruments is less than that by the ITF system. Further, the 10 July 1996 STERAO storm is likely not a typical thunderstorm whose flash, physical structure, and dynamics characteristics are necessarily meaningful for extrapolation of 25 LNO x production to the global scale.

Initialization
The simulation performed is similar to those described by Skamarock et al. (2000Skamarock et al. ( , 2003 and Barth et al. (2001. The environment was assumed to be homogeneous, thus a single profile was used for initialization. The initial profiles of the mete-5 orological data were obtained from sonde and aircraft data (Skamarock et al., 2000).
The convection was initiated with three warm (3 • C perturbation) bubbles oriented in a NW to SE line. WRF is configured to a 160×160×20 km 3 domain with 161 grid points in each horizontal direction (1 km resolution) and 51 grid points in the vertical direction with a variable resolution beginning at 50 m at the surface and stretching to 1200 m at 10 the top of the domain. The simulation was integrated at a 10 s time step. To keep the convection near the center of the model domain, the grid is moved at 1.5 m s −1 eastward and 5.5 m s −1 southward. The simulation was integrated for a 3-h period. While the observed storm lasted from 21:30 to 03:00, only the multicell and supercell stages are simulated. These stages correspond to 23:15-02:15 UTC in the observations.

15
Initial mixing ratios of the chemical species are the same as those in . Of interest to this study are NO, NO 2 and O 3 mixing ratios. NO and NO x initial mixing ratios are relatively high near the surface, low in the mid-troposphere, and moderately high in the UTLS region ( Fig. 2). O 3 mixing ratios are 60 ppbv near the surface and are fairly constant with height to 12 km m.s.l. where mixing ratios increase 20 into the stratosphere. Because of the short integration time and small domain, there is very little effect of lightning-produced NO on O 3 mixing ratios in these simulations. The initiation process is the same in all the simulations.
The details of the control experiment (REF) are summarized in Table 1. The total flash rate is computed from the non-precipitation and precipitation ice mass flux prod- Interactive Discussion bimodal distribution correspond to the −15 • C and the −50 • C isotherms, respectively. The flash length of an individual storm is assumed to be lognormal in the range 1 to 400 km. As in the observations (Defer et al., 2001, the percentage of the short flashes (<1 km) is set to 47%. Among these short flashes, 36% are considered as short duration flashes (<1 ms). In this simulation, it is assumed that the short duration

10
The dynamics and the microphysics structure of this storm has been studied previously (Skamarock et al., 2000, BDB07). These studies showed that the WRF simulations compare well to the radar data. Based on the storm structure and dynamics, three different stages have been identified in this storm: a multicell ( Table 3). The mean flash rate during the multicell stage is 9.2 fl. min −1 for the simulated flash rate, which is similar to the 10.5 fl. min −1 observed. During the transition stage, the ice mass Interactive Discussion flux parameterization underestimates the mean flash rate (15.7 vs. 20.4 fl. min −1 ) and overestimates it in the supercellular stage (43.8 vs. 33.8 fl. min −1 ). BDB07 have shown that the supercellular stage begins earlier in the simulation compared to observations and the simulated mass flux is larger than observed in this stage of the storm. The differences in the simulated and observed microphysics and dynamics features cause 5 the differences in the observed and simulated flash rate. Figure 4 shows horizontal and vertical cross sections of the NO mixing ratio during the different stages of the storm. During the multicell stage, NO production by lightning occurs mainly in the convective cores. At 3600 s, five different cells are present in the domain: three of them are located along a NW-SE axis as initialized and the two others  shows that high values of the NO mixing ratio are mainly located in the range 6-7.5 km altitude and 11-14 km altitude due to the bimodal distribution of the flash segments.
At 5400 s, the flash rate is lower than in the multicell stage (see Fig. 3) which can explain why peak NO values decrease (Fig. 4c). NO starts to spread over a large region in the anvil. The vertical cross section during the transition stage ( Fig. 4d) confirms that 20 less LNO x is produced at this time of the simulation.
In the supercell stage (t=7200 s), the flash rate is higher than in the two other stages with values in the range 19-61 fl. min −1 leading to high values of the NO mixing ratio. In the anvil, values higher than 4000 pmol mol −1 extend horizontally for 25 km from the convective core. NO mixing ratios up to 1000 pmol mol −1 can be found 100 km 25 downwind of the updraft maximum. Even though the LNO x parameterization produces LNO x only in a small region in and downwind of the convective core, NO molecules are transported and diluted in the whole cloud. Lightning flashes are then responsible for a large amount of NO in the whole system during the supercell stage (Fig. 4f). The NO mixing ratios from the model are next compared to the University of North Dakota's (UND) Citation aircraft measurements. The observed NO vertical crosssection (Fig. 5) is a result of projecting the NO aircraft measurements collected between 23:16 and 00:36 UTC onto the across-anvil plane, which is ∼60 km downwind of the convective core. Several regions of NO mixing ratio higher than 540 pmol mol −1 5 can be seen at the altitude of 11.5 km and up to 13 km which is in agreement with observations (Fig. 5).
NO transects across the anvil during the multicell and the transition stage of the storm have been plotted in Fig. 6. After 1 h of simulation, the transect is 10 km downwind of the southeastern cell. The simulated transect compares well with observations. 10 A peak of 1800 pmol mol −1 is simulated and NO>500 pmol mol −1 extends over a distance of 20 km in the simulation. The distance over which the observed NO mixing ratio is higher than 500 pmol mol −1 is 30 km. Thus, WRF with its new LNO x parameterization is able to simulate peaks of NO higher than 1000 pmol mol −1 in the region where lightning flashes are mostly triggered and propagate. After 1 h 30 min of simulation, 15 the transect is located 50 km downwind of the main convective core. The trends of the simulated and observed transects are the same, however the simulated values are lower in magnitude compared to the observations. The lightning activity in this storm started at 21:52 UTC, i.e. 108 min before the first flash was triggered in the simulation (25 min after the beginning of the simulation or 23:40 UTC). NO may have accumulated 20 in the environment of the storm leading to larger values than simulated. In summary, the WRF model coupled with the new LNO x scheme gives results in good agreement with observations both near the convective cores and in the anvil.
The effect of the transport and lightning production of NO can be seen by comparing NO and NO x final mixing ratios with their initial values. Figure  The NO x mixing ratio at the ground is reduced from 600 pmol mol −1 (Fig. 2a) to 350 pmol mol −1 (Fig. 2b). When comparing the NO and LNO x curves, it can be seen that the two peaks at 13.0 km m.s.l. and 6.5 km m.s.l. are caused by NO production by lightning flashes. The intense vertical motions and the high electrical activity in the 10 July 1996 STERAO storm causes the modeled NO mixing ratio to increase by 120% at 5 13 km m.s.l. and by 400% at 6.5 km m.s.l.
The simulated NO x flux through the anvil is 5.72×10 −8 moles m −2 s −1 , which is similar to the value derived from observations reported by Skamarock et al. (2003) Table 2). In this simulation, a typical flash is 21.7 km long. The mean production rate per flash calculated by dividing the total production of NO by the total number of flashes is 121.3 moles of NO (7.3×10 25 molecules NO). This production rate is in the lower range of the estimates from the last 5 years as reported by Schumann and Huntrieser (2007). Furthermore, our estimate is lower than most of the 3-D CRM studies.

15
To examine the key processes in the LNO x parameterization for CRMs, the sensitivity of the NO mixing ratio to several parameters is investigated. In particular, the impact of the total flash rate, the CG rate, the flash length, the short duration flashes, the spatial distribution of the NO molecules, and the NO production rate per flash are studied. A summary of the name and conditions of the experiments is given in Table 1. 20

Total flash rate
When simulating the production of LNO x in thunderstorms, the total flash rate is fundamental. Two different parameterizations for the total flash rate have been tested in this section and compared to observations from the ONERA interferometer (Defer et al., 2001). First, the total flash rate can be deduced from the non-precipitation and pre- cipitation ice mass flux product as described in Sect. 2 (Eq. 1) and used in the REF simulation. Secondly, a sensitivity simulation is performed using the flash data from the ONERA interferometer (FR OBS, Table 1). To determine the flash rate in each convective cell i , the cell flash rate (F (i ) ) is assumed to be in proportion to the cell ice mass flux product (f (i ) ). That is, where F is the total flash rate. Thirdly, the total flash rate can be determined from the maximum vertical velocity w max (m s −1 ) following Price and Rind (1992): Two different simulations use the Price and Rind (1992) approach for the flash rate. 10 In the first simulation, w max is the maximum vertical velocity in the whole domain: the total flash rate is then computed for the whole domain (FR PR92). The distribution of the flash rate in each cell is in proportion to the cell ice mass flux product (Eq. 3). Because the Price and Rind (1992) parameterization has been used several times in CRMs (Pickering et al., 1998;Fehr et al., 2004;Barth et al., 2007b) where several 15 cells can be identified, a second simulation in which the total flash rate per cell is computed from the maximum vertical velocity in each individual cell (FR PR92 CELL) is performed. That is, Eq. (4) is applied for each cell. Equation (4)  Interactive Discussion collision between more or less rimed ice particles (Takahashi, 1978;Jayaratne et al., 1983;Saunders et al., 1991, e.g.). Using an explicit electrical scheme in a CRM, Barthe and Pinty (2007b) showed that the first flash was triggered 20 min after ice particles appeared in the cloud. Since the maximum vertical velocity in their simulated storm did not exceed 20 m s −1 , a shorter delay can be expected in the 10 July storm which is 5 more intense. The simulated 10 July STERAO storm starts to produce ice particles at 16 min and the flux hypothesis allows the first flash to be triggered at 25 min, i.e. 9 min later. In the FR OBS simulation, the observed flash rate is available since 21:52 UTC, but it is only taken into account starting at 25 min (23:40 UTC) to allow the cloud to develop, to electrify and to trigger lightning. 10 As some differences exist between the time of the observed and simulated stages of this storm, the lightning activity is compared in each stage. Table 3 shows highly variable results for FR PR92 and FR PR92 CELL. An important point is that the flash rate in FR PR92 and FR PR92 CELL is in advance compared to REF and to observations (Fig. 3). For example, at the end of the transition stage the total flash rates in the 15 FR PR92 and FR PR92 CELL simulations increase 10 min before the REF simulation and the observed flash rate increase. In the Price and Rind (1992) parameterization for the total flash rate, only the maximum vertical velocity is considered. There is no information about the ice content or about conditions favorable for the non-inductive separation mechanism which can lead to this lag. Thus, contrary to the total flash rate 20 deduced from the flux hypothesis, the total flash rate from Price and Rind (1992) seems to be not as suitable for use in CRM as it does not take into account the microphysical development of the storm which is of primary importance for the cloud electrification. Using an explicit electrical scheme to simulate a STEPS supercellular storm, Kuhlman et al. (2005) concluded that there is no correlation between the maximum vertical ve-25 locity and the total flash rate. Deierling (2006) confirmed this conclusion using radar and lightning data.
To examine the effect of using different flash rate parameterizations, the NO mixing ratio downwind of the convective core is analyzed. Figure 5 (Fig. 3) since 1193 flashes are produced in FR PR92 and 1184 in REF, 5 but instead are attributed to the lower flash rate during the multicell stage (487 flashes in FR PR92 and 728 flashes in REF). Indeed, NO molecules are produced mostly near the convective core and it takes ∼40 min to be transported in the anvil to 50 km downwind. Because the number of flashes in the multicell stage is higher for FR OBS and FR PR92 CELL than for the REF simulation, the NO mixing ratios for FR OBS 10 and FR PR92 CELL are also higher than those for REF in the vertical cross section at 6000 s. Along the transects across the convective core during the multicell stage ( Fig. 6), all the simulations exhibit peaks higher than 1000 pmol mol −1 characteristic of local production of LNO x . The peak of NO reaches 2700 pmol mol −1 in FR PR92 CELL, which is the simulation with the most flashes in the first hour, and 1200 pmol mol

Cloud-to-ground flash rate
Boccippio et al. (2001) analyzed four years of OTD and NLDN flash rate data over the United States and found that the CG to IC ratio may be dominated by storm type, morphology and level of organization instead of environmental parameters. Here, we 5 examine the influence of the CG to IC ratio on the NO mixing ratios. In this study, the cloud-to-ground ratio α is defined by: Price and Rind (1993) estimated the cloud-to-ground ratio α P R from the depth Z of the layer from the freezing layer to the cloud top (simulation CG PR93). is computed taking the average altitude where the total hydrometeor mixing ratio decreases to 10 −5 kg kg −1 . Two other simulations (Table 1) have been performed with the observed α OBS ratio (CG OBS), and with a constant value (CG CST). The constant value is α CST =0.02 which corresponds to α OBS averaged during the simulated storm. These simulations are compared to REF in which CG flashes are not considered. 20 Figure 3 displays the observed and simulated N CG /(N IC +N CG ) ratio. The α P R ratio from Price and Rind (1993)  from an elevated charge region (MacGorman et al., 1989). MacGorman et al. (1989) hypothesized that strong updrafts can suspend the negative charge center to higher altitudes than in ordinary storms. This elevated charge would favor IC flashes over CG flashes. Here, we can investigate the impact of the high α ratio on the LNO x production compared to the other simulations. 10 Pickering et al. (1998) used Eq. (7) to estimate the CG ratio in seven different storms. In the two mid-latitude continental events in which CG flash data were available, the simulated CG rate was in reasonable agreement with observations. For the total flash rate in the simulated 21 July 1998 EULINOX storm, Fehr et al. (2004) rescaled Eq. (7) by a factor of 1.10. As for the total flash rate estimated by the Price and Rind (1992) parameterization, it is not clear if Eq. (7) needs to be rescaled depending on the storm or on the model. More tests should be done on several convective cases with the same model and with available observations. When comparing the NO mixing ratios along the transects and in the vertical cross section across the anvil (Figs. 5 and 6), the CG OBS results are in better agreement 20 with the REF simulation and with observations than CG PR93 and CG CST results. The CG PR93 and CG CST results have NO mixing ratios lower than the observations both in the convective region and in the anvil. Consequently, it is important to have the right number and temporal distribution for the CG flashes. Because the α P R ratio is ∼10 times greater than α OBS , more NO is produced in the mid troposphere at the expense of 25 the upper troposphere resulting in a 33% increase of NO in the mid-troposphere and a 20% decrease in the upper troposphere (Fig. 8). However, even if the NO mixing ratio is low compared to the reference run and to the observations, it reaches 1000 pmol mol  Despite an increase in the CG rate by a factor 10 between CG CST and CG PR93 (the total number of flashes being constant), NO produced by lightning is not significantly impacted. It causes a decrease of 3.5% in the flux through the anvil and an increase of 2.8% in the total amount of nitrogen produced during the storm lifetime. The flux through the anvil is decreased when there are more CG flashes because even 5 if more NO molecules are produced, they are produced at lower altitudes and it takes more time to be transported in the anvil.

Altitude of the upper and lower modes for the bimodal distribution
Typically the IC discharge has a bilevel structure (Shao and Krehbiel, 1996), which is correlated with the main negative and upper positive charge regions of the storm. Rison 10 et al. (1999), Thomas et al. (2001) and Krehbiel et al. (2000) deduced from the use of the Lightning Mapping Array in New Mexico and Oklahoma that the main negative charge is located in the middle level at 5-6 km m.s.l. while the upper positive charge is centered at about 10-12 km altitude. By examining electric field soundings in different thunderstorms, Stolzenburg et al. (1998) found that the temperature at the center of the 15 main negative charge region varies from −4 • C to −32 • C with a mean value of −15.7 • C. Thus, the production of NO x and its subsequent distribution by updrafts and downdrafts may be sensitive to the altitudes chosen for the bimodal distribution.
Two sensitivity tests have been performed to investigate the impact of a change in the altitude of the two levels where flash segments are preferentially distributed on the  (Fig. 6). However, it is slightly higher in ISO UP 45 and ISO LOW 20 than in REF.
The LNO x produced in the lower mode for ISO LOW 20 is more readily available at high altitude than in the REF simulation.
The impact of moving the altitude of the upper and lower modes in the vertical distribution of the flash segments has an impact on the flux of NO x through the anvil 5 (Table 2). Lowering the upper mode and raising the lower mode both result in an increase of the NO x flux through the anvil. DeCaria et al. (2000) tested different temperatures (−40 • C and −30 • C) for the upper isotherm of their bimodal distribution in the 12 July 1996 STERAO storm. They found that by lowering the upper mode of the flash distribution the maximum of the NO plume occurs at a lower altitude. 10

Flash length
Previous studies have estimated the flash length to be in the range 20-50 km . For the 10 July 1996 STERAO storm analysis, Defer et al. (2001) have shown that the flash length in this storm is not constant and varies from 0.02 and 474 km. They have concluded that the average value of the flash length is 19 km but 15 is 34 km if short duration flashes (flashes <1 ms) are not considered. Pinty and Barthe (2008) conducted an ensemble of simulations of two idealized electrified storms with a cloud-resolving model coupled to an explicit electrical scheme. They showed that the number of segments in an individual flash is highly variable from flash to flash but looks like a lognormal distribution. 20 To explore the importance of the flash length to the production of NO from lightning, a simulation is performed using a constant flash length of 21 km (LENGTH 21, Table 1). The value of 21 km corresponds to the mean flash length simulated when a lognormal distribution for the flash length is used and when short duration flashes are taken into account. The aim of these sensitivity tests is to investigate the impact of using a 25 constant or a varying flash length.
Using the more realistic lognormal distribution for the flash length (REF) leads to results fairly similar to a constant value (LENGTH 21) Fig. 8 where the LENGTH 21 and REF profiles differ between 5.5 and 10.5 km m.s.l. These mean profiles show that LENGTH 21 tends to produce more NO molecules in the lower part of the cloud. This difference can be due to the treatment of the short flashes <1 km since they are made up of only one point. Due to the bimodal vertical distribution that is used, the short flashes tend to be produced in the upper part of the cloud (∼12 km m.s.l.). This is slightly higher than the analysis of Defer et al. (2001) who showed that the VHF sources with strong radiation were mostly located in altitude between 7.5 and 10.5 km m.s.l. The increase of sources at lower altitude increases the NO production per flash from 121.3 moles fl. −1 in REF to 125.0 moles fl. −1 in LENGTH 21 (=3.1%) because of the pressure dependence in the NO production (Eq. 2).

Short duration flashes
In the 10 July 1996 STERAO storm, Defer et al. (2001) defined the short duration flashes as discharges lasting less than 1 ms while the mean flash duration is 240 ms in this storm. They reported that between 22:50 and 23:30 UTC up to 45% of the 15 recorded flashes were short duration flashes. Analyzing data from the lightning mapping array, Harlin et al. (2003) identified three different types of short duration flashes lasting less than 80 ms: the isolated events, the precursor events, and the high source power events. Several studies have shown that taking into account or not the short duration flashes does not modify the flash rate 20 trend but can significantly change its magnitude (Lang et al., 2000;Wiens et al., 2005). Whether these short duration flashes should be considered or not as part of the total lightning flash rate and as a source of nitrogen oxides is a matter of debate.
To investigate the possible impact of short duration flashes in the LNO x production, a sensitivity simulation (NO SDF) has been performed. In this simulation, the lognormal 25 distribution is used but it is hypothesized that the short duration flashes do not produce NO molecules. It is worth recalling that the reference run (REF) assumes that short duration flashes produce as many NO molecules per meter of flash as "normal" flashes. It is assumed that the flash length of a short duration flash is 1 km, and that 36% of the flashes that are 1 km in length are short duration flashes.
In the NO SDF simulation 4253 flashes have been triggered among which 699 were short duration flashes. If short duration flashes are not considered in the total flash rate, 146.2 moles of NO are produced per flash on average. If the total amount of NO 5 produced by lightning is divided by the total flash rate (4253), one flash produces on average 120.6 moles of NO.
If short duration flashes are assumed to produce as many molecules per meter of flash as a normal flash, their impact on the NO mixing ratio is not significant. In Figs. 5, 7 and 8, the NO SDF and REF simulations give similar results. As noted in the previous 10 section, the simulated short flashes are mostly produced at altitudes between 11 and 12 km m.s.l., which is higher than in the observations.
Even though this sensitivity test concludes that short duration flashes are not very important for the LNO x production, it must be kept in mind that little is known about the physics of these discharges. Thus the impact of the short duration flashes on the NO 15 budget should be studied further when more measurements of short duration flashes are available in conjunction with chemistry data.

Spatial distribution of the NO molecules
Previous parameterizations (Pickering et al., 1998;DeCaria et al., 2000DeCaria et al., , 2005 assumed that the NO molecules produced by the lightning flashes are instantly diluted 20 over the whole volume of the cloud or in the 20 dBZ contour. This instant dilution does not produce the NO peaks observed by instruments onboard airplanes (Barth et al., 2007b;Ott et al., 2007). The goal of these sensitivity tests is to investigate how the NO mixing ratio is affected by an instantaneous dilution of the LNO x source.
Four are distributed uniformly in the 20 dBZ volume above the −15 • C isotherm (VOL 20). Third, the simulation VOL CLD follows Pickering et al. (1998) where the LNO x is distributed over the entire cloud above the −15 • C isotherm. In these first three simulations the LNO x source is only distributed in the detected electrified cell in proportion to the flash rate of each cell. Finally, the last simulation (VOL ALLCLD) follows Pickering et al.

5
(1998) like VOL CLD but there is no cell identification. In these four simulations, the flash length is held constant (21 km) since the sensitivity test in Sect. 5.4 has shown no significant impact of the flash length distribution. In the VOL 20 DC, VOL 20, VOL CLD and VOL ALLCLD, the NO molecules produced by a single flash are distributed over 9400, 4300, 17 700 and 23 400 grid points on average, respectively. In contrast, the REF simulation places NO by a single flash on ∼20 grid points.
To identify the impact of the bimodal distribution on the NO mixing ratio, the results of the VOL 20 DC and VOL 20 simulations are compared. In Fig. 5, the cross section of NO mixing ratio across the anvil for VOL 20 DC is fairly similar to REF.
The region of high NO mixing ratio is much larger in VOL 20 than in VOL 20 DC. The 15 study of VOL CLD and VOL ALLCLD results gives some insight on the impact of the cell identification. Figure 5 shows that three different spots of NO mixing ratio higher than 540 pmol mol −1 are visible for VOL CLD. In VOL ALLCLD there is only one large region where the NO mixing ratio exceeds 540 pmol mol −1 . When the vertical distribution of the NO molecules is uniform (VOL 20, VOL CLD, VOL ALLCLD), the region 20 with high NO mixing ratio across the anvil (Fig. 5) increases with the volume in which NO molecules are instantly diluted. Moreover this volumetric distribution of the NO molecules over a large area does not allow the peaks in the transects to be reproduced near the convective core where an enhancement in the NO mixing ratio is observed (Fig. 7). In these transects (Fig. 7), the NO mixing ratios near the convective core for 25 VOL 20, VOL 20 DC, VOL CLD and VOL ALLCLD are lower than observations, but are similar in magnitude to the observations 50 km downwind of the convective core.
The NO x flux in all these simulations is larger than that determined from the observations and in the REF simulation (Table 2) VOL 20 DC,VOL 20, VOL CLD, VOL ALLCLD simulations are placed at thousands of grid points, many of which are near and within the anvil (see 20 dBZ contour in Fig. 1). The VOL 20 DC NO x flux is similar to that predicted by other models simulating this storm (Barth et al., 2007b) indicating that their overprediction of the NO x flux may be due to the larger region of the NO lightning source compared to the filamentary re-5 gion in the REF simulation. Despite the significant differences in the NO x flux, the total amount of NO produced from lightning is fairly similar among these 4 sensitivity simulations and the REF simulation.

Production of NO per flash
A recent review of lightning production of NO reported that the LNO x production rate per 10 flash length varies from 1×10 21 to 13×10 21 molecules m −1 . This large range of values has been obtained in different storms by different ways: laboratory experiments, modeling studies with different models and different parameterizations or airborne measurements. Barth et al. (2007b) have also reported a large range of NO moles produced per flash (between 36 and 465 moles of NO per 15 IC flash) for the same storm simulated here but using different models and different parameterizations.
The impact of the number of NO molecules produced per flash length unit is investigated in this section. All the simulations use the Wang et al. (1998)  This sensitivity analysis is the one that impacts most the NO mixing ratio both in the convective core and in the anvil. Figure 5 shows that the NO mixing ratio in the anvil in the WANG 10 simulation is far too high compared to observations. In WANG 10,25 there is a large region with NO mixing ratio higher than 540 pmol mol −1 that extends over 50 km horizontally. Conversely, the WANG 1 simulation exhibits NO mixing ratio in the anvil less than 300 pmol mol −1 . The large difference in the NO mixing ratio between the WANG 1, REF and WANG 10 simulations can be also seen in Fig. 7. The NO moles per flash, the NO peak value, and the total amount of LNO x produced during the storm increase quasilinearly as the a and b parameters increase ( Table 2). The NO mixing ratio increase is not exactly linear since it is impacted by the NO x chemistry. The values of the LNO x 5 average profiles are affected by changing the a and b parameters, but the trends are the same for the three simulations.
Two additional simulations are performed to study the impact of a larger NO production rate by CG flashes. Pickering et al. (1998) used the values given by Price et al. (1997) (111 mol(NO) IC −1 and 1113 mol(NO) CG −1 ) to simulate seven case studies 10 representing different environments. Following the method of Price et al. (1997) to determine the amount of NO produced per flash in the 21 July EULINOX storm, Fehr et al. (2004) used the value of 350 mol(NO) CG −1 . The 480 moles of NO produced per IC were obtained by a fit to observations. DeCaria et al. (2005) who simulated the 12 July 1996 STERAO storm using the GCE model with the Price et al. (1997) method 15 determined that the IC mean production rate is 460 moles of NO per IC flash. They also concluded from their sensitivity tests that IC and CG flashes produce the same amount of NO per flash. However, they assumed an instantaneous dilution of the NO molecules and only compared with the column NO x mass from observations. Ott et al. (2007) deduced from their study that 360 moles of NO per flash matches best with ob-20 servations. Thus, the value of 121 moles of NO produced per flash derived here for the 10 July 1996 STERAO storm is in the lower range of values used or deduced from previous CRM studies except when an explicit electrical scheme is used. In previous parameterizations, the number of NO moles per flash was mostly deduced from Price et al. (1997) who suggested a ratio of 10 between the production rate per CG and IC,25 or it was deduced from observations, but only in one part of the storm. Few of the previous studies have used a geometric approach for the flash segments distribution or compared the NO mixing ratio to observations in different regions of the storm.
The impact of the NO production rate per IC and CG is investigated with simulations Introduction PROD CG 2 and PROD CG 10. In these two sensitivity tests, the CG flash rate is given by the NLDN observations. In PROD CG 2, the a and b coefficients of Eq. (2)  sitivity tests, the number of NO molecules produced per meter of flash in a CG flash for a particular altitude is multiplied by 2 (PROD CG 2) and 10 (PROD CG 10) compared to an IC flash. Figures 5,7 and 8 show no significant influence of increasing the NO production rate per CG flash length by 2 or 10. The NO x flux through the anvil in PROD CG 2 and PROD CG 10 is higher than in REF where CG flashes are not considered (Table 2). Comparing CG OBS, PROD CG 2 and PROD CG 10 shows an increase in the average number of NO moles per flash and in the NO x flux through the anvil. The average number of NO moles per flash is increased by 0.5% and 5.4% between CG OBS and PROD CG 2, and between CG OBS and PROD CG 10, respectively. However, be-15 cause only a small fraction of flashes are CG in this storm, it is difficult to conclude whether the IC and CG NO production ratio is important or not.

Conclusions
A new lightning-produced NO x parameterization has been developed for use in cloudresolving models. This parameterization is based on three factors: flash rate, spatial 20 distribution and number of NO molecules produced per flash. The flash type is assumed to be intracloud based on observations of the storm simulated in this study. The new lightning-produced NO x parameterization has three unique characteristics. First, a vertical velocity threshold is used to identify the cells that can produce lightning. Second, the flash rate is estimated from the non-precipitation and precipitation ice mass flux product, which has the benefit of containing information on both the dynamical and microphysical state of the storm. Third, the source location of the NO is filamentary  Ott et al. (2007). This parameterization has been tested on the 10 July 1996 STERAO storm. The predicted flash rate is in good agreement with observations for both the magnitude and trend of the flash rate. The distribution of the NO mixing ratio in the anvil, and the NO mixing ratio values near the convective core and in the anvil agree well with the aircraft measurements taken in the 10 July 1996 5 STERAO storm. Several sensitivity tests have been performed to determine the parameters that the lightning-produced NO x parameterization is most dependent upon. The flash rate is an important factor. For these simulations, using the predicted flash rate based on the maximum vertical velocity (Price and Rind, 1992) produced significantly different 10 results from the observations for both flash rate and NO mixing ratios. Using observed flash rate data when available is desirable, but the modeled storm dynamics and physics must reproduce the observations well. For the 10 July 1996 STERAO storm, the duration of each stage of the storm (multicell, transition, supercell) was correlated with the lightning flash rate. Thus, for this case, it is also important to simulate the duration of each stage accurately to have corresponding microphysics and dynamics with observed flash rate.
The spatial placement of the lightning-NO x source is also an important factor. Placement of NO in large regions of the storm such as within the 20 dBZ contour (DeCaria et al., 2000) or the cloud boundary at temperatures <−15 • C (Pickering et al., 1998) re-20 sulted in high NO mixing ratios in the anvil region of the storm. Further these previous studies give much higher magnitudes of NO produced per flash than estimated in the current study using a virtual flash path (i.e. small NO source region). The sensitivity of NO to the vertical distribution of the source of NO from lightning was substantial for NO mixing ratios in the anvil. The resulting vertical profiles of NO x from lightning 25 always had two altitudes of peak NO even when the NO source was uniformally distributed vertically. Thus, lightning-NO x parameterizations in large scale models (where the convection is parameterized) should represent this NO source in a similar bimodal fashion. Introduction The magnitude of the NO produced per length of flash strongly influences the NO mixing ratio near the convective cores and in the anvil region. Tests with twice the NO produced and with one-fifth the NO produced show a quasi-linear response of NO mixing ratio.
The lightning-produced NO x parameterization was shown to not be very sensitive to 5 a number of parameters including the cloud-to-ground to intra-cloud lightning ratio, the altitudes (within ∼1 km) of the peaks of the NO source, the lightning flash length, the presence of short duration flashes, and the ratio of NO produced by CG lightning to NO produced by IC lightning. Some of these parameters that are sensitive to flash type will likely have a greater influence on NO produced from lightning in other storms because 10 the 10 July 1996 STERAO storm primarily contained intracloud flashes. These sensitivity simulations can provide guidance for both future modeling and measurement strategies. The storm simulated for this study represented a severe storm occurring in a regime of high cloud bases, high wind shear and convective available potential energy. Storms from other regimes e.g. the southern Great Plains with 15 warmer cloud bases, the airmass thunderstorms typically produced by sea breezes, and tropical storms will be simulated. Despite a number of recent field experiments (STEPS, TROCCINOX, SCOUT-O3, and AMMA), there is still a need for concomitant observations of cloud dynamics, microphysics, electrical activity, and chemistry. In particular, measurements of NO x near the convective core should be done to pinpoint the 20 spatial extent of the lightning source of NO. The modeling and experimental studies of a variety of storms should then be able to elucidate the amount of NO typically produced from a thunderstorm. Barthe, C., Pinty, J.-P., and Mari, C.: Lightning-produced NO x (10), 2005. 6620 Lang, T. J., Rutledge, S. A., Dye, J. E., Venticinque, M., Laroche, P., and Defer, E.: Anomalously low negative cloud-to-ground lightning flash rates in intense convective storms observed during STERAO-A, Mon. Weather Rev., 128, 160-173, 2000. 6623, 6624, 6626 Table 1. Summary of the sensitivity tests. CG ratio is the fraction of total flashes that are cloudto-ground. MF: derived from the non-precipitation and precipitation ice mass flux product; Obs.: Observations; PR92: Price and Rind (1992) for the global system; PR92/cell: Price and Rind (1992) Table 1.