Preliminary Study on the Generating Mechanism of the Atmospheric Vertical Electric Field before Earthquakes

Abstract

Precursor signals for earthquakes, such as radon anomalies, thermal anomalies, and water level changes, have been studied in earthquake prediction over several centuries. The atmospheric vertical electric field anomaly has been observed in recent years as a new and valuable signal for short-term earthquake prediction. In this paper, a physical mechanism of the atmospheric vertical electric field anomaly before the earthquake was proposed, based on which the Wenchuan earthquake verified the correctness of the model. Using Monte Carlo simulations, the variation of the radon concentration with height before the earthquake was used to simulate and calculate the ionization rates of radioactive radon decay products at different heights. We derived the atmospheric vertical electric field from −593 to −285 V/m from the surface to 10 m before the earthquake by solving the system of convection-diffusion partial equations for positive and negative particles. Moreover, negative atmospheric electric field anomalies were observed in both Wenjiang and Pixian before the Wenchuan earthquake on 12 May, with peaks of −600 V/m in Pixian and −200 V/m in Wenjiang. The atmospheric electric field data obtained from the simulation were shown to be in excellent concordance with the observed data of the Wenchuan earthquake. The physical mechanism can provide theoretical support for the atmospheric electric field anomaly as an earthquake precursor.

1. Introduction

Earthquakes(EQs) are one of the most serious natural disasters. There have been hundreds of Ms ⩾ 7 EQs worldwide from 2000 to now [1]. EQs cause direct damage to the buildings in the residential areas in the center of the EQ and cause different levels of damage to buildings of different materials and structures [2,3]. Secondary disasters such as tsunamis and landslides triggered by EQs cause enormous human casualties, construction damages, and property losses [4,5,6,7]. Increasingly vital research on pre-EQ disaster phenomena was conducted by researchers around the world to avoid significant losses caused by EQ disasters. For instance, Namgaladze, Karpov, and Knyazeva have found that low-frequency electromagnetic radiation is emitted in the preparation period for EQs in the ionosphere [8]. The abnormal phenomenon has been observed by ground-based stations and satellite detectors [8,9]. Moreover, the movement of the crust during the preparation period of many EQs also caused anomalous increases in radon concentrations, including water radon and gas radon [10,11,12,13,14]. In addition, there were also other observable anomalies that the researchers studied before the EQ, such as thermal anomalies and ground light phenomena [15,16]. However, the pre-EQ anomalies mentioned above generally occur on a time scale of a few days to a few months during the preparation period of an EQ, and it is impossible to predict the occurrence time and location of an EQ accurately. Furthermore, the measurement of radon concentration and low-frequency electromagnetic waves are associated with high uncertainties.

In recent decades, EQ forecasters have repeatedly detected negative abnormal signals in the atmospheric vertical electric field in front of the EQ. During fair weather and excluding weather conditions such as haze, thunder, and dust, this signal has occurred at the hourly level ahead of an EQ, and numerous observations have demonstrated the accuracy of this anomalous signal [17,18,19]. Although the atmospheric electric field anomaly was proven to be an important precursor before EQs, the cause of its formation was not clearly understood. In previous work [20], researchers claimed that atmospheric vertical electric field anomalies were caused by radon ionization of the air; however, the authors failed to propose a specific mechanism of interaction from radon concentration to atmospheric electric field. A mechanism for the formation of anomalous electrode effects by near-surface air ionization was proposed by Boyarchuk, Lomonosov, and Pulinets in 1997 [21]. They obtained an anomalous reduction effect of the atmospheric electric field within 1 m of the surface. Existing studies have confirmed that radon ionization of air contributes to the atmospheric vertical electric field anomaly, and some explorations of radon ionization atmospheric processes have also been conducted. However, up to now, there is no quantitative study on the signal of atmospheric vertical electric field anomaly caused by radon concentration.

