Assessing and improving cloud-height-based parameterisations of global lightning flash rate, and their impact on lightning-produced NO_x and tropospheric composition in a chemistry–climate model

Although lightning-generated oxides of nitrogen (LNO_x) account for only approximately 10 % of the global NO_x source, they have 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 LNO_x 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 conduct a critical assessment of flash-rate parameterisations that are based on cloud-top height and validate them within the Australian Community Climate and Earth System Simulator – United Kingdom Chemistry and Aerosol (ACCESS-UKCA) global chemistry–climate model using the Lightning Imaging Sensor and Optical Transient Detector (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⁻¹ compared to the observed 9.16 flashes s⁻¹ and leading to LNO_x being underestimated proportionally. We formulate new flash-rate parameterisations following Boccippio's (2002) scaling relationships between thunderstorm electrical generator power and storm geometry coupled with available data. The new parameterisation for land performs very similarly to the corresponding PR92 one, as would be expected, whereas the new oceanic parameterisation simulates the flash-rate observations much more accurately, giving a global average over the ocean of 8.84 flashes s⁻¹. The use of the improved flash-rate parameterisations in ACCESS-UKCA changes the modelled tropospheric composition – global LNO_x increases from 4.8 to 6.6 Tg N yr⁻¹; the ozone (O₃) burden increases by 8.5 %; there is an increase in the mid- to upper-tropospheric NO_x by as much as 40 pptv, a 13 % increase in the global hydroxyl radical (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 NO_x and ozone in the Southern Hemisphere and over the ocean are improved by this new flash-rate parameterisation.

1 Introduction

Oxides of nitrogen (NO_x≡ NO (nitric oxide) + NO₂ (nitrogen dioxide)) play an important role in tropospheric chemistry by acting as a precursor to ozone (O₃) and the hydroxyl radical (OH), which are the principal tropospheric oxidants (Labrador et al., 2005). As a greenhouse gas, O₃ 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. OH plays a critical role in the chemical cycles of many trace gases, including methane (CH₄) 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 (Schumann and Huntrieser, 2007; 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₂) and oxygen (O₂) molecules into free N and O within the flash channel. These in turn react with ambient N₂ and O₂ to produce NO, which remains after the lightning channel cools. There is a conversion between NO and NO₂, during which ozone is generated in the presence of HO₂ and organic peroxy radicals, collectively called RO₂ (Bucsela et al., 2019). 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 (Schumann and Huntrieser, 2007).

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₃ and OH (Murray, 2016). For example, although LNO_x emissions are of similar magnitude to those from biomass burning or soils, their contribution to the total tropospheric ozone column is about 3 times larger (Dahlmann et al., 2011). This is because in the middle to upper troposphere where LNO_x is directly released, the O₃ production efficiency per unit NO_x is much higher (∼100 molecules O₃ per molecule NO_x) than that near the surface (∼10–30 molecules O₃ 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₃, OH, nitric acid (HNO₃), and peroxyacetyl nitrate (PAN) to the magnitude and vertical distribution of LNO_x. Modelling studies by Grewe (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 deep convection as a result of future climatic change have 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., 2014, 2018; Iglesias-Suarez et al., 2018). In addition, lightning intensity and distribution, and its uncertainty, in a future climate have implications 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 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

$$\begin{array}{}\text{(1)}& {\text{LNO}}_x=F\times P_{\mathrm{NO}},\end{array}$$

where F is the lightning flash rate and P_NO 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, 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 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 (Finney et al., 2014, 2016), product of CAPE and precipitation rate (Romps et al., 2014; Zhu et al., 2019), and multi-parameter 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 the Goddard Earth Observing System with chemistry (GEOS-Chem) (Hudman et al., 2007), Model for OZone And Related chemical Tracers (MOZART) (Emmons et al., 2020), Community Atmosphere Model with chemistry (CAM-Chem) (Lamarque et al., 2012), European Centre Hamburg Modular Earth Submodel System (ECHAM5/MESSy) (Jöckel et al., 2006), and the UK Met Office Unified Model coupled with the United Kingdom Chemistry and Aerosol global atmospheric composition model (UM-UKCA) (Archibald et al., 2020), 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 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 remain relatively constant compared to CCD values, meaning that it is largely the 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 parameterisations had considerable biases compared to satellite observations of spatial flash density but performed second best (behind their parameterisation based on ice flux) out of the five parameterisations tested. A recent comparison by Gordillo-Vázquez et al. (2019) of six flash-rate parameterisation schemes showed a relatively good performance by those based on cloud-top height. However, a number of global modelling studies have 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., 2014, 2016; Clark 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 flash-rate parameterisations following Boccippio's (2002) (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 (Australian Community Climate and Earth System Simulator – United Kingdom Chemistry and Aerosol; 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) 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 available as global climatological, interannual, and seasonal distributions. The veracity of the new modelled LNO_x estimates is examined through comparison of modelled and reanalysis tropospheric NO₂ 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.

2 ACCESS-UKCA global chemistry–climate model

We use the UKCA global atmospheric composition model (Abraham et al., 2012; http://www.ukca.ac.uk, last access: 4 May 2021) coupled to ACCESS (Bi et al., 2013; Woodhouse et al., 2015). In our simulations, ACCESS is essentially the same as the UK 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-only mode with prescribed monthly mean sea surface temperature (SST) and sea ice fields, and the UM's original land-surface scheme (Joint UK Land Environment Simulator; 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., 2013) is used, and ozone is coupled interactively between chemistry and radiation. The aerosol component includes sulfur 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 terrain-following 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 stratosphere. The model dynamical time step is 20 min, 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 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).

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.

The ACCESS-UKCA setup used here incorporates some additional modifications compared to the base UM-UKCA version 8.4, and these include the following:

• 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 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₂ (cm³ molecule⁻¹ s⁻¹) for the chemical reaction HO₂+ NO → HNO₃:

  $$\begin{array}{}\text{(2)}& \begin{array}{rl}{k}_{\mathrm{2}}/k& =\mathit{\left\{}\right[(\mathrm{530}\pm \mathrm{10})/T]+\mathrm{4.8}\times {\mathrm{10}}^{-\mathrm{4}}\times P\\ & -(\mathrm{1.73}\pm \mathrm{0.07})\mathit{\}}/\mathrm{100},\end{array}\end{array}$$

  where pressure P is in hPa and temperature T is in K, and k is the rate constant (cm³ molecule⁻¹ s⁻¹) for the included chemical reaction HO₂+ NO → NO₂+ HO. The coefficient 8.53 is correct when P is expressed in Torr (Butkovskaya et al., 2007), but it should be 4.8 when P is in hPa¹. Additionally, in the parameterisation $k=c_{\mathrm{0}}\times {\mathrm{10}}^{-\mathrm{12}}\exp (\mathrm{270}/T)$, the model uses c₀=3.3 for this reaction, but for k in the branching ratio (Eq. 2) it uses c₀=3.6, which we change to 3.3 for consistency.

The above changes led 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₃, and this increase is towards a reported global modelling average (see Sect. 4.2).

In the model, the tropopause is calculated every time step. In the extratropics (latitudes $\ge \left|\mathrm{28}\right|$^∘), the tropopause is the pressure level of the 2 PVU (potential vorticity unit) surface, and in the tropics (latitudes $\le \left|\mathrm{13}\right|$^∘) it is the pressure level of the 380 K potential temperature isentropic surface. Between the two regions, a weighted average of the two definitions is used following the method of Hoerling et al. (1993).

2.1 Implementation of lightning flash-rate parameterisation in the model

In ACCESS-UKCA, the convective cloud bottom level (H_b) and top level (H) are diagnosed on a time step basis from the UM convection scheme. This scheme represents the subgrid-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 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 time step, 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).

A minimum convective cloud scale needs to be specified for it to constitute a thunderstorm. In our model, a minimum convective cloud thickness (H−H_b) 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 (Price and 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 implicitly assumed that a convective cloud always translates to a thunderstorm, and this would likely lead to unrealistically 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 Sect. 3.2).

The PR92 flash-rate expressions (discussed later in Sect. 3.2) were developed based on observations of individual thunderstorms. Price and Rind (1994) (see their Fig. 1) developed a spatial calibration factor (c) to adjust these expressions for varying model resolutions:

$$\begin{array}{}\text{(3)}& c=\mathrm{0.97241}\exp \left(\mathrm{0.048203}\mathrm{\Delta }x\mathrm{\Delta }y\right),\end{array}$$

where Δx is the longitudinal resolution and Δy 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

$$\begin{array}{}\text{(4)}& f_{\mathrm{L},\mathrm{O}}={\displaystyle \frac{c{F}_{\mathrm{L},\mathrm{O}}A}{A_{\mathrm{c}}}},\end{array}$$

where A is the model grid box area (which is a function of latitude), A_c is the area of model grid box centred at 30^∘ N (Allen 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 f_L,O which is then apportioned into CG and IC flash rates using an empirical parameterisation for the ratio $z_{\mathrm{R}}=\mathrm{IC}/\mathrm{CG}$ developed by Price and Rind (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 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 and obtained z_R=3.76 and a CG to total flash ratio of 0.21 over the tropics (35^∘ N–35^∘ S).

2.2 Production of NO per flash

Both the moles of NO produced per flash, P_NO, and the variation of this parameter between CG and IC flashes are poorly constrained by atmospheric observations. The overall rate, P_NO, regulates the amount of nitrogen oxides produced by lightning, whereas the variation of P_NO 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, P_NO is set at S_f×10²⁶ molecules NO per flash, where the scaling factor S_f=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¹⁶ molecules NO J⁻¹. 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¹⁶ molecules NO J⁻¹.

In bottom-up models, in addition to flash rate, P_NO 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 (Bucsela et al., 2019). 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, 2016).

The value P_NO=330 moles NO per flash used in our model lies close to the middle of the range of current literature. Recent estimates include a global average value of 310 moles NO per flash obtained by Miyazaki et al. (2014) using an assimilation of multiple satellite measurements of atmospheric composition and the LIS/OTD lightning flash data into a global CTM; 665 moles NO per flash estimated by Nault et al. (2017) using airborne observations of atmospheric composition, satellite-based Ozone Monitoring Instrument (OMI) NO₂ columns and the GEOS-Chem model; 280±80 moles NO per flash by Marais et al. (2018) using the OMI NO₂ columns and satellite-based lightning data together with GEOS-Chem; 180±100 moles NO per flash by Bucsela et al. (2019) for three northern midlatitude regions that were primarily continental; and 170±100 moles NO per flash by Allen et al. (2019) for the tropics. The last two stem from the same OMI NO₂ columns and ground-based lightning measurements. Values used in calculating global estimates of LNO_x include 360 moles NO per flash by Ott et al. (2007), 500 moles NO per flash for selected extratropical regions, and 260 moles NO per flash for the rest of the globe by Murray et al. (2012).

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. (2007, 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 P_NO may not be constant over the globe, with higher production rates in 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 Sect. 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 vertical in log-pressure coordinate from 500 hPa to the cloud top for IC flashes, and from 500 hPa to surface for CG flashes (see Sect. 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_{\mathrm{R}}=\mathrm{IC}/\mathrm{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 which the CG-generated LNO_x is distributed as the observed minimum freezing level plus half of the minimum cold-cloud thickness, i.e. $(\mathrm{2.7}+\mathrm{5.5}/\mathrm{2})\approx \mathrm{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.

3 Flash-rate parameterisations based on convective cloud-top height and LNO_x

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 (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 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 that convective velocity does. (However, as Williams (1985) writes, this does not establish whether convective motion is directly responsible for the charge separation.)

3.1 Boccippio's (2002) extension of Vonnegut's (1963) model and derivation of scaling relationships for land and ocean

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 conceptualised 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 conceptualised 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 (the product of charge density multiplied by 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 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

$$\begin{array}{}\text{(5)}& F\sim k_{\mathrm{1}}w{H}^{\mathrm{4}},\end{array}$$

where the coefficient $k_{\mathrm{1}}\sim \mathit{\gamma }{\mathit{\epsilon }}^{-\mathrm{1}}{\mathit{\rho }}^{\mathrm{2}}$, H is cloud-top height (m), w is updraft velocity (m s⁻¹), ρ is charge density (coulomb m⁻³), ε is permittivity constant (coulomb² J⁻¹ m⁻¹), and the coefficient γ has units of J⁻¹. By assuming that each flash is responsible for a constant electrostatic energy and that the variability in the charge density of the dominant dipole regions is small, Boccippio (2002) treated k₁ 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⁻¹.

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 w(z) values for land and ocean are considered to be power-law fits of H,

$$\begin{array}{}\text{(6)}& w_{\mathrm{L},\mathrm{O}}=k_{\mathrm{L},\mathrm{O}}{H}^{a_{\mathrm{L},\mathrm{O}}}.\end{array}$$

The above approach based on Eq. (5) leads to the following self-consistent scaling relationships, where a_L,O and k_L,O are coefficients for land and ocean:

$$\begin{array}{}\text{(7)}& {\displaystyle }{F}_{\mathrm{L},\mathrm{O}}=k_{\mathrm{1}}{k}_{\mathrm{L},\mathrm{O}}{H}^{a_{\mathrm{L},\mathrm{O}}+\mathrm{4}}\text{(8)}& {\displaystyle }{F}_{\mathrm{L},\mathrm{O}}=k_{\mathrm{1}}{k}_{\mathrm{L},\mathrm{O}}^{-\mathrm{4}/a_{\mathrm{L},\mathrm{O}}}{w}^{\mathrm{1}+(\mathrm{4}/a_{\mathrm{L},\mathrm{O}})}.\end{array}$$

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) derived

$$\begin{array}{}\text{(9)}& F_{\mathrm{L}}=\mathrm{2.13}\times {\mathrm{10}}^{-\mathrm{5}}{H}^{\mathrm{5.09}}\end{array}$$

(note there is a negative sign missing in the first exponent of this equation in Boccippio, 2002), and

$$\begin{array}{}\text{(10)}& F_{\mathrm{O}}=\mathrm{4.09}\times {\mathrm{10}}^{-\mathrm{5}}{H}^{\mathrm{4.38}},\end{array}$$

where $k_{\mathrm{1}}=\mathrm{1.4314}\times {\mathrm{10}}^{-\mathrm{5}}$ was used in deriving the relationships in Eqs. (9) and (10), which are plotted in Fig. 1.

Figure 1 Various relationships, including those from this study (TS), between lightning frequency F (flashes per minute) and convective cloud-top height H (km). The land data points are from PR92, and the ocean data are from Molinié and Pontikis (1995) (MP95) and Ushio et al. (2001) (Us01); Bo02 – Boccippio (2002), Mi99 – Michalon et al. (1999), and PR92 – Price and Rind (1992).

3.2 Price and Rind's (1992) (PR92) parameterisations

PR92, following Vonnegut (1963) and Williams (1985), present a simple lightning parameterisation for calculating global lightning distributions. Using continental storm observations used by Williams (1985) to establish the fifth-power dependency, PR92 derived an empirical relationship between continental lightning flash rate and H:

$$\begin{array}{}\text{(11)}& F_{\mathrm{L}}=\mathrm{3.44}\times {\mathrm{10}}^{-\mathrm{5}}{H}^{\mathrm{4.9}}.\end{array}$$

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 intensity of convective updrafts enhances cloud electrification, and derived empirical relations based on observations between maximum updraft velocity (w, m s⁻¹) and H (km) for continental and marine clouds:

$$\begin{array}{}\text{(12)}& {\displaystyle }{w}_{\mathrm{L}}=\mathrm{1.49}{H}^{\mathrm{1.09}}\text{(13)}& {\displaystyle }{w}_{\mathrm{O}}=\mathrm{2.86}{H}^{\mathrm{0.38}}.\end{array}$$

Eliminating H from Eqs. (11) and (12) yields

$$\begin{array}{}\text{(14)}& w_{\mathrm{L}}=\mathrm{14.66}{F}_{\mathrm{L}}^{\mathrm{0.22}}.\end{array}$$

Now, PR92 made the crucial assumption that Eq. (14) is independent of location, so it is valid for the marine environment too. Thus, by equating Eqs. (14) and (13), they obtained the following expression for marine environment:

$$\begin{array}{}\text{(15)}& F_{\mathrm{O}}=\mathrm{6.4}\times {\mathrm{10}}^{-\mathrm{4}}{H}^{\mathrm{1.73}}.\end{array}$$

(Our calculation suggests that the coefficient 6.4 in the above equation should actually be 5.94.) The relationships in Eqs. (11) and (15) are plotted in Fig. 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. an inverse relationship between updrafts and cloud heights.

PR92 parameterisations (11) and (15) are widely used in global chemistry transport models and coupled chemistry–climate models.

3.3 Michalon et al.'s (1999) ocean parameterisation

Michalon et al. (1999) identified that the PR92 parameterisation, when used over water in global models, was producing 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 $F=A{N}^{\mathrm{2}/\mathrm{3}}{H}^{\mathrm{5}}$, 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 $(N_{\mathrm{O}}/N_{\mathrm{L}}{)}^{\mathrm{2}/\mathrm{3}}$ assuming that cloud droplet concentrations for continental (N_L) and marine (N_O) clouds are different:

$$\begin{array}{}\text{(16)}& F_{\mathrm{O}}=\mathrm{3.44}\times {\mathrm{10}}^{-\mathrm{5}}{\left({\displaystyle \frac{N_{\mathrm{O}}}{N_{\mathrm{L}}}}\right)}^{\mathrm{2}/\mathrm{3}}{H}^{\mathrm{4.9}}=\mathrm{6.57}\times {\mathrm{10}}^{-\mathrm{6}}{H}^{\mathrm{4.9}},\end{array}$$

where two “standard” continental and maritime cloud droplet concentrations of N_L=600 and N_O=50 per mg of air, respectively, were used. Equation (16) is plotted in Fig. 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 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.

3.4 This study: alternative flash-rate parameterisations

This study includes (a) a reanalysis of the Williams (1985) and PR92 data for lightning flashes vs. 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 in Eq. (6) with the observed maximum updraft velocity from PR92 given in Eq. (12) gives k_L=1.49 and a_L=1.09. Substituting these values into Eq. (7) yields

$$\begin{array}{}\text{(17)}& F_{\mathrm{L}}=\mathrm{1.49}{k}_{\mathrm{1}}{H}^{\mathrm{5.09}}.\end{array}$$

At this stage, k₁ 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 Fig. 1 as solid blue circles). This gives $k_{\mathrm{1}}=\mathrm{1.612}\times {\mathrm{10}}^{-\mathrm{5}}$ and

$$\begin{array}{}\text{(18)}& F_{\mathrm{L}}=\mathrm{2.40}\times {\mathrm{10}}^{-\mathrm{5}}{H}^{\mathrm{5.09}}.\end{array}$$

The relationships of PR92 (Eq. 11), Boccippio (2002) (Eq. 9), and this study (TS) (Eq. 18) are presented in Fig. 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, 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.

Concerning the oceanic parameterisation, there are limited data on flash rate vs. cloud-top height for marine environments. Figure 1 presents Ushio et al.'s (2001) (Us01) data (triangles) over the ocean in the tropics and extratropics (their Fig. 3b and d), which they obtained by averaging flash rates every 1 km in cloud height. Also shown are the Molinié and Pontikis' (1995) (MP95) flash-rate data (solid black circles) over the 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 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 Fig. 1.

Applying Boccippio's (2002) scaling relationships to obtain an equivalent relationship for marine clouds involves equating Eq. (6) with the observed oceanic convective updraft velocity vs. cloud-top height from PR92 given in Eq. (13). This yields k_O=2.86 and a_O=0.38. Substituting these into the scaling relationship (Eq. 7) gives

$$\begin{array}{}\text{(19)}& F_{\mathrm{O}}=k_{\mathrm{1}}\times \mathrm{2.86}\times H^{\mathrm{4.38}}=A_{\mathrm{O}}{H}^{\mathrm{4.38}}.\end{array}$$

For this study, fitting Eq. (19) to the ocean data in Fig. 1 leads to $A_{\mathrm{O}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ (which gives $k_{\mathrm{1}}=\mathrm{0.7}\times {\mathrm{10}}^{-\mathrm{5}}$) and

$$\begin{array}{}\text{(20)}& F_{\mathrm{O}}=\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{5}}{H}^{\mathrm{4.38}}.\end{array}$$

This marine parameterisation yields flash rates that are approximately an order of magnitude smaller than the PR92 continental formula and roughly 2 orders of magnitude larger than the PR92 marine parameterisation. In Fig. 1, the relationships for lightning flash rates and associated cloud-top height over the ocean for PR92, Boccippio's (2002) Eq. (10), Michalon et al. (1999), and TS Eq. (20) are shown. 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 k₁ 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 one or both assumptions in the theory 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 k₁ in Eqs. (18) and (20) for land and ocean can be interpreted in terms of N_O and N_L being different from N_O<N_L.

3.5 Global model runs with various flash-rate parameterisations

In order to assess the above flash-rate parameterisations against global lightning flash observations and to investigate how they influence tropospheric composition via their impact on LNO_x generation, we conduct the following five runs of the ACCESS-UKCA global chemistry–climate model incorporating seven specified flash-rate parameterisations:

• Run 1 (PR92): the default PR92 parameterisations: continental Eq. (11) and marine Eq. (15);

• Run 2 (this study – TS1): the new parameterisations: continental Eq. (18) and marine Eq. (20);

• Run 3 (this study – TS2): the PR92 continental Eq. (11) and new marine Eq. (20);

• Run 4 (Mi99): the PR92 continental Eq. (11) and Michalon et al.'s (1999) marine Eq. (16);

• Run 5 (Bo02): Boccippio's (2002) Eqs. (9) and (10).

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 set up 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 Newtonian relaxation nudging (Uhe and Thatcher, 2015) within model levels 20–45 (between altitudes ∼3 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 did Runs 1 and 2 with meteorological nudging for the years 2005–2006 with the same initial conditions as for the free-running simulations, and the results are summarised in Sect. 5.)

3.6 Comparison of the modelled lightning flash rates with satellite observations

We analyse the global gridded lightning flash data from OTD on the OrbView-1 (formerly MicroLab-1) satellite and the Lightning Imaging Sensor (LIS) on the Tropical Rainfall Measuring Mission (TRMM) satellite, which are described by Cecil et al. (2014) and available from https://lightning.nsstc.nasa.gov/data/data_lis-otd-climatology.html (V2.3.2015, last access: 5 May 2021). 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⁻² yr⁻¹ and flashes km⁻² d⁻¹, respectively. The OTD data available from July 1995 to January 2000 cover all 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.

Figure 2 Observed vs. modelled flash rates from the five model runs (corresponding to PR92, TS1 and TS2, Mi99, and Bo02) conducted for the year 2006. The vertical clusters of points are the values over the globe, ocean, land, Southern Hemisphere (SH), Northern Hemisphere (NH). (a) The observations are the LIS/OTD climatological data, and (b) the observations are the LIS data for the year 2006 that are available for the latitude range ±42.5^∘ (with the modelled values also corresponding to this range).

Table 1 Observed and simulated lightning flash frequencies for the various model runs. The values in the parentheses are for the latitudinal range ±42.5^∘ (which correspond to the latitudinal coverage of the LIS data).

Table 1 gives the observed and modelled lightning flash frequencies (flashes s⁻¹) averaged over the globe, Northern Hemisphere (NH), Southern Hemisphere (SH), land, and ocean for the five model runs, and these are plotted in Fig. 2. In Fig. 2a, the observations are the combined LIS and OTD climatological data, whereas in Fig. 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 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 parameterisations lead to an almost equal partitioning of flashes in the Northern Hemisphere and Southern Hemisphere, 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 predicts flash rate 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 that from Run 2.

Figure 3 Monthly variation of the observed and modelled flash rates from the five model runs for the year 2006: (a) NH, (b) SH, (c) globe, (d) land, and (e) ocean. The observations are the LIS/OTD climatological data (with 1 standard deviation variability) and the LIS data for the year 2006.

The monthly variation of the observed and modelled flash rates is shown in Fig. 3. The observed variation for the year 2006 is very similar to the climatological variation, mostly within the 1 standard deviation climatological variability shown. The model runs underpredict in spring in both hemispheres and overpredict in autumn in the Southern Hemisphere. In the Northern Hemisphere (Fig. 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–May and do well for the other months. In the Southern Hemisphere (Fig. 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 (Fig. 3d). In the global plot (Fig. 3c), while Run 1 always underestimates the observed flash rate, which is mostly because of its underprediction of the flash rate over the ocean (Fig. 3e), the other runs underestimate it in the first 3 months of the year and overestimate it during September–November, and these differences can be explained in terms of the hemispheric differences shown above. The nature of monthly variation in Fig. 3d for land-based flash rates is very similar to that in Fig. 3c, indicating that continental flashes dominate the global total. The large underprediction by the PR92 scheme in Run 1 and an overprediction in Run 4 over the ocean can be seen in Fig. 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 Normalised mean square error (NMSE), fractional bias (FB), and correlation coefficient (r) for the monthly varying observed climatology and modelled flash rates shown in Fig. 3.

Normalised mean square errors (NMSEs) ($\stackrel{\mathrm{‾}}{(M-O{)}^{\mathrm{2}}}/\left(\overline{M}\overline{O}\right)$) calculated from the monthly varying observed climatology (O) and modelled (M) flash-rate time series shown in Fig. 3 are given in Table 2 for the globe, land, and ocean. Also given are the values of fractional bias (FB) $\mathrm{2}(\overline{O}-\overline{M})/(\overline{O}+\overline{M})$ which varies between −2 (overestimation) and +2 (underestimation), and the (Pearson) correlation coefficient (r). Considering these performance statistics together suggests that the Run 2 flash-rate parameterisations from this study yield the best comparison with the data.

Figure 4 Observed and modelled zonal mean flash densities (flashes km⁻² yr⁻¹) from the five model runs over the (a) globe, (b) land, and (c) ocean. The observations are the LIS/OTD climatology and the LIS data for the year 2006.

Figure 4 presents the observed and modelled zonal mean flash density (flashes km⁻² yr⁻¹) over the globe, land, and ocean. All modelled global distributions and the data agree that the flash density is largely concentrated in the tropics. The observed peak for the year 2006 in Fig. 4a is better simulated by Runs 2–4 than by Run 1 (underestimation) or Run 5 (overestimation). The results over land (Fig. 4b) are very similar for all runs. Over the ocean (Fig. 4c), while the default oceanic parameterisation (Run 1) yields a near-zero flash density distribution and Run 5 overestimates the data considerably, the new flash parameterisation (in Runs 2 and 3) performs much better. There are significant distributional differences compared to the data. It is clear in these plots that the observed latitudinal distributions of flash density are wider than the modelled ones, with larger observed flash densities in the subtropics stretching into the midlatitudes (roughly 20–40^∘ in both hemispheres) than modelled. The reason for this may be the inherent limitation of the simple flash parameterisation approach based on convective cloud-top height or uncertainty and/or 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 5 Global distribution of the mean annual lightning flash density (flashes km⁻² yr⁻¹): (a) LIS/OTD satellite climatology, (b) 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).

Figure 5 compares the various global distributions of the mean annual lightning flash density at N96 resolution. The LIS data for the year 2006 in Fig. 5b are only available for the latitudinal range ±42.5^∘ and have no sampling for some regions within this range. Where there are data, there is a good agreement between the observed distribution in Fig. 5b and the LIS/OTD climatology in Fig. 5b, showing high flash density over land in the tropics and subtropics, and also in the lower midlatitudes. There is also some significant flash density over the ocean at these latitudes, particularly over the Pacific, western 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 Fig. 5c is very similar to other global modelling studies that use the same PR92 scheme, e.g. Allen and Pickering (2002) using the Goddard Earth Observing System – Stratospheric Tracers of Atmospheric Transport Data Assimilation System (GEOS-STRAT DAS), Tost et al. (2007) using ECHAM5/MESSy, Murray et al. (2012) using GEOS-Chem, Finney et al. (2014) using ERA-Interim reanalyses, Finney et al. (2016) using UM-UKCA, and Clark et al. (2017) using CAM5. It is remarkable that the simple PR92 scheme based on the convective cloud-top height is able to simulate the broad observed global distribution of flash density over land at low latitudes (except parts of India) but it does not properly reproduce the extension of lightning flash density into the temperate latitudes, particularly in the Northern Hemisphere. Over the ocean, in contrast to the observations, the PR92 scheme predicts almost zero flash density. However, as shown in Fig. 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 Fig. 5c and d are nearly the same.

Table 3 NMSE, FB, and r for the spatially varying annual-mean modelled flash rates and observed climatology shown in Fig. 5a.

The area-weighted NMSE, FB, and r comparing the spatial patterns of the observed climatology presented in Fig. 5a with the annually averaged modelled patterns of flash rate are given in Table 3 for the five model 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 (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 is 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).

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⁻¹), 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 of the observed global average flash rate to the modelled global average flash rate and is calculated by doing a model pre-run. 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 overadjust the flash rate (and hence LNO_x) over land to compensate for the deficiency in the oceanic parameterisation. (Scaling can also be applied to tune the amount of NO produced per flash to get a desired total global 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.

3.7 Modelled LNO_x and comparison

3.7.1 Global LNO_x

The modelled global lightning-generated NO_x using the various lightning flash-rate schemes is presented in Table 4. With the new flash-rate parameterisations (Run 2), the modelled global LNO_x increases to 6.6 Tg N yr⁻¹ from 4.8 Tg N yr⁻¹ 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, 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 the NH, SH, 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.

Table 4 Modelled global lightning-generated NO_x using various lightning flash-rate schemes (Tg N yr⁻¹).

The amount of global LNO_x produced, L_G (Tg N yr⁻¹), is a function of the global average flash rate, f_s (flash s⁻¹), and the moles of NO produced per flash, P_NO:

$$\begin{array}{}\text{(21)}& L_{\mathrm{G}}=\mathrm{441.5}\times {\mathrm{10}}^{-\mathrm{6}}{f}_{\mathrm{s}}{P}_{\mathrm{NO}}.\end{array}$$

If the climatological average f_s=46.5 flash s⁻¹ based on the LIS/OTD satellite data is used, then

$$\begin{array}{}\text{(22)}& L_{\mathrm{G}}=\mathrm{20.5}\times {\mathrm{10}}^{-\mathrm{3}}{P}_{\mathrm{NO}}.\end{array}$$

If the NO production per flash differs for IC and CG flashes, then P_NO 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 P_NO=330 moles NO per flash as used in ACCESS-UKCA and the modelled f_s from Table 1 are substituted. These values can be compared with a global estimate of $L_{\mathrm{G}}=\mathrm{5}\pm \mathrm{3}$ Tg N yr⁻¹ based on Schumann and Huntrieser's (2007) review. More recently, there have been estimates of LNO_x incorporating top-down approaches, which we divide into (a) verification studies using other constraints and (b) model experiments.

Verification studies include the following. Using a CTM representing the LIS/OTD flash data, Martin et al. (2007) obtained an estimate of 6±2 Tg N yr⁻¹ that best reproduced satellite observations of tropospheric NO₂, O₃ and HNO₃. Miyazaki et al. (2014) obtained $L_{\mathrm{G}}=\mathrm{6.3}\pm \mathrm{1.4}$ Tg N yr⁻¹ and a global average production of 310 moles NO per flash using an assimilation of multiple satellite measurements of NO₂, O₃, HNO₃ and CO, and the LIS/OTD flash data into a global CTM. The global values for Runs 2–4 in Table 4 compare very well with the above constrained estimates, considering their differences in NO per flash used (the ratio $L_{\mathrm{G}}/P_{\mathrm{NO}}\approx \mathrm{0.02}$ in these estimates as per Eq. 22). Using upper-tropospheric airborne observations of LNO_x and global satellite-retrieved tropospheric NO₂ column densities along with GEOS-Chem, Nault et al. (2017) estimated 665 moles NO per flash with global LNO_x at ∼9 Tg N yr⁻¹. Using global satellite data of NO₂ columns and lightning flashes together with the GEOS-Chem model, Marais et al. (2018) derived a global production rate of 280±80 moles NO per flash, with a global LNO_x emission of 5.9±1.7 Tg N yr⁻¹.

Boersma et al. (2005) analysed above-cloud tropospheric NO₂ column retrievals from the Global Ozone Monitoring Experiment (GOME) satellite observations for the year 1997 for cloudy scenes over tropical oceans and continents, and found that the above-cloud annual-mean NO₂ column increases sharply with convective cloud-top height (H) – as H^5.1 for continents and H^4.6 for oceans, where H>6.5 km. Considering that these above-cloud NO₂ columns primarily consist of contributions from lightning-generated NO_x, which is a direct function of flash rate, there is very good agreement between these power-law exponents and those in the flash-rate relationships of Eqs. (11) and (18) for continents and (20) for oceans. The analysis of Boersma et al. (2005) demonstrates that the exponent of 1.73 in the PR92 oceanic flash-rate relationship (15) is unrealistic.

Model experiments include the following. Using a 3-D cloud-resolving model coupled with observations from a thunderstorm and assuming P_NO=360 moles NO per flash, Ott et al. (2007) estimated a global LNO_x of 7 Tg N yr⁻¹. Using a CTM constrained by the LIS/OTD flash together with a production of 500 moles NO per flash for all extratropical lightning north of 23^∘ N in America and 35^∘ N in Eurasia, and 260 moles NO per flash for the rest of the globe, Murray et al. (2012)² determined a global LNO_x of 6±0.5 Tg N yr⁻¹.

Figure 6 Global distribution of the annual-averaged LNO_x (10⁻¹³ kg N m⁻² s⁻¹) for the year 2006: (a) model simulation with the default PR92 flash-rate parameterisations (Run 1), (b) model simulation with the new flash-rate parameterisations from this study (Run 2), and (c) the distribution obtained by Miyazaki et al. (2014, plot redrawn) using assimilation of multiple satellite datasets into a global CTM for the year 2007. The respective global LNO_x totals are 4.8, 6.6, and 6.3 Tg N yr⁻¹.

The modelled mean global distributions of LNO_x from the two runs presented in Fig. 6a and b are essentially in proportion to the flash density distributions given in Fig. 5c and d, 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 Fig. 6c the annual LNO_x distribution obtained by Miyazaki et al. (2014) using an assimilation of satellite measurements of atmospheric composition and the LIS/OTD lightning flash data into a global CTM for the year 2007. This plot is a reproduction of their Fig. 6 (middle-left plot) based the data³ supplied by Kazuyuki Miyazaki (personal communication, 2020) at a horizontal resolution of 2.8^∘ × 2.8^∘. Over the ocean, the new flash-rate scheme (Fig. 6b) agrees much better with the assimilated field than does the PR92 scheme, but it 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 Fig. 6b nor indicated by the observed flash-rate distributions in Fig. 5a and b (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 Fig. 6c at low latitudes (except parts of India) but do not properly 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⁻¹ 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 more extensive measurements would be necessary for a better evaluation of the predicted LNO_x distribution.

Figure 7 Average vertical distribution of percentage of LNO_x mass per kilometre for (a) tropical continental, (b) tropical marine, (c) midlatitude continental, and (d) subtropical regimes. The total LNO_x for these regimes calculated using model Run 1 (PR92) is 3.69, 0.035, 1.09, and 0.74 Tg N yr⁻¹, respectively, whereas that calculated using Run 2 (TS1) is 4.10, 1.09, 1.16, and 0.92 Tg N yr⁻¹, respectively. The vertical profiles from Pickering et al. (1998) and Ott et al. (2010) are also shown (where available).

3.7.2 Vertical distribution

The vertical distribution of LNO_x in the model at a grid point location is a parameter (see Sect. 2.2) that is essentially unconstrained by observations. Figure 7 presents the modelled average vertical distribution of percentage of LNO_x mass in each 1 km layer for (a) tropical continental, (b) tropical marine, (c) midlatitude continental, and (d) subtropical regimes. The non-uniform shape of the averaged modelled vertical distributions is largely caused by the averaging of the LNO_x profile from every time step over spatial and temporal variations in the cloud-top height. Also shown for comparison are the average profiles based on thunderstorm cases simulated by Pickering et al. (1998) using a 2-D convective cloud-resolving tracer transport model and those by Ott et al. (2010) using a 3-D convective cloud-resolving chemical transport model, with both studies using parameterised lightning. The profiles of Pickering et al. (1998) show peaks near the surface, as significant mass is transported to the boundary layer by downdrafts, and in the upper troposphere (the so-called “C-shaped” profile), whereas those by Ott et al. (2010) show very little LNO_x mass in the boundary layer with the majority of LNO_x remaining in the middle and upper troposphere (the so-called “backward C-shaped” profile) where it is originally produced. Our model profiles match better with the Ott et al. (2010) profiles, but the model gives an almost uniform distribution of LNO_x mass below 5 km, whereas the latter decrease to almost zero. There is not a large difference between the modelled profiles for the various regimes, except that the tropical ones are almost uniformly distributed between 5 and 12 km, whereas the midlatitude and 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.

Figure 8 Zonal-averaged tropospheric total NO₂ column (N_v,trop) (in units 10¹⁵ molecules cm⁻²) obtained from Run 1 (with the default PR92 flash-rate parameterisations) and Run 2 (with the new flash-rate parameterisations from this study), and the OMI satellite data (green diamonds), for the year 2006. The orange data points are the OMI satellite data ($N_{\mathrm{v},\mathrm{trop},\mathrm{180}}$) over the longitude 180^∘ W.

3.7.3 Modelled tropospheric total column NO₂ and validation

As lightning impacts atmospheric NO_x directly, any changes in the modelled total tropospheric NO_x can be examined and compared with available observations. The modelled variations of the zonal-averaged tropospheric total NO₂ column (N_v,trop) for the globe presented in Fig. 8 show a broad peak in the industrialised Northern Hemisphere dominated by surface NO₂ emissions. Within the latitudes ±30^∘, the new lightning flash-rate parameterisations yield 12 % larger values of tropospheric total NO₂ column compared to the default PR92 scheme and the difference between the two is $\approx \mathrm{0.1}\times {\mathrm{10}}^{\mathrm{15}}$ molecules cm⁻².

Figure 8 also presents tropospheric NO₂ column retrievals from the Ozone Monitoring Instrument (OMI) satellite overpasses at ∼13:30 local time (LT) for the year 2006. These zonal averages were derived from the OMI monthly mean tropospheric NO₂ columns (QA4ECV, version 1.1) given at a horizontal resolution of 0.125^∘ × 0.125^∘ (https://www.temis.nl/airpollution/no2.php, last access: 6 May 2021; Boersma et al., 2018). While the model simulations qualitatively agree with the satellite observations in Fig. 8, it is apparent that except for the latitudes 30^∘ S–60^∘ S the modelled values are generally higher than the observations. Notwithstanding any model shortcomings in predicting the global NO₂ distribution, there are limitations of the OMI satellite data used. Firstly, there is limited sensitivity of the OMI sensor to NO₂ below the cloud level, where most NO₂ is situated, and thus cloudy tropospheric NO₂ 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, 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 Copernicus Atmosphere Monitoring Service (CAMS) (Inness et al., 2019; https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-global-reanalysis-eac4-monthly, last access: 6 May 2021) for comparison with the monthly averaged modelled NO₂. 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₂, the model assimilated the tropospheric column retrievals from the SCIAMACHY, OMI, and GOME-2 satellite overpasses at ∼10:00, 13:30, and 09:30 LT, respectively (for the year 2006, only SCIAMACHY and OMI data were available for assimilation). Monthly averaged total vertical column NO₂ reanalysis data (version – ECMWF Atmospheric Composition Reanalysis 4) are available and used here.

We obtain the tropospheric NO₂ column (N_v,trop) from the CAMS total NO₂ column (N_v) as follows:

$$\begin{array}{}\text{(23)}& N_{\mathrm{v},\mathrm{trop}}=N_{\mathrm{v}}-N_{\mathrm{v},\mathrm{180}}+N_{\mathrm{v},\mathrm{trop},\mathrm{180}},\end{array}$$

where N_v,180 is the CAMS total NO₂ column over the Pacific (180^∘ W) and $N_{\mathrm{v},\mathrm{trop},\mathrm{180}}$ is the tropospheric NO₂ 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₂ column amounts and a constant and negligibly small $N_{\mathrm{v},\mathrm{trop},\mathrm{180}}$ (Lauer et al., 2002). $N_{\mathrm{v},\mathrm{trop},\mathrm{180}}$ is not available from the CAMS data, but that obtained from the OMI data discussed above is shown in Fig. 8 (orange data points) for the year 2006. While these $N_{\mathrm{v},\mathrm{trop},\mathrm{180}}$ values are small, they are neither constant nor negligibly small compared to the modelled or OMI N_v,trop plotted in Fig. 8 and would reflect contributions from sources such as lightning, aviation, shipping, and possibly regional transport over the ocean. The $N_{\mathrm{v},\mathrm{trop},\mathrm{180}}$ values are also greater and possibly more uncertain than the differences between the two model simulations.

Figure 9 Zonal-averaged total tropospheric NO₂ column (N_v,trop) (in units 10¹⁵ molecules cm⁻²) modelled using the default PR92 parameterisations (Run 1), the new lightning flash-rate parameterisations from this study (Run 2), and the CAMS reanalysis data for the year 2006 over (a) globe, (b) ocean, and (c) land. The model variation without LNO_x is also shown.

Since $N_{\mathrm{v},\mathrm{trop},\mathrm{180}}$ is not available from the CAMS data, we take the OMI-derived $N_{\mathrm{v},\mathrm{trop},\mathrm{180}}$ shown in Fig. 8 and use it in Eq. (23) to determine N_v,trop from the CAMS data. (While these OMI NO₂ 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 N_v,trop as a function of latitude for the globe is presented in Fig. 9a, together with the corresponding modelled variations. The modelled variations agree well with the CAMS data, better than with the direct OMI data in Fig. 8. A comparison with the model variation without LNO_x in Fig. 9a suggests that the modelled increase in NO_x due to lightning is largely confined to ±35^∘.

In Fig. 9b for the ocean, the agreement with the data is again good, except for considerable model underestimation for latitudes greater than 35^∘ which, presumably not related to LNO_x, may be due to factors such as possible underestimation of the calculated $N_{\mathrm{v},\mathrm{trop},\mathrm{180}}$ and overestimation of shipping emissions of NO_x for these latitudes. Within ±35^∘, the new parameterisation yields 22 % larger NO₂ 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₂ columns are somewhat better simulated by the model with the new oceanic flash rate than with the PR92 parameterisation, particularly in the northern tropics.

For land (Fig. 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.

Table 5 Mean (10¹⁵ molecules cm⁻²), NMSE, FB, and r for the modelled and CAMS reanalysis NO₂ columns shown in Fig. 9 for latitudes between ±30^∘.

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 somewhat better prediction. Obviously, there are other sources of uncertainties in the NO₂ 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.

4 Impact on tropospheric composition

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₃, OH, methane lifetime, and CO.

Figure 10 Zonal annual-mean total tropospheric NO_x (as NO₂, pptv) modelled using (a) the default PR92 parameterisations (Run 1) and (b) the new lightning flash-rate parameterisations from this study (Run 2). The difference between Run 2 and Run 1 is shown in panel (c).

4.1 Oxides of nitrogen

The modelled tropospheric NO₂ columns and their comparison with observations have already been presented in Sect. 3.7.3. Figure 10 presents the zonal distribution of total NO_x (as NO₂) 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₂) 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 lightning-generated NO_x. The new lightning parameterisations cause an increase in the mid- to upper-tropospheric NO_x (Fig. 10c), particularly within the tropics and subtropics, and this increase is by as much as 40 pptv 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 pptv, and it is 35.2 pptv over the ocean and 94.0 pptv over land. With the new flash-rate scheme, these values increase by 8.7 pptv (15.7 %), 9.9 pptv (28.0 %), and 6.3 pptv (6.7 %), respectively, and can be compared with the values 36.9, 20.7, and 68.5 pptv, respectively, obtained from the model simulation with no LNO_x emissions. To some extent, the modelled tropospheric averages also depend on how the tropopause is defined.

4.2 Ozone

Tropospheric ozone is a byproduct of the oxidation of CO, CH₄, 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₃ burden increases from 284 to 308 Tg O₃, a rise of 8.5 % over the default PR92 scheme (Run 1) (cf. 219 Tg O₃ with no LNO_x emissions in the model). The new burden is closer to the ACCMIP (Atmospheric Chemistry and Climate Model Intercomparison Project) multi-model mean of 337±23 Tg O₃ 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.

Figure 11 Mean relative difference (%) between the global annual-mean ozone mixing ratios predicted by Run 2 using the new lightning flash-rate parameterisations from this study and the default PR92 parameterisations (Run 1): (a) at 20 m (the lowest model level) and (b) at 6400 m (∼450 hPa).

The mean relative difference (%) between the global ozone mixing ratios predicted using the new lightning flash-rate parameterisations and the default PR92 parameterisations is shown in Fig. 11. Near the surface (Fig. 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. southeastern US and northern Australia) where ozone has increased, and there are a few regions in the midlatitudes to high latitudes in the Northern Hemisphere where ozone has decreased very slightly. Tropospheric ozone chemistry is complex, but broadly speaking the O₃ increases in the Southern Hemisphere are influenced by low ambient NO_x concentrations where the O₃ production increases with NO concentration. O₃ is produced through photodissociation of NO₂, which is produced through oxidation of NO by HO₂ and RO₂ radicals (e.g. NO + HO₂ → NO₂ + OH). In the Northern Hemisphere, the increase in O₃ is lower beyond the tropics, 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 more responsive to the availability of NO than that of the (more polluted) Northern Hemisphere.

At an altitude of 6400 m (∼450 hPa) (Fig. 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.

Figure 12 Comparison of the modelled monthly averaged ozone concentrations with observations at five oceanic/coastal ground stations for the year 2006. The two model runs are with the default PR92 parameterisations (Run 1) and new lightning flash-rate parameterisations (Run 2, this study). The values of r and NMSE are also shown.

In Fig. 12, we compare the modelled monthly averaged ozone with ground-based in situ observations from the Global Atmosphere Watch – World Data Centre for Reactive Gases (GAW-WDCRG, http://ebas.nilu.no, last access: 6 May 2021; https://www.gaw-wdcrg.org, last access: 6 May 2021) 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 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 used the modelled concentrations from the bottom model level for all sites. The two model simulations describe the observed monthly variations reasonably well, except at Mauna Loa and Mace Head (the relatively large disagreement at Mauna Loa is 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₃ 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 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 at Ushuaia, where the differences are larger during winter to spring.

Figure 13 Profiles of annual-mean ozone constructed using the SHADOZ ozonesonde measurements at eight sites for the year 2006, and the corresponding modelled profiles obtained using the default PR92 parameterisations (Run 1) and the new lightning flash-rate parameterisations from this study (Run 2). The modelled profiles without LNO_x are also shown as a reference.

In Fig. 13, we compare the modelled profiles of annual-mean ozone with those obtained from the Southern Hemisphere ADditional OZonesondes (SHADOZ) ozonesonde measurements at eight tropical/subtropical sites (Witte et al., 2017; https://tropo.gsfc.nasa.gov/shadoz, last access: 6 May 2021) for the year 2006. These sites are Hilo (19.40^∘ N, 155.0^∘ W), Paramaribo (5.81^∘ N, 55.21^∘ W), Costa Rica (9.94^∘ N, 84.04^∘ W), Natal (Brazil) (5.42^∘ S, 35.38^∘ E), Ascension Island (7.98^∘ S, 14.42^∘ W), Irene (25.9^∘ S, 28.22^∘ E), Nairobi (1.27^∘ S, 36.80^∘ E), Réunion Island (21.1^∘ S, 55.48^∘ E). The above eight sites were selected based on data availability and to have a mix of Northern Hemisphere and Southern Hemisphere locations where O₃ could be more impacted by LNO_x. The profile data are given at a greater resolution in height than the vertical model resolution and were binned in the model vertical levels and averaged. The number of observed profiles at a location typically varied between one and five per month. At Irene, Nairobi, and Réunion 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.

The observed O₃ profiles are simulated better by the new lightning flash-rate parameterisations (Run 2) than by the default PR92 parameterisations (Run 1) at Costa Rica (Fig. 13c) and for all the southern hemispheric sites (Fig. 13d–h). At northern latitudes beyond Costa Rica, the PR92 scheme performs better at Hilo and Paramaribo (Fig. 13a and b). The model with the PR92 scheme has a general tendency to overestimate ozone at most levels within 15–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 Fig. 15). The model describes the general shape of the observed profiles reasonably well. The model profiles with no LNO_x included highlight the importance of LNO_x on tropospheric ozone, particularly in the middle to upper troposphere.

Figure 14 Zonal annual-mean tropospheric ozone (ppbv) modelled using (a) the default PR92 parameterisations (Run 1) and (b) the new lightning flash-rate parameterisations from this study (Run 2). The difference between Run 2 and Run 1 is shown in panel (c).

Figure 15 Zonal distribution of tropospheric ozone concentration (ppbv) for the year 2006: (a) observed distribution based on the global monthly mean vertical ozone profile database available from Bodeker Scientific; (b) the relative difference (%) between the concentration modelled using the default PR92 lightning parameterisations (Run 1) and the observations; (c) the relative difference (%) between the concentration modelled using the new lightning parameterisations (Run 2) and the observations; and (d) the relative difference (%) between the concentration modelled without any LNO_x and the observations.

The modelled zonal annual-mean tropospheric O₃ from the two runs and the difference between the two are presented in Fig. 14. In the lower troposphere, the modelled ozone is smaller over the Southern Hemisphere than over the Northern Hemisphere (Fig. 14a and b). The new flash-rate parameterisations result in O₃ increases everywhere (Fig. 14c). Closer to the surface, the increase is approximately 2 ppbv in the Southern Hemisphere and 0.5–2 ppbv in the Northern Hemisphere. The largest increases are nearly 8 ppbv in the tropics at altitudes ∼9 km.

The modelled zonal ozone distribution can be compared with observations. We use the monthly mean vertical ozone profile data for the year 2006, given as zonal means, from the Bodeker Scientific database (Hassler et al., 2009; http://www.bodekerscientific.com/data/monthly-mean-global-vertically-resolved-ozone, last access: 6 May 2021) which combines measurements from several satellite-based instruments and ozone profile measurements from the global ozone-sonde network. The database spans the period 1979 to 2016 with 5^∘ latitude resolution and 70 altitude levels (1 to 70 km). Different “tiers” of data are provided, and we used the highest tier (Tier 1.4 vn1.0) data. For comparison with the model predictions, these monthly profile data were regridded to the model resolution. The modelled monthly tropospheric mask was averaged zonally and then applied to the monthly regridded data.

There is an agreement, both in magnitude and distribution, between the modelled zonal ozone in Fig. 14a–b and the data plotted in Fig. 15a. The model reproduces the observed lower levels of tropospheric ozone in the Southern 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 ($=[\overline{M}-\overline{O})/\overline{O}]\times \mathrm{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 been 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 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 (Fig. 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 overestimation within 10–50^∘ N. Additional factors that influence tropospheric ozone distribution in the model include dynamics, including interhemispheric 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 the upper troposphere (Worden et al., 2011).

The volume-weighted tropospheric O₃ obtained from the PR92 scheme is 51.5 ppbv over the globe, 48.8 ppbv over the ocean, and 56.6 ppbv over land. With the new flash-rate scheme, these values increase by 4.1 ppbv (8.0 %), 4.4 ppbv (9.1 %), and 3.6 ppbv (6.3 %), respectively, and can be compared with the values 38.7, 36.2, and 43.4 ppbv, respectively, obtained from the model simulation with no LNO_x emissions.

4.3 Hydroxyl radical

OH is the dominant oxidising (and cleansing) agent in the global troposphere and controls the atmospheric abundance and chemical lifetime of most natural and anthropogenic gases, such as methane (CH₄). The tropospheric abundance of OH is determined by a complex series of chemical reactions involving tropospheric ozone, methane, 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.

Figure 16 Zonal annual-mean tropospheric OH (×10⁵ molecules cm⁻³) modelled using (a) the default PR92 parameterisations (Run 1) and (b) the new lightning flash-rate parameterisations from this study (Run 2). The difference between Run 2 and Run 1 is shown in panel (c).

The modelled zonal total annual-mean tropospheric OH in Fig. 16 shows highest OH concentrations near the surface in the tropics, with values as high as (25–30) × 10⁵ molecules cm⁻³ 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 (Fig. 16c), particularly in the upper troposphere in the tropics, by as much as 5×10⁵ molecules cm⁻³.

Figure 17 Mean relative difference (%) between the hydroxyl radical (OH) predicted using the new lightning flash-rate parameterisations (Run 2) and the default PR92 parameterisations (Run 1): (a) at 20 m (the lowest model level) and (b) at 6400 m (∼450 hPa).

The annual-mean relative OH difference (%) between Run 2 and Run 1 near the surface (Fig. 17a) shows an increase in OH over the Southern Hemisphere oceans and Antarctica, and pockets of increase and slight decrease in the Northern Hemisphere. In the mid-troposphere, at a model height of 6.4 km (∼450 hPa) (Fig. 17b), there are large areas showing an increase in OH by up to 20 %–25 % with the new flash-rate parameterisations, particularly in the tropics and Southern Hemisphere. The broad hemispheric differences in OH are qualitatively similar to those for O₃. With an increase in NO due to the new flash-rate parameterisation, OH increases (e.g. via the recycling of HO₂ by reaction with NO). In highly polluted air, NO₂ can be an OH sink (Lelieveld et al., 2016). Of course, transport would also influence these patterns, both horizontally and vertically.

Overall, we find that, with the new flash-rate parameterisations, there is a 13 % increase in the annual-average volume-weighted global tropospheric OH, from 10.6×10⁵ to 12.0×10⁵ molecules cm⁻³. The increase over the ocean is by 1.6×10⁵ (16.3 %) and that over land by 0.9×10⁵ molecules cm⁻³ (7.6 %). For comparison, the respective annual-average values obtained for the globe, ocean, and land from the model simulation with zero LNO_x emissions are 7.6×10⁵, 7.3×10⁵, and 8.1×10⁵ molecules cm⁻³.

The global amount can be compared with the ACCMIP multi-model mean of $\mathrm{11.1}\pm \mathrm{1.6}\times {\mathrm{10}}^{\mathrm{5}}$ molecules cm⁻³ derived by Naik et al. (2013) for the year 2000. Recent observationally based values reported by Wolfe et al. (2019) for August 2016 are $\mathrm{12.6}\pm \mathrm{2.9}\times {\mathrm{10}}^{\mathrm{5}}$ for the Northern Hemisphere and $\mathrm{8.1}\pm \mathrm{1.9}\times {\mathrm{10}}^{\mathrm{5}}$ molecules cm⁻³ for the Southern Hemisphere, and these for February 2017 are $\mathrm{8.8}\pm \mathrm{2.1}\times {\mathrm{10}}^{\mathrm{5}}$ and $\mathrm{11.4}\pm \mathrm{2.8}\times {\mathrm{10}}^{\mathrm{5}}$ molecules cm⁻³, respectively. These can be compared with the corresponding modelled values 16.9×10⁵, 7.3×10⁵, 8.0×10⁵, and 11.6×10⁵ using the PR92 scheme, and 18.7×10⁵, 8.5×10⁵, 9.2×10⁵, and 13.5×10⁵ molecules cm⁻³ using the new scheme. The LNO_x-induced increase in OH due to the new scheme adds to the model 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 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₄ against loss by tropospheric OH (${\mathit{\tau }}_{{\mathrm{CH}}_{\mathrm{4}}\mathrm{\_}\mathrm{OH}}$, defined as the division of the global total atmospheric CH₄ burden and the globally integrated CH₄ 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 burden of non-lightening-related NO_x in ACCESS-UKCA and/or a more intense photolysis.

Figure 18 The difference between the zonal total annual-mean tropospheric CO (ppbv) modelled using the new lightning flash-rate parameterisations (Run 2) and the default PR92 parameterisations (Run 1).

4.4 Carbon monoxide

There is a decrease in the modelled total tropospheric CO with the use of the new lightning parameterisations, as evident from the zonal annual-mean difference plot in Fig. 18 (this CO reduction is coupled to the OH increase via the reaction OH + CO → CO₂ + H). In the lower troposphere, the decrease is by approximately 4–6 ppbv (∼7 %) in the Southern Hemisphere and 2–4 ppbv (∼3 %) in the Northern Hemisphere. This reduction gets a little larger in the middle to upper troposphere in the tropics. Overall, the reduction in the volume-weighted global CO is 4.5 ppbv (5.6 %). Over the ocean, it is 4.7 ppbv (6.2 %) and over land it is 4.0 ppbv (4.5 %). The volume-weighted tropospheric CO obtained from the PR92 scheme is 80.3 ppbv over the globe, 75.4 ppbv over the ocean, and 89.9 ppbv over land. With the new flash-rate scheme, these values decrease by 4.5 ppbv (5.6 %), 4.7 ppbv (6.2 %), and 4.0 ppbv (4.5 %), respectively, and can be compared with the values 96.4, 91.8, and 105.4 ppbv, respectively, obtained from the model simulation with no LNO_x emissions.

Figure 19 Comparison of the modelled monthly averaged CO concentrations (ppbv) with observations at four oceanic/coastal ground stations for the year 2006. The two model runs are with the default PR92 parameterisations (Run 1) and new lightning flash-rate parameterisations (Run 2). The values of r and NMSE are also shown.

In Fig. 19, we compare the modelled monthly averaged CO with surface flask observations from the same GAW-WDCRG sites as in Fig. 12 (data from Minamitorishima were missing for 4 months, so they 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₃, the lightning changes alone do not reconcile the differences between the modelled CO and observations.

5 Model simulations with meteorological nudging

As mentioned in Sect. 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⁻¹, respectively, which are very similar to 32.92, 32.56, and 0.36 flashes s⁻¹, 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⁻¹, respectively, which are on average 5 % higher than 44.96, 35.88, and 9.08 flashes s⁻¹, 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 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 Fig. 5c and d, with the overall spatial correlation with the LIS/OTD satellite climatology (Fig. 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 %.

6 Conclusions

We have critically examined parameterisations of lightning flash rate that are based on the cloud-top-height approach. Testing of the widely used Price and Rind (1992) (PR92) parameterisations within the ACCESS-UKCA global chemistry–climate 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⁻¹ compared to the observed 34.92 flashes s⁻¹, the oceanic parameterisation severely underestimates the observed flash rate, yielding on average 0.33 flashes s⁻¹ compared to the observed 9.16 flashes s⁻¹. This leads to lightning-generated NO_x (LNO_x) being underestimated 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⁻¹ compared to the observed 34.92 flashes s⁻¹, the new parameterisation for ocean performed more accurately, giving a globally averaged flash rate of 8.84 flashes s⁻¹ compared to the observed 9.16 flashes s⁻¹. We also tested an oceanic parameterisation by Michalon et al. (1999), which gives a global oceanic average of 11.31 flashes s⁻¹. With the new parameterisations, there was an increase in global LNO_x from 4.8 to 6.6 Tg N yr⁻¹, with the new estimate comparable to 6.3±1.4 Tg N yr⁻¹ 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 NO produced per flash in the scientific literature. The model's use of 330 moles NO produced per flash is relatively 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 default PR92 parameterisations, mainly due to the change in the oceanic flash-rate component. In particular, the following impacts were observed.

There was an increase in the mid- to upper-tropospheric NO_x by as much as 40 pptv (as NO₂) in the northern tropics. There was an overall increase in the global NO_x by 8.7 pptv (15.7 %) and by 9.9 pptv (28.0 %) over the ocean. A better agreement was yielded of the modelled tropospheric NO₂ columns with the CAMS reanalysis data over the ocean.

The tropospheric O₃ burden increased by 8.5 %, from 284 to 308 Tg O₃, closer to a multi-model estimate of 337±23 Tg O₃ (Young et al., 2013), the latter supported by measurement climatology. Overall, the distribution of the 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.

There was a 13 % increase in the annual-average volume-weighted global tropospheric OH, from 10.6×10⁵ to 12.0×10⁵ molecules cm⁻³, and a decrease in the global annual-mean methane lifetime against loss by tropospheric OH by 6.7 %.

There was an overall reduction in the global CO by 4.5 ppbv (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 due to an almost-fifth-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; Iglesias-Suarez et al., 2018), primarily as a result of increases in the depth of convection. This is also the case with CAPE and precipitation-rate-based parameterisations (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 proposed here which drastically alters the flash rate over the ocean. In addition, 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.