Natural Sources of Ionization and Their Impact on Atmospheric Electricity

We present a study of atmospheric electricity using the chemistry‐climate model SOCOL considering ionization by solar energetic particles during an extreme solar proton event (SPE), galactic cosmic rays (GCR), and terrestrial radon (Rn‐222). We calculate the global distribution of the atmospheric conductivity and fair‐weather downward current density (Jz) using atmospheric ionization rates from all sources. We found that Jz is enhanced (by more than 3.5 pA/m2) in radon source and polar regions. Contribution of Rn‐222 is essential at middle and low latitudes/altitudes where GCR‐induced air conductivity is reduced. The model results are in good agreement with the available observations. We also studied the effects of an extreme SPE, corresponding to the 774 AD event, on the atmospheric electricity and found that it would lead to a large increase of Jz on a global scale. The magnitude of the effects depends on location and can exceed background value more than 30 times over the high latitudes (a conservative upper bound). Such an assessment has been performed for the first time.


Introduction
The main sources of atmospheric ionization are continuous solar UV radiation, highly energetic galactic cosmic rays (GCR), and energetic electron precipitation (EEP) as well as sporadic solar proton events (SPEs). The ionization produced by these sources affects the chemistry and dynamics of the atmosphere via enhanced formation of active hydrogen (HOx) and nitrogen (NOx) oxides, which further participate in catalytic ozone loss cycles (e.g., Mironova et al., 2015). These chemical processes thereby affect the entire climate system (Rozanov, 2018). The lower atmosphere over land is mostly ionized by natural ground sources (Eisenbud & Gesell, 1997), since energetic particles, except for high-energy GCR, do not reach these low altitudes, especially at middle and low latitudes. An important source of ionization near the ground is the radioactive isotope radon-222 (Pearson & Moses, 1966), which is a product of the uranium-238 decay. Rn-222 decays by emitting -particles, which have a very short range (several cm), and thus ionize air very locally. Continuous monitoring of Rn-222 since the late 1950s shows the importance of Rn-222 for human life due to a hazard for health. Rn-222 is injected into the atmosphere via exhalation because of changes in pressure and air temperature as well as evaporation of soil moisture from the surface where Rn-222 is decomposed (Israel, 1970(Israel, , 1973Mohnen, 1977). The half-life time of Rn-222 is 3.826 days (Pearson & Moses, 1966), which is too short to reach the middle atmosphere, but it is sufficient to affect near-ground ionization (until to 3 km) and the global electric circuit (GEC). Rn-222 has a diurnal cycle with a maximum concentration during the night and a seasonal variation as it strongly depends on dry convective activity, soil freezing that reduces Rn-222 emissions and moist convection in oceanic areas (Jacob & Prather, 1990). Chen et al. (2016) observed the highest radon concentration at 6 a.m. and lowest at 4 p.m. and found the lowest concentration of Rn-222 during a spring season using measurements at the SMEAR II station in southern Finland. In the outdoor air, Rn-222 concentration is estimated to be between 0.1 and 10 Bq m −3 with the averaged value equal to 3 Bq m −3 over the land (UNSCEAR, 1982). However, in some countries, for example, Sweden, higher values of up to 60 Bq m −3 can be observed over land (UNSCEAR, 1982). In the seas, the concentration of Rn-222 is assumed to be a few orders of magnitude smaller and estimated to be ≈1 Bq m −3 on average over the oceans with the exhalation rate of about 10 −9 Bq·cm −2 ·s −1 (UNSCEAR, 1982). Considering the importance of Rn-222 for atmospheric electricity, it should be noted that released Rn-222 creates a so-called reversible electrode effect or changing the sign of the electric field as mentioned in Shuleikin (2017).
The ionization impact on the atmospheric electricity was partly studied earlier by Rycroft et al. (2008). The ionization rate induced by GCR, EEP, SPE, and Rn-222 depends on the altitude and latitude differently for each source (e.g., Mironova et al., 2015;Roffman, 1972). The energetic-electron-driven ionization increases with height to the maximum in mesosphere/thermosphere, while for GCR, the maximum is at about 18 km and then decreases with heights (e.g., Bazilevskaya et al., 2008;Phillips et al., 2016). The Rn-222-induced ionization maximizes in the boundary layer irrespective of the latitude but having a geographical pattern related to the radon emission. Since the electric current in a circuit is largely defined (for a given voltage) by a region with the lowest conductivity, the ionization at lower heights is important for the GEC of the Earth, that controls the fair-weather downward current density (J z ) and ground-level electric fields and can further impact cloud microphysics and close-to-the-ground meteorology (e.g., Burns et al., 2007;Lam & Tinsley, 2016). While the influence of energetic particles on the ozone layer and the climate was studied earlier (e.g., Rozanov et al., 2012), the role of Rn-222 has been poorly analyzed yet. Baumgaertner et al. (2013) showed that ionization induced by Rn-222 resides primarily in the planetary boundary layer within 3 km of the surface layer (1,000 hPa) reaching 28 ions·cm −3 ·s −1 . The main objective of this study is to properly model the ionization rates produced by Rn-222 and to estimate the contribution of all ionization sources to the atmospheric conductivity and the fair-weather downward current density (J z ), which characterize the GEC state. Special attention is paid to a comparison of the obtained results with the available observations and other model results to estimate the model performance in the GEC calculations. Despite having some global-scale measurements of the ionospheric potential (IP) (e.g., Markson, 2007), measurements of the current density are generally not available. The GEC measurements from surface potential gradient measurements are difficult because of local influences such as aerosol variations (Nicoll et al., 2019). It is the main reason why we can rely mostly on the model intercomparison for the model evaluation.
The possible effect of an extreme solar proton event on the GEC was also studied. For a better understanding of this process, it is necessary to separate the effects of SPE and the effects of interplanetary and geomagnetic disturbances (such as geomagnetic storms) (e.g., Reiter, 1992). The effects of these events on atmospheric electricity have not been studied in detail previously and are based on separate empirical studies (e.g., Rycroft et al., 2000). Some models can explain the modulation of the global atmospheric electrical circuit by SPEs (e.g., Farrell & Desch, 2002), but these models do not distinguish SPE impact from the influence of CME (coronal mass ejection) and Forbush decrease. Tacza et al. (2018) reported that usually, SPE enhances the electrical conductivity by about 10%. This conclusion is based on the atmospheric electric fair-weather field measurements at ground level and is supported qualitatively by models results (Cobb, 1967;Markson, 1978). However, Willett (1979) disagreed with the conclusions presented in Markson (1978) showing a smaller effect of SPE on the air conductivity and the GEC, but also concluded that differences in the results may lie with different approaches used to represent thunderstorm clouds since neither models used was able to represent such clouds properly due to lack of knowledge of physics behind them. Also, it should be mentioned that these results were obtained for a relatively weak SPE, which inspires studies of much more extreme event similar to SPE of 774 AD (Sukhodolov et al., 2017). Tacza et al. (2018) analyzed variations of the fair-weather electric field considering two types of events: solar flares with X-ray fluxes higher than GOES-Class C1 without any solar energetic particles (SEPs), and strong SPEs with fluxes of energetic (>100 MeV) protons. Here, we develop this approach by analyzing extreme solar events and analyzing their impact on the GEC.
Sometimes, sporadic events of strong gamma-ray bursts originating from our galaxy or beyond may affect the ionosphere (e.g., Fishman & Inan, 1988;Nickolaenko et al., 2011;Williams & Mareev, 2014), but these events are left beyond the current work, which is focused on the effect of ionization of the lower atmosphere on the atmospheric vertical current.

Model Description and Verification
The chemistry-climate model (CCM) SOCOL (SOlar Climate Ozone Links) was developed in the Physical-Meteorological Observatory and World Radiation Center (PMOD/WRC) in cooperation with the Institute for Atmospheric and Climate Science, ETH Zurich, Switzerland. The CCM SOCOL v2 was designed to study the impact of different external factors such as GCR and SPEs on the Earth climate system and the ozone layer. It couples the middle-atmosphere general circulation model MA-ECHAM4 and the chemistry-transport module (CTM) MEZON (Model for the Evaluation of Ozone Trends). The MA-ECHAM4 model was developed at the Max Planck Institute for Meteorology in Hamburg, Germany. It has a horizontal resolution of about 3.75 • × 3.75 • and 39 vertical layers that cover the heights ranging from the Earth's surface to about 80 km. The dynamical core of the CCM SOCOL v2 is based on the semi-implicit time-stepping scheme and exploits a time resolution of 15 min. However, the chemistry, transport, and full radiation calculations in the model are performed every two model hours. The CTM MEZON shares the MA-ECHAM4 horizontal and vertical grids and treats 41 different atmospheric chemical species (Schraner et al., 2008).
Atmospheric ionization rates utilized in CCM SOCOL v2 were calculated using the cosmic-ray-induced ionization model CRAC:CRII (Usoskin & Kovaltsov, 2006;Usoskin et al., 2010), which can compute the ionization rates, due to GCR and SEPs, in the atmosphere at any spatial-temporal location for the known energy spectrum of incoming particles. The model is based on the Monte Carlo codes CORSIKA (v.6.617 August 2007) and FLUKA (v.2006(v. .3b March 2007 to simulate lower energy nuclear processes. The model output is the number of ion pairs produced per one gram of air per second at a given atmospheric depth and geomagnetic latitude (Usoskin et al., 2011).
The basic CCM SOCOL v2 version considers only ionization from GCR, energetic particle precipitation (EPP), and SEPs, while the Rn-222 has not been taken into account. To fill this gap and to consider sources, we developed a module simulating Rn-222 emission (Schery & Wasiolek, 1998), transport and decay. The latter was treated using the 3.8-day half-life time. The geographical distribution of the emission was prescribed according to Kazil et al. (2010) after conversion to units of g/(m 2 s), as required for the transport module. The convective transport was treated using the method described by Jacob and Prather (1990), which suggests that dry convection should be applied for regions with instability and moist convection-for regions where convective clouds are present. The Rn-222 concentration in the atmosphere depends on transport and is not solely defined by its emissions. Rn-222 in the atmosphere can accumulate in certain locations depending on the preferential wind patterns and decay. In this study, we performed numerical experiments corresponding to the period time of 2004-2005. In the case of excess of the signal over the noise, ensemble calculations can be omitted. The obtained distribution of the radon concentration is illustrated in Figure 1a. There are more sources of atmospheric radon over the equatorial regions in South Asia, south part of Northern America, South America, Africa, and Australia, due to geological features. Unfortunately, there is no information about Rn-222 sources over oceans, and accordingly, it was neglected here.
While equatorial regions have more intensive sources of radon emission (see Figure 1a), the maximum of Rn-222 concentration appears between 20 • N and 40 • N latitudes due to the advective transport pattern in the northern hemisphere. Figure 1b shows Rn-222-induced ionization rates (IR) calculated from the Rn-222 mass-mixing ratio (C Rn222 ) applying the unit conversion between mBq/(m 2 s) and g/(m 2 s), so that 1 g of Rn-222 corresponds to 5.69 · 10 15 Bq, and the air density ( , kg/m 3 ): The radon-related ionization rate reaches a maximum near the Earth's surface, where the Rn-222 concentration is enhanced over the middle-and low-latitude areas. During the boreal summer, high radon-related ionization rates are located primarily in the northern hemisphere between 40 • N and 70 • N, caused by a strong and stable Rn-222 emission. A reason for the high level of Rn-222 during summer might be related to variations of the soil dryness, which allows the gas to travel more easily. High radon-related ionization rates often occur during boreal winter in regions with the lowest temperature, for example, northern Eurasia. A reason for the maxima in boreal winter is probably related to the prevalence of low-level temperature inversions in winter, which can also maintain high concentrations of smoke/aerosol. It also may be used in the simulations of aerosol size distribution, cloud properties, and the related climate effect (Zhang et al., 2011). The obtained results were compared with the results of Rn-222 simulations using Community Earth System Model (CESM1) Whole Atmosphere Community Climate Model (WACCM) (Lucas, 2017;Lucas et al., 2015). Our simulations with the CCM SOCOLv2 yield higher ionization rates (≈12 Bq/cm 3 vs. ≈9 Bq/cm 3 by CESM1). In general, the distribution of ionization rates in the SOCOLv2 is similar and agrees with the results of the CESM1 (WACCM) (Baumgaertner et al., 2013). Figure 2a shows GCR-induced ionization rates utilized by the CCM SOCOLv2. For details on the GCR ionization rates, see section 3. It is worth noting that in boundary layer, the Rn-222-induced ionization rates are higher than GCR-induced ionization rates by a factor of 2-6 depending on the region. In particular, the ionization rate may reach a value of 12 ion pairs cm −3 s −1 in Rn-222 active regions, while it is less then 1.7 ion pairs cm −3 s −1 for GCR even in polar regions. Thus, radon dominates in ionization near the surface above ground, but its effect quickly fades with altitude. As mentioned above, the calculated ionization rates resulted from Rn-222 and GCR distributions are used to calculate the air conductivity in the atmosphere. Calculation of air conductivity ( ) includes ion concentration (n) and ion mobility ( ± ), with e being the elementary charge: = n · e · ( + + − ). (2) The computed conductivity was compared with the results of the WACCM model (Lucas, 2017). Figure 3 illustrates the latitude-pressure distribution of the simulated zonal mean atmospheric conductivity. While the results of WACCM yield slightly lower conductivity values, the spatial distributions of conductivity are similar, namely, ≈10 −11 -10 −10 S/m as obtained by CCM SOCOLv2 versus ≈10 −10.5 -10 −10 S/m by WACCM. An increase of GCR-induced ionization rates toward the geomagnetic poles results in enhanced air conductivity at high latitudes. Temperature gradients also lead to latitudinal dependence of the air conductivity because of their impact on ion mobility from latitudinal temperature variations. It is important to note that local sources of aerosol can affect the total air conductivity of the atmosphere leading to high spatial and temporal variabilities. The conductivity, columnar (R c ), and global resistance can be calculated using the model output data. The columnar resistance is defined as the vertical integral of the reciprocal of conductivity i , where dz i is the ith layer thickness: (3) Finally, the downward electric current can be calculated as where IP is the ionospheric potential. The description and verification of IP calculations with the CCM SOCOLv2 are presented by Karagodin et al. (2019). The ionospheric potential varies in time but is the same in space, as it is the global sum. For each model step, it is represented by one number, which may be changed at the next step. Atmospheric electric parameters can be sensitive to changes in the local weather conditions, particularly in precipitation intensity (in the convective area) and wind speed. The calculated downward current J z for the fair-weather conditions along the Reading (UK) longitude is presented in Figure 4. For comparison, the J z is presented for two cases: ionization only due to GCR; and ionization induced by both Rn-222 and GCR. The difference between the two cases is about 0.2-0.6 pA/m 2 and appears in the Rn-222 active regions. In the GCR case (green curve in Figure 4), the line is everywhere below or equal to the blue curve (GCR + Rn222 case); however, the largest departures occur in zones with a huge landmass in Northern Hemisphere between 20 • N-60 • N. A comparison with the vertical current density map constructed by the Tropical Rainfall Mission Measurement (TRMM) satellite precipitation radar measurements (Lucas, 2017) and observations conducted by geophysical observatories in the United Kingdom (Nicoll & Harrison, 2014) showed that the CCM SOCOLv2 can reproduce the pattern of the atmospheric electrical parameters distribution. The values of downward current calculated with the SOCOLv.2 lie roughly between 1.4 and 2.4 pA/m 2 , being close to those from the TRMM (1.5 to 3 pA/m 2 ).

Experiment Setups
We considered the background ionization from GCR corresponding to a moderate phase of the solar cycle. The modulation potential, characterizing the energy spectrum of GCR fluxes as modulated by solar magnetic activity (e.g., Usoskin et al., 2005), was set to 600 MV (Koldobskiy et al., 2019). The flux and composition of GCR were modeled according to Gil et al. (2015). The geomagnetic field was set up using the IGRF (International Geomagnetic Reference Field) for epoch 2005. For the extreme SPE, we considered a scenario for the strongest known event of the year 774 AD (Sukhodolov et al., 2017). According to this scenario, the energy spectrum of SEPs was similar to that of the strongest directly observed SPE of 23 February 1956, but the fluence was scaled up with a factor of 100, which is an upper conservative limit of the extreme SPE integrated flux (Kovaltsov & Usoskin, 2014;Usoskin, 2017). Atmospheric ionization profiles were calculated for the selected scenario using the CRAC:CRII model (Usoskin & Kovaltsov, 2006;Usoskin et al., 2010). We have computed 3D distributions of the atmospheric ionization due to external factors for two scenarios: (1) the quiet conditions using only GCR and Rn-222, and (2) the extreme case, by adding also an extreme SPE to Scenario 1, as an input for the CRAC:CRII model. Figure 2 shows near-surface ionization rates for these two scenarios. For Scenario 2, the ionization rate can reach 30-100 ion pairs cm −3 s −1 (see panel b) exceeding those due to Rn-222 ( Figure 1b) and GCR (Figure 2a). The ionization rates are higher over the polar areas with reduced geomagnetic shielding.

Experiment Results
Our results imply that during a quiet period (Scenario 1; see Figure 5a), the atmospheric current density is affected by several competing processes of roughly equal importance: the geomagnetic shielding of GCR (visible as the latitudinal gradient); orography (enhanced spots on the Antarctic Plateau and Himalaya) due to the thinner atmosphere and thus higher GCR-related ionization rate; Rn-222 emission (spots in regions of enhanced volcanic and geological activity) (Iskandar et al., 2018). The overall spatial variability of the fair-weather vertical current is within a factor of two roughly between 1.5 and 3 pA/m 2 .
On the other hand, the vertical current density during the extremely disturbed period (Scenario 2 of an extreme SPE case; see Figure 5b) reaches much higher (up to a factor of 25) values exhibiting a different geographical distribution, where the most important factors are the latitudinal gradient and orography, both related to SEPs while near-ground Rn-222 does not play an important role. The variability over the globe is up to 1.5-order of magnitude for this scenario (from the maximum of 3 to about 90 pA/m 2 in the Antarctic plateau region). The strongest enhancement of the vertical current is expected in Central Antarctica and Greenland, while it is only a few tens of percent in the equatorial region, because of the geomagnetic shielding and a relatively soft SEP energy spectrum. We note that because of the implicit assumption that IP is generated independently of the atmospheric conductivity, a possible feedback mechanism of the IP reduction by the enhanced return current during an extreme SPE is not modeled. However, since SPE events are short, global ionospheric capacity is big, and the greatly enhanced return current is located mostly in spatially small polar regions, the effect is not expected to be very strong. In this sense, we provide a conservative upper bound of the SPE effect. However, this effect, obtained for the worst-case extreme event, becomes small when scaled down to the scale of events observed during the modern era: the strongest of directly observed SPEs of 23 February 1956 was a factor ≈100 weaker.

Conclusion
We have performed and analyzed model experiments corresponding to the time period of 2004-2005. The ionization rates from Rn-222 and GCRs were used to calculate the global distribution of the atmospheric conductivity and the fair-weather vertical current density J z . For the first time, we estimated the influence of an extreme SPE of 774 AD on atmospheric electricity. For the extreme scenario with the ionization induced by GCR, extreme SPE, and Rn-222, the atmospheric electricity response is substantial. The atmospheric vertical current is dramatically enhanced after the SPE, particularly in high-elevated polar regions (Greenland and Antarctica). We note that only the effects of the increased ionization rate due to SEPs were considered here, while other effects, such as the influence of a geomagnetic storm on the ionospheric potential (see Equation 4), were not treated, making this estimate conservative. It is concluded that an extreme SPE (the worst-case SPE scenario, corresponding to the event of 774 AD) may lead to a dramatic (a factor of up to 30) increase of the vertical atmospheric current density, reaching ≈90 pA/m 2 in high-elevated polar regions, with only a small enhancement in equatorial regions. Such extreme SPEs are rare events, and the effect of more regular SPEs are several orders of magnitude smaller. The quiet scenario with the ionization induced only by GCR and Rn-222 is observed during most of the time. In this case, the atmospheric vertical current may be up to ≈3.8 pA/m 2 in radon source regions and in the polar regions. Contribution of Rn-222 is important in middle-and low-latitude/altitude regions where the GCR-induced air conductivity is significantly lower than that at greater heights/latitudes. This indicates the importance of the Rn-222-induced ionization for the state of the GEC at quiet geomagnetic conditions.