In this paper, the physical mechanism of the formation of negative atmospheric vertical electric field anomalies during the EQ preparation period is proposed by studying the pre-EQ radon gas concentration data. This paper aims to obtain detailed information on the variation of the atmospheric electric field with height obtained by solving the proposed physical model using radon gas and atmospheric electric field data before the Wenchuan EQ. The first part of this paper briefly analyzes the formation mechanism of the atmospheric vertical electric field before the earthquake and describes in detail the formation process of anomalous atmospheric electric fields in a progressive way, including the variation of radon concentration with height, simulating the ionization rate of radon stable decay daughter at different heights and solving the convective diffusion equation of positive and negative particles to obtain the atmospheric electric field. In the second part, the 2008 Wenchuan EQ in China is used as an example for simulation and calculation, and the correctness of the model is verified by taking the radon concentration before the Wenchuan EQ as input. In the model of this paper, a Monte Carlo simulation of the radon ionization rate at different heights is proposed for the first time, and the model is verified by using EQ cases. The electrode effect proposed by Boyarchuk, Lomonosov, and Pulinets only considered the effect of radon ionization in the atmosphere within 10 cm of the surface, and it was a small-scale model [21]. Moreover, Boyarchuk, Lomonosov, and Pulinets only pointed out that the electrode effect can cause atmospheric electric field anomalies, and they did not validate it for a particular EQ. This mechanism provides a theoretical explanation for the pre-EQ atmospheric vertical electric field anomaly as an EQ precursor.

2. Methodology

In the build-up to a strong EQ, the crushing and tearing of the Earth’s crust occur during the crushing process. Cracks in the Earth’s crust have been extended to the surface, resulting in numerous microcracks at the surface [22]. In this context, radon, a radioactive gas in the Earth’s crust, diffuses into the atmosphere along with the cracks. As shown in Figure 1, the negative signal of the atmospheric vertical electric field due to the radon ionization process in the atmosphere is demonstrated. Radon diffuses from the Earth’s crust to the atmosphere through microcracks and then undergoes alpha decay, beta decay, and gamma decay in the air. A high amount of positive and negative ion pairs are produced by ionizing alpha particles, electrons, and gamma-rays in the air. The positive and negative ion pairs diffuse and migrate in the air. Then eventually, the distribution of positive and negative particles reverses the normal downward positive atmospheric vertical electric field in fair weather.

2.1. Radon Concentration

2.1.1. Vertical Distribution

In the troposphere over the mainland, the radioisotope of primary interest for small ion production is ${}^{222}$Rn, with a half-life of 3.8 days. Radioactive radon is transported upward with atmospheric movement, resulting in radon concentration variation with height. We reviewed the ${}^{222}$Rn in fair weather relative to altitude profile radioactive decay product distribution. For the fair weather mentioned above, Harrison and Nicoll (2018) defined fair weather in terms of the following detailed descriptions: (1) the visual range should be above 2 to 5 km, and the relative humidity should be less than 95%, (2) there are no low-level clouds and negligible cumulonimbus clouds, and (3) the wind speed at 2 m should be less than 7 m/s [23]. The average vertical distribution of ${}^{222}$Rn in different parts of the continent in fair weather was demonstrated by the experimental results [24]. As shown in Figure 2, the experimental results and best-fit curves for radon gas concentrations at different altitudes are presented. The best-fit trend line is logarithmically described by the following Equation:

$$y=-2.351ln\left(x\right)+13.377$$

(1)

where x (dpm/m${}^3$) is the Radon concentration and y (km) is the height. The fitting accuracy is $R^2=0.70453$.

2.1.2. Decay Products

Radioactive radon can undergo alpha decay, beta decay, and gamma decay in the air, producing many alpha particles, electrons, and gamma rays. A quantitative relationship between radon concentration and the corresponding decay products is first required to obtain the radon ionization rate at various heights. The decay paths of ${}^{222}$Rn and its daughter are depicted in Figure 3. We simulated the decay of 100,000 radon in a vacuum using Monte Carlo simulation software and recorded the number of stable decay products. As shown in

Table 1. Stabilization products of radon decay.

Particle Type	Number	Average Energy	Range in Air
${}^{222}$Rn	100,000	0	0
alpha	400,000	6.123 MeV	4.94 cm
e	520,866	263.4 keV	63.16 cm

, the alpha to beta particles ratio in the decay products is 4:5.2.

