Parameterizations for global thundercloud corona discharge distributions

Four parameterizations, distinguishing between land and ocean, have been developed to simulate global distributions of thundercloud streamer corona discharges (also known as Blue LUminous Events or BLUEs) mainly producing bluish optical emissions associated with the second positive system of N₂ accompanied by no (or hardly detectable) 777.4 nm light emission. BLUEs occur globally about 12 times less frequently (Soler et al., 2022) than lightning flashes. The four schemes are based on non-linear functions of the cloud-top height (CTH), the product of the convective available potential energy (CAPE) and total precipitation (TP), the product of CAPE and specific cloud liquid water content (CLWC), and the product of CAPE and specific cloud snow water content (CSWC). Considering that thunderstorms occur on hourly timescales, these parameterizations have been tested using hourly ERA5 data (except for CTH, not available in ERA5) for the meteorological variables considered, finding that the proposed BLUE schemes work fine and are consistent with observations by the Atmosphere–Space Interactions Monitor (ASIM). Moreover, the parameterizations have been implemented in a global chemistry–climate model that generates annual and seasonal global distributions for present-day and end of 21st century climate scenarios. Present-day predictions are in reasonable agreement with recent observations by the ASIM. Predictions for the end of the 21st century suggest BLUE occurrence rates that range between 13 % higher (∼ 3 % K⁻¹) and 52 % higher (∼ 13 % K⁻¹) than present-day average occurrences of BLUEs.

1 Introduction

The availability of regular space observations of total (intracloud and cloud-to-ground) lightning since 1995 has generated large datasets that have allowed researchers to derive annual and seasonal geographical distributions of total lightning, resulting in an annual average flash rate of ∼ 45 ± 2 flashes s⁻¹ between ±52° latitude (Christian et al., 2003; Cecil et al., 2014; Blakeslee et al., 2020). Predictions of global total lightning flash rate and geographical distribution are increasingly important, since lightning is a frequent natural hazard, considered a proxy for severe weather, a cause of large wildfires (Komarek, 1964; Pyne et al., 1998; Latham and Williams, 2001; Pérez-Invernón et al., 2021, 2022a, 2023), and a direct source of nitric oxide (NO) (Huntrieser et al., 2002; Schumann and Huntrieser, 2007; Pérez-Invernón et al., 2022b) in the troposphere that impacts the balance of important upper troposphere and lower stratosphere (UTLS) chemical species such as nitrogen dioxide (NO₂), ozone (O₃) and key oxidizing radicals such as hydroxyl (OH) and hydroperoxyl (HO₂), (Schumann and Huntrieser, 2007; Finney et al., 2016; Gordillo-Vázquez et al., 2019). Besides this, recent studies also suggest a direct production of OH and HO₂ by lightning strokes (Brune et al., 2021). All these reasons supported the need to incorporate lightning into chemistry–climate models.

The sub-grid spatial dimensions of lightning require their parameterization using different input meteorological variables and functional forms (Price and Rind, 1992; Grewe et al., 2001; Allen and Pickering, 2002; Finney et al., 2014; Luhar et al., 2021). The implementation of lightning parameterizations in different global chemistry–climate models (Tost et al., 2007; Romps et al., 2014a; Finney et al., 2014; Gordillo-Vázquez et al., 2019) have been tested against total lightning observations from low Earth orbit (Christian et al., 2003; Cecil et al., 2014; Blakeslee et al., 2020) and, very recently, also using data from geostationary satellites (Zhang et al., 2021).

Corona discharges, occurring both in the lab and in thunderclouds, are characterized by cold ionization waves known as streamers. Corona discharges are formed by numerous streamers. The electromagnetic counterpart of thundercloud corona discharges are narrow bipolar events (NBEs) (Rison et al., 2016; Soler et al., 2020). They produce bluish optical emissions (250–450 nm), leading to the adoption of the term Blue LUminous Events (BLUEs) for their optical counterpart.

While the hot and thermal air plasma in lightning stroke channels mostly excites atomic species like oxygen atoms released from thermal dissociation of O₂, leading to 777.4 nm optical emissions typical of lightning flashes, streamer corona discharges are cold non-thermal plasmas where only heavy particles are cold and electrons are very hot (up to 10 eV). Thus, corona discharges are able to activate (excite) molecular species like N₂, O₂ and H₂O by non-thermal equilibrium electron-impact collisions (Gordillo-Vázquez and Pérez-Invernón, 2021), which cause distinct bluish optical emissions mostly associated with the second positive system of N₂ radiative de-excitations.

Research results since the early 1970s indicate that, in addition to lightning, thundercloud leaderless kilometre-scale corona electrical discharges formed by hundreds of millions of streamers (Liu et al., 2019; Cooray et al., 2020; Li et al., 2021; Soler et al., 2022) are relatively common, ∼ 45 ± 2 lightning flashes s⁻¹ vs. 3.5 BLUE flashes s⁻¹, i.e. ∼ 12 times less frequent than global average number of lightning flashes in thunderclouds around the globe. In particular, recent laboratory experiments (Jenkins et al., 2021) suggest that observations during thunderstorms reporting sudden local enhancements of O₃ (Shlanta and Moore, 1972; Brandvold et al., 1996; Zahn et al., 2002; Minschwaner et al., 2008; Brune et al., 2021) and OH and HO₂ (Brune et al., 2021) could be associated with dim leaderless corona discharges (BLUEs) in storm clouds (Brune et al., 2021). These episodes suggest a probable regional atmospheric chemistry impact of thundercloud coronas, a subject which is still poorly quantified (Gordillo-Vázquez and Pérez-Invernón, 2021).

