Modeling of the Lower Ionosphere During Solar X‐Ray Flares of Different Classes

This study presents the results obtained from modeling the lower ionosphere response to C‐, M‐ and X‐class solar X‐ray flares. This model is based on a 5‐component scheme for the ionization‐recombination cycle of the ionospheric D‐region. Input parameters for the plasma‐chemical model under different heliogeophysical conditions corresponding to selected X‐ray flares were determined by using data received from the AURA, SDO, and GOES satellites. Verification of the obtained results was carried out with use of ground‐based radiophysical measurements taken at the geophysical observatory Mikhnevo. The results obtained from a comparison of the calculated and experimental radio wave amplitude variations along six European very low frequency (VLF) paths show that the average normalized root mean square error is ∼7%, 14%, and 18% for C‐, M‐, and X‐class flares, respectively. Qualitative and quantitative analysis of the verification results for the VLF signal amplitude show good predictive capability of the model built for describing weak and moderate ionospheric disturbances.

capability at solving problems for VLF signal propagation. Moreover, the quantity of radiophysical data is many orders of magnitude greater than the electron concentration measurements taken for the D-region, which makes it possible to carry out verification under almost any heliogeophysical conditions including solar flares.
This work is devoted to modeling the lower ionosphere during solar X-ray flares of different classes and analysis of the results obtained by verification with ground-based radiophysical measurements taken at the Mikhnevo geophysical observatory of the Sadovsky Institute of Geospheres Dynamics (Gavrilov et al., , 2019a).

Experimental Data Obtained at the Mikhnevo Geophysical Observatory
Since 2014, Mikhnevo has been continuously monitoring the amplitude and phase characteristics of electromagnetic signals from VLF transmitters located throughout the world (Gavrilov, Zetser, et al., 2019).
To verify the results of modeling, we used measurements of the amplitude of the signals received at the Moscow region observatory Mikhnevo (55°N 38°E) from six transmitters located at the European midlatitudes ( Figure 1). The main characteristics of the transmitters are listed in Table 1.
The amplitude of the signal passing through the simulated medium was calculated by using the LWPC program (Ferguson, 1998). The configured nominal power of the transmitters was 1 kW. Amplitude values were taken in (dB) for the purpose of quantitative comparison of theoretical and experimental daily variations. The amplitude shift due to the difference between the unknown actual transmitter power and the power used for the calculations was determined from the difference between the theoretical amplitude value and the experimental data obtained at a quiet heliogeophysical day that preceded each considered flare (Palit et al., 2013).

Plasma-Chemical Model of the Lower Ionosphere
As stated in the introduction, at present there are a significant number of empirical and theoretical models of the lower ionosphere, describing its state with some accuracy. For solving applied problems of radio wave propagation, the two-parameter Wait-Ferguson model (Ferguson, 1995;Wait & Spies, 1964) is most often used. On the one hand, this model makes it possible to simulate the radioequivalent ionosphere to successfully calculate radio wave propagation. On the other hand, the modeled vertical electron concentration profiles do not correspond to the real profiles, since the model is based on an exponential Ne profile.
When choosing the number of plasma-chemical processes on which the model under development is based, it is necessary to take into account not only the expected accuracy of the results but also the possibility BEKKER ET AL.    (Krivolutsky et al., 2015;Turunen et al, 1992Turunen et al, , 1996Verronen et al., 2016). However, such models require extremely high computational power and time resources, and, therefore, the applied tasks of forecasting radio wave propagation cannot be promptly solved using these models. Moreover, the number of unknown reaction rate constants increases with the number of considered plasma-chemical processes, which can lead to higher errors.
The system of differential equations for the ionization-recombination cycle describes the behavior of charged and neutral components, the dynamics of which are most important at D-region heights. Analysis of the results obtained from calculation of the electron concentration in four-component (Glukhov et al., 1992), five-component (Egoshin et al., 2012), and eight-component (Kudryavtsev & Romanyukha, 1995) plasma-chemical models of the lower ionosphere shows no significant difference between the results for the latter two models. However, not accounting for the variations of concentration of the negative cluster ions that are absent in the four-component model leads to a noticeable decrease in the Ne concentration at altitudes less than 70 km, which fundamentally affects the results of the calculation of the VLF radio wave propagation both under calm conditions and at X-ray flares.
Therefore, the five-component ionization-recombination cycle of the ionospheric D-region (1) was selected as the best model for this work. This system describes the behavior of the concentration of electrons Ne and four ion types: NO + , O 2 -, positive and negative cluster ions XY + , XY -. It takes into account almost all major photochemical processes occurring in the lower ionosphere.