The radon concentration is the number of alpha particles resulting from decay per cubic meter per minute. As a result, we deduced the variation in the concentration of radon’s main decay products, alpha and beta particles, with the height from the radon distribution.

2.2. Ionization Rate at Different Heights

To obtain the ionization rates caused by radon at different heights, we used Monte Carlo simulation software GEANT4 (Geometry And Tracking) for simulation [25]. The ionization rates of alpha and beta particles in the air were simulated from 0 to 100 m. A series of isotropic surface sources of alpha and beta particles, corresponding to the concentration of the corresponding particles at different heights, were set up separately. Following that, the ionization rates at different heights were calculated by combining the average ionization energy of the particles in the air with the energy deposition of detectors set at different heights.

2.2.1. Geometric Models and Particle Source

As shown in Figure 4, a geometric model was built that involved 38 identical detectors, recording energy deposition at different altitudes. Each detector is a thin slice with a radius of 5 m and a height of 1 mm. We configured one detector per 1 cm interval for 0–9 cm, one detector per 10 cm interval for 10–90 cm, one detector per 1 m interval for 1–9 m, and one detector per 10 m interval for 10–100 m. Since the atmospheric density varies very little within 0–100 m, we assumed sea level air for the world and the detector material. Sea level air is comprised of 0.0000124 C, 0.755267 N, 0.231781 O and 0.012827 Ar with a density of 1.20479 × 10${}^{-3}$ g/cm${}^3$.

Whereas alpha particles have a range of 4.94 cm in air, all isotropic sources were arranged ±5 cm around the detector; if this distance was exceeded, alpha particles could never access the detector. Moreover, surface sources were established at 1 cm intervals to measure the alpha emitted within a volume of 5 m radius and 1 cm height at that location. For instance, the detector located at 1 m was surrounded by an alpha source of 5 m radius at 0.95, 0.96, 0.97, 0.98, 0.99, 1, 1.01, 1.02, 1.03, 1.04, and 1.05 m. The range for beta particles is 63.16 cm. Due to the longer range of the beta particles, we did not utilize all 38 detectors as we did for alpha particles. Only 10–100 m detectors were used, and a beta isotropic surface source was designed with a radius of 5 m and spacing of 10 cm. The beta surface sources around the ten detectors were positioned at ±60 cm. As with alpha particles, surface sources were placed at 10 cm intervals to measure the beta emitted within a volume of 5 m radius by 10 cm height at that location. The number of particles emitted per second from the source at different heights was fixed from the distribution of alpha and beta particles with height.

2.2.2. Physical Lists and Data Statistics

The stable products of ${}^{222}$Rn decay are alpha with an average energy of about 6.123 MeV and beta with an average energy of about 263.4 keV. Alpha and beta particles entering the air undergo complex physical processes that produce many secondary charged particles. Secondary particles are tracked to the lowest possible energy to improve the accuracy of the simulated process. Different physical processes occur when particles of different types and energies interact with air. This paper considered the following physical processes: low-energy electromagnetic processes, hadron processes, decay processes, and photonuclear processes. A secondary particle generation threshold of 250 eV was imposed within the entire model.

Deposition of energy at different height locations was recorded by distinct detectors in the model. First, the total energy deposition for each detector was calculated by multiplying the average energy deposition by the number of simulated particles identified and the duration. Jesse discussed the generation of ions in the air by ionizing radiation [26]. He reported the mean energy required to generate an ionization in the air for alpha and beta particles, 35 eV for alpha particles and 33.8 eV for beta particles. Next, the number of positive and negative particle pairs deposited within each detector was obtained by dividing the total energy deposited by the average ionization energy of alpha and beta particles. Finally, the ionization rates of detectors at various heights were evaluated.

The following case is illustrated with an example of a detector. In this paper, we assumed that the average energy deposition of the detector is E eV, the number of simulated events is N pcs, the ionization duration event is t s, the total energy deposition is $N_{tot}$ eV, the average ionization energy of particles in the air is W eV, the number of positive and negative particle pairs in the detector is $N_i$ ion parts, and the ionization rate is Q cm${}^3·$s. The technique outlined above was presented as Equation (2), and then Equation (3) was obtained.