In this study, we present four parameterizations, distinguishing between land and ocean, to simulate global distributions of thundercloud corona discharges producing BLUEs (Soler et al., 2020, 2021, 2022; Li et al., 2021). The proposed storm cloud corona schemes are based on a non-linear dependence of cloud-top height (CTH) and on non-linear combinations of pairs of meteorological parameters, such as convective available potential energy (CAPE) and total precipitation (TP), CAPE and cloud liquid water content (CLWC), or CAPE and cloud snow water content (CSWC), which are all available from satellite data and atmospheric reanalyses (used to build the parameterizations) as well as in global chemistry–climate models. Some of these meteorological variables (CTH; Price and Rind, 1992; TP; Allen and Pickering, 2002; Romps et al., 2014b, 2018; and CAPE; Romps et al., 2014b, 2018) have been previously used to build different lightning parameterizations (Price and Rind, 1992; Allen and Pickering, 2002; Romps et al., 2014b, 2018). Other variables, like CSWC and CLWC, have not been used before (to the best of our knowledge). Both CSWC and CLWC can contribute to the electrification of the thundercloud since collision of graupel and ice water crystals at temperatures less than 253 K results in a negative charge transfer to the graupel that falls to lower regions of the cloud. The lighter, positive charged ice crystals stay in the higher regions of the cloud.

The physics behind CAPE (used in three of the four developed corona discharge schemes) relies on findings shown in (1) Soler et al. (2021) and (2) Husbjerg et al. (2022). Figure S12 by Soler et al. (2021) first showed that the seasonal CAPE is, in general, stronger in regions with more BLUEs. On the other hand, Fig. 5b by Husbjerg et al. (2022) showed that the CAPE associated with thunderstorms producing lightning flashes have a median value of 1000 J kg⁻¹, while thunderstorms producing BLUEs require stronger convection than needed for lightning alone. The CAPE found in the scenarios of thunderstorms that produce BLUEs range (median values) between 1280 J kg⁻¹ (slow BLUES, i.e, those buried in the thunderclouds) and 1570 J kg⁻¹ (fast BLUES, i.e., those appearing in the top of thunderclouds). As indicated by Husbjerg et al. (2022), a CAPE larger than 2000 J kg⁻¹ usually indicates deep convection. Cells generating fast BLUEs occur 25 % of the time in the region of deep convection. For cells generating only slow BLUE discharges, it is 17 %, while for regular lightning only 10 % of the events have a CAPE greater than 2000 J kg⁻¹. Therefore, there is a strong link between deep convection and the generation of BLUE discharges. Another consequence is that, in general, the occurrence of lightning is more probable since BLUEs do not require, so much, the presence of deep convection. These results lead us to consider CAPE as a plausible meteorological variable to track the occurrence of BLUEs.

Regarding the terms “fast” and “slow” mentioned above, please note that both of them underlie the scattering of the light emitted by BLUEs in thunderclouds. Since fast BLUEs are located in the cloud top, the scattering of their light emission is smaller (than that of slow BLUEs located in the bottom of the cloud). Consequently, the rise and decay times of the light curves (as seen by ASIM photometers) are faster than the rise and decay times of the light curve associated with slow BLUEs that last longer.

Most previous lightning parameterizations have been tested in a number of global atmospheric circulation models. This was done to explore how the different lightning schemes compare with available lightning observations in the present, to establish correlations with meteorological/climatic patterns and to predict possible future global lightning occurrence rates and geographical distributions in the context of a variety of future climatic scenarios. Our goal here is to proceed similarly using BLUE parameterizations since models allow for looking into the future (end of the 21st century) to reach preliminary answers to how BLUEs geographical distribution and the global occurrence rate will change in a warmer atmosphere. Therefore, the parameterizations are first tested on reanalysis data. After that, a global chemistry–climate model is used to (i) test the corona schemes against present-day climatic scenarios (both annual and seasonal), when observations are available by the Atmosphere–Space Interactions Monitor (ASIM), and to (ii) approximately predict the occurrence rate and annual geographical distribution of thundercloud coronas in future (2091–2095) climate scenarios.

The next section describes the data, observations, and modelling employed to build and test the proposed thundercloud corona parameterizations. Section 3 explains the procedures followed to develop the schemes for cloud corona discharges. Section 4 evaluates the climatological performance of the storm cloud corona parameterizations at the present day and for an end of the 21st century climate scenario. The last section of the paper presents the main conclusions.

2 Data description, observations and modelling

2.1 ECMWF ERA5 and Copernicus CLARA datasets

The European Centre for Medium-Range Weather Forecasts (ECMWF) provides the ERA5 global atmospheric reanalysis data product (Hersbach et al., 2020). The ERA5 data cover the Earth on a spatial resolution of 0.25° latitude and longitude and resolve the atmosphere using 137 levels from the surface up to a pressure of 0.01 hPa (∼ 80 km height). Single or surface level data are also available. ERA5 combines large numbers of historical observations into global estimates using advanced modelling and data assimilation systems. ERA5 provides hourly (also sub-daily and monthly) estimates of a large number of atmospheric, land and oceanic climate variables.

