{'fraction_non_alphanumeric': 0.03987438583480115, 'fraction_numerical': 0.031658477797348664, 'mean_word_length': 4.546792646869777, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 4, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Assessing and improving cloud-height based parameterisations of global lightning flash rate, and their impact on lightning-produced NOx and tropospheric composition in a chemistry-climate model

Although lightning-generated oxides of nitrogen (LNOx) account for only approximately 10% of the global NOx source, it has a disproportionately large impact on tropospheric photochemistry due to the conducive conditions in the tropical upper troposphere where lightning is mostly discharged. In most global composition models, lightning flash rates used to calculate LNOx are expressed in terms of convective cloud-top height via the Price and Rind (1992) (PR92) parameterisations for land and ocean, where the oceanic parameterisation is known to greatly underestimate flash rates. We 15 conduct a critical assessment of flash-rate parameterisations that are based on cloud-top height and validate them within the ACCESS-UKCA global chemistry-climate model using the LIS/OTD satellite data. While the PR92 parameterisation for land yields satisfactory predictions, the oceanic parameterisation, as expected, underestimates the observed flash-rate density severely, yielding a global average over the ocean of 0.33 flashes s-1 compared to the observed 9.16 flashes s-1 and leading to LNOx being underestimated proportionally. We formulate new/alternative flash-rate parameterisations following Boccippio’s 20 (2002) scaling relationships between thunderstorm electrical generator power and storm geometry coupled with available data. The new parameterisation for land performs very similar to the corresponding PR92 one, as would be expected, and the new oceanic parameterisation simulates the flash-rate observations more accurately, giving a global average over the ocean of 8.84 flashes s-1. The use of the improved flash-rate parameterisations in ACCESS-UKCA changes the modelled tropospheric composition—global LNOx increases from 4.8 to 6.6 Tg N yr-1; the ozone (O3) burden increases by 8.5%; there 25 is an increase in the midto upper-tropospheric NOx by as much as 40 ppt (by volume); a 13% increase in the global hydroxyl (OH); a decrease in the methane lifetime by 6.7%; and a decrease in the lower tropospheric carbon monoxide (CO) by 3–7%. Compared to observations, the modelled tropospheric NOx and ozone in the Southern Hemisphere and over the ocean are improved by this new flash rate parameterization.


Introduction
Oxides of nitrogen (NO x ≡ NO (nitric oxide) + NO 2 (nitrogen dioxide)) play an important role in tropospheric chemistry by acting as a precursor to ozone (O 3 ) and the hydroxyl (OH) radical, which are the principal tropospheric oxidants (Labrador et al., 2005). As a greenhouse gas, O 3 is most potent in the upper troposphere, whereas near the Earth's surface it is an air pollutant, adversely impacting human health and plant productivity. The OH radical plays a critical role in the chemical 5 cycles of many trace gases, including methane (CH 4 ) and carbon monoxide (CO), in the atmosphere.
Lightning is the dominant source of NO x in the middle to upper troposphere, and the only direct natural source remote from the Earth's surface Banerjee et al., 2014). Lightning predominantly occurs over land in the tropics associated with deep atmospheric convection. The extreme heat in a lightning flash channel in the atmosphere, which can extend over tens of kilometres, allows for the dissociation of nitrogen (N 2 ) and oxygen (O 2 ) molecules into free N and O 10 within the flash channel. These in turn react with ambient N 2 and O 2 to produce NO, which remains after the lightning channel cools. There is a conversion between NO and NO 2 , during which ozone is generated in the presence of HO 2 and organic peroxy radicals, collectively called RO 2 ( . A large uncertainty in the amount of NO x produced by lightning has been reported, with most global estimates ranging between 2 and 8 Tg nitrogen (N) per year . 15 Although lightning accounts for only approximately 10% of the global NO x source, the lightning-generated NO x (referred to as LNO x ) has a disproportionately large contribution to the tropospheric burdens of O 3 and OH (Murray, 2016). For example, although LNO x emissions are of similar magnitude as those from biomass burning or soils, their contribution to the total tropospheric ozone column is about three times larger (Dahlmann et al., 2011). This is because in the middle to upper troposphere where LNO x is directly released, the O 3 production efficiency per unit NO x is much higher (∼ 100 molecules O 3 20 per molecule NO x ) than that near the surface (∼ 10-30 molecules O 3 per molecule NO x ) as a result of the higher amount of UV radiance, lower concentrations and longer lifetimes of NO x (a few days rather than hours), and lower temperatures affecting ozone loss chemistry at such altitudes (e.g. Williams, 2005;Dahlmann et al., 2011). Of any major emission source, variability in global mean OH is most sensitive to LNO x (Murray, 2016). Using a global chemistry transport model, Labrador et al. (2005) observed marked sensitivity of NO x , O 3 , OH, nitric acid (HNO 3 ) and peroxyacetyl-nitrate (PAN) to 25 the magnitude and vertical distribution of LNO x . Modelling studies by Grewe et al. (2007) and Dahlmann et al. (2011) found that of all the major sources of NO x , LNO x is the dominant source for tropospheric ozone (up to 40%) in the tropics and Southern Hemisphere. LNO x is also important when studying natural variability of tropospheric composition because lightning occurrence is influenced by natural climate variability drivers such as El Niño and La Niña in the tropics. Similarly, potential changes in 30 deep convection as a result of future climatic change has a bearing on LNO x and thus tropospheric ozone and associated radiative feedbacks, and some of these have been explored through modelling (e.g. Banerjee at al., 2014Banerjee at al., , 2018Iglesias-for projections of lightning-induced fire activity (Krause et al., 2014).
A realistic representation of LNO x source strength and its global distribution is thus of vital importance for understanding tropospheric chemistry and its impacts. In most global circulation models used for climate applications, convection is diagnosed/parameterised (i.e. clouds are not explicitly resolved) with limited cloud microphysics. A pragmatic way to 5 predict lightning flash rates in these models is to use parameterisations based on simple physics of electrical charge and correlations between the flash rate and appropriate parameters representing convection.
The most common conceptual model used to compute the amount of LNO x in global models is where F is the lightning flash rate and is the amount of NO produced per flash. This calculation is carried out in atmospheric models by first calculating the lightning flash rate within a model grid column at every model time step, 10 partitioning it into intracloud (IC) and cloud-to-ground (CG) flash-rate components, applying an amount of NO mass produced per flash, and then vertically distributing the calculated NO mass in the column.
Various approaches to estimate lightning flash rate in global models have been followed in the past. Price and Rind (1992, hereafter PR92) developed simple empirical parameterisations for calculating lightning flash rate in terms of convective cloud-top height over land and ocean. The model of PR92 was based on the assumption of an electric dipole structure for a 15 thunderstorm with two equal but opposite charge volumes, separated by a distance of the order the vertical cloud dimension as developed by Vonnegut (1963) and Williams (1985). Flash-rate parameterisations based on convection parameters other than cloud-top height have also been developed, e.g. convective precipitation and upward mass flux (Allen and Pickering, 2002); convective available potential energy (CAPE) (Choi et al., 2005), cold-cloud depth (Futyan and Del Genio, 2007;Yoshida et al., 2009), maximum vertical velocity and updraft volume (Deierling and Petersen, 2008), upward cloud ice flux 20 (Finney et al., 2014;2016), product of CAPE and precipitation rate (Romps et al., 2014;Zhu et al., 2019), and multiparameter regression fits (Luo et al., 2017).
The PR92 parameterisations are by far the most widely used ones by default in global chemistry transport models (CTMs) and coupled chemistry-climate models, such as GEOS-Chem (Hudman et al., 2007), MOZART (Emmons et al., 2020), CAM-Chem (Lamarque et al., 2012), ECHAM5/MESSy (Jöckel et al., 2006), and UM-UKCA (Archibald et al., 2020), 25 perhaps primarily because they are based on convective cloud-top height which can be easily diagnosed from a model's convection scheme. Additionally, with its use of an electric dipole structure for a thunderstorm, the framework underlying the PR92 parameterisations has some theoretical support. The PR92 parameterisations also perform reasonably well. For example, Clark et al. (2017) tested flash-rate parameterisations based on cloud-top height, cold cloud depth (CCD), mass flux, convective precipitation rate, and cloud-top height with column-integrated cloud droplet number concentration, in a 30 global model, and found that the PR92 parameterisations had the best correlation with the observations, closely followed by the CCD based parameterisation of Yoshida et al. (2009). The PR92 scheme had a higher value of the spatial standard deviation compared to observations due to a large land-ocean contrast in this parameterisation. The quantity CCD is defined as the convective cloud-top height minus the freezing level. The thunderstorm data analysis presented by Price and Rind (1993) indicates that freezing levels remains relatively constant compared to CCD values, meaning that it is largely the 5 cloud-top height that provides the variation in the lightning flash rate in the CCD based scheme, which suggests that the cloud-top height and the CCD based schemes would perform very similarly. Finney et al. (2014) found that the PR92 lightning flash parameterisation had considerable biases compared to satellite observations of spatial flash density but performed 2 nd best (behind their parameterisation based on ice flux) out of the 5 parameterisations tested. A recent comparison by Gordillo-Vázquez et al. (2019) of six flash-rate parameterisation schemes showed a relatively good 10 performance by those based on cloud-top height. However, a number of global modelling studies has demonstrated that the PR92 parameterisations underestimate the lightning flash density over the ocean compared to satellite observations (Allen and Pickering, 2002;Tost et al., 2007;Finney et al., 2014Finney et al., , 2016Clark et al. (2017).
In this paper, we critically assess the performance of currently used flash-rate parameterisations for land and ocean based on convective cloud-top height. We also derive new/alternative flash-rate parameterisations following Boccippio's (2002) 15 (Bo02) scaling relationships founded on theory linking thunderstorm electrical generator power and storm geometry applied to the available data. We implement in a global chemistry-climate model, viz. ACCESS-UKCA (Woodhouse et al., 2015), (a) these new parameterisations, (b) an oceanic flash-rate parameterisation suggested by Michalon et al. (1999) (Mi99), and (c) use the standard PR92 flash-rate parameterisation as a default. The model results are tested using global satellite data of flash density from the Optical Transient Detector (OTD) and Lightning Imaging Sensor (LIS) (Cecil et al., 2014), which are 20 available as global climatological, inter-annual and seasonal distributions. The veracity of the new modelled LNO x estimates are examined through comparison of modelled and reanalysis tropospheric NO 2 column amounts. The impacts of the modelled LNO x based on the new flash-rate parameterisations on tropospheric composition involving NO x , ozone, the hydroxyl radical, methane lifetime, and carbon monoxide are examined, including comparison with observations where available or appropriate. 25

The ACCESS-UKCA global chemistry-climate model
We use the United Kingdom Chemistry and Aerosol (UKCA) global atmospheric composition model (Abraham et al., 2012; http://www.ukca.ac.uk) coupled to ACCESS (Australian Community Climate and Earth System Simulator) (Bi et al., 2013;Woodhouse et al., 2015). In our simulations, ACCESS is essentially the same as the U.K. Met Office Unified Model (UM) (vn8.4) since the ACCESS specific ocean and land-surface components are not invoked as the model is run in atmosphere-30 only mode with prescribed monthly-mean sea surface temperature (SST) and sea ice fields, and the UM's original landsurface scheme (viz. JULES) is used. The atmosphere component of the UM vn8.4 is the Global Atmosphere (GA 4.0) (Walters et al., 2014). The UKCA configuration used here is the so-called StratTrop (or Chemistry of the Stratosphere and Troposphere (CheST)) (Archibald et al., 2020), which, in essence, combines the tropospheric chemistry scheme described by O'Connor et al. (2014) and the stratospheric chemistry as described by Morgenstern et al. (2009). The tropospheric chemistry scheme includes the O x , HO x and NO x chemical cycles, and the oxidation of CO, methane, and other volatile organic carbon species (e.g. ethane, propane and isoprene). The Fast-JX photolysis scheme (Neu et al. 2007;Telford et al. 5 2013) is used, and ozone is coupled interactively between chemistry and radiation. The aerosol component includes sulphur chemistry. The total number of reactions, including aerosol chemistry, is 306 across 86 species.
The atmospheric model has a horizontal resolution of 1.875° in longitude and 1.25° in latitude, and 85 staggered terrainfollowing hybrid-height levels extending from the surface to 85 km (the so-called N96L85 configuration). The vertical resolution decreases with height, with the lowest 65 levels (up to ∼ 30 km) lying within the troposphere and lower 10 stratosphere. The model dynamical timestep is 20 minutes, the UKCA chemical solver is called every 60 min. It is a symbolic backward Euler solver with Newton-Raphson iteration and runs to convergence, halving the step when required (see Esentürk et al. (2018) for more details). A global monthly-varying emissions database for reactive gases and aerosols is used, which includes both anthropogenic, biomass burning and natural components, whereas for carbon dioxide, methane, nitrous oxide and ozone depleting substances, concentrations are prescribed instead of emissions (Woodhouse et al., 2015). 15 We have backported the LNO x subroutines from a more recent version of the model (at UM vn11.0) to vn8.4. This is to ensure that any refinements that may have occurred in the new version are used in our study. However, it is found that there are no major LNO x parameterisation differences between the two versions, with the new version continuing to use the original PR92 flash-rate parameterisations, except that the amount of NO produced per flash is taken to be the same for both IC or CG flashes in the new version, which leads to small changes in the LNO x production compared to the old version. 20 The ACCESS-UKCA setup used here incorporates some additional modifications compared to the base UM-UKCA version 8.4, and these include: • Dry deposition of ozone to the ocean is now based on a process-based scheme developed by Luhar et al. (2018) (the Ranking 1 configuration in their Table 1).
• Dry deposition of all relevant species is applied at the lowest model level, instead of it being distributed through the 25 vertical extent of the atmospheric boundary layer (Luhar et al., 2017).
• The coefficient 8.53 is corrected to 4.8 in the following branching ratio expression used to compute the rate constant k 2 (cm 3 molecule -1 s -1 ) for the chemical reaction HO 2 + NO → HNO 3 : / = {[(530 ± 10)/ ] + 4.8 x 10 × -(1.73 ± 0.07)}/100 , where pressure P is in hPa and temperature, T is in K, and k is the rate constant (cm 3 molecule -1 s -1 ) for the included chemical reaction HO 2 + NO → NO 2 + HO. The coefficient 8.53 is correct when P is expressed in Torr 30 (Butkovskaya et al., 2007), but it should be 4.8 when P is in hPa 1 . Additionally, in the parameterisation = × 10 exp(270/ ), the model uses = 3.3 for this reaction, but for k in the branching ratio Eq. (2) it uses = 3.6, which we change to 3.3 for consistency.
The above changes lead to an increase in the modelled tropospheric ozone burden by about 12% (the first two changes by ∼ 7% and the last by ∼ 5%) to 284 Tg O 3 and this increase is towards the global modelling average (see Section 4.2). 5 In the model, the tropopause is calculated every timestep. In the extratropics (latitude ≥ |28|°), the tropopause is the pressure level of the 2-PVU (potential vorticity units) surface, and in the tropics (latitude ≤ |13|°) it is the pressure level of the 380 K potential temperature isentropic surface. Between the two latitudes, a weighted average of the two definitions is used following the method of Hoerling et al. (1993).

Implementation of lightning flash-rate parameterisation in the model 10
In ACCESS-UKCA, the convective cloud bottom level (H b ) and top level (H) are diagnosed on a timestep basis from the UM convection scheme. This scheme represents the sub-grid scale transport of heat, moisture and momentum associated with cumulus clouds within a grid box. The scheme (at GA4.0) is summarised by Walters et al. (2014); it uses a mass flux convection scheme based on Gregory and Rowntree (1990) with various extensions to include downdraughts and convective momentum transport. It consists of three stages: (a) diagnosis to determine whether convection is possible from the boundary 15 layer, (b) a call to the shallow or deep convection scheme for all points diagnosed deep or shallow by the first step, and (c) a call to the mid-level convection scheme for all grid points. H b is taken to be the air parcel ascent start level and H is set to be the top of the ascent. On each timestep, the flash rate is then determined by using the calculated values of H b and H at each grid point. Examples of evaluation of the distribution of cloud depths simulated by the UM include those by Klein et al. (2013) and Hardiman et al. (2015). 20 A minimum convective cloud scale needs to be specified for it to constitute a thunderstorm. In our model, a minimum convective cloud thickness ( − ) of 5 km is required for the lightning NO x to be activated. The selected threshold of 5 km is very similar to observations of the smallest vertical scale of thunderstorms presented by several researchers, viz. Rind (1992, 1993), Molinié and Pontikis (1995), and Ushio et al. (2001). While prescribing a minimum convective cloud thickness of 5 km for lightning is somewhat arbitrary, having no such threshold value is unrealistic because then it would be 25 implicitly assumed that a convective cloud always translates to a thunderstorm, and this would likely lead to unrealistically 1 This inconsistency stems from a typo in the document "IUPAC Task Group on Atmospheric Chemical Kinetic Data Evaluation -Data Sheet I.A3.45 NOx15" (http://iupac.pole-ether.fr/htdocs/datasheets/pdf/NOx15_HO2_NO.pdf) in which a constant of 8.53 instead of 4.8 is mistakenly specified in the branching ratio expression in the "Rate coefficient data" table.
high flash rates. (We found that increasing or decreasing the minimum cloud thickness value by 1 km from 5 km resulted in a change of -3.2% and 1.7%, respectively, in the modelled global flash rate using the PR92 scheme.) The model diagnosed H is used in the flash-rate parameterisations, which calculate the lightning flash rate (F, flashes per minute) as a function of H (km) (see Section 3.2).
The PR92 flash-rate expressions (discussed later in Section 3.2) were developed based on observations of individual 5 thunderstorms. Price and Rind (1994) (see their Figure 1) developed a spatial calibration factor (c) to adjust these expressions for varying model resolutions: where is the longitudinal resolution and is the latitudinal resolution of the model (in degrees).
The total flash rate (f) within a grid cell (flashes per minute per grid box) is then calculated as where A is the model grid box area (which is a function of latitude), is the area of model grid box centred at 30°N (Allen 10 and Pickering, 2002), and the subscripts L and O to refer to land/continental and ocean/marine, respectively. Any model grid cell that has a non-zero land surface fraction is considered land for the purposes of lightning NO x calculation. Conversely, only grid cells with 100% water surface coverage are considered ocean.
We use Eq. (4) along with Eq. (3) to calculate the total flash rate , which is then apportioned into cloud-to-ground (CG) and intracloud (IC) flash rates using an empirical parameterisation for the ratio z R = IC/CG developed by Price and Rind 15 (1993) (PR93) based on thunderstorm observations in the western United States. In this parameterisation, z R increases as a function of the thickness (dH) of the cold cloud region in thunderstorms (from 0°C to cloud top), and dH is parameterised as a decreasing function of latitude. The PR93 parameterisation has been used frequently, with further validation for case studies reported by Pickering et al. (1998) and Fehr et al. (2004). Allen and Pickering (2002) and Grewe et al. (2001) used it in global atmospheric chemistry models, with the former evaluating it for cases in the US. The averaged values of z R and the 20 CG to total flash ratio obtained from the PR93 parameterisation in the present study are 3.14 and 0.24, respectively. These values are comparable to z R ∼ 4 and the CG to total flash ratio ∼ 0.2 obtained by Barthe and Barth (2008) using dH calculated directly from modelled cloud temperature and total hydrometeor mixing ratio in the PR93 parameterisation, and to z R = 2.64-2.94 obtained by Boccippio et al. (2001) using satellite-and ground-based lightning observations over the continental US.
Using IC/CG measurements, Bond et al. (2002) derived a parameterisation for z R as a linearly decreasing function of latitude 25 and obtained z R = 3.76 and the CG to total flash ratio = 0.21 over the tropics (35°N-35°S).

Production of NO per flash
Both the moles of NO produced per flash, , and the variation of this parameter between CG and IC flashes are poorly constrained by atmospheric observations. The overall rate, , regulates the amount of nitrogen oxides produced by lightning, whereas the variation of between CG and IC flashes regulates the level at which lightning nitrogen oxides are introduced into the atmosphere, both are critical variables. In this study, is set at × 10 26 molecules NO per flash where 5 the scaling factor = 2 by default irrespective of whether a flash is IC or CG, which is equivalent to 330 moles NO per flash. Assuming a mean energy release of 0.67 GJ per IC flash and 6.7 GJ per CG flash (Price et al., 1997), with 24% of the total modelled flashes being CG, the production of 330 moles NO per flash corresponds to 9.4 × 10 16 molecules NO J -1 . If we use a mean energy release of 0.9 GJ per IC flash and 3.0 GJ per CG flash based on Schumann and Huntrieser (2007), then the NO production is calculated to be 14.2 × 10 16 molecules NO J -1 . 10 In bottom-up models, in addition to flash rate, is a key source of uncertainty, with a review by Schumann and Huntrieser (2007) suggesting a range of ∼ 33-660 with an average of 250 moles NO per flash, and similarly 70-700 moles . This large range, in part, reflects spatial variation in the frequency and uncertainty in the yield of CG and IC flashes, which may involve a varying level of dependencies on environmental variables such as peak current, rate of energy dissipation, channel length, air density, and strokes per flash (Murray et al., 2016). We assume that both CG and IC flashes yield the same amount of NO, which follows studies such as DeCaria et al. (2005), Ridley et al. (2005), Ott et al. (2007Ott et al. ( , 2010 and Cummings et al. (2013). On the other hand, some studies consider or find that the less frequent CG flashes yield a greater amount of NO per flash than IC flashes (Price et al., 1997;Koshak et al., 2014;Luo et al., 2017), A few studies suggest that may not be constant over the globe, with higher production rates in 30 extratropics than tropics (Huntrieser et al., 2008;Murray et al., 2012) and globally variable production rates (Miyazaki et al., 2014). Differences in land and ocean production rates have also been noted. Boersma et al. (2005) found that land flashes were ~1.6 times more productive than those over the ocean, and conversely Allen et al. (2019) estimated marine flashes to be twice as productive than those over land. Clearly, further measurements and process understanding are needed to reconcile differences in LNO x production.
Details about global LNO x reported in previous studies are given in Section 3.7.1.
In our model, the calculated amount of LNO x at a grid point location at a given time step (60 min) is distributed evenly in the 5 vertical in log-pressure coordinate from 500 hPa to the cloud top for intra-cloud (IC) flashes, and from 500 hPa to surface for cloud-to-ground (CG) flashes (see Section 3.7.2). The method is motivated by the data analysis of Price and Rind (1993).
Their observations from 139 thunderstorms cover cold cloud thickness (i.e., the cloud top height minus the freezing level) values ranging between 5.5-15 km and freezing level values between 2.7-5 km. The ratio z R = IC/CG increases from 0 to 4.6 with cold cloud thickness from 5.5 to 15 km but remains relatively constant with freezing level. We take the level below 10 which the CG generated LNO x is distributed as the observed minimum freezing level plus half of the minimum cold cloud thickness, i.e. (2.7+5.5/2) ≈ 5.5 km. The selected 500 hPa level is closest to this 5.5 km value. Since the amount of NO produced per flash is taken to be the same for both IC and CG flashes, the partitioning of the total flash rate into the CG and IC flash rates only influences the shape of the vertical distribution, with the total LNO x released remaining independent of the partitioning. 15

Flash-rate parameterisations based on convective cloud-top height and LNOx
The approach of parameterising lightning flash frequency in terms of cloud-top height has its origins in the simple scaling relationships suggested by Vonnegut (1963) for the electrical power output of a thundercloud with the cloud size. Vonnegut's model assumes an electric dipole structure with two equal but opposite charge volumes, separated by a distance on the order of the vertical cloud dimension, and a cloud aspect ratio of approximately unity. The charge transport velocity 20 (charging current) that supports the dipole, flows in the vertical and is assumed to exist across the horizontal cross-section of the cloud. The electrical power generated is proportional to the fifth power of cloud dimension and the lightning flash rate is proportional to the electrical power generated. Williams (1985) extended this work using Vonnegut's (1963) model and available atmospheric observations over land to demonstrate from observations that (a) the convective velocity for clouds ranging in size from a few kilometres to 17 km is proportional to the size of the cloud and (b) cloud height and flash-rate 25 data from three US locations (New England, Florida and New Mexico) are in good agreement with the predicted slope of 5 from Vonnegut (1963). An important consequence of the last finding is that charge transport velocity must scale with cloud dimension, in the same way as convective velocity does. However, as Williams (1985) writes, this does not establish whether convective motion or falling precipitation is directly responsible for the charge separation. Boccippio (2002) presents a systematic physical account of Vonnegut's (1963) model and the subsequent work on this issue by Williams (1985), PR92 and others. Boccippio (2002), using the laws of electricity, derived a fundamental scaling relationship between thunderstorm electrical generator power and storm geometry, which provides a possible theoretical basis for linking lightning to thunderstorm dynamics, microphysics, and geometry. Boccippio (2002) takes Vonnegut's 5 conceptualized thunderstorm of a quasi-steady state electrical dipole, with the horizontal and vertical scales of the two dipole charge centres comparable and varying with storm scale. The storm generator current is conceptualized as a net charge transport velocity, which maintains the dipoles. The electrical generator power (watts) is calculated as the generator current multiplied by the potential drop between the dipoles. The generator power is further expressed by assuming tangential spherical charge centres with volumetric charge density (±ρ) and radius R, maintained by a generator current density (= 10 charge density × charge transport velocity). Scale similarity between horizontal and vertical dimensions is invoked, and the dipole separation (= 2R) is assumed to vary linearly with cloud-top height, a more readily observable parameter. This assumption is based on observations that in many storms the lower negative charge region remains relatively constant in height and that most upper positive charge is carried on small ice crystals with negligible terminal velocity. The lightning flash rate is taken to vary linearly with lightning power dissipation, the latter is assumed to vary monotonically with 15 generator power. A further significant and explicit simplification replaces charge transport velocity with storm updraft velocity w. With the assumptions above, the flash rate F scaling relation is

Boccippio's (2002) extension of Vonnegut's model and derivation of scaling relationships for land and ocean
where the coefficient ~ , ε is permittivity constant (coulomb 2 joule -1 m -1 ), and the coefficient γ has units of joule -1 . By assuming that that each flash is responsible for a constant electrostatic energy and that the variability in the charge density of the dominant dipole regions 20 is small, Boccippio (2002) treated as a single constant, irrespective of whether the flash is over land or ocean.
In all the following parameterisations, F L,O is in flashes per minute, H is in km, and w is in m s -1 .
It is apparent that in the above approach, flash rate is set by storm dimension, with a small direct contribution from generator current density (which is taken to vary linearly with updraft velocity). In the case where ( ) for land and ocean are considered to be power law fits of H, i.e. 25 , = , , , the above approach based on Eq. (5) leads to the following self-consistent scaling relationships, where , and , are coefficients for land and ocean: , = , / , ( / , ) .
These self-consistent scaling relationships are central to this study.
Using these scaling relationships, the data of Williams (1985), PR92 and LIS satellite data, Boccippio (2002) (note there is a negative sign missing in the first exponent of this equation in Boccippio (2002)), and = 4.09 × 10 . , where = 1.4314 × 10 was used in deriving the relationships (9) and (10), which are plotted in Figure 1.
There were no direct F vs. H data to fit a relationship similar to Eq. (11) for the marine environment. Therefore, PR92 made the assumption that the charge separation velocity and the convective updraft velocity are equal, noted that an increase in the 10 intensity of convective updrafts enhances cloud electrification, and derived empirical relations based on observations between maximum updraft velocity (w, m s -1 ) and H (km) for continental and marine clouds = 1.49 . , = 2.86 . .
Eliminating H from Eqs. (11) and (12) Now, PR92 made the crucial assumption that Eq. (14) is independent of location, so is valid for the marine environment too.
Thus, by equating Eq. (14) and Eq. (13) they obtained the following expression for marine environment 15 (Our calculation suggests that the coefficient 6.4 in the above equation should actually be 5.94.) The relationships (11) and (15) are plotted in Figure 1, which show that the oceanic flash frequencies that are roughly 2 to 3 orders of magnitude smaller than those obtained for continental clouds.
According to Boccippio (2002), the derivation of Eq. (15) includes significant formal inconsistency and yields non-physical cloud-height predictions upon inversion and other non-physical behaviours, e.g. inverse relationship between updrafts and 5 cloud heights.
The PR92 parameterisations (11) and (15) are widely used in global chemistry transport models and coupled chemistryclimate models.  Michalon et al. (1999) identified that the PR92 parameterisation, when used over water in global models, was producing 15 lightning flash frequencies that did not agree with observations. To address this, they proposed that cloud electrification is directly related to cloud droplet concentration (N) and droplet size, and derived = / , where A is a proportionality constant. They retained the continental PR92 expression (11) by considering that it has been directly calibrated by using observed F and H values, and used it for the ocean by multiplying it with a factor of ( / ) / assuming that cloud

Michalon et al.'s (1999) ocean parameterisation
where two "standard" continental and maritime cloud droplet concentrations of N L = 600 and N O = 50 per mg of air, respectively, were used. Eq. (16) is plotted in Figure 1. One factor for a smaller droplet concentration in marine clouds is suggested to be more intense droplet coalescence in such clouds (Rosenfeld and Lensky, 1998).
In the above approach, the values of N L and N O are not well constrained. There is a large variability in cloud droplet concentration over land and ocean, as reflected in values observed in field experiments as well as those prescribed in cloud 5 microphysics schemes in global models (convective cloud droplet concentrations are not usually predicted in global climate models) (Rosenfeld and Lensky, 1998;Gultepe and Isaac, 2004). Boccippio (2002) cautions against this approach, and we suggest that given the uncertainty in the mean droplet concentrations, the approach of Michalon et al. (1999) is essentially empirical.

This study: Alternative flash-rate parameterisations 10
This study includes (a) a reanalysis of the Williams (1985) and PR92 data for lightning flashes versus cloud top height over land into the self-consistent scaling relationships framework of Boccippio (2002) and (b) the derivation of a new relationship for the oceanic environment using these scaling relationships.
For land, considering initially the relationship of updraft velocity with cloud-top height: equating the scaling relationship Eq.
At this stage is undetermined. We proceed by fitting Eq. (17) directly to the F vs. H data for land reported by Williams (1985) and compiled by PR92 (which are reproduced in Figure 1 as blue solid circles). This gives = 1.612 × 10 and = 2.40 × 10 . .
The relationships PR92 Eq. (11), Boccippio (2002) Eq. (9) and this study (TS) Eq. (18) are presented in Figure 1 along with the F vs. H data just discussed. Although these behaviours look almost identical and are well within the scatter of the data, 20 Eq. (18) shows a slightly higher flash rate for higher H than the other two. The almost fifth power dependence on H makes F L very susceptible to even small changes in the F L -H relationship and in how cloud-top height is calculated in the model.
Note that variation of flash rate at higher H is important for both global distributions of flash rate and how it may change with changing convective activity, e.g. climate change.
environments. Figure 1 presents Ushio et al.'s (2001) (Us01) data over the ocean in the tropics and extratropics (their Figure   3b and 3d), which they obtained by averaging flash rates every 1 km in cloud height. These are shown as triangles in Figure   1. Also shown are the Molinié and Pontikis' (1995) (MP95) flash-rate data over French Guyana coastal zone which we averaged over every 1 km in cloud height for heights greater than 10 km (below this height, the number of data points is not 5 sufficient for averaging). Because these are coastal observations, it is possible that the air masses in which the thunderstorm clouds developed could be mixed air masses (i.e. both continental and marine). However, these data do show a qualitative agreement with the Ushio et al. (2001) data in Figure 1.
Applying the scaling relationships (Boccippio 2002) to obtain an equivalent relationship for marine clouds involves equating Eq. (6) with the observed oceanic convective updraft velocity versus cloud top height from PR92 given in Eq. (13). This 10 yields = 2.86 and = 0.38. Substituting these into the scaling relationship Eq. (7) gives For this study, fitting Eq. (19) to the ocean data leads to = 2 × 10 (which gives = 0.7 × 10 ) and This marine parameterisation yields flash rates that are approximately an order of magnitude smaller than the PR92 continental formula and roughly two orders of magnitude larger than the PR92 marine parameterisation. In Figure 1 are shown the relationships for lightning flash rates and associated cloud-top height over the ocean for PR92, Boccippio's (2002)  15 Eq. (10), Michalon et al. (1999) and this study (TS) Eq. (20). Clearly the PR92 ocean equation is unrealistic, and the relationship of Boccippio (2002) gives marine flash rates that are twice as large as Eq. (20) and are not supported by the data plotted. The relationships of Michalon et al (1999) and this study group together around the oceanic data.
The values of implicit in Eqs. (18) and (20) are different for land and ocean, although as per Boccippio's (2002) theory, they should be the same (as used in deriving Eqs. (9) and (10)). This implies that either or both assumptions in the theory 20 that each flash is responsible for a constant electrostatic energy and that the variability in ρ is small are not true, possibly due to differences in cloud microphysics between land and ocean. If we use the same logic concerning aerosol microphysics as Michalon et al. (1999) used in deriving Eq. (16), then the different values of in Eqs. (18) and (20) for land and ocean can be interpreted in terms of and being different with < .
• Run 5 (Bo02): Boccippio's (2002) formulae Eqs. (9) and (10). 5 Note that Runs 2 and 5 only differ in the values of the linear coefficients used in the flash-rate parameterisations. We also did a model run with no LNO x emissions for a broad comparison, where appropriate, with the modelled changes in tropospheric composition resulting from the changes in the flash-rate parameterisations.
ACCESS-UKCA was setup as a free running simulation for 2 years (2005)(2006) for each of the above runs, and the simulation was started using model initial conditions taken from a previously spun-up, nudged model run that used a 10 Newtonian relaxation nudging (Uhe and Thatcher, 2015) within model levels 20-45 (between altitudes ∼ 3 km to 14 km) and the default lightning scheme. The variables nudged were the horizontal wind components and potential temperature by using ECMWF's ERA-Interim reanalyses (Dee et al., 2011) on pressure levels. The idea was to start the simulation with meteorological/transport errors minimised in the free troposphere to the extent possible. The first year of the free running simulation was used as a spin-up period and the model output for the year 2006 used for analysis reported below. (We also 15 did Runs 1 and 2 with nudging for the years 2005-2006 with the same initial conditions as for the free running simulations, and the results are summarised in Section 5).

Comparison of the modelled lightning flash rates with satellite observations
We analyse the global gridded lightning flash data from the Optical Transient Detector (OTD) on the OrbView-1 (formerly MicroLab-1) satellite and the Lightning Imaging Sensor (LIS) on the Tropical Rainfall Measuring Mission (TRMM) 20 satellite, which are described by Cecil et al. (2014) and available from https://lightning.nsstc.nasa.gov/data/data_lis-otdclimatology.html (V2.3.2015). The high resolution (0.5° × 0.5°) mean annual flash climatology (HRFC), low resolution (2.5° × 2.5°) mean annual flash climatology (LRFC), and low resolution (2.5° × 2.5°) monthly time series (LRMTS) data products are useful for our analysis. The climatology and time series data are flash density values with the units of flashes km -2 yr -1 and flashes km -2 day -1 , respectively. The OTD data available from July 1995 to January 2000 cover all 25 latitudes, whereas the LIS data available from February 2000 to February 2014 cover ± 42.5°. The climatology data cover all latitudes. For comparison with the modelled flash parameters, the satellite data were spatially regridded to the model N96 resolution (1.875° longitude × 1.25° latitude) using the Climate Data Operators (CDO) software. Table 1 gives the observed and modelled lightning flash frequencies (flashes s -1 ) averaged over the globe, Northern Hemisphere (NH), Southern Hemisphere (SH), land, and ocean for the five model runs, and these are plotted in Figure 2. In 30 Figure 2a, the observations are the combined LIS and OTD climatological data, whereas in Figure 2b, the observations are the LIS data for the year 2006 which are only available for the latitudinal range ± 42.5º with the modelled values also given for this range. The observations for the year 2006 are very well correlated to the climatology, and the two are very similar in magnitude as well since most lightning activity would fall within ± 42.5º. The data show that ∼ 80% of the global lightning flashes are over land, and ∼ 55% are in the Northern Hemisphere. On average, the default PR92 parameterisations (Run 1) underestimate the flash frequency data by 28% for the globe, by ∼ 13% for land, and by ∼ 96% for the ocean. Clearly, the 5 oceanic PR92 flash-rate parameterisation Eq. (15) does not work well at all over the ocean, yielding an almost zero flash frequency compared to the data. In contrast, the new flash-rate parameterisations (Run 2, Eqs. (18) and (20)) greatly improve the estimation of flash frequency over the ocean, with some overestimation (by ∼ 15%) of the climatology and giving nearly the same value as the year 2006 observational data. There is an improvement over land too compared to the PR92 formula.
Globally, the new parameterisations yield almost the same flash frequency as the data. Both the PR92 and the new 10 parameterisations lead to an almost equal partitioning of flashes in the Northern and Southern hemispheres, compared to ∼ 55% in the Northern Hemisphere indicated by the data. As expected, the Run 3 flash rate is nearly the same as that for Run 1 for land, and that for Run 2 for the ocean. For Run 4, with the PR92 continental Eq. (11) and Michalon et al.'s (1999) marine Eq. (16), the flash rate over land is almost the same as Run 1, whereas for the ocean Run 4 gives a flash rate that is about 50% higher than the climatology and 23% higher than the data for 2006. Over the ocean, Run 5 (Bo02) predicts flash rate 15 nearly twice as large as that observed and that from Run 2, which leads to an overestimation of the observed global total.
Over land the Run 5 estimate is about 10% lower than Run 2.

Modelled flashes per second
Modelled flashes per second The monthly variation of the observed and modelled flash rates is shown in Figure 3. In the Northern Hemisphere (Figure 3a), the model simulates the observed variation qualitatively with a peak in July, but while Run 1 underestimates the observed flash rate for all months, the other Runs mostly underestimate during February -10 May and do well for the other months. In the Southern Hemisphere (Figure 3b), again the model is able to simulate the observed variation well qualitatively, but a significant overprediction for January -April and an underprediction for August -October is apparent. The underprediction in spring in the Northern Hemisphere and overprediction in autumn in the Southern Hemisphere could be due to a displacement of lightning activity across the equator. The underprediction in spring in the Southern Hemisphere appears to be due to model deficiency over land (Figure 3d). In the global plot (Figure 3c), while Run 1 always underestimates, which is mostly because of its underprediction of the flash rate over the ocean (Figure   3e), the other Runs underestimate in the first three months of the year and overestimates during September -November, and these differences can be explained in terms of the hemispheric differences shown above. The nature of monthly variation in Figure 3d for land-based flash rates is very similar to Figure 3c, indicating that continental flashes dominate the global total. 5 The large underprediction by the PR92 scheme in Run 1 and an overprediction in Run 4 over the ocean can be seen in Figure   3e; there are also differences in the monthly variation. Run 5 overestimates the observed oceanic variation the most out of all Runs, and this leads to an overprediction of the observed global variation except for September and October.  Table 2 for the globe, land and ocean. Also given are the values of fractional bias FB = 2( − )/( + ) , which varies between -2 (overestimation) and +2 (underestimation), and the (Pearson) correlation coefficient (r). Considering these performance statistics together suggests 5 that the Run 2 flash-rate parameterisations from this study yield the best comparison with the data.  based on convective cloud-top height or uncertainty/biases in the modelled convection (e.g. Allen and Pickering, 2002;Tost et al., 2007). Another potential factor could be greater vertical wind shear outside the tropics which extends the horizontal lightning channel length (Huntrieser et al., 2008), which is not accounted for in the cloud-top height-based approaches. The LIS/OTD observations have some limitations too, such as a short sampling duration (just minutes) for a particular global location and lightning detection efficiencies not being perfect (Clark et al., 2017). Hereafter, we only present plots from Run 1 (default) and Run 2 (new), but the results from Runs 3-5 are included in all comparison Tables except Table 5.  Figure 5b, showing high flash density over land in the tropics and subtropics, and also lower mid latitudes. There is also some significant flash density over the ocean at these latitudes, particularly over the Pacific, western 15 Atlantic, western Indian Ocean near southern Africa, and the seas around the maritime continent (i.e., largely Indonesia, the Philippines and Papua New Guinea). The distribution modelled using ACCESS-UKCA with the default PR92 flash-rate scheme (Run 1) shown in Figure 5c is very similar to other global modelling studies that use the same PR92 scheme, e.g. Allen and Pickering (2002)  Hemisphere. Over the ocean, in contrast to the observations, the PR92 scheme predicts almost zero flash density. However, as shown in Figure 5d, ACCESS-UKCA with the new flash-rate parameterisations (Run 2) simulates the oceanic distribution of flash density much better than the PR92 scheme, although it is clear that there are some significant spatial differences (e.g. low bias over western Indian Ocean near southern Africa, and high bias over equatorial Indian Ocean and the Pacific) compared to the corresponding observations and climatology. The modelled flash-density distributions over land in Figure  5 5c and Figure 5d are nearly the same.

LIS satellite data for the year 2006 (available only for ± 42.5° latitudes), (c) model simulation with the default PR92 flash-rate parameterisations (Run 1), and (d) model simulation with the new flash-rate parameterisations from this study (Run 2).
The area-weighted NMSE, FB and correlation coefficient (r) comparing the spatial patterns of the observed climatology presented in Figure 5a with the annually-averaged modelled patterns of flash rate are given in Table 3 for the five model  15 runs. The NMSE and FB values clearly show that Run 1 performs the worst, a result dominated by the oceanic component.
While the NMSE values are nearly the same for Runs 2-4, Run 2 has the best FB values for both land and ocean. While Run 5 has the best FB value for the globe, this is fortuitous because the model underestimation and overestimation for land and ocean, respectively, counteract in the calculation of global FB. Thus, these statistics should be considered separately for land and ocean in examining the model performance. The spatial pattern correlation stays essentially the same for all Runs 20 (presumably because the underlying independent model variable, the cloud-top height, is the same in all model runs), and it is lower for the ocean -suggesting that further understanding of convection and lightning processes and their parameterisations are needed. The global correlations in Table 3 are very similar to those reported by Gordillo-Vazquez et al.
(2019) for cloud-height based schemes, but those for the ocean are lower in our study.
Based on the above flash-rate comparisons, Run 2 (TS1) performs the best, followed by Run 3 (TS2), Run 4 (Mi99) and Run 1 (PR92). 5 Modelled flash rates depend critically on modelled convection parameters (e.g. the cloud-top height) used by flash-rate parameterisations and on the representativeness of these parameterisations themselves. Thus, it is common in a global model to match the globally averaged modelled lightning flash rate to the observed value, e.g. that based on the LIS/OTD climatology (∼ 46 flashes s -1 ), by applying a constant scaling factor to the modelled global flash-rate spatial distribution (e.g., Tost et al., 2007;Finney et al., 2016;Clark et al., 2017;Gordillo-Vázquez et al., 2019), where the scaling factor is the ratio 10 of the observed global average flash rate to the modelled global average flash rate and is calculated by doing a model prerun. However, such a scaling would be misleading when there are large differences in the spatial representativeness of the flash rate computed by the parameterisation used in a model. For example, scaling the PR92 derived global flash-rate distribution would over-adjust the flash rate (and hence LNO x ) over land to compensate for the deficiency in the oceanic parameterisation. (Scaling can also be also applied to tune the amount of NO produced per flash to get a desired total global 15 LNO x amount, as per Eq. (1).). In the present study, no scaling factor was applied to the modelled flash rate, nor was it necessary. Figure 5a.

Global LNOx
The modelled global lightning-generated NO x using the various lightning flash-rate schemes are presented in Table 4. With the new flash-rate parameterisations (Run 2), the modelled global LNO x increases to 6.6 Tg N yr -1 from 4.8 Tg N yr -1 for Run 1, an increase of 38%, most of which is due to the change in the oceanic flash-rate component. Of the total global LNO x , 5 about 20% is generated over the ocean in Run 2, compared to ∼ 1% in Run 1. The partitioning into NH and SH is almost equal for both schemes. The Run 3 and Run 4 total LNO x emissions are similar to the Run 2 value, whereas the Run 5 value is 14% greater than that for Run 2 due to the higher value of the oceanic component. Given the same value of NO emitted per flash used for both IC and CG flashes, the partitioning of the global LNO x into NH, SH, and Land and Ocean for all Runs in Table 4 is very similar to the partitioning of flash rate for the corresponding Runs in Table 1. 10 The amount of global LNO x produced, (Tg N yr -1 ), is a function of the global average flash rate, (flash s -1 ), and the 15 moles of NO produced per flash, : If the climatological average = 46.5 flash s -1 based on the LIS/OTD satellite data is used, then = 20.5 × 10 .
If the NO production per flash differs for IC and CG flashes then can be taken as a weighted average over mean IC and CG flash fractions. The values in Table 4 are consistent with Eq. (21) when = 330 moles NO per flash as used in ACCESS-UKCA and the modelled from Verification studies include: Using a CTM representing the LIS/OTD flash data, Martin et al. (2007)  moles NO per flash, with a global LNO x emission of 5.9 ± 1.7 Tg N yr -1 . Boersma et al. (2005) analysed above-cloud tropospheric NO 2 column retrievals from the GOME satellite observations for the year 1997 for cloudy scenes over tropical oceans and continents, and found that the above-cloud annual-mean NO 2 column increases sharply with convective cloud-top height (H)-as . for continents and . for oceans, where H > 6.5 km. Considering that these above-cloud NO 2 columns primarily consist of contributions from lightning-generated NO x , 20 which is a direct function of flash rate, there is a very good agreement between these power-law exponents and those in the flash-rate relationships (11) and (18)  The modelled mean global distributions of LNO x from the two runs presented in Figure 6a and b are essentially in proportion to the flash density distributions given in Figure 5c and Figure 5d, respectively. The new flash-rate scheme (Run 2) leads to a larger and broader distribution of LNO x over the ocean compared to the PR92 scheme, while over land they are very similar.
In the absence of any direct measurements of global spatial distribution of LNO x for comparison we present in Figure 6c the annual LNO x distribution obtained by Miyazaki et al. (2014) using an assimilation of satellite measurements of atmospheric 5 composition and the LIS/OTD lightning flash data into a global CTM for the year 2007. This plot is a reproduction of their Figure 6 (middle-left plot) based the data 3 supplied by K. Miyazaki (personal communication, 2020) at a horizontal resolution of 2.8° × 2.8°. Over the ocean, the new flash-rate scheme (Figure 6b) agrees much better with the assimilated field than does the PR92 scheme, but is clear that the oceanic LNO x distribution in the plot with assimilation is broader, more diluted in the tropics, and even extends to high latitudes which is not seen in Figure 6b nor indicated by the observed flash-10 rate distributions in Figure 5a and Figure 5b (this could be due to limitations of the data assimilation used). Over land, the LNO x distributions predicted by both PR92 and the new scheme are similar and broadly agree with Figure 6c at low latitudes (except parts of India), but do properly not describe the extension of LNO x into the temperate latitudes, particularly in the Northern Hemisphere. Figure 6c yields a total LNO x of 6.36, 3.67, 2.69, 5.58 and 0.78 Tg N yr -1 for the globe, NH, SH, land and ocean, respectively, which except for SH are closer to the Run 2 values than to the Run 1 values in Table 4. Direct and 15 more extensive measurements would be necessary for a better evaluation of the predicted LNO x distribution.

Vertical distribution
The vertical distribution of LNO X in the model at a grid point location is a parameter (see Section 2.2) that is essentially unconstrained by observations. Figure 7  subtropical ones show a peak at 6.5 km. There are no direct measurements to verify any of the LNO x profiles and we believe further work is needed to constrain them.

Modelled tropospheric total column NO2 and validation
As lightning impacts atmospheric NO x directly, any changes in the modelled total tropospheric NO x can be examined and 10 compared with available observations. The modelled variations of the zonal averaged tropospheric total NO 2 column ( , ) for the globe presented in Figure 8 show (https://www.temis.nl/airpollution/no2.php; Boersma et al., 2018). While the model simulations qualitatively agree with the satellite observations in Figure 8, it is apparent that except for the latitudes 30°S-60°S the modelled values are generally 20 higher than the observations. Notwithstanding any model shortcomings in predicting the global NO 2 distribution, there are limitations of the OMI satellite data used. Firstly, there is limited sensitivity of the OMI sensor to NO 2 below the cloud level, where most NO 2 is situated, and thus cloudy tropospheric NO 2 retrievals (with cloud radiance fraction > 0.5 or cloud fraction > 0.2) cannot be interpreted as valid down to the Earth's surface (Boersma et al., 2017). Secondly, given that the OMI satellite data are representative of overpass time and thus comparing them with the mean model fields averaged over full diurnal periods introduces uncertainty.
As an alternate to the OMI data, we use data from the global reanalysis of atmospheric composition produced by the https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-reanalysis-eac4-monthly) for comparison with the 5 monthly-averaged modelled NO 2 . The reanalysis was produced by assimilating space observations of aerosols and reactive gases using a 4D-Var method in an ECMWF global atmospheric model with 60 pressure levels (from 1000 to 1 hPa) and a horizontal resolution of 0.75° × 0.75°. For NO 2 , the model assimilated the tropospheric column retrievals from the SCIAMACHY, OMI, and GOME-2 satellite overpasses at ∼ 10:00 LT, 13:30 LT and 09:30 LT, respectively (for the year 2006, only SCIAMACHY and OMI data were available for assimilation). Monthly-averaged total vertical column NO 2 10 reanalysis data (version -ECMWF Atmospheric Composition Reanalysis 4) are available and used here.
We obtain the tropospheric NO 2 column ( , ) from the CAMS total NO 2 column ( ) as follows: where , is the CAMS total NO 2 column over the Pacific (180°W) and , , is the tropospheric NO 2 column over 180°W. This is one approach used with satellite data to separate the tropospheric and stratospheric amounts (Inness et al., 2019), which assumes a longitudinal homogeneity of the stratospheric NO 2 column amounts and a constant and negligibly 15 small , , (Lauer et al., 2002). , , is not available from the CAMS data, but this obtained from the OMI data discussed above is shown in Figure 8  contributions from sources such as lightning, aviation, shipping, and possibly regional transport over the ocean. The , , values are also greater and possibly more uncertain than the differences between the two model simulations. 20  Since , , is not available from the CAMS data, we take the OMI derived , , shown in Figure 8 and use it in Eq. (23) to determine , from the CAMS data. (While these OMI NO 2 column data at 180°W have the same limitations as the global OMI data mentioned earlier, we expect that the impact of these limitations would be much smaller at 180°W with almost no surface source contributions compared to terrestrial locations.) The resulting CAMS zonal annual-mean 10 , as a function of latitude for the globe is presented in Figure 9a, together with the corresponding modelled variations.
The modelled variations agree well with the CAMS data, better than with the direct OMI data in Figure 8. A comparison with the model variation without LNO x in Figure 9a suggests that the modelled increase in NO x due to lightning is largely confined to ±35°.
In Figure 9b for the ocean, the agreement with the data is again good, except for considerable model underestimation for 15 latitudes greater than 35° which, presumably not related to LNO x , may be due to factors such as possible underestimation of the calculated , , and overestimation of shipping emissions of NO x for these latitudes. Within ±35°, the new parameterisation yields 22% larger NO 2 column values compared to the default model setup, whereas as expected the two model curves are virtually the same for the other latitudes. The CAMS NO 2 columns are somewhat better simulated by the model with the new oceanic flash rate than with the PR92 parameterisation, particularly in the northern tropics. 20 For land (Figure 9c), the two model curves are nearly identical and overall compare very well with the CAMS reanalysis data. An overestimation of the reanalysis data in the southern tropical region is evident.
The statistics in Table 5 show that PR92 and TS1 have similar agreement with CAMS, except for a noticeably better FB over the ocean with TS1. We expect no significant changes over land because the LNO x is almost unchanged in the two runs.
Clearly, the comparison also depends on the selected value of NO produced per flash. We have used the model default value of 330 moles NO per flash. However, if we were to match the average CAMS column value in Table 5, the new parameterisation with 310 moles NO per flash, the value suggested by Miyazaki et al. (2014), would probably yield a 5 somewhat better prediction. Obviously, there are other sources of uncertainties in the NO 2 comparison, such as those associated with the CAMS reanalysis and the assimilated satellite columns, which are documented in the appropriate references cited above, as well as those to do with model inputs (e.g. emissions) and processes.   Figure 9 for latitudes between ± 30º. We present the impact of LNO x determined from the flash-rate parameterisations from Runs 1 and 2 on tropospheric composition, namely total NO x, O 3 , OH and CO.

Oxides of nitrogen (NOx)
The modelled tropospheric NO 2 columns and their comparison with observations have already been presented in Section 5 3.7.3. Figure 10 presents the zonal distribution of total NO x (as NO 2 ) from the two model simulations and the difference between the two. In the lower troposphere, the two modelled distributions of the zonal annual-mean NO x (as NO 2 ) are virtually identical, with highest levels predicted within latitudes 20-60º N. These levels are governed by surface emissions of NO x which are the same in both simulations. The secondary concentration maximum at ∼15 km is due to the lightninggenerated NO x . The new lightning parameterisations cause an increase in the mid-to upper-tropospheric NO x (Figure 10c), 10 particularly within the tropics and subtropics, and this increase is by as much as 40 ppt (by volume) in the northern tropics.
There are some localised decreases in concentration in the lower troposphere over the Northern Hemisphere.
The volume-weighted global tropospheric NO x obtained from the PR92 scheme is 55.1 ppt, and it is 35.2 ppt over the ocean and 94.0 ppt over land. With the new flash-rate scheme, these values increase by 8.7 ppt (15.7%), 9.9 ppt (28.0%) and 6.3 ppt (6.7%), respectively, and can be compared with the values 36.9, 20.7 and 68.5 ppt, respectively, obtained from the model 15 simulation with no LNO x emissions. To some extent, the modelled tropospheric averages also depend on how the tropopause is defined.

Ozone (O3)
Tropospheric ozone is a by-product of the oxidation of carbon monoxide (CO), CH 4 and other volatile organic compounds in the presence of NO x and is thus impacted by LNO x .
With the new flash-rate parameterisations (Run 2), the modelled tropospheric O 3 burden increases from 284 to 308 Tg O 3 , a rise of 8.5% over the default PR92 scheme (Run 1) (cf. 219 Tg O 3 with no LNO x emissions in the model). The new burden is 5 closer to a the ACCMIP (Atmospheric Chemistry and Climate Model Intercomparison Project) multi-model mean of 337 ± 23 Tg O 3 reported by Young et al. (2013), the latter value is consistent with measurement climatologies (this, however, does not necessarily mean that LNO x in these models is represented correctly). The Run 3 and Run 4 ozone burdens are 306 and 308 Tg, respectively.
The mean relative difference (%) between the global ozone mixing ratios predicted using the new lightning flash-rate 10 parameterisations and the default PR92 parameterisations is shown in Figure 11. Near the surface (Figure 11a), there are significant increases in ozone over the tropical oceans, especially in the Pacific and western Indian Ocean, and in most of the Southern Hemisphere (roughly by 8% on average). Over land, there are regions (e.g. south-eastern U.S.A. and northern Australia) where ozone has increased, and there are a few regions in the mid to high latitudes in the Northern Hemisphere where ozone has decreased very slightly. Tropospheric ozone chemistry is complex, but broadly speaking the O 3 increases in 15 the Southern Hemisphere are influenced by low ambient NO x concentrations where the O 3 production increases with NO concentration. O 3 is produced through photodissociation of NO 2 which is produced through oxidation of NO by HO 2 and RO 2 radicals (e.g., NO + HO 2 → NO 2 + OH). In the Northern Hemisphere, the increase in O 3 is less beyond the tropic, partly because the smaller oceanic area results in a smaller increase in LNO x through the use of the new oceanic flash-rate parameterisation. Carpenter et al. (1997) suggest that the tropospheric production potential of the Southern Hemisphere is 20 more responsive to the availability of NO than that of the (more polluted) Northern Hemisphere.
At an altitude of 6400 m (∼ 450 hPa) (Figure 11b), there are even bigger increases in global ozone using the new flash-rate parameterisations, particularly in the tropics, by as much as 25%, and in the Southern Hemisphere. This is because most LNO x emissions occur in the middle to upper tropical troposphere where the photochemical production of ozone is most efficient. 25

5
In Figure 12, we compare the modelled monthly averaged ozone with ground-based in-situ observations from the World Data Centre for Reactive Gases (GAW-WDCRG, http://ebas.nilu.no; https://www.gaw-wdcrg.org;) for the year 2006 at five stations: Ushuaia (54.85ºS,68.31ºW),Cape Grim (40.68ºS,144.69ºE),Mauna Loa (19.54ºN,155.58ºW),Minamitorishima (24.29ºN,153.98ºE),and Mace Head (53.33ºN,9.90ºW). Apart from data availability and covering a range of latitudes, the site selection was based on these sites being either oceanic or coastal so that the relatively large difference between the PR92 10 oceanic flash rate parameterisation and the new one could be examined against the observations. The hourly data were averaged to monthly values, and only those observational months were considered for which there were more than 75% valid hourly data points. Mauna Loa is located at an elevation of 3397 m on an island which is smaller in size than the grid resolution of the model and therefore it is difficult to correspond the sampling height to a particular vertical model level. We likely due to the model resolution issues). There are small, but noticeable, differences in the modelled ozone from the two simulations. The relative change in the modelled yearly averaged O 3 at these ground-based sites with the use of the new lightning parameterisation is small, at 5.9%, 1.3%, -1.9%, 5.9% and 0%, respectively. There is some improvement in the modelled seasonal variation at Ushuaia, Cape Grim and Minamitorishima with the new LNO x scheme, but for the other two 20 sites the model-data differences are much larger than those due to the LNO x changes. Generally, factors such as model's transport and chemical mechanisms, and input precursor emissions and their distributions are probably more influential in governing ozone model-data differences than LNO x near the Earth's surface. There is no clear indication if the differences in ozone between the two models are larger in the winter or summer, except for Ushuaia where the differences are larger during winter to spring. 25 Ozone (ppb) Ozone (ppb) Ozone (ppb) binned in the model vertical levels and averaged. The number of observed profiles at a location typically varied between 1 to 5 per month. At Irene, Nairobi and La Reunion Island, there were no data for 3-4 months. The observed profiles were averaged over the year. The modelled monthly-averaged profiles were averaged over the months for which profile observations were available. Clearly, the number of observed profiles is not sufficient for averaging over the year, and, therefore, the model-data comparison is essentially qualitative. 5 The observed O 3 profiles are simulated better by the new lightning flash-rate parameterisations (Run 2) than by the default PR92 parameterisations (Run 1) at Costa Rica ( Figure 13c) and for all the southern hemispheric sites (Figure 13d-Figure   13h). At northern latitudes beyond Costa Rica, the PR92 scheme performs better at Hilo and Paramaribo (Figure 13a and b).
The model with the PR92 scheme has a general tendency to overestimate ozone at most levels within 15°N-50°N, which is probably due to reasons not related to LNO x , and the use of the new lightning flash-rate parameterisations worsens it (see 10 Figure 15). The model describes the general shape of the observed profiles reasonably well. The model profiles with no LNO x included highlights the importance of LNO x on tropospheric ozone, particularly in the mid to upper troposphere. The modelled zonal annual-mean tropospheric O 3 from the two Runs and the difference between the two are presented in Figure 14. In the lower troposphere, the modelled ozone is smaller over the Southern Hemisphere than over the Northern 20 Hemisphere (Figure 14a and Figure 14b). The new flash-rate parameterisations result in O 3 increases everywhere ( Figure   14c). Closer to the surface, the increase is approximately 2 ppb (by volume) in the Southern Hemisphere and 0.5-2 ppb in the Northern Hemisphere. The largest increases are nearly 8 ppb in the tropics at altitudes ∼ 9 km. There is an agreement, both in magnitude and distribution, between the modelled zonal ozone in Figure 14a and b, and the data plotted in Figure 15a. The model reproduces the observed lower levels of tropospheric ozone in the Southern 15 Hemisphere. The observed high levels just below the tropopause within 10-40ºN are somewhat better reproduced by the new flash-rate parameterisations. Figure 15b and c represent the relative differences (=[ − )/ ] × 100%) between the annual averaged modelled and observed ozone for the PR92 scheme and the new parameterisations, respectively. On average, the model underestimation in the Southern Hemisphere has reduced with the new parameterisations, but there are areas such as that within 10-20 ºS below 5 km where there is some overprediction. There is a clear improvement in the predicted ozone 20 between 10ºS -10ºN throughout the troposphere using the new scheme while that in the northern high latitudes below 6 km remains unaffected. In the Northern Hemisphere, the new scheme tends to overestimate ozone within 10-50 ºN below ∼ 9 km; this is where there was already some overprediction by the PR92 scheme. In this region, even when LNO x is not included in the model (Figure 15d), the model-data differences are either small or there is some overestimation (near the surface) of the observed ozone. This suggests that there are likely to be factors other than LNO x responsible for the predicted 25 overestimation within 10-50 ºN. Additional factors that influence tropospheric ozone distribution in the model include dynamics, including inter-hemispheric mixing and stratosphere-to-troposphere exchange, precursor emissions, and how chemical mechanisms are represented. Considering the above, we can say that the new flash-rate scheme leads the modelled tropospheric ozone in the right direction, which is also supported by the fact that it causes the tropospheric ozone burden to improve, as stated above. Getting ozone in the upper troposphere correct is climatically important as surface temperature is more sensitive to changes in ozone in the upper troposphere and near the tropopause than those in the lower atmosphere (Forster and Shine, 1997) and, similarly, radiative forcing due to tropospheric ozone is more sensitive to ozone abundance in 5 the upper troposphere (Worden et al., 2011).
The volume-weighted tropospheric O 3 obtained from the PR92 scheme is 51.5 ppb over the globe, 48.8 ppb over the ocean and 56.6 ppb over land. With the new flash-rate scheme, these values increase by 4.1 ppb (8.0%), 4.4 ppb (9.1%) and 3.6 ppb (6.3%), respectively, and can be compared with the values 38.7, 36.2 and 43.4 ppb, respectively, obtained from the model simulation with no LNO x emissions. 10

Hydroxyl (OH) radical
The hydroxyl radical (OH) is the dominant oxidizing (and cleansing) agent in the global troposphere and controls the 20 atmospheric abundance and chemical lifetime of most natural and anthropogenic gases, such as methane (CH 4 ). The tropospheric abundance of OH is determined by a complex series of chemical reactions involving tropospheric ozone, methane, carbon monoxide (CO), non-methane volatile organic compounds (NMVOCs), and NO x , and also the amount of solar radiation and humidity (Naik et al., 2013). Through these reactions, the amount of LNO x produced also impacts OH.
The modelled zonal total annual mean tropospheric OH in Figure 16 shows highest OH concentrations near the surface in the tropics, with values as high as (25 -30) × 10 5 molecules cm -3 at ∼ 20º N. The concentrations decrease with altitude, but there is a secondary maximum in the upper troposphere at around 13 km. There is an increase in the OH concentration using the new flash-rate parameterisations (Figure 16c), particularly in the upper troposphere in the tropics, by as much as 5 × 10 5 molecules cm -3 . 5  Hemisphere. The broad hemispheric differences in OH are qualitatively similar to those for O 3 . With an increase in NO due to the new flash-rate parameterisation, OH increases (e.g. via the recycling of HO 2 by reaction with NO). In highly polluted air, NO 2 can be an OH sink (Lelieveld et al., 2016). Of course, transport would also influence these patterns, both horizontally and vertically.

5
Overall, we find that, with the new flash-rate parameterisations, there is a 13% increase in the annual-average volumeweighted global tropospheric OH, from 10.6 × 10 5 to 12.0 × 10 5 molecules cm -3 . The increase over the ocean is by 1.6 × 10 5 (16.3%) and that over land by 0.9 × 10 5 molecules cm -3 (7.6%). For comparison, the respective values obtained from the model simulation with zero LNO x emissions are 7.6 ×, 7.3 × and 8.1 × 10 5 molecules cm -3 .
The global amount can be compared with the ACCMIP multi-model mean of 11.1 ± 1.6 × 10 5 molecules cm -3 derived by 10 Naik et al. (2013) for the year 2000. Recent observationally based values reported by Wolfe et al. (2019) for August 2016 are 12.6 ± 2.9 × 10 5 for the Northern Hemisphere and 8.1 ± 1.9 × 10 5 molecules cm -3 for the Southern Hemisphere, and these for February 2017 are 8.8 ± 2.1 × 10 5 and 11.4 ± 2.8 × 10 5 molecules cm -3 , respectively. These can be compared with the corresponding modelled values 16.9 ×, 7.3 ×, 8.0 × and 11.6 × 10 5 using the PR92 scheme, and 18.7 ×, 8.5 ×, 9.2 × and 13.5 × 10 5 molecules cm -3 using the new scheme. The LNO x -induced increase in OH due to the new scheme adds to the model 15 high bias in the OH burden in summer, whereas it reduces the magnitude of the bias in winter with the bias shifting from low to high. The model value in the Northern Hemisphere in August is considerably larger than the observation even with the PR92 scheme. It is known that the UKCA StratTrop configuration yields substantially larger OH in the Northern Tropics at low altitudes compared to observations and to the ACCMIP multi-model estimates (Archibald et al., 2020).
The surface methane concentrations were prescribed in the model and methane was allowed to undergo loss processes in the 20 rest of the atmosphere. With an overall increase in OH using the new flash-rate parameterisations in ACCESS-UKCA, the global annual mean lifetime of CH 4 against loss by tropospheric OH ( _ , defined as the division of the global total atmospheric CH 4 burden and the globally integrated CH 4 loss rate by reaction with tropospheric OH) decreases by 6.7%, from 7.5 to 7.0 years. This value without the LNO emissions in the model is 9.2 years. The modelled methane lifetime is lower than the multi-model mean 9.7 ± 1.5 years reported by Naik et al. (2013), which could be due to a higher tropospheric 25 burden of non-lightening related NO x in ACCESS-UKCA and/or a more intense photolysis.

Carbon monoxide (CO)
There is a decrease in the modelled total tropospheric carbon monoxide (CO) with the use of the new lightning parameterisations, as evident from the zonal annual-mean difference plot in Figure 18 (this CO reduction is coupled to the OH increase, via the reaction OH + CO → CO 2 + H). In the lower troposphere, the decrease is by approximately 4-6 ppb (∼ 7%) in the Southern Hemisphere and 2-4 ppb (∼ 3%) in the Northern Hemisphere. This reduction gets a little larger in the 5 mid to upper troposphere in the tropics. Overall, the reduction in the volume-weighted global CO is 4.5 ppb (5.6%). Over the ocean is 4.7 ppb (6.2%) and that over land 4.0 ppb (4.5%). The volume-weighted tropospheric CO obtained from the PR92 scheme is 80.3 ppb over the globe, 75.4 ppb over the ocean and 89.9 ppb over land. With the new flash-rate scheme, these values decrease by 4.5 ppb (5.6%), 4.7 ppb (6.2%) and 4.0 ppb (4.5%), respectively, and can be compared with the values 96.4, 91.8 and 105.4 ppb, respectively, obtained from the model simulation with no LNO x emissions. 10 Figure 18: The difference between the zonal total annual-mean tropospheric CO (ppb, by volume) modelled using the new lightning flash-rate parameterisations (Run 2) and the default PR92 parameterisations (Run 1).
In Figure 19, we compare the modelled monthly-averaged CO with surface flask observations from the same GAW-15 WDCRG sites as in Figure 12 (data from Minamitorishima were missing for 4 months, so are not presented). With the use of the new lightning parameterisation, the relative change in the modelled yearly averaged CO at Ushuaia, Cape Grim, Mauna Loa and Mace Head is -8.1%, -9.8%, -3.8%, and -0.3%, respectively. The modelled ground-level CO is affected only very marginally by the flash-rate modification compared to the magnitude of the model-data differences, except at Ushuaia and at Cape Grim during the austral summer. Clearly, as in the case of ground-level O 3 , the lightning changes alone do not 20 reconcile the differences between the modelled CO and observations.

Model simulations with meteorological nudging
As mentioned in Section 3.5, we also did Run 1 (PR92) and Run 2 (TS1) with meteorological nudging (which could impact convection indirectly in the model) with the same initial conditions as for the free running simulations. Broadly speaking, the results were not too different from the respective free running simulation results for the year 2006 reported above. The averaged global, land and ocean flash rates obtained from Run 1 with nudging were 32.87, 32.50 and 0.37 flashes s -1 , 10 respectively, which are very similar to 32.92, 32.56 and 0.36 flashes s -1 , respectively, obtained from the free running Run 1 in Table 1. These obtained from Run 2 with nudging were 46.88, 36.80 and 10.08 flashes s -1 , respectively, which are on average 5% higher than 44.96, 35.88 and 9.08 flashes s -1 , respectively, obtained from the free running Run 2 in Table 1 (so a scaling factor of 0.95 would make the nudged model flash rates approximately match the free-running model flash rates).
The total LNO x produced over the globe, land and ocean changed in the same proportions since the NO produced per flash 15 was the same.
The annual-averaged global spatial distributions of the modelled flash rate with nudging for both Runs were very similar to the respective free running model plots shown in Figure 5c and Figure 5d, with the overall spatial correlation with the LIS/OTD satellite climatology (Figure 5a) slightly improved from 0.72 to 0.75, and some improvement in the model performance over the Indian subcontinent but a slight deterioration over the southern US. Nudging would have an impact on tropospheric composition, which we have not presently investigated, but we estimate that with flash rate constrained the modelled tropospheric composition with and without nudging would, on average, be within ∼5%.

Conclusions
We have critically examined parameterisations of lightning flash rate that are based on the cloud-top height approach. 5 Testing of the widely used Price and Rind (1992) (PR92) parameterisations within the ACCESS-UKCA global chemistryclimate model for the year 2006 using the global LIS/OTD satellite data has revealed that while the parameterisation for land yields satisfactory predictions, with a globally averaged flash rate of 31.03 flashes s -1 compared to the observed 34.92 flashes s -1 , the oceanic parameterisation severely underestimates the observed flash rate, yielding on average 0.33 flashes s -1 compared to the observed 9.16 flashes s -1 . This leads to lightning-generated NO x (LNO x ) being underestimated 10 proportionally over the ocean and thus influencing tropospheric composition. Any interannual variability in lightning was not investigated.
Following Boccippio's (2002) scaling relationships between thunderstorm electrical generator power and storm geometry as the basis, we derived alternative flash-rate parameterisations. While the new parameterisation for land performed slightly better than the corresponding PR92 one, giving a globally averaged flash rate of 34.23 flashes s -1 compared to the observed 15 34.92 flashes s -1 , the new parameterisation for ocean performed more accurately, giving a globally averaged flash rate of 8.84 flashes s -1 compared to the observed 9.16 flashes s -1 . We also tested an oceanic parameterisation by Michalon et al. (1999), which gives a global oceanic average of 11.31 flashes s -1 . With the new parameterisations, there was an increase in global LNO x from 4.8 to 6.6 Tg N yr -1 , with the new estimate comparable to 6.3 ± 1.4 Tg N yr −1 obtained by Miyazaki et al. (2014) using an assimilation of multiple satellite datasets into a global CTM. There is a large uncertainty in the amount of 20 NO produced per flash in the scientific literature. The model's use of 330 moles NO produced per flash is close to the average value 310 moles NO per flash determined by Miyazaki et al. (2014) using data assimilation, but requires better constraining.
The use of the new flash-rate parameterisations in ACCESS-UKCA demonstrated a considerable impact on the modelled tropospheric composition compared to the defaults PR92 parameterisations, mainly due to the change in the oceanic flash-25 rate component. In particular, the following impacts were observed: • An increase in the mid-to upper-tropospheric NO x by as much as 40 ppt by volume (as NO 2 ) in the northern tropics. An overall increase in the global NO x by 8.7 ppt (15.7%) and by 9.9 ppt (28.0%) over the ocean. A better agreement of the modelled tropospheric NO 2 columns with the CAMS reanalysis data over the ocean.
• The tropospheric O 3 burden increased by 8.5%, from 284 to 308 Tg O 3 , closer to a multi-model estimate of 337 ± 30 modelled ozone in the troposphere improved somewhat compared to global observations in the Southern Hemisphere. There are considerable ozone biases in the model in the Northern Hemisphere that are not related to LNO x .
• A 13% increase in the annual-average volume-weighted global tropospheric OH, from 10.6 × 10 5 to 12.0 × 10 5 molecules cm -3 . 5 • A decrease in the global annual mean methane lifetime against loss by tropospheric OH by 6.7%.
• An overall reduction in the global CO by 4.5 ppb by volume (5.6%).
The approach of parameterising lightning flash rate in terms of convective cloud-top height works well given its simplicity, and continues to be useful in accounting for LNO x in global models, although there were some significant spatial distributional differences in the modelled flash-rate density compared to the satellite data. The approach is also very sensitive 10 due to an almost 5 th power dependence on cloud-top height. With increased computational power in the future, it may be possible to understand and represent global LNO x in a better process-based manner through a cloud-resolving modelling framework with an explicit prediction of the electrical activity in storms.
Recent global chemistry-climate modelling studies using flash-rate parameterisations based on the convective cloud-top height show an increase in LNO x emissions in a future warming climate (e.g., Banerjee et al. 2014;Clark et al., 2017;15 Iglesias-Suarez et al., 2018), primarily as a result of increases in the depth of convection. This is also the case with a CAPE and precipitation rate based parameterisation (Romps et al., 2014). Conversely, it is found that a flash-rate scheme based on convective mass flux (Clark et al., 2017) and that based on upward cloud ice flux (Finney et al., 2018) predict a global decrease in future lightning flash density (under the RCP8.5 scenario). The study of Finney et al. (2018) argues that the ice flux method performs better than the PR92 cloud-top height method, and this should be reanalysed using the scheme 20 proposed here which drastically alters the flash rate over the ocean. As well, there is an existing uncertainty as to which physical parameterisation approach best represents the reality and the feedbacks that are important for lightning under climate change. This will be difficult as parameters such as cloud ice content and/or updraught mass flux used in flux-based lightning schemes are poorly constrained by available observations.

Author contributions 5
AKL devised the study, performed the flash-rate parameterisation formulation and model runs, analysed model output and data, and wrote the paper, with contributions and comments from all co-authors. IEG assisted with the flash-rate parameterisations and advised on various components of the paper, MTW performed some of the early model-data comparison for atmospheric composition and advised on the model setup, and NLA advised on the model configuration and provided relevant technical details on UM-UKCA.