$$E\times N=N_{\mathit{tot}};\frac{N_{\mathit{tot}}}{W}=N_i;\frac{N_i\times {10}^{-6}}{\pi \times 5^2\times t}=Q$$

(2)

$$\frac{E\times N}{W}\times \frac{{10}^{-6}}{\pi \times 5^2\times t}=Q$$

(3)

2.3. Physical Lists and Data Statistics

Positive and negative particles are ionized by radon decay particles at different heights. However, the different mobility of differently charged ions leads to the inability of positive and negative particles to fully compensate for each other. Then, the reverse atmospheric vertical electric field is produced in the case of strong turbulent conditions diffusion. For simulating the electrical state of the turbulent ground layer, a steady-state kinetic equation was applied in this paper.

$$\begin{array}{c}\hfill -\frac{\partial }{\partial z}\left[D_t\frac{\partial n_{+}}{\partial z}\right]+b_{+}\frac{\partial }{\partial z}\left(E{n}_{+}\right)=Q-\alpha n_{+}{n}_{-}\\ \hfill -\frac{\partial }{\partial z}\left[D_t\frac{\partial n_{-}}{\partial z}\right]-b_{-}\frac{\partial }{\partial z}\left(E{n}_{-}\right)=Q-\alpha n_{+}{n}_{-}\\ \hfill \frac{\partial E}{\partial z}=4\pi e\left(n_{+}-n_{-}\right)\end{array}$$

(4)

Here, n± is the concentration of positive and negative ions; $D_t\left(z\right)=(Kz+\gamma )/(z+\beta )$ is the turbulent diffusion coefficient; b± is the mobility of positive and negative ions; Q is the ionization rate, which is mainly caused by radon ionization of air near the ground; $\alpha$ is the ion recombination coefficient; e is the electron charge; and ${\epsilon }_0$ is the vacuum dielectric constant.

The fluctuations in positive and negative ion concentrations and the atmospheric vertical electric field with height are provided by solving Equation (4). It proves that the formation of the atmospheric vertical anomaly is caused by the ionization of radioactive radon decay products into the atmosphere.

3. Wenchuan EQ Analysis and Discussion

3.1. Radon Concentration before the Wenchuan EQ

The Wenchuan Ms8.0 EQ in 2008 in China was used as an example for our study to verify the validity of the mechanism. The distances between the epicenter of the Wenchuan EQ and the detection stations at Pixian, Wenjiang, and Guzan are 49, 52, and 202 km, respectively. As shown in Figure 5a, before the Wenchuan EQ, Guzan shows a significant increase in the daily median radon concentration from May 7 to the moment of the EQ [19]. As shown in Figure 5b, the visibility and low cloud cover at the Wenjiang and Pixian stations on the day of the EQ. As shown in Figure 5c, the wind velocity at Wenjiang and Pixian stations throughout the day on 12 May. Therefore, it can be determined that it was fair weather in the epicenter area of the Wenchuan EQ on 12 May. The atmospheric vertical electric field is disturbed only a little by meteorological effects.

It can be seen from Figure 5a that the radon concentration before the Wenchuan EQ peaked at 8.98 Bq/L on May 10. For receiving the variation of radon concentration with height before the Wenchuan EQ, it was considered that the peak radon concentration before the Wenchuan EQ was the ground, y = 0 km, radon concentration, i.e., (8.98 Bq/L, 0 km). The vertical distribution of radon concentration varies with the changes in meteorological parameters [27]. For the simulation in this paper, we considered the constant migration and turbulent diffusion coefficients. Thus the slope of the natural logarithm of radon concentration with height was considered to be the same as the curve of Equation (1). The intercept of the curve is determined from the radon concentration data before Wenchuan earthquake (8.98 Bq/L, 0 km). The radon concentration variation with height was obtained as:

$$y=-2.351ln\left(x\right)+31.026$$

(5)

where x (dpm/m${}^3$) is the Radon concentration and y (km) is the height.

The concentration of alpha particles at different heights was obtained from Equation (5), and the variation of beta particles with height was obtained from the alpha to beta ratio in Section 2.1.2. Since we configured the Monte Carlo simulation alpha particle source equivalent to the number of alpha particles generated in a cylinder with a radius of 5 m and a height of 1 cm in 1 h, it was possible to obtain the number of alpha particles simulated at different heights from the alpha particle concentration at different heights.

