Cloud ice fraction governs lightning rate at a global scale

Lightning flash rate is strongly influenced by cloud microphysics, such as cloud ice properties, but this relationship is poorly constrained. Here we analyze 20 years of satellite-derived lightning flash rate data and cloud water data from the ERA-Interim reanalysis above continental and ocean regions at a global scale. We find a robust modified gamma function relationship between cloud ice fraction and lightning rate. Lightning rate increases initially with increasing cloud ice fraction in stratocumulus, liquid clouds. Maximum flash rates are reached at a critical cloud ice fraction value that is associated with high top, large optical thickness, deep convective clouds. Beyond the critical value, lightning rate decreases as the ice fraction increases to values representative of cirrus, ice clouds. We find consistent critical ice fraction values over continental and oceanic regions, respectively, with a lower value over the continent due to greater cloud thickness at similar cloud top height. We suggest that our findings may help improve the accuracy of lightning forecast and hazard prediction. Lightning rate varies non-linearly with cloud ice fraction and between ocean and continent regions due to differences in convection depth, according to analyses of 20-years of lightning observations and cloud water data.

L ightning is caused by the electrification of deep convective storms, which have been measured globally by satellite sensors and land-based networks 1 . Lightning strikes can result in great loss of personal property 2 and the disturbance of ecosystems through wildfire ignition 3 . The importance of lightning is self-evident, and previous exploratory efforts have been contributed to the lightning prediction based on parameterized charging 4,5 . But up to now, it is still a challenge for weather and climate models to accurately simulate the processes that lead to the electrification of clouds 6 . Cloud microphysics, especially cloud water, are widely acknowledged to have a strong impact on thunderstorms and lightning activities [7][8][9] . Hence, in order to improve the lightning forecast, it is critical to determine the relationship between the lightning flash rate and the cloud water content on a global to regional scale.
The liquid water path (LWP) and the ice water path (IWP), which are two key parameters for describing cloud phase and water content, are the integrals of liquid and ice water from the surface to the top of the atmosphere, respectively. The understanding of the relationships between cloud water and electrification of thunderstorms is one of the key points of theoretical and applied physics of clouds. Previous studies have verified the positive correlation between IWP and lightning 10,11 , which is consistent with the non-inductive charging (NIC) mechanism. The NIC mechanism is generally regarded as the most plausible explanation for the rapid initial electrification of the thunderstorm, and is based on the charge transferred by the adhesion of water particles to the surface of small ice crystals, which does not require any external electric field 12 . Indeed, the mixed-phase clouds are required to facilitate cloud electrification and generate strong electric fields required to initiate lightning, as well as provide regions where NIC occurs efficiently through collisions between graupel/hail and ice crystals with the coexistence of supercooled liquid water [13][14][15] . Lightning produced from thunderstorms is associated with charge properties. The charge transfer occurs via the NIC mechanism when ice particles collide with and rebound from riming graupel particles, and the polarity of charge transfer tends to be related to the cloud liquid water content and ambient temperature 16,17 . Under conditions that charge transfer typically occurs in or around the mixed-phase region of the thunderstorm, the graupel particles appear to acquire negative charge, whereas cloud ice particles acquire positive charge 18 . Lightning is initiated when and where the electric field has a sufficiently large magnitude between oppositely charged regions. The generation of lightning is affected by the combined effects of various phase particles in the cloud, therefore, not only IWP but also LWP has an impact on lightning.
To quantify the influence of cloud water path (CWP, including LWP and IWP) on lightning, we take into account the ratio of IWP to (IWP+LWP), referred to as ice fraction 19,20 , in this study. This work intends to establish a climatological relationship between the observed lightning flash rate and the ice fraction from long-term (during 1995-2014) measurement datasets. Due to the difference of atmospheric electricity caused by the cloud thermodynamic and microphysical characteristics between continental and oceanic regions 21 , we separate them to analyze the relationship between ice fraction and lightning on a global scale. Although the properties of clouds over the continent and ocean are different, the relationships between ice fraction and lightning are consistent, and both conform to the modified gamma function. Moreover, it is noted that the characteristics of lightning activities are highly related to the dynamic and microphysical processes in the clouds 22 . We further study the characteristics of different cloud types corresponding to different positions on the relationship curve, so the physical implications of the curve are revealed. Then, a conceptual model is proposed to illustrate the effect of cloud ice fraction on lightning activity. In addition to the climatological relation mechanism between ice fraction and lightning on the global scale, this uniform law also exists in regions with various climate properties.
Through our analysis, we define a relationship between ice fraction and lightning that holds at regional and global scales, which is a modified gamma function. The discovery of climatological linkage between cloud and lightning will promote our understanding of the physical mechanism of lightning and the prediction of its occurrence.