Input Data
The system input parameters are the ionization rate q, concentration of neutrals It is known that ionization rate is the key parameter responsible for changes in the Ne concentration during solar flares, which is why its calculation must be carried out particularly carefully. At the same time, for a correct modeling of the lower ionosphere during disturbances of various natures, we need to know the state of the medium before and after the perturbation. In other words, we must correctly calculate the background Ne values. There are at least two reasons for this. First, the accuracy for calculating the background concentration of Ne significantly affects the quality of the modeling of small solar flares (Palit et al., 2013). Second, vertical profiles for Ne during a calm Sun are used to normalize the radio wave amplitude during verification of the results (see paragraph 2). At the same time, the accuracy of the calculation of the background electron concentration values in addition to the ionization rate significantly depends on the temperature and concentration of neutrals (Bekker, 2018). Therefore, great attention was paid to the accuracy of determining not only q but also other parameters of the system 1.  (Livesey et al., 2013) was performed. This period includes a maximum and two minimums for the solar activity.
It is obvious that when solving the system of ionization-recombination cycle under some specific heliogeophysical conditions, use of averaged values rather than separate measurements from a satellite as input parameters is more correct. In addition, it is necessary to approach such averaging carefully, because the results of the calculation of the electron concentration significantly depend on the selected ranges for the heliogeophysical conditions -latitude, longitude, zenith angle, season, index F10.7. On the one hand, the selected ranges for the heliogeophysical conditions should be wide enough to include a representative sample of measurements, and, on the other hand, a considered parameter should not change significantly within the given limits.  sity function of temperature P(T) normalized to its maximum. As shown, the temperature has a latitudinal dependence. Therefore, the latitude step ∆φ must be selected in such a way that the temperature does not change significantly. For example, during the summer months, when choosing ∆φ = 10°, the difference between the nearby latitudinal T values does not exceed 3.5%. The final value of ∆φ was chosen so that the temperature change did not exceed 5%.
In addition, one should pay attention to how the function P(T) varies with season. Not only does the average value profile differ but also its dispersion changes: the spread of values for the winter period is much larger than that for the summer. This point was also taken into account for the simulation, and temperature values were considered for each month separately. For solar activity, as a rule, the dispersion of the values increases with its growth, but the median remains practically unchanged, especially within the latitudinal range we are modeling (30-60°). That is why separation of data according to levels of solar activity was not carried out.
The next parameter essentially influencing Ne at the D-region altitudes is the concentration of neutrals M. It is found that the average concentration does not depend on latitude. However, an increased spread in values is clearly observed on approach to the poles. A change in seasonal value is seen only at h ≤ 80 km and not for the entire range of latitudes. Figure 3 shows the seasonal variation of the neutrals concentration M for different latitude ranges at an altitude of 75 km. The colors in the figure correspond to the values of the probability density function P(M) normalized to its maximum. As follows from Figure 3, the seasonal M profile is observed at latitudes of φ > 40°, and, within the equatorial part, its value remains almost constant.
Since the concentration at low latitudes does not depend on the season, it is natural that the dispersion of concentrations here is lower.
The daily variation of temperature and concentration of neutrals was not detected at any of the altitudes.
The behavior of small neutral components is nearly the same as that of the concentration of neutrals concentration required for the calculation of q are given in (ppmv) (Anderson et al., 1986;Bekker, 2018;Brunelli & Namgaladze, 1988). Other regularities registered make an insignificant contribution to the already discovered dependences, so they are not discussed here.