3.2. Results and Discussion

3.2.1. Monte Carlo Simulation of Ionization Processes

Based on the variation of alpha and beta ion concentrations with height, the Monte Carlo simulations were used to simulate their ionization rates at different heights. For the Wenchuan EQ alpha and beta particle concentrations given in Section 3.1, the variation of the ionization rate with height was derived from the simulation method in the methodology, as shown in

Table 2. Alpha particle simulation parameters and results.

Altitude/m	Average Energy Deposition/keV	Simulation Number	Total Energy Deposition/eV	Ion Pairs	Ionization Rate/cm${}^3·$s
0.1	60.8249	2.7924 × 10${}^8$	1.6985 × 10${}^{13}$	4.8527 × 10${}^{11}$	1716.3
1	60.8241	2.7913 × 10${}^8$	1.6978 × 10${}^{13}$	4.8508 × 10${}^{11}$	1715.62
2	60.824	2.7901 × 10${}^8$	1.6971 × 10${}^{13}$	4.8487 × 10${}^{11}$	1714.89
3	60.8242	2.7889 × 10${}^8$	1.6963 × 10${}^{13}$	4.8467 × 10${}^{11}$	1714.17
4	60.8238	2.7877 × 10${}^8$	1.6956 × 10${}^{13}$	4.8446 × 10${}^{11}$	1713.43
5	60.824	2.7866 × 10${}^8$	1.6949 × 10${}^{13}$	4.8426 × 10${}^{11}$	1712.7
6	60.8241	2.7854 × 10${}^8$	1.6942 × 10${}^{13}$	4.8405 × 10${}^{11}$	1711.98
7	60.8241	2.7842 × 10${}^8$	1.6935 × 10${}^{13}$	4.8384 × 10${}^{11}$	1711.25
8	60.8239	2.7830 × 10${}^8$	1.6927 × 10${}^{13}$	4.8364 × 10${}^{11}$	1710.52
9	60.8239	2.7818 × 10${}^8$	1.6920 × 10${}^{13}$	4.8343 × 10${}^{11}$	1709.79
10	60.8237	2.7806 × 10${}^8$	1.6913 × 10${}^{13}$	4.8322 × 10${}^{11}$	1709.06
20	60.823	2.7688 × 10${}^8$	1.6841 × 10${}^{13}$	4.8117 × 10${}^{11}$	1701.78
30	60.8226	2.7571 × 10${}^8$	1.6769 × 10${}^{13}$	4.7912 × 10${}^{11}$	1694.55
40	60.8215	2.7454 × 10${}^8$	1.6698 × 10${}^{13}$	4.7708 × 10${}^{11}$	1687.32
50	60.8206	2.7337 × 10${}^8$	1.6627 × 10${}^{13}$	4.7505 × 10${}^{11}$	1680.14
60	60.8196	2.7221 × 10${}^8$	1.6556 × 10${}^{13}$	4.7302 × 10${}^{11}$	1672.98
70	60.8209	2.7106 × 10${}^8$	1.6486 × 10${}^{13}$	4.7103 × 10${}^{11}$	1665.91
80	60.8216	2.6991 × 10${}^8$	1.6416 × 10${}^{13}$	4.6903 × 10${}^{11}$	1658.86
90	60.8207	2.6876 × 10${}^8$	1.6346 × 10${}^{13}$	4.6703 × 10${}^{11}$	1651.8
100	60.8228	2.6762 × 10${}^8$	1.6277 × 10${}^{13}$	4.6507 × 10${}^{11}$	1644.84

and

Table 3. Beta particle simulation parameters and results.