In order to build our thundercloud corona parameterizations, we have selected as input variables annual averages of the cloud-top height (CTH), convective available potential energy (CAPE) and total precipitation (TP), which are single (surface) level variables, as well as the annual averages of the specific cloud liquid water content (CLWC) and specific snow liquid water content (CSWC) at 450 hPa. CAPE, TP, CLWC and CSWC are taken from ERA5 (Hersbach et al., 2018a, b, 2020). CTH (not available in ERA5) is taken from the CLARA product family (Karlsson et al., 2017) of the essential climate variable (ECV) cloud properties of Copernicus (the European Union's Earth observation programme).

Note that hourly data have been averaged to obtain daily values and then averaged again into monthly values and yearly values. We have tested the BLUE parameterizations with hourly data (except for CTH for which only monthly data are available), and the result is shown in Fig. 1 (see Sect. 3) based on hourly data of the meteorological parameters used.

Figure 1 Global annual average occurrence rate (coronas s⁻¹) of nighttime thundercloud coronas observed by ASIM between 1 April 2019 and 31 March 2021 over land (red dots) and over the ocean (blue crosses) (a) vs. corona discharge schemes based on the cloud-top height (CTH) (b), CAPE × TP (c), the CAPE × specific CLWC (d) and CAPE × specific CSWC (d). The solid red (land) and blue (ocean) lines correspond to the functional forms for the nighttime corona occurrence rates of the four proposed schemes: $C_{\mathrm{F}\mathrm{1}}^{\mathrm{l},\mathrm{o}}$ = ${\mathit{\alpha }}_{\mathrm{1}}^{\mathrm{l},\mathrm{o}}$ × (CTH)${}^{{\mathit{\beta }}_{\mathrm{1}}^{\mathrm{l},\mathrm{o}}}$, $C_{\mathrm{F}\mathrm{2}}^{\mathrm{l},\mathrm{o}}$ = ${\mathit{\alpha }}_{\mathrm{2}}^{\mathrm{l},\mathrm{o}}$ × (CAPE × TP)${}^{{\mathit{\beta }}_{\mathrm{2}}^{\mathrm{l},\mathrm{o}}}$, $C_{\mathrm{F}\mathrm{3}}^{\mathrm{l},\mathrm{o}}$ = ${\mathit{\alpha }}_{\mathrm{3}}^{\mathrm{l},\mathrm{o}}$ × (CAPE × CLWC)${}^{{\mathit{\beta }}_{\mathrm{3}}^{\mathrm{l},\mathrm{o}}}$ and $C_{\mathrm{F}\mathrm{4}}^{\mathrm{l},\mathrm{o}}$ = ${\mathit{\alpha }}_{\mathrm{4}}^{\mathrm{l},\mathrm{o}}$ × (CAPE × CSWC)${}^{{\mathit{\beta }}_{\mathrm{4}}^{\mathrm{l},\mathrm{o}}}$, where the fitting parameters (listed in Sect. 3) are obtained from the best approximations of nighttime corona occurrence rate observed by ASIM over land and ocean as a function of values of CTH taken from the essential climate variable (ECV) cloud properties of Copernicus (the European Union's Earth observation programme), from CAPE × TP, CAPE × CLWC and CAPE × CSWC taken from ERA5 reanalysis. The quality of the approximations covering the 2 years (1 April 2019 to 31 March 2021) is evaluated with $R^{{\mathrm{2}}^{\mathrm{l},\mathrm{o}}}$ metrics, whose values are listed in Sect. 3. The dashed red and blue lines in panels (a)–(d) show the upper and lower limits of the fitting curves associated with the upper/lower errors in the fitting coefficients over land and ocean, respectively.

2.2 Observations used

We use global observations of nighttime thundercloud corona discharges (also known as Blue LUminous Events or BLUEs) recorded by the high sampling rate (100 kHz) photometer in the near UV (337 nm/4 nm) of the Modular Multispectral Imaging Array (MMIA) that is part of ASIM (Chanrion et al., 2019; Soler et al., 2020, 2021). ASIM is on board the International Space Station (ISS), and, due to the inclination (∼ 52°) of the ISS orbit, locations near the Equator are observed less frequently than those at higher latitudes.

The worldwide corona observations used here have been recently published (Soler et al., 2022) and span a period of 2 years of MMIA level 1 (calibrated) data from 1 April 2019 to 31 March 2021. In particular, we use the annual and seasonal averages associated with the global distribution of BLUEs obtained by the algorithm described in Soler et al. (2021), adding an extra step for filtering high-energy (and cosmic ray) candidates. The resultant distribution (used here) is named GD-2 (Soler et al., 2022).

MMIA observations of storm cloud coronas exhibit strong features in the 337 nm/4 nm photometer with negligible (barely above the noise level 0.4 µW m⁻²) signal in the 777.4 nm/5 nm photometer, which is also continuously monitored (Soler et al., 2020, 2021; Li et al., 2021).

2.3 Modelling

As an illustration of their applicability, the developed BLUE parameterizations have been incorporated into a chemistry–climate model. We use the ECHAM/MESSy Atmospheric Chemistry (EMAC) model, which is a chemistry–climate model that couples the fifth-generation European Centre HAMburg general circulation model (ECHAM5) and the second version of Modular Earth Submodel System (MESSy) to link multi-institutional computer codes, known as MESSy submodels (Jöckel et al., 2010, 2016). Such submodels are used to describe tropospheric and middle-atmosphere processes as well as their interactions with oceans, land and influences coming from anthropogenic emissions.

The thundercloud corona parameterizations described below are used in MESSy for usage within the ECHAM/MESSy Atmospheric Chemistry (EMAC) model. Cloud corona schemes are implemented as a new element of EMAC to account for atmospheric electricity phenomena that, so far, only includes lightning parameterizations. The present-day simulations without interactive chemistry are performed without lightning chemical emissions. For the climate simulations, for the end of the 21st century, we have used the lightning parameterization by Grewe et al. (2001) based on the updraught flux of mass. The variables CAPE, convective precipitation, and specific cloud liquid water content (CLWC) are calculated by the submodel CONVECT (Tiedtke, 1989; Nordeng, 1994), while the large-scale precipitation is imported from the CLOUD submodel (Roeckner et al., 2006). The total precipitation is calculated as the sum of the convective and large-scale precipitation. Following the same approach as the LNOX (lightning NO_x) submodel for the calculation of lightning (Tost et al., 2007), corona frequencies are ignored below a certain cutoff, in order to avoid introducing artefacts in the simulations. We set the corona frequency to zero if the cloud thickness is lower than 3 km.

Two sets of simulations are performed: one that covers the present-day climatic state and another one for the end of the 21st century under Representative Concentration Pathway 6.0 (RCP6.0). RCP6.0 simulations are used to estimate future occurrence and geographical patterns of thundercloud coronas. We do not have enough computational resources to run all the possible future scenarios. Therefore, we had to choose only one of them. We chose the RCP6.0 scenario that is one of the two intermediate stabilization pathways (higher medium). The EMAC simulations are performed with 720 s time step length and in the T42L90MA resolution, i.e. with a 2.8° × 2.8° quadratic Gaussian grid in latitude and longitude with 90 vertical levels starting at the surface and reaching up to the 0.01 hPa pressure level (Jöckel et al., 2016). Present-day simulations are set up by using the namelist set-up for purely dynamical simulations (referred to as the E5 set-up, no chemistry) in the mode of free-running simulations. Present-day simulations were started in January 2000 using ERA-Interim reanalysis meteorological fields (ECMWF, 2011) as initial conditions. Please note that ERA-Interim starting fields have no impact on the simulation (since we are adopting ERA5 for the meteorological parameters). The RCP6.0 simulation is set up following the simulation RC2-base-04 of Jöckel et al. (2016) and Pérez-Invernón et al. (2023). The sea surface temperatures (SSTs) and the sea-ice concentrations (SICs) are prescribed from simulations with the Hadley Centre Global Environment Model version 2 – Earth System (HadGEM2-ES) model (Collins et al., 2011; The HadGEM2 Development Team, 2011). Projected mixing ratios of the greenhouse gases and SF₆ are incorporated from Eyring et al. (2013). Anthropogenic emissions are taken from monthly values provided by Fujino et al. (2006) for the RCP6.0 scenario. We refer to Jöckel et al. (2016) for more details about the simulation set-up.

We start the 6 years for RCP6.0 in January 2090 and consider 1 year of spin-up to reach equilibrium. Section 4 below shows a comparison of the results of implementing thunderstorm corona schemes in EMAC for the years 2091 to 2095 with present-day observations (recorded by ASIM) and with simulations for the years 2000 to 2009. We are aware of the limitations of this comparison in terms of the short period of available space observations. However, it is important to highlight that this is the first global-scale continuous observation of BLUEs, and, though limited, we consider that it is worth showing such a comparison between the present day and the warmer atmosphere expected for the end of the 21st century. The present-day simulations cover 10 years. The projection (for the end of the 21st century) simulations span 5 years. It is important to note that both simulations cover more than 2 years, which are the total number of years used to develop the parameterizations. We do not simulate more years as considerable computational resources are needed. The projection simulations are initialized by using the prescribed conditions of year 2090, previously obtained from the transient climate simulations RC2-base-04 by Jöckel et al. (2016). With this approach, the climate state of 2090 (as projected for the RCP6.0 scenario) is already established, and the production of BLUEs, since based on the meteorological parameters, adapts nearly immediately. Finally, please note that despite the fact that the projection simulation spans only 5 years instead of 10 years, the mean and the standard deviation of the obtained global rate of BLUEs for the end of the 21st century are significantly larger than those obtained from the present-day simulations.

3 Thundercloud corona parameterizations

Electrical activity in thunderclouds in the form of lightning flashes has been previously correlated with CAPE (Williams et al., 1992; Pawar et al., 2012), with precipitation (Battan, 1965; Petersen and Rutledge, 1998; Allen and Pickering, 2002), and even with a linear combination of both of them (Romps et al., 2014a, 2018).

The experimental finding that thunderstorm electrification (in terms of substantial charge transfer) needs the presence of water droplets (Saunders et al., 1991) indicates that the liquid water content in thunderclouds can be a good proxy for electrical activity.

We have chosen the convective available potential energy (CAPE) as a proxy of deep convection as shown by Ukkonen and Mäkelä (2019). Also, as illustrated in Fig. S12 of the Supplement (Soler et al., 2021), the seasonal CAPE shows that, in general, there are stronger CAPE values in regions with more BLUEs. However, CAPE can also be very high in the ocean where it is not that common to find many BLUEs, so it is important to distinguish between land and ocean. Accompanying Fig. S12 in Soler et al. (2021), the seasonal CAPE vs. BLUEs per second relationships were quantified using the Pearson linear correlation coefficient (R) (varying between −1 and +1, with +1 being perfect linear correlation, 0 null linear correlation, and −1 perfect linear anti-correlation). The correlations resulted better for the zonal distributions (0.77 < R < 0.89, with DJF and MAM being the best) than for the meridional distributions (0.42 < R < 0.64, with DJF and JJA being the best).

Following the above results, Fig. 5b in Husbjerg et al. (2022) showed that by clustering the BLUE discharge dataset, cells which generate fast (close to cloud top) BLUE discharges have a median CAPE of 1390 J kg⁻¹ compared to 1128 J kg⁻¹ for cells generating only slow (deep in the cloud) BLUE discharges, further indicating that stronger cells are more likely to generate fast BLUE discharges. For comparison, the median CAPE for lightning flashes was 816 J kg⁻¹ in Husbjerg et al. (2022).

The above results lead us to use CAPE as a plausible meteorological variable to track the occurrence of BLUEs. Thus, three of the proposed parameterizations are based on CAPE multiplied by another meteorological variable like the total precipitation (previously used in lightning schemes like the one proposed in Romps et al. (2014b) or the cloud liquid (or snow) water contents (CLWC or CSWC)), whose presence in thunderclouds can contribute to cloud electrification. He et al. (2022) developed a parameterization based on the product of CAPE times the charging rate of collisions between graupel and other types of hydrometeors. Finally, we used the cloud-top height (CTH) variable to somehow test its quality, as was proposed in the popular lightning scheme presented by Price and Rind (1992).

We propose four parameterizations for corona discharges in thunderclouds based on (i) a non-linear dependence of the cloud-top height (CTH (km)) and non-linear combinations of (ii) the convective available potential energy (CAPE (J kg_air ⁻¹)) and the total precipitation (TP (m d⁻¹)), (iii) CAPE and the specific cloud liquid water content (CLWC (kg_liquid kg_moist−air ⁻¹)), and (iv) CAPE and the specific cloud snow water content (CSWC (kg_ice kg_moist−air ⁻¹)), where the CLWC and CSWC are defined as the mass of cloud liquid and snow water droplets per kilogram of the total mass of moist air.

We build the thundercloud corona schemes by using them in a 2-year data period (1 April 2019 to 31 March 2021). We consider the annual average number of thundercloud coronas per second given in each grid cell by the GD-2 distribution provided by Soler et al. (2022), as well as the annual average of the meteorological variables. The CTH, CAPE, TP, and the specific CLWC and CSWC are taken within 2° × 2° grid cells where cloud coronas take place according to GD-2.

The detailed procedure to build the corona parameterizations follows a number of steps. First, we take both 2° × 2° maps – observed global annually averaged cloud corona occurrence rate (BLUEs s⁻¹) and annually averaged of monthly averaged meteorological variables chosen as proxies – and transform them into 1D arrays. Second, we sort the meteorological dataset array following ascending values. Third, we sort the corona array by doing the same index changes we made in the second step. Fourth, we divide each array into 40 chunks, cutting each interval where the percentile distribution is a divisor of 2.5 % (neglecting the geographical position). We choose 40 chunks as an optimum. If we take more chunks, we may have too many repeated values between one chunk and the next, and if we take less values, it may not be enough to characterize the curve. Fifth, we calculate the mean of each chunk; this will result in 40 points for the predictor (meteorological data) and 40 points for the predictand (corona data). Note that the previous procedure is performed at the global scale. Sixth, we split the grid into land and ocean cells; for that purpose, we create a mask of land and ocean cells and then split each 2° × 2° cell into 400 cells of 0.1° × 0.1° size, and we calculate in the centre of that cell if it is over land (1) or over ocean (0). To finish, we assign to each 2° × 2° cell the average of the 400 cells of 0.1° × 0.1° (the average is always in the interval [0,1]). Then, if the value is larger than 0.5 it is land; otherwise, it is ocean. Finally, we approximate, separating between land (l) and ocean (o), the data with the equation y^l,o = α^l,o × $x^{{\mathit{\beta }}^{\mathrm{l},\mathrm{o}}}$ by a least squares fit, where x represents the values of the selected meteorological variables used as proxies and y^l,o corresponds to the occurrence rate (in events s⁻¹) of thundercloud coronas over land (y^l) and ocean (y^o), respectively. Note that the chosen mathematical formula (y^l,o = α^l,o × $x^{{\mathit{\beta }}^{\mathrm{l},\mathrm{o}}}$) for the parameterization of BLUEs is formally similar to previous ones for lightning flashes (Price and Rind, 1992; Michalon et al., 1999; Luhar et al., 2021; Romps et al., 2014b), including the product of two magnitudes, a constant α multiplied by a variable x raised to a power (β), where x can be a single meteorological variable (Price and Rind, 1992; Michalon et al., 1999; Luhar et al., 2021) or the product of two meteorological variables (Romps et al., 2014b). However, as indicated in the Introduction, the physics behind the corona schemes is not exactly the same as that underlying lightning parameterizations.

Figure 1a–d represent the global annual average (based on hourly data for CAPE, TP, CLWC and CSWC, and monthly data for CTH) of the nighttime corona occurrence rate (coronas s⁻¹) versus CTH, CAPE × TP, CAPE × CLWC and CAPE × CSWC, which result in the four proposed schemes: $C_{\mathrm{F}\mathrm{1}}^{\mathrm{l},\mathrm{o}}={\mathit{\alpha }}_{\mathrm{1}}^{\mathrm{l},\mathrm{o}}\times$ (CTH)${}^{{\mathit{\beta }}_{\mathrm{1}}^{\mathrm{l},\mathrm{o}}}$, $C_{\mathrm{F}\mathrm{2}}^{\mathrm{l},\mathrm{o}}={\mathit{\alpha }}_{\mathrm{2}}^{\mathrm{l},\mathrm{o}}\times$ (CAPE × TP)${}^{{\mathit{\beta }}_{\mathrm{2}}^{\mathrm{l},\mathrm{o}}}$, $C_{\mathrm{F}\mathrm{3}}^{\mathrm{l},\mathrm{o}}={\mathit{\alpha }}_{\mathrm{3}}^{\mathrm{l},\mathrm{o}}\times$ (CAPE × CLWC)${}^{{\mathit{\beta }}_{\mathrm{3}}^{\mathrm{l},\mathrm{o}}}$ and $C_{\mathrm{F}\mathrm{4}}^{\mathrm{l},\mathrm{o}}={\mathit{\alpha }}_{\mathrm{4}}^{\mathrm{l},\mathrm{o}}$ × (CAPE × CSWC)${}^{{\mathit{\beta }}_{\mathrm{4}}^{\mathrm{l},\mathrm{o}}}$, where the parameters are listed in Table 1 for land and Table 2 for ocean. These parameters are obtained from the best approximations of nighttime corona occurrence rate observed by ASIM as a function of values of CTH, CAPE × TP, CAPE × CLWC and CAPE × CSWC. The quality of the approximations covering the 2 years (1 April 2019 to 31 March 2021) is evaluated with R² metrics that results in $R^{{\mathrm{2}}^{\mathrm{l}}}$ and $R^{{\mathrm{2}}^{\mathrm{o}}}$ listed in Table 1 for land and Table 2 for ocean. In all cases, the fitting produces a $R^{{\mathrm{2}}^{\mathrm{l}}}$ score equal or above 0.91 and a $R^{{\mathrm{2}}^{\mathrm{o}}}$ score equal or above 0.94, which indicates a strong correlation between BLUEs and the meteorological variables used. The dashed red (for land) and dashed blue (for ocean) lines in Fig. 1a–d show the upper and lower limits of the fitting curve associated with the upper/lower errors in the fitting coefficients.

Table 1 Parameters for the thundercloud corona discharge parameterizations over land.

Table 2 Parameters for the thundercloud corona discharge parameterizations over ocean.

The parameterizations are implemented in EMAC using scaling factors as proposed in the lightning parameterization schemes (Tost et al., 2007). The applied scaling factors ensure a yearly occurrence rate of 3.5 BLUEs s⁻¹ during the first year of present-day simulations. Also, please note that the parameterizations are based on data of BLUEs, without including lightning data. Therefore, the developed corona schemes are independent of lightning. Regarding the sensitivity to corona discharges (BLUEs), we have calculated the spatial correlation coefficients (ranging between 0, no agreement, and 1, maximum agreement) between the simulated and the observed climatology of BLUEs (see Sect. 4.1).

4 Results

4.1 ERA5, CLARA and model simulations vs. observations

Figure 2 shows the 2-year average (1 April 2019 through 31 March 2021) nighttime geographical distribution of global corona (BLUE) electrical activity in thunderclouds according to the GD-2 distribution derived from ASIM observations (Soler et al., 2022) (panel a) and annual global predictions for BLUE occurrence rate based on hourly ERA5 data introduced in the corona (BLUE) parameterizations C_F2 (panel b) and C_F3 (panel c). Figure S1 shows the BLUE occurrence rate based on hourly ERA5 data for C_F4. Note that the colour bars have the same scale. Considering that thunderstorms occur on hourly timescales, these plots show that when ERA5 hourly data are considered for the meteorological variables, the proposed BLUE parameterizations are reasonably consistent with observations by ASIM. The global annual average BLUE rates obtained from the four adopted parameterizations using ERA5 hourly data (for CAPE, TP, CLWC and CSWC) and CLARA monthly data (for CTH) produce 3.50 BLUEs s⁻¹ for the four schemes.

Figure 2 Two-year average (1 April 2019 through 31 March 2021) nighttime geographical distribution of global corona (BLUE) electrical activity in thunderclouds according to the GD-2 distribution derived from ASIM observations (Soler et al., 2022) (a), annual global predictions for BLUE occurrence rate based on hourly ERA5 data introduced in the corona (BLUE) parameterizations C_F2 (b) and C_F3 (c). Note that the colour bars have the same scale. Considering that thunderstorms occur on hourly timescales, these plots show that when hourly ERA5 data for the meteorological variables are considered, the proposed BLUE parameterizations work fine.

Figures S2 and S3 show the 2-year average (1 April 2019 through 31 March 2021) nighttime geographical distribution of global corona (BLUE) electrical activity in thunderclouds according to the GD-2 distribution derived from ASIM observations (Soler et al., 2022) (panel a) and annual global predictions for BLUE occurrence rate based on CLARA monthly data introduced in the corona (BLUE) parameterizations C_F1, as well as monthly ERA5 data introduced in the corona (BLUE) parameterizations C_F2, C_F3 and C_F4.

Figure 3 illustrates the comparison between the 2-year average (1 April 2019 through 31 March 2021) nighttime corona GD-2 distribution observed by ASIM (panel a) and synthetic annual global average distributions obtained for the present-day climate state model simulations applying the C_F1 (panel b), C_F2 (panel c) and C_F3 (panel d) cloud corona schemes, respectively. Figure S4 shows synthetic annual global average distributions of cloud coronas obtained using C_F4. We calculated the synthetic annual global average by accounting for all time steps throughout the diurnal cycle, assuming that daytime coronas in thunderclouds causing BLUEs are equally probable as those occurring at nighttime. Both observations and simulation results show four thundercloud corona chimneys clearly distinguishable over the Pacific Ocean, the Americas, Africa/Europe and Asia/Australia. While the observed maximum in the Americas is located north of Colombia, model predictions locate it in a region partly covering the north of Peru, southern Colombia and eastern Brazil. Note that the spatial correlations (comparing observations in panel a with model simulations in panels b–d in Fig. 3) are 0.4689 (for C_F2), 0.4542 (for C_F4), 0.4226 (for C_F1) and 0.3620 (for C_F3). Although the global annual average rate for observations and predictions is the same (3.5 coronas s⁻¹), the geographical distribution of coronas in thunderclouds predicted by the C_F2 scheme resembles ASIM observations more faithfully. For instance, according to available ASIM observations, the west of North America and the south of Australia hardly exhibit thundercloud corona activity. Tornado Alley shows the highest activity in North America. These features are better reproduced by the C_F2 scheme than by the others (see Fig. 4). In addition, the predicted cloud corona discharge geographical distribution in Africa/Europe given by the C_F2 scheme follows observations better than those predicted by the C_F3 scheme or the C_F4 scheme (shown in Figs. S4 and S5). Figure 4 shows the zonal (left panels) and meridional (right panels) nighttime geographical distributions of global corona (BLUEs) electrical activity in thunderclouds according to the GD-2 distribution derived from ASIM observations (orange line) (Soler et al., 2022) and zonal/meridional distributions of the annual global chemistry–climate model predictions (using 10-year simulations, blue line) for BLUE occurrence rate according to corona parameterizations C_F1 (top panel), C_F2 (middle panel) and C_F3 (bottom panel). The predicted zonal distributions show a reasonable agreement with observations. However, only the longitudinal distribution of C_F2 (middle panel) follows the observed longitudinal distribution in terms of the relative importance of the four BLUE chimneys.

Figure 3 Two-year average (1 April 2019 through 31 March 2021) nighttime geographical distribution of global corona (BLUE) electrical activity in thunderclouds according to the GD-2 distribution derived from ASIM observations (Soler et al., 2022) (a), annual global chemistry–climate model predictions (using 10-year simulations) for BLUE occurrence rate according to corona parameterizations C_F1 (b), C_F2 (c) and according to C_F3 (d). Note that the colour bars have the same scale.

Figure 4 Zonal (latitudinal, left panels) and meridional (right panels) nighttime geographical distributions of global corona (BLUEs) electrical activity in thundercloud according to the GD-2 distribution derived from ASIM observations (orange line) (Soler et al., 2022); zonal/meridional distributions of the annual global chemistry–climate model predictions (using 10-year simulations, blue line) for BLUE occurrence rate according to corona parameterizations C_F1 (top row), C_F2 (middle row) and C_F3 (bottom row).

Thundercloud corona recordings by ASIM are conditioned by nighttime-only observations, the South Atlantic Anomaly (SAA) and the relatively reduced number of detections due to the short period of observations (2 years). These circumstances constrain the parameterizations developed in this study. However, when implemented in global atmospheric chemistry–climate models, simulations are able to predict the BLUE occurrence rate and geographical distribution in space–time regions where observations were not completely available. Regarding the SAA, please note that ASIM was not shutoff over the South Atlantic Anomaly (SAA); this was already discussed by Soler et al. (2022). Note that on March 2019 there was an update to the ASIM-MMIA cosmic ray rejection algorithm software (“ON” only over the SAA before March 2019; ON everywhere after March 2019) that could have influenced the originally obtained BLUE global distribution, the so-called GD-1 (see Soler et al., 2021, and Fig. 1 by Soler et al., 2022). That was the main reason that moved us (in Soler et al., 2022) to consider a new BLUE dataset between 1 April 2019 and 31 March 2021, generating the so-called GD-2 BLUE global distribution. The GD-2 distribution (used in this paper) already shown in Fig. 2 of Soler et al. (2022) was generated by our modified BLUE search algorithm. This included a new condition (with respect to the algorithm originally presented by Soler et al., 2021) that the 337 nm events (the associated 337 nm photometer light curve) are removed over the entire planet (not only in the SAA) when their rise times (τ_rise) are ≤ 40 µs and their total duration times (τ_total) are ≤ 150 µs (see Fig. 2 by Soler et al., 2022). Comparing Fig. 1 (GD-1) and Fig. 2 (GD-2, used in this paper) of Soler et al. (2022), it is clear that in GD-2 the radiation belt particles (RBPs) and cosmic rays (CRs) are overall removed (not only in the SAA), but we think that, most probably, GD-2 underestimates the number of BLUEs. It should be clear that GD-2 is the ASIM-observed distribution of BLUEs that we are adopting in this paper (see Fig. 2a).

Figure 5 presents a comparison of the seasonal behaviour of in-cloud coronas resulting from the GD-2 distribution of ASIM observations (left column) and the model-predicted (right column) seasonal distribution according to the C_F2 corona scheme for the present-day climate scenario. We notice that the C_F2 scheme correctly reproduces the seasonal global average occurrence rate, placing the maximum during the boreal summer (JJA) and the minimum in the boreal winter (DJF). In particular, ASIM nighttime observations indicate that the seasonal global average occurrence rates of in-cloud coronas are the following: 3.73 (SON), 2.60 (DJF), 3.75 (MAM) and 4.01 (JJA) events s⁻¹. Model predictions based on the C_F2 corona scheme result in 3.37 (SON), 3.26 (DJF), 3.39 (MAM) and 3.89 (JJA) events s⁻¹, respectively. The seasonal distributions of the C_F1, C_F3 and C_F4 schemes are shown in the Supplement (Figs. S6, S7 and S8), and all of them exhibit their maxima in JJA (3.81 for C_F1, 3.96 for C_F3 and 3.93 events s⁻¹ for C_F4).

Figure 5 Two-year average (1 April 2019 through 31 March 2021) nighttime seasonal climatology of global corona (BLUE) electrical activity in thunderclouds according to the GD-2 distribution derived from ASIM observations (Soler et al., 2022), resulting in 2.60 (DJF), 3.75 (MAM), 4.01 (JJA) and 3.73 (SON) coronas (or BLUEs) s⁻¹ (left column), and global annual average chemistry–climate model predictions (using 10-year simulations) for seasonal BLUE occurrence rate and geographical distribution according to the corona parameterization C_F2, resulting in 3.26 (DJF), 3.39 (MAM), 3.89 (JJA) and 3.37 (SON) coronas (or BLUEs) s⁻¹ (right column). Note that the colour bars have the same scale.

Model simulations based on the C_F2 scheme indicate that during the boreal winter (DJF) electrical activity in the form of thundercloud coronas is more important in the north of South America. However, this is not evident from DJF observations in South America, which exhibit disperse in-cloud corona activity probably due to limited ASIM observations in this region due to the presence of the SAA. The model-predicted in-cloud corona discharge activity in the central and southern parts of Africa, its surrounding seas and the Asia/Australia region during the boreal winter is similar to available ASIM observations, including the peak of thundercloud corona activity in the north of Australia.

The in-cloud corona observations by ASIM during the boreal summer season is in good agreement with model predictions. In particular, model simulations show that thundercloud corona discharges in Africa are mainly restricted to its central region; also, significant in-cloud corona activity is predicted in India and surrounding seas, including Indonesia and the eastern region of China. Model simulations also predict the largest number density of in-cloud coronas in North America during boreal summer, which is in agreement with ASIM seasonal observations of thundercloud corona discharges.

During the boreal autumn and spring seasons, the observed in-cloud corona maxima in America appear between the north of Colombia and central America, reaching up to the west coast of Mexico (in the boreal spring). Model predictions of in-cloud activity maxima during these seasons point to southern Colombia and the maritime region next to the east coast of Mexico (in the boreal spring). ASIM observations indicate that corona activity in central Africa is more important during the boreal spring season than during the autumn, the same is predicted by the global chemistry–climate model using the C_F2 corona scheme.

Figures 6 and 7 illustrate the end of the 21st century climate scenario and its comparison with the present day. Figure 6 shows chemistry–climate simulations that result in the geographical distribution of global annual corona (BLUE) occurrence according to different BLUE parameterizations for the end of the 21st century (global annual average of years 2091 to 2095). Panels (a), (b), (c) and (d) correspond to end of 21st century BLUE occurrence rates of 4.02 ± 0.04 coronas s⁻¹ for C_F1, 4.04 ± 0.07 coronas s⁻¹ for C_F2, 5.33 ± 0.03 coronas s⁻¹ for C_F3 and 4.02 ± 0.04 coronas s⁻¹ for C_F4, respectively.

Figure 6 Chemistry–climate simulations showing the geographical distribution for the end of the 21st century (global annual average of years 2091 to 2095) of global annual corona (BLUE) occurrence according to different BLUE parameterizations. Panels (a), (b), (c) and (d) correspond to end of 21st century global annual average BLUE occurrence rate of 4.02, 4.04, 5.33 and 4.02 coronas s⁻¹, respectively. Note that the colour bars have been scaled.

Figure 7 Chemistry–climate simulations showing the variation (percentage) between the geographical distribution for the end of the 21st century and that of the present day (global annual average of years 2000 to 2009) of global annual corona (BLUE) occurrence according to different BLUE parameterizations. Note that the colour bars have been deliberately saturated at the upper ends due to the high variability of the plotted risk.

Figure 7 shows a comparison between predicted global annual average geographical distributions of in-cloud coronas under the present-day conditions and by the end of the 21st century climate scenarios for the C_F1 (panel a), C_F2 (panel b), C_F3 (panel c) and C_F4 (panel d) thundercloud corona schemes. We find that, in general, the C_F3 corona scheme predicts a large (5.33 coronas ± 0.03 s⁻¹) global annual average occurrence rate for thundercloud coronas by the end of the 21st century when compared to present-day simulations, leading to 3.50 ± 0.01 events s⁻¹ (in agreement with ASIM nighttime observations). In particular, for the end of the 21st century, the C_F2 and C_F3 schemes predict global annual averages of 4.04 and 5.33 events s⁻¹, respectively. The standard deviation of the occurrence rates for thundercloud coronas in present-day simulations range between 0.006 (for C_F1) and 0.008 for (C_F3), so present-day occurrence rates are significantly lower than those of the end of the 21st century projections. Note that the p values from a t test comparing the periods 2001–2009 with 2091–2095 are 3.51 × 10⁻¹² (for C_F1), 2.41 × 10⁻¹⁰ (for C_F2), 1.77 × 10⁻¹⁹ (for C_F3) and 3.74 × 10⁻¹³ (for C_F4).

As said above, one of the most noticeable features in the predicted RCP6.0 future scenario for the end of the century is the 5.33 coronas ± 0.03 s⁻¹ for the scheme C_F3. The underlying reason for this can be seen in the changes of the CLWC in Fig. S9 (panel c). The liquid content at 440 hPa increases considerably over the continents because of the larger convection over the continents at the end of the 21st century. Moreover, there is less ice and snow at the 440 hPa pressure level, so Fig. 8 shows that the C_F4 scheme produces less BLUEs in many continental regions.

Figure 8 Curves show different 1-hourly BLUE flash densities in the period 1 April 2019 to 31 March 2021 for the C_F2 (blue line), C_F3 (green line) and C_F4 (orange line) BLUE parameterizations.

Finally, the end of the 21st century climate scenario results in BLUE occurrence rates between 13 % higher (∼ 3 % K⁻¹) (for the C_F1, C_F2 and C_F4 schemes) and 52 % higher (∼ 13 % K⁻¹) (for the C_F3 scheme based on CAPE × CLWC) than the present day. According to chemistry–climate simulations (see Fig. 6) for the end of the 21st century, the globally averaged temperature at the surface increases by about 4 K (Pérez-Invernón et al., 2023) by the end of the 21st century (RCP6.0 scenario, 2091–2095) compared with the present-day scenario (2000–2009).

4.2 Robustness of BLUE parameterizations on the 1-hourly timescale

We have developed the four BLUE parameterizations using ERA5 hourly (for CAPE, TP, CLWC and CSWC) and Copernicus CLARA monthly (for CTH) data since the temporal resolution in chemistry transport models is of the order of minutes. To check that the BLUE parameterizations behave reasonably well at model timescales, 1-hourly BLUE flash densities are built for the period 1 April 2019 to 31 March 2021 for the three parameterizations (C_F2, C_F3 and C_F4) as shown in Fig. 8. All the tested parameterizations exhibited approximately 95 % of cells with values less than 10⁻⁴ BLUEs km⁻² h⁻¹. It can be seen that C_F4 produces the most homogeneous occurrence of BLUEs in time and space, while C_F2 exhibits the most inhomogeneous distribution, with more cells experiencing a high rate of BLUEs. The simulated distributions of BLUE density in a chemistry–climate model can be compared with the distributions of flash density obtained by using lightning parameterizations (Finney et al., 2014, Fig. 5). The distributions of BLUEs obtained in this work show decreasing trends similar to those of flashes concerning the number of cells versus the density value. Therefore, we can conclude that the parameterizations of BLUEs obtained in this study are as applicable to chemistry–climate models as those of lightning.

5 Conclusions

Corona discharge activity in thunderclouds is found to be positively but non-linearly correlated with the following: CTH, CAPE × TP, CAPE × CLWC and CAPE × CSWC. These findings allowed us to develop four parameterizations for global chemistry–climate model simulations of in-cloud corona activity distinguishing between land and ocean.

The four corona schemes were tested against a 2-year dataset of worldwide nighttime in-cloud corona observations recorded by ASIM between 1 April 2019 and 31 March 2021. Projections to the end of the 21st century are also performed.

The global annual average BLUE rates obtained from the four adopted parameterizations using ERA5 annual-averaged data (for CAPE, TP, CLWC and CSWC) and CLARA annual-averaged data (for CTH) produce ∼ 3.50 BLUEs s⁻¹ for the four schemes. The obtained spatial correlations are 0.4689 (for C_F2), 0.4542 (for C_F4), 0.4226 (for C_F1), and 0.3620 (for C_F3). Present-day model predictions are in reasonable agreement with recent observations by ASIM.

Predictions for the end of the 21st century suggest BLUE occurrence rates range between 13 % higher (∼ 3 % K⁻¹) (for the C_F1, C_F2 and C_F4 schemes) and 52 % higher (∼ 13 % K⁻¹) (for the C_F3 scheme based on CAPE × CLWC) than present-day global average occurrence rates of BLUEs.

In-cloud corona schemes can be helpful to test global and/or regional chemical impacts of corona discharges in thunderstorms since in-cloud coronas are known to directly produce not only greenhouse gases such as ozone (O₃) and nitrous oxide (N₂O) but also oxidant species such as hydroxyl (OH) and hydroperoxyl (HO₂).