Results and discussion
Global scale analysis. The global spatial distribution of the annual average lightning flash rate from 1995 to 2014 is presented in Fig. 1 using lightning data from Optical Transient Detector (OTD) and Lightning Imaging Sensor (LIS) instruments, which are onboard the MicroLab-1 and the Tropical Rainfall Measuring Mission (TRMM) satellites, respectively. Lightning activities over the continent occur more frequently than those over the ocean. Although part of the reasons has been explained in the previous study [23][24][25] , the mechanism underlying this difference has not been well confirmed. The average lightning flash rate between 38°S and 38°N over continental and oceanic regions is 12.13 (fl km −2 year −1 ) and 1.26 (fl km −2 year −1 ), respectively. This result is related to the difference in convective vitality between continental clouds and oceanic clouds 26 .
The cloud water path can modulate lightning activity by affecting the microphysical electrification processes. From the perspective of the time duration and space coverage, cloud water data from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis dataset are used in this study, which were evaluated against the TRMM observation data over the ocean (Supplementary Fig. 1, Supplementary Discussion). As shown in Fig. 2, the response curve of the lightning flash rate to the cloud ice fraction first increases and then decreases, with a peak turning point at which the ice fraction is 0.22 and 0.34 over the continent and ocean, respectively. That is to say, there exists a specific ice fraction value in the mixed-phase cloud, which maximizes the occurrence of lightning activity. When the ice fraction of the clouds decreases (increases) and they tend to be the liquid (ice) clouds, the lightning flash rates will decrease. However, this threshold value varies with the cloud characteristics in the different regions. Compared with the continental clouds, the ice fraction of the oceanic clouds is larger when the lightning flash rate reaches its maximum. Bang and Zipser 24,25 noted that oceanic precipitation features with lightning tend to be larger and more mature than those over land, which is consistent with our findings.
Based on Eq. (1) (see "Methods"), the corresponding curves over the continental and oceanic regions are well fitted with R 2 ¼ 0:99 and R 2 ¼ 0:95, respectively (the solid lines in Fig. 2). The coefficients are empirical values via long-term data statistics. According to the coefficients of the fitting curves, the value of cloud ice fraction (x c ) and lightning flash rate (L max ) at the critical threshold can be deduced from the Eqs. (2a) and (2b) (see "Methods"), respectively, in the continental and oceanic regions. For the clouds over the continental regions, x c ¼ 0:22 and L max ¼ 0:044 (fl km −2 d −1 ). For the clouds over the oceanic regions, x c ¼ 0:34 and L max ¼ 0:008 (fl km −2 d −1 ). The difference of x c in the continental and oceanic regions indicates that the oceanic convective clouds need a larger proportion of IWP to generate strong lightning activities, which is consistent with previous findings 24, 25 . Here, the macro characteristics of clouds, which include cloud top pressure (CTP) and cloud optical thickness (COT), are compared between the continent and the ocean during the occurrence of strong lightning density ( Supplementary  Fig. 2). The results show that when the lightning flash rate near its peak (the threshold values are presented in Supplementary Table 2), the difference in CTP between the continent and ocean are similar but the COT of the continent is higher than that of the ocean, especially focusing on the difference in maximum COT. In other words, when a strong lightning density occurs, the top of the thunderstorms over the continent and ocean can achieve a similar height, but the cloud thickness of the continent is much higher with a lower cloud base height. Ice fraction is smaller when the thunderstorm has a fixed cloud top height but a lower cloud base height. In addition, it suggests that given similar cloud top heights, the indication of greater COT over the continent points to a deeper mixed-phase region, which is more favorable for extensive NIC. Therefore, the values of x c are influenced by the macro characteristics of the deep convective thunderstorm clouds in each region. In addition, the values of L max indicate that the lightning flash rate produced by deep convective clouds is 5.5 times higher over the continent than that over the ocean. This difference is smaller than the annual average lightning difference between the continent and ocean, which is~9.6 times (Supplementary Table 1). To identify why it differs, the probability distribution functions (PDFs) of ice fraction over continent and ocean are discussed, and the result is shown in Supplementary Figure 3. The peak of continental ice fraction distribution is close to its x c (0.22), but the oceanic peak is lower than its x c (0.34). It is apparent that, on average, thunderstorms with the highest lightning flash rate occur more frequently over the continent than that over the ocean. Stratocumulus with a smaller ice fraction than the oceanic x c occur more frequently over the ocean, which reduces the oceanic average lightning flash rate.
To further understand the response of lightning activities to various ice fractions, we select three sets of continental data based on the threshold values of lightning rate and ice fraction. The specific grouping conditions are listed in Table 1. Group-2 indicates the data with the highest incidence of lightning when the value of ice fraction is equal to or near x c . Group-1 (Group-3) represents the datasets with smaller lightning rate and smaller (larger) ice fraction on the left (right) side of the curve maximum. The Group-1, 2, and 3 represent the thunderstorms with mostly liquid clouds, the thunderstorms with strong lightning density, and the thunderstorms with mostly ice clouds, respectively. Figure 3 displays the PDFs of LWP, IWP, convective available potential energy (CAPE), and the content of total column water vapor (TCWV) of each group (see "Methods" for the datasets information). From the perspective of CWP ( Fig. 3a and b), the peaks of Group-2 occur at similar values to the LWP of Group-1 and the IWP of Group-3, but the probability densities at the peak CWP values of Group-2 are more prominent than those of Group-1 and Group-3. This suggests that thunderstorms with high lightning rates have higher concentrations of ice (liquid) water than ice (liquid) clouds with low lightning rate, and Group-2 storms have more balanced ice fractions, which is more conducive to NIC. In addition, Fig. 3c and d point out the descending order of the peaks of CAPE and TCWV is Group-2, Group-1, and Group-3. The occurrence of thunderstorms with high lightning rates needs vertical atmospheric instability with high values of CAPE 25,27,28 and sufficient water vapors [29][30][31] . Noting the peaks of Group-1 and Group-3, under the same lightning rate circumstances, ice clouds do not require larger values of CAPE and TCWV than liquid clouds. The abundant ice crystals make lightning more frequent in ice clouds than that in liquid clouds at the same values of CAPE 25 and TCWV 32 . Liquid clouds require greater CAPE and TCWV to achieve the same charging when compared with the ice clouds. This might be the reason why in Fig. 2 the lightning rate increases faster with small ice fraction (0-0.22) and then decreases slowly with large ice fraction (0.22-1). Moreover, from the perspective of storm morphology, developing cells tend to electrify quickly once ice forms (rapid increase in lightning with small ice fraction), but lingering charged ice particles aloft in a decaying cell maintain electrification and continue to participate in lightning discharges  when dynamic updraft support wanes and warm precipitation processes diminish (reduced liquid water content, increasing ice fraction, and slow decrease in lightning). Furthermore, based on cloud-type definitions of the International Satellite Cloud Climatology Project (ISCCP), cloud types and macro characteristics of the three groups are summarized. As presented in Fig. 4, the clouds of Group-2 are featured by high cloud top and large optical thickness (e.g., deep convections). However, the clouds of Group-1 and Group-3 are dominated by the stratocumulus and cirrus, respectively. The macrophysical characteristics of the three groups are consistent with the properties of the parameters in Fig. 3.
Based on the above understanding, a conceptual model is shown in Fig. 5. This model suggests that the lightning activities first increase when cloud ice fraction increases from 0 to x c and then decrease with the increase of cloud ice fraction from x c to 1. When the cloud ice fraction is smaller (greater) than x c , the clouds tend to be liquid (ice) clouds, producing less lightning. However, when cloud ice fraction is equal to or near x c , the clouds tend to be thunderstorms with sufficient ice crystals, graupels, hails, and supercooled liquid cloud droplets. These cloud hydrometeors provide the areas where the NIC occurs efficiently through collisions between them, then produce the most lightning activities. In other words, small ice fractions correspond to shallow and warm clouds, and large ice fractions correspond to cirrus, cirrostratus, and cirrocumulus. These liquid or ice clouds do not necessarily generate charges, though large ice fraction conditions still may be charged as in decaying convections. Lightning flash rates are highest where there are deep convections, which are known to generate a lot of lightning. Therefore, creatively proposing and accurately calculating the value of x c contributes to our understanding of the lightning processes in thunderstorm clouds.
Adaptability of the function on the regional scale. To determine whether this global scale model also holds at regional scales, we select four regions where lightning strikes frequently occur. The four regions (North America, Asia, South America, and South Africa) we selected are marked with red rectangles in Fig. 1. Due to the conspicuousness of lightning over the continent, we only extracted the data over land when analyzing regional characteristics. In line with the global scale research, the same data processing methods are used to study the relationships between lightning flash rate and cloud ice fraction in the four regions.
As shown in Fig. 6, the relationships between lightning flash rate and ice fraction in the four regions are consistently shaped like the modified gamma curve. We fit their correlations using Eq. (1) (see "Methods"), and the values of R 2 displayed in Fig. 6 are all >0.9 in the four regions. This result indicates that the empirical model is valid on both global and regional scales. In addition, x c and L max calculated by the Eqs. (2a) and (2b) (see "Methods") are shown in Fig. 6. The values of x c in the four regions are all~0.25, which is slightly higher than the global average due to the fitting deviation. The comparison of the L max values among the four regions is consistent with the level of average lightning flash rate in each region (Fig. 1). Figure 1 shows that the average lightning flash rate in South Africa is the highest among the four regions, thus the value of L max is highest in South Africa.
In addition, we select a region with a weak lightning rate located in North Africa, which is marked with a blue rectangle in Fig. 1. The same data processing method is used as the global scale research. As shown in Supplementary Fig. 4, the relationship between lightning flash rate and ice fraction in North Africa is consistently fitted by the modified gamma function. The value of x c is much higher than the global analysis, which is more similar to the features over the ocean. In other words, the strong lightning regions may have dominated the signal in the global analysis.
Applicability of regions with different cloud properties in Asia. Cloud properties are affected by multiple elements, such as  thermodynamics, microphysics and aerosols. As one of the impact factors, the atmospheric aerosol is shown to play an important role in suppressing or invigorating the development of convective cloud 33 , thus affecting the occurrence of lightning 34,35 . Hypotheses have been proposed for decades about the role of aerosols as activated cloud condensation nuclei on delaying the warm rain process, enhancing the mixed-phase process and invigorating deep convective cloud vertical development 36 . However, different aerosol types have varying effects on convective clouds 37 .
We select three continental regions, including Southeast Asia (SEA), South Asia (SA) and West Asia (WA), dominated by different aerosol types and cloud properties in Asia, to see whether or not this relationship applies to them. The locations of these three regions are outlined by red rectangles in Supplementary Fig. 5. We only extracted the data over land due to the conspicuousness of lightning over the continent.
As shown in Supplementary Fig. 5, high loadings of smoke and dust aerosols are prevalent in SEA and WA, respectively. In addition, SA is an area dominated by both dust and smoke aerosols. Biomass burning is known to be an important contributor to smoke aerosols in SA and SEA 38,39 . WA and part of SA are considered as the important dust aerosols source areas in the world 40,41 . Aerosols affect clouds and lightning activities via both radiative and microphysical effects 34 . The differences of aerosol types and cloud properties in the three regions can help us to study the relationship between lightning flash rate and cloud ice fraction from a more detailed regional view. Finally, we can find out whether the characteristics of each region are consistent with the results obtained on the global scale.
Regional differences in lightning activity and cloud ice fraction are presented in Supplementary Fig. 6. Wang et al. 34 demonstrated the roles played by the different aerosol types and the relative importance of dynamics-thermodynamics versus aerosols in affecting lightning. Among the three regions, the lightning flash rate is the lowest in WA and the highest in SEA ( Supplementary Fig. 6a), which is consistent with the result analyzed in the dust and smoke regions 34 . The drier air in the dusty region tends to evaporate cloud droplets and suppresses the development of convective clouds, thus inhibiting lightning activity in combination with the dust radiative effect which causes surface cooling and stabilizes the atmosphere. On the contrary, in the smoke dominant regions, the air might be moist relative to the dust areas which is less susceptible to the suppression of clouds. The absorbing smoke aerosols also heat the lower atmosphere and lead to unstable stratification which invigorates convection and lightning activity 34 . Therefore, the lightning rate shown in Supplementary Fig. 6a is the highest in SEA and the lowest in WA. Dust aerosols serve as ice nuclei in mixed-phase clouds, and thereby invigorate ice production and alter the properties of clouds 42 . Therefore, cloud ice fraction is highest in WA and lowest in SEA. Beyond these differences in lightning activities and clouds, it is worth studying whether there are commonalities among the three regions.
The same data processing methods are used to study the relationships between lightning flash rate and cloud ice fraction in the three regions (WA, SA, and SEA), except that the data are averaged every 100-point. The scatter plots are displayed in Supplementary Fig. 7. Clearly, the correlations of lightning flash rate and ice fraction in the three regions are also shaped like the modified gamma curve. Equation (1) (see "Methods") is then used to fit their correlations, showing the R 2 = 0.87, 0.93, and 0.96 in WA, SA, and SEA, respectively. The larger correlation coefficients further illustrate the applicability of this function on the regional scale, regardless of the characteristics of aerosol, cloud and lightning in each region.
In addition, x c and L max are calculated using the Eqs. (2a) and (2b) (see "Methods"), and the results are presented in Supplementary Fig. 7. The value of x c is the highest in WA, indicating the noticeable ice production associated with the role of dust aerosols as ice nuclei. Furthermore, the macro characteristics of clouds in the three regions during the occurrence of strong lightning density are compared in Supplementary Fig. 8. The results indicate that when the lightning flash rate near its peak (the threshold values are presented in Supplementary Table 3), the difference in average CTP among the three regions are similar, which are centered similarly near 500 hPa, but the COT of SA is the highest while the COT of WA is the lowest. The result is somewhat consistent with the discussions about the comparison between continental and oceanic thunderstorms. Therefore, the values of x c are affected by the macro characteristics of the deep convective thunderstorm clouds in each region. The value of L max reflects the ability of deep convective clouds to generate lightning. The difference in L max among the three regions agrees with the averaged lightning flash rate shown in Supplementary Fig. 6a.
It seems that the aerosol characteristics and ice fractions are relatively similar in WA and SA, yet their COTs and lightning rates are different. We should be aware that not only aerosols but also dynamic-thermodynamic conditions (e.g., CAPE, wind shear, divergence, convective organization, temperature, pressure, humidity, etc.) can affect clouds and lightning activities 34,43 . More aspects that constrain lightning need to be discussed in future research to have a better grasp of the empirical modified gamma function in a specific region.
Conclusions. Using 20-year datasets, this study presents the relationship between the lightning flash rate and the ice fraction, and explains the intrinsic physical mechanism of their relationship. The response curve of the lightning rate to cloud ice fraction first increases and then decreases. Lightning rate increases initially with increasing cloud ice fraction in shallow and warm clouds (e.g., stratocumulus). Maximum flash rates are reached at a critical cloud ice fraction value that is associated with deep convective clouds. Beyond the critical value, lightning rate decreases as the ice fraction increases to values representative of ice clouds (e.g., cirrus, cirrostratus, and cirrocumulus). By conducting the theoretical analysis, the results indicate a modified gamma function (LðxÞ ¼ ax α e Àbx β ; 0 ≤ x ≤ 1; see "Methods" for the specific parameter information), which reflects the relationship between the lightning rate and the ice fraction. It fits well in the continental and oceanic regions covered by different lightning intensities and cloud properties. In addition, the parameters of the function can readily interpret the physical implications of the relationship between the intensity of lightning and the macro and micro characteristics of clouds in each region. Based on the robust relationship between lightning flash rate and ice fraction, this study can contribute to the recognition and prediction of lightning, and will assist the communities or agencies with the lightning disaster prevention and mitigation. Future work should concentrate on the precise regional parameters that influence the relationship, and provide a better knowledge of the regional parameters for the empirical investigations and prognostic use.