Altitude/m	Average Energy Deposition/keV	Simulation Number	Total Energy Deposition/eV	Ion Pairs	Ionization Rate/cm${}^3·$s
10	462.855	4.27 × 10${}^9$	9.89 × 10${}^{11}$	2.93 × 10${}^{10}$	103.45
20	462.958	4.25 × 10${}^9$	9.85 × 10${}^{11}$	2.91 × 10${}^{10}$	103.04
30	462.936	4.24 × 10${}^9$	9.80 × 10${}^{11}$	2.90 × 10${}^{10}$	102.59
40	462.963	4.22 × 10${}^9$	9.76 × 10${}^{11}$	2.89 × 10${}^{10}$	102.17
50	462.893	4.20 × 10${}^9$	9.72 × 10${}^{11}$	2.88 × 10${}^{10}$	101.72
60	462.878	4.18 × 10${}^9$	9.68 × 10${}^{11}$	2.86 × 10${}^{10}$	101.28
70	462.954	4.16 × 10${}^9$	9.64 × 10${}^{11}$	2.85 × 10${}^{10}$	100.87
80	462.938	4.15 × 10${}^9$	9.60 × 10${}^{11}$	2.84 × 10${}^{10}$	100.44
90	462.892	4.13 × 10${}^9$	9.56 × 10${}^{11}$	2.83 × 10${}^{10}$	100
100	462.978	4.11 × 10${}^9$	9.52 × 10${}^{11}$	2.82 × 10${}^{10}$	99.59

and Figure 5. Table 2 and Table 3 show that the heights and simulation number are the simulation input conditions, and the average energy deposition is the output from the GEANT4 calculation. It was noteworthy that the total energy deposition is the case of one-hour accumulated radon concentration. Next, the ionization rate variation was obtained by Equation (3).

As shown in Figure 6, the black dots are the results obtained from the simulation, and the red curve is the best-fit result. An exponential fit was used as the fitting function to match the actual ionization rate variation more closely with height. As a result, the ionization rate fitted alpha and beta ionization rate curves were Equations (6) and (7).

$$Q_{\alpha }=e^{21.263-4.25\times {10}^{-4}z}$$

(6)

$$Q_{\beta }=e^{18.459-4.245\times {10}^{-4}z}$$

(7)

here the ionization rate fit is close to perfect, where the $Q_{\alpha }$ fitting accuracy $R^2=1$ and the $\beta$ fitting accuracy $R^2=0.99995$. Therefore, the variation of the ionization rate with height was $Q_{\alpha }+Q_{\beta }$ caused by radon.

3.2.2. Simulation Results and Discussion

From the previous work, the variation of the ionization rate with height due to radon was given as $Q=Q_{\alpha }+Q_{\beta }$. In this model, the mobility $b_{\pm }$ of positive and negative ions is set to 2.5 × 10${}^{-4}$ and 2.2 × 10${}^{-4}$ m${}^2·$(V·s)${}^{-1}$[28], respectively, and the ion recombination coefficient $\alpha$ is set to 1.6 × 10${}^{-12}$. For the turbulent diffusion coefficient $D_t\left(z\right)=(Kz+\gamma )/(z+\beta )$, where $\beta$ is set to 10 m, $\gamma$ is set to 5 × 10${}^{-5}$ m${}^3·$s${}^{-1}$ and the strong turbulence coefficient K is set to 5 m${}^2·$s${}^{-1}$ [21].

For the set of partial differential Equation (4), the boundary conditions near the Earth’s surface were considered as:

$$n_{+}\left(z=0\right)=n_{-}\left(z=0\right)=450{\mathrm{cm}}^{-3},n_1\left(10\right)=n_2\left(10\right)={(Q\left(10\right)/\alpha )}^{1/2}$$

(8)

$$E\left(10\right)=\frac{j_0}{{\lambda }_{10}},{\lambda }_{10}=e\times \left(b_{+}{n}_{+}\left(10\right)+b_{-}{n}_{-}\left(10\right)\right)$$

(9)

In Equation (8), $j_0$ is the current density at $z=0$. For the Wenchuan EQ, Kuo et al. developed a model for the piezoelectric effect in rocks and have calculated a current density of $I_{rock}$ = 500 nA·m${}^{-2}$ for an EQ fault area of $A_{rock}$ = 200 × 30 km${}^2$ [29]. The total current input to the atmosphere equals the current output from rocks.

$$I_{air}=A_{air}\times J_0=I_{rock}=A_{rock}\times J_{rock}$$

(10)