Calculation of the Ionization Rate
As noted above, one of the key parameters of the system of Equation 1 is the ionization rate q. The main source of the Earth's atmosphere ionization is the electromagnetic radiation of the Sun with wavelengths of λ ≤ 134 nm.
Under calm conditions, the main contribution to the ionization of the medium by photons with λ ≤ 102.57 nm arises from the ionization of the main atmospheric components. Most of this radiation is absorbed above 90 km (Solomon & Qian, 2005). Softer radiation ionizes small components, mainly the NO molecule, with an ionization threshold of 134 nm. The loss of the photon flux with λ > 101 nm occurs due to interaction with O 2 molecules. Particle ionization in this range occurs by the photoelectric effect mechanism, which leads to the production of only one electron with an energy that is not sufficient to knock out additional electrons by electron impact.
As shown by satellite measurements, solar radiation with λ > 101 nm has weak temporal variations. Measurements of the EUVE channel (115-127 nm) from the GOES-13, -14, -15 satellites (Machol & Viereck, 2016;Machol et al., 2016) show that the maximum of the radiation flux W EUVE (W/m 2 ) is approximately two times higher than its minimum. Moreover, a sharp change in W EUVE values is observed during some solar flares.
Taking into account the above, an algorithm for calculating the ionizing effect of UV radiation on the D-layer of the ionosphere was developed based on the LASP Reference Spectrum (No. 1 March 25-29, 2008) (https://lasp.colorado.edu/lisird/data/whi_ref_spectra/).
The following functions are determined for α particles based on the LASP spectrum: where   O is the absorption cross section of O 2 molecules, and   ion, is the ionization cross section of the α particles. The cross sections were taken from previous papers (Fennelly & Torr, 1992;Keller-Rudek et al., 2013) and an electronic resource (http://satellite.mpic.de/ spectral_atlas/cross_sections/). The electron production rate at a point with a particle concentration of N α can be calculated by using the following formula: where W LASP = 6.79·10 −3 W/m 2 is the radiation flux of the reference spectrum in the 115-127 nm range, and W EUVE is the measured flux in the channel.
This work takes into account the ionization of NO molecules in the 102.57-134 nm range and O 2 in the 101-102.57 nm range.
Ionization by galactic cosmic rays Q CR plays an important role below ∼70 km under calm conditions. This source of ionization for mid-latitudes is based on previous papers by Thomas & Bowman (1985) and Ermakov et al. (1997) and does not change over time.
The X-ray ionization Q XR (100 eV-100 keV) plays a decisive role in the ionization of the lower ionosphere during an X-ray flare. Ionization occurs via two different ways: photoelectric effect and Compton scattering, and photoelectrons lead to additional ionization by electron impact in the general case. In the framework of this work, the electron production rate inducted by X-ray radiation was determined as: where ε = 33 eV/electron is the electron energy price, hν is the photon energy, E(hν) is the average energy transferred from a photon to the medium in one interaction, E c (hν) is the average energy of a Compton electron, N β is the concentration of particles β, F(hν) is the spectral density of photons arriving at the considered point in space, F XR (hν) is the spectral density of photons in the Earth's orbit, τ(hν) is the optical depth, S β are quantities similar to 3 for particles β,   ion, ,   com, , and   abs, are cross sections for photoionization, Compton scattering, and total absorption of photons by particles, respectively. In Equation 6, the particles N 2 , O 2 , and O were taken into account, and, in Equation 8, Ar was added to them. Cross sections were determined according to the EPDL97 model (Cullen et al., 1997). Molecular cross sections were calculated as double cross sections of the corresponding atoms.
The spectral density of the photons was calculated from measurements in the XL (0.1-0.8 nm) and XS (0.05-0.4 nm) channels of the GOES-13 and -15 satellites (with correction) and the QD channel (0.1-7 nm) of the SDO satellite (Woods et al., 2012) according to the HMSXS model (Korsunskaja, 2019). Satellite data averaged with a step of 1 min were used.
Thus, the total ionization rate was calculated by using the following formula: The solar gamma-ray is not taken into account in this approximation.

Calculation of Electron Concentration
Altitude profiles for the electron concentration were calculated for separate VLF paths, for which further verification of the results was carried out. We calculated the profiles along the paths with steps of ∼300 km that could vary slightly so that the profiles were distributed evenly.
For the simulation, we selected several X-ray flares of different classes (C-, M-and X-class), meeting the following conditions: 1. Available X-ray solar flux measurements from the GOES and SDO satellites to be used for the ionization rate calculation 2. Full illumination of the European VLF paths during the flares 3. Operable condition of VLF signal transmitters and a receiver located at the Mikhnevo geophysical observatory, allowing for verification of the results According to these criteria, we selected several X-ray flares of different classes that occurred in October 2013 and June 2014.  Table 2 were used for normalization of the amplitude values.
We selected these two series of consecutive flares of different classes because we also wanted to check if the temporal dynamics of small neutral atmospheric components during the X-radiation growth should be  , that fell within the above ranges, were used to construct the probability density functions with an altitude step of ∆h = 5 km. These functions were used for calculation of the most likely values, which we applied to solve the system of differential equations for the ionization-recombination cycle. The ionization rate q was determined as a function of time with a step size of ∆t = 1 min. The system of Equation 1 was solved for each set of the obtained vertical profiles for the input parameters. At the output, we obtained the electron concentrations under calm conditions and during solar flares of different classes.   (Hedin, 1991). The blue error bars represent the standard deviation of the AURA data selected for the calculation.
The curves obtained show a difference between the T and M values obtained from the MSIS model and the most likely values obtained from the satellite. In addition, in Figure 5, one can see a significant contribution to the accuracy of the Ne concentration calculation made by not only the ionization rate but also by the temperature and concentration of the neutral components. Therefore, we assume that the statistical analysis of the input parameters based on the satellite data should fundamentally improve the D-region simulation quality at least for calm conditions and C-and M-class flares when the concentration of Ne is the most sensitive to the considered variations.
BEKKER ET AL.

Model Verification Based on the Experimental Data Obtained for VLF Radio Wave Propagation during X-Ray Flares and Discussion of the Results
As mentioned above, we used the amplitude characteristics of the VLF transmitter signals (Table 1) received at the Mikhnevo geophysical observatory for verification of the results obtained for modeling of the lower ionosphere at the flares listed in Table 2. BEKKER ET AL.
10.1029/2020JA028767 9 of 15   It is found out that the constructed plasma-chemical model allows us to predict the amplitude response to weak ionospheric disturbances with sufficiently high accuracy, which is evident in the bottom panels of Figure 6. Moreover, three recurrent C-class flares that occurred on June 10, 2014 are correctly simulated, which indicates a minor influence of the dynamics of small neutral components.
As for the M-class flares, verification of the lower ionosphere model built according to the radiophysical data showed that the model allows describing quantitative variation of the experimental value of the amplitude received from the GBZ transmitter. The difference between the experimental and calculated values reaches a maximum during the highest ionization rate, with values lying within a range of 1.5 dB.
The top right panel of Figure 6 presents the results obtained by simulating the ionospheric response to a X-class flare that occurred on October 25, 2013. The amplitude dynamics simulation is not correct in the part following the flare start at 08:00 UT. The current heuristic model of the radiation spectrum and ionization of the lower ionosphere (Korsunskaja, 2019) does not describe ultrahard X-radiation, which could be the reason for the difference between the theoretical and experimental results. Analysis of the Sun's radiation spectrum based on RHESSI observatory data (http://sprg.ssl.berkeley.edu/∼tohban/browser/) has proven the presence of ultrahard X-radiation in this flare.
It is obvious that X-class flares have the greatest influence on the amplitude variations, which can lead to a loss of signal at these frequencies. Therefore, the ability of the model to predict the qualitative and quantitative behavior of experimental values for the amplitude A exp during high energy flares is of great interest in this work. As shown in Figure 7, the constructed model describes quite well the response of the VLF signal amplitude to the initial sharp increase of the solar X-radiation flux. It should be noted that the starting amplitude response to the flare that occurred on October 25, 2013 was described qualitatively and quantitatively quite correctly on 5 paths out of 6, but, as follows from the behavior of the radio wave amplitude, the subsequent relaxation of the medium differed markedly from the real one. At the same time, within 2 h after the disturbance, the theoretical value of the amplitude reaches the experimental value.
As for the right panels in Figure 7, the X1.5 flare (12:36 UT), which occurred on June 10, 2014, is described worse than the X2.2 flare (11:36 UT) on all paths. It is most likely that this behavior is associated with incorrect accounting for the relaxation processes of the medium after the first perturbation. Thus, as a consequence, the D-region ionization during the second flare is calculated with inaccurate initial conditions.
The following work will include further development of the ionization model aimed at determining the correct account of ultrahard X-radiation, rectification of the relaxation processes and rate constants and advanced verification for all available paths and transmitter frequencies.
In addition to qualitative evaluation of the results obtained, we analyzed the quantitative difference in [dB] between the simulated and experimental values for the radio wave amplitude during the flares considered. For this purpose, a function D(t) = A theor (t) -A exp (t) was analyzed for 4 days including the days with solar flares listed in Table 2.
As quality criteria for the obtained results, we selected two parameters: the root mean square error (RMSE) and the D(t) module integral rationed by the integration time interval (from start to end of the particular flare), as given in formulas 10, 11: It is obvious that the smaller the value for these parameters, the better the model describes the measured values for the signal amplitude.  Table 1 and the frequencies at which they operate. The red, yellow, and green colors show the lines corresponding to the X-, M-, and C-class flares, respectively. RMSE, root mean square error.
The obtained criteria provided in Figure 8 confirm the conclusion made for qualitative evaluation of the agreement between the theoretical and experimental results; generally, the smaller the class of simulated flare, the better it is described with use of the constructed model. There are some exceptions. For example, the values for the lowest criteria correspond to M-class flares on the DHO-Mikhnevo path.
In addition to the dependence on flare energy, one can notice some correlation of the results with the transmitter frequency. For example, simulated data for the paths from lower frequency transmitters: GBZ, ICV, FTA, and GQD have a better agreement with the experimental data for the X-class flares. At the same time, there is a significant difference in the path lengths for the transmitters: GBZ, ICV, FTA, GQD and those for DHO and TBB. The latter two show the worse agreement with the experiment. Thus, in terms of modeling high-energy flares, the issue of correct calculation of the electron concentration profile defining the geometry of the waveguide remains undecided.
As a result of the verification using data from ground-based radiophysical measurements taken at the Mikhnevo geophysical observatory, it is found that the average normalized RMSE is ∼7%, 14%, and 18% for C-, M-, and X-class flares, respectively. Qualitative and quantitative analysis of the verification results for the VLF signal amplitude shows good predictive capability for the model built for describing weak and moderate ionospheric disturbances generated by solar flares.

Conclusions
The verification results obtained for a model simulating the ionospheric D-region carried out according to ground-based radiophysical measurements taken at the Mikhnevo geophysical observatory prove that application of satellite data obtained for the neutral atmosphere and ionization rates calculated on the BEKKER ET AL.
10.1029/2020JA028767 12 of 15 basis of actual radiation flux values within the X-ray and ultraviolet ranges allows one to describe the lower ionosphere during solar flares with adequate accuracy. Additionally, it should be emphasized that a year of high solar activity is characterized by a C-class flare level for the background X-ray flux. This means that additional ionization of the lower ionosphere by X-ray radiation must be taken into account at all times rather than only at the time of an actual large flare.
When simulating the ionosphere under the condition of X-class flares, we obtained contradicting results. In our opinion, the main reason for errors occurring for modeling of the D-region during a high-level of ionization is the unaccounted for ultrahard X-radiation. This leads to inaccurate calculation of the electron concentration profile, and, consequently, errors in determining the absolute value of the amplitude response to a sharp growth in ionization. Consequently, this affects the accuracy in describing environmental relaxation after the disturbance. At the same time, we note that the model recovers the background amplitude, varying by ∼0.5 dB from the experimental value within 1-2 h after a X-class flare event.

Data Availability Statement
The VLF data set used in the present work was collected by the Mikhnevo observatory of Sadovsky Institute of Geospheres Dynamics of Russian Academy of Sciences (http://geospheres-dynamics.com/). The experimental atmospheric data obtained by the AURA satellite are available at https://disc.gsfc.nasa.gov/. The satellite data for the radiation flux obtained by the GOES and SDO satellites are available at http:// satdat.ngdc.noaa.gov/sem/goes/data/ and http://lasp.colorado.edu/eve/data_access/evewebdataproducts/. The MSIS model profiles can be obtained under https://ccmc.gsfc.nasa.gov/modelweb/models/msis_vitmo. php. All the other data necessary to reproduce the reported findings are available at https://doi.org/10.5281/ zenodo.4055274.