Methods
Lightning data. The LIS/OTD lightning datasets observed by the space-borne OTD and LIS instruments are used to obtain the global distribution and variation of total lightning (cloud-to-cloud, intra-cloud, and cloud-to-ground) between 38°S and 38°N 44 . The low-resolution monthly time series (LRMTS) product during the study period from July 1995 to January 2014, which provide the daily-lightning flash rate (fl km −2 d −1 ) per month with a 2.5°× 2.5°gridded climatology dataset, are used in this study.
Cloud data. The data of LWP and IWP are taken from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis dataset, corresponding to the reanalysis products of total column cloud liquid water and total column cloud ice water, respectively 45 . Here, the monthly average data during the study period with a grid of 2.5°× 2.5°have been used.
The monthly average ISCCP H-series cloud product with a spatial resolution of 10 km is used during the study period, which includes the cloud-type, CTP and COT 46 . In this study, we interpolate the ISCCP data into a grid of 2.5°× 2.5°, which matches with the grid of lightning data spatially.
Atmospheric stability and water vapor data. The data of CAPE and TCWV are taken from the ECMWF ERA-Interim reanalysis dataset during the study period 45 , and the same preprocessing method as cloud water data are used.
Aerosol data. The Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) provides monthly average data of dust, organic carbon (OC), black carbon (BC), and total extinction aerosol optical depth (AOD) at a spatial resolution of 0.5°× 0.625°4 7 . In this work, we interpolate these data to match with the grid of lightning data spatially.
Determining the cloud ice fraction-lightning relationship. In order to study the relationship between lightning and cloud ice fraction, lightning data are arranged in ascending order of ice fraction, and then the arranged data are averaged every 1000-point.
Given the shape of the relation curve between lightning flash rate and cloud ice fraction, an ideal function is proposed. The function shows well fittings with the measurement data in both continental and oceanic regions, which is a modified gamma function, written as follows: where, the LðxÞ is the lightning flash rate at the ice fraction x, and a, α, b, and β are positive constants. Then, we let the derivative of Eq. (1) be 0, and the critical point with the maximal lightning flash rate can be obtained. The value of cloud ice fraction (x c ) and lightning flash rate (L max ) at the critical point are expressed as follows: where ε ¼ α β is a positive constant.

Code availability
Computer code used in the analysis is available on request from the corresponding author.