The electromagnetic field parameters of a large area of China were abnormal before the Wenchuan EQ. Abnormal electromagnetic signals were observed at the Gao Bei station in Hebei Province, 1440 km from the epicenter, on May 9, three days before the main EQ. Electromagnetic anomalies were also observed at the Qingxian station in Hebei province, 1439 km from the epicenter, and at Ningjin county in Hebei province, 1241 km from the epicenter. However, there was no electromagnetic signal anomaly observed in Sanhe County, Hebei Province, 1541 km from the epicenter [30]. Thus, it was assumed that the area of the pre-EQ air anomaly region A = 1440 × 1440 km${}^2$ calculated by Equation (9) yields $J_0\approx$ 1.5 nA.

The partial differential Equation (4) was solved mathematically by the boundary conditions given for the Wenchuan EQ. As shown in Figure 7, we solve the variation of the positive and negative particle concentrations and the atmospheric vertical electric field with height. From the data in Figure 7, the separation of positive and negative particles can be seen.

As shown in Figure 7, we calculate that the atmospheric vertical electric field from the surface to 10 m is −593 to −285 V/m. Negative atmospheric electric field anomalies were observed at both Wenjiang and Pixian on 12 May, with a peak of −600 V/m at Pixian and −200 V/m at Wenjiang [19]. The data of atmospheric electric field sounding at Wenjiang station were consistent with the data obtained from the model calculation in this paper. Moreover, there were some deviations in the data of the Pixian station. The results were consistent with those calculated by Boyarchuk, Lomonosov, and Pulinets, which showed an anomaly of reduction in the atmospheric electric field [21]. The difference was that the calculations used as input conditions were the atmospheric electric field on fair weather and the radon ionization source within 10 cm of the ground. An analogous process was observed in the case of measuring the electric field in the underground nuclear weapons test area: the electric field drastically decreased at the time of the explosion [31]. In this case, the near-surface layer was ionized by the fission fragments emerging from the soil. This experiment was similar to the mechanism of this paper and was a direct demonstration of this paper. Since we could not obtain the detailed radon concentration detection data before the Wenchuan EQ, the radon concentration used in the calculation was based on the gradient of the natural logarithm of radon concentration with height in fair weather. The variation of radon concentration with height was obtained using the radon concentration at the ground level y = 0 km before the Wenchuan EQ. Moreover, the parameters in the equation of motion were according to the average parameters on fair weather, which might have some errors with the actual situation. Therefore, the parameters used in the simulation deviated from the actual parameters. If the actual meteorological parameters before the EQ and the variation of radon concentration with height were available in the simulation, the accuracy of the simulation would be improved. In general, these results still indicate that the mechanism of our initial exploration matches well with the measured values and indicates a good validity of our model.

4. Conclusions

A research method based on a combination of Monte Carlo simulations and the numerical solution was proposed to study the mechanism of the formation of the atmospheric vertical electric field of the short-term precursor signals of EQs. An innovative Monte Carlo simulation method was established to calculate the ionization rate at different heights from data on the variation of radon concentration with height before an EQ. Next, the ionization rate data were combined with the partial differential equations describing the motion of positive and negative particles to numerically solve for the variation of positive and negative particles and radon concentration with height. Finally, the model was validated with the 2008 Wenchuan EQ in China, which had both radon concentration data and atmospheric electric field data before the EQ. In addition, the ionization of alpha and beta particles at different heights, which were the main decay stabilization products of radon, was considered in the model. Alpha particles were the best ionized of the radon decay products, and their fitted ionization rates decayed exponentially with height. The results showed that the atmospheric electric field results obtained from the simulations compounded well with the atmospheric electric field data detected before the Wenchuan EQ. The atmospheric electric field measured by the Wenchuan EQ is direct evidence of the correctness of the model. Moreover, our atmospheric electric field simulations’ results were consistent with those obtained by Boyarchuk, Lomonosov, and Pulinets. In this paper, Holzer’s experiments on the atmospheric electric field effects of nuclear explosions in 1972 fundamentally proved the feasibility of the atmospheric electric field anomalies generated by radon ionization. Therefore, the mechanism model can effectively explain the mechanism of atmospheric electric field anomalies before EQ.