Low Ionosphere under Inﬂuence of Strong Solar Radiation: Diagnostics and Modeling

Featured Application: The presented data can be used in practice in different areas of science and in several possible ways: for an analysis of the Earth atmosphere; for the investigation of coupling atmospheric electricity with biological systems; for the identiﬁcation and classiﬁcation of solar X-ray ﬂares; for obtaining the daytime atmosphere parameters induced by extreme solar radiation; and for understanding and preventing extreme space weather. Abstract: Solar ﬂares (SFs) and intense radiation can generate additional ionization in the Earth’s atmosphere and affect its structure. These types of solar radiation and activity create sudden ionospheric disturbances (SIDs), affect electronic equipment on the ground along with signals from space, and potentially induce various natural disasters. Focus of this work is on the study of SIDs induced by X-ray SFs using very low frequency (VLF) radio signals in order to predict the impact of SFs on Earth and analyze ionosphere plasmas and its parameters. All data are recorded by VLF BEL stations and the model computation is used to obtain the daytime atmosphere parameters induced by this extreme radiation. The obtained ionospheric parameters are compared with results of other authors. For the ﬁrst time we analyzed physics of the D-region—during consecutive huge SFs which continuously perturbed this layer for a few hours—in detail. We have developed an empirical model of the D-region plasma density and gave a simple approximative formula for electron density. propagates circle a mostly overland great circle investigated VLF signal changes and observed that the shape of the signal curve changes dramatically with time. Characteristic examples of SID event that happened during occurrence of strong SF is presented Figure 2. SF October 2005 from 09:42–10:08 UT with a peak at 09:59 The dramatical changes of GQD amplitude and phase values


Introduction
In today's science, special attention is given to the extreme weather events, climate change, preservation, and protection, because they have been identified as crucial for sustainable development in our century. Consequently, a very important question nowadays in modern society is: can we predict the magnitude of impact of explosive solar events (such as solar flares) on Earth, humans, electronic equipment, and nature in general? By analyzing past instances of these phenomena, we cannot predict the occurrence of each new event itself with certainty, but we can statistically estimate the consequences of phenomena on ionospheric parameters, perhaps predict damage to electronic equipment, predict disruption of GPS, etc.
The ionosphere is a part of the atmosphere that contains ionized gases with various kinds of particles, ions, and tens of different species [1][2][3][4] with the degree of ionization which depends mainly on the incident radiation that is solar extreme ultraviolet and X-ray radiation [5,6]. This ionizing solar radiation depends on the Sun activity during the solar cycle and varies around order of magnitude during the solar minimum and the solar maximum. Furthermore, ionosphere as a huge segment of atmosphere has a tendency to be constantly separated in different regions D, E, and F, with different physical characteristics and chemistry [7][8][9][10][11] which depend on the incident radiation. During stronger solar activity, the violent solar radiation and activity can induce sudden ionospheric disturbances (SIDs) as well as affect electronic equipment on the ground and signals from space, disrupt GPS, and bring damage to power grid components as a consequence [12][13][14][15]. Furthermore, many investigations indicate potential connections between this extreme activity and natural disasters such as tropical cyclones or forest fires (see e.g., [16][17][18]).
Signals of very low frequency (VLF) in the narrow band (3-30 kHz) of electromagnetic (EM) spectrum are generated in Earth's atmosphere by natural sources such as lightning discharges and meteor showers, and likewise VLF EM can be created by man-made sources, e.g., radio emitters [19,20]. The sensitivity of propagation of man-made VLF signals in the lower ionosphere makes them an ideal tool for remote monitoring of this medium under unperturbed/perturbed conditions [9,21]. During solar flares (SFs) and consequently SID events, an increase in ionospheric electron concentration and other ionospheric parameters at all altitudes is noticeable. As a result of radiation effects, the solar-induced SID and plasma irregularities cause perturbations in the received amplitude and phase of VLF EM signals mainly in the D-region which is part of the ionosphere.
In this contribution, we focus on amplitude and phase data of radio signals recorded by BEL system (Belgrade, Serbia) [22][23][24]. Emitted VLF signals [25] are constantly recorded by the equipment at BEL site since 2003. Intensity of X-ray flux during SFs, which are monitored by Geostationary Operational Environmental Satellite (GOES), are simultaneously analyzed with recorded data by BEL system of receivers. For these events, amplitude and phase perturbations are measured and analyzed during the daytime.
The aim of this study is to accurately model the perturbed D-region. We determine ionospheric parameters during huge SFs which may continuously perturb this region for several hours. Furthermore, for the perturbed ionospheric D-region, we give a simple approximative formula for altitude-dependent electron density which makes results easy for use. The results presented here are ready for further use with a particular emphasis on the astrophysical applications, more precisely in the rapid investigation of the low ionosphere under the influence of strong solar radiation.
The paper is structured as follows: in Section 2 we describe techniques for monitoring the D-region and models for determination of electron densities; in Section 3 we present and discuss the results of the investigation; and in Section 4 conclusions are presented.

VLF Technique for Monitoring Low Altitude Ionosphere
At ionosphere altitudes from 50 km up to 90 km (D-region), measurements mainly rely on EM wave propagation techniques [26]. In our study, synchronal monitoring of phase and amplitude of VLF narrow band radio signal data by the receivers during SFs occurrences were used to obtain daytime ionospheric parameters. For this reason, the perturbation of the VLF phase and amplitude was estimated as a difference between values of the perturbed signal induced by solar X-ray flare and signal in the normal condition.
The analysis of VLF data was carried out simultaneously with the examination of the correlative solar X-ray fluxes collected from GOES satellite (see e.g., Figure S1). The intensity of X-ray irradiance was recorded by GOES satellite by X-ray sensors in the band 0.1-0.8 nm. We based our study on investigation of series of stronger flares, i.e., flares of C-, M-, and X-classes. For our study, the effects of A-and B-class flares on VLF signals are rarely recorded and were non-essential.
During ionospheric disturbances, a standard numerical procedure for obtaining ionospheric daytime parameters and electron densities was used. It was based on a comparison of the recorded changes of amplitude and phase with the corresponding values obtained in simulations by the LWPC numerical software package [27,28] and Wait's two-component exponential model, as presented in [26,29].

Simulated VLF Signals
The LWPC program package [30] can be used for simulation of VLF propagation along any particular great circle path (GCP) under various diurnal, seasonal, and solar cycle variations in the ionosphere. This program package commonly executes the calculations for all possible modes and was tested with experimental data. In the published papers, there are few models for simulation of VLF propagation in the Earth-ionosphere waveguide (EIWG), such as the FDTD method, modefinder, etc. [31,32]. These models have been applied in many studies for calculations of various ionospheric parameters when characteristics and the ionospheric state allow applying some approximations.
In this contribution, by applying the LWPC code, the propagation of VLF signal was simulated for normal ionospheric condition, with the aim to estimate the best fitting pairs of Wait's parameters H' nor and β nor (reflection height and slope of sharpness, respectively; nor denote normal ionosphere condition) to obtain values of VLF amplitude and phase as close as possible to the recorded data for a selected day [33].
The next step was to simulate propagation of VLF signal through EIWG in the perturbed lower ionosphere during the flare event induced by an increased X-ray flux. For this step, we needed monitored, i.e., observed, VLF data to examine the signal perturbations during the solar X-ray flare.
We selected the pair of H' and β and used it as input parameter in the LWPC code to obtain the VLF signals, i.e., the simulated data of amplitude and phase at the receiver site. The process was rerun until differences between the simulated signal (∆A sim and ∆P sim ) for the disturbed and normal condition matched the measured ones ∆A and ∆P. After accomplishing a fair matching: ∆A sim ≈ ∆A and ∆P sim ≈ ∆P, the corresponding pair of H' per and β per was finally obtained (sim means simulated, per denote perturbed condition). The corresponding pair of H' per and β per was used for calculation of the electron density at a function of height for moments during the solar flare event (see schematic presentation in Figure S2).
Electron densities were resolved from the observed VLF amplitude and phase perturbations by a trial-and-error method where density profile was adjusted until the simulated amplitude and phase paired with observed disturbed data, as explained in the above paragraph. In this way, by processing data, the acquired Wait's parameters H' per and β per were used for our further calculations.

Procedure for Determination of Electron Density: Two-Component Exponential Model
The procedure was based on the Wait's two-component exponential model [34], which uses time-dependent parameters H' with values 74-50 km and β with values 0.3-0.6 km −1 , to describe the electron density height profile in the lower ionosphere. We note that there is a lack of agreement for sharpness values obtained by different models and techniques.
Although the values of the β parameter are in the range between 0.3-0.6 km −1 , in some papers β values are higher i.e., 1.0 km −1 at daytime and~3.0 km −1 at nighttime conditions (see e.g., [35]). Such a deviation illustrates the general deficiency of precision in the research of the ionosphere D-region and the need for new and more precise methods and techniques for its investigation.
This two-parameter system provides a model of a vertically stratified ionosphere with an electron density profile Ne that increases exponentially with altitude: Here, Ne and h are given in m −3 and km, respectively. The model [34] was successfully used in previous comparisons of VLF measurements and theory (e.g., [9,33]) and we also use it in our work.
Important plasma quantities used for describing the characteristics of the D-region are height dependent conductivity parameter; ω r (h) defined as where ω p [s −1 ] is defined as the plasma frequency and ν(h) is the effective electronneutral collision frequency. Quantity H' is the reflection height, meaning height at which ω r (h) = 2.5 × 10 5 rad/s. The positive ion density can be approximated to be equal to the electron density at high altitudes where electron density has values larger than 10 8 m −3 , while at lower altitudes the positive ion density is around 10 8 m −3 or can be calculated based on charge neutrality.
As it was noted, the electron densities can be obtained from the observed VLF amplitude and phase perturbations by a trial-and-error method (see Sections 2.1 and 2.2) where density profile is adjusted until the simulated amplitude and phase (using LWPC code) match with observed data. In this way, the obtained Wait's parameters H' per and β per can be used in Equation (1) for further calculations.

Results
In this research we studied the amplitude (A) and phase (P) data, obtained by monitoring VLF signals emitted by worldwide transmitters during solar-induced SIDs. The VLF signals which we examined are so-called sky waves, which travel via Earth-ionosphere duct and fluctuate in real-time due to the physical conditions through duct [34]. On the contrary, so-called ground waves closely follow and interact with terrestrial surface, and are attenuating. This category of EM waves is not investigated here.
All the data were registered by two receiver systems at the Belgrade site: AbsPAL and AWESOME. The stations can monitor VLF signals under different time resolutions. These kind of monitored data must be independently methodized and differently processed. The receivers can simultaneously collect several signals emitted by different emitters (located at different countries and territories) at the fixed frequencies. This makes it possible to monitor huge regions of the Earth and enables detection of local and global perturbations. The quality of the monitored data depends on the characteristics of the emitters, i.e., the geophysical location and the emitted power. Furthermore, the quality of the recorded data depends on the length of the VLF route, physical properties of the Earth-ionosphere duct, electromagnetic pollution, etc.
The technicalities and description of the Belgrade site are presented in [22]. Location of emitters and the receiver site and GCPs of subionospherically propagating VLF signals are provided in Figure 1 and also in Table S1 in online material. Here, we present our results of amplitude and phase perturbations on the GQD/22.10 kHz radio signal observed over a long time period. Behaviors of amplitude and phase data on VLF radio signals induced by different intensity of solar X-ray fluxes and observed on Belgrade site were examined. These results are presented. The signal propagates from the GQD transmitter to the receiver site over the great circle distance of 1982 km along a short and mostly overland great circle path that covers two time zones (propagating from Anthorn (54.72° N, 2.88° W), UK to Belgrade (44.85° N, 20.38° E)-see Figure 1. The GQD transmitter is located in the UK, but a great segment of the radio signal propagates over path in the same time zone where the receiver is.  The signal propagates from the GQD transmitter to the receiver site over the great circle distance of 1982 km along a short and mostly overland great circle path that covers two time zones (propagating from Anthorn (54.72 • N, 2.88 • W), UK to Belgrade (44.85 • N, 20.38 • E)-see Figure 1. The GQD transmitter is located in the UK, but a great segment of the radio signal propagates over path in the same time zone where the receiver is.
We investigated VLF signal changes and observed that the shape of the signal curve changes dramatically with time. Characteristic examples of SID event that happened during occurrence of strong SF is presented in Figure 2. SF X3.6-class on 9 October 2005 lasted from 09:42-10:08 UT with a peak at 09:59 UT. The dramatical changes of GQD amplitude and phase values at peak time can be observed. Figure 1. The geographic location of Belgrade system of VLF receivers (blue circle) and the GQD transmitter (red circles) (54.73° N, 2.88° W) Anthorn UK with GCP of sub-ionospheric propagating VLF signal.
The signal propagates from the GQD transmitter to the receiver site over the grea circle distance of 1982 km along a short and mostly overland great circle path that cover two time zones (propagating from Anthorn (54.72° N, 2.88° W), UK to Belgrade (44.85° N 20.38° E)-see Figure 1. The GQD transmitter is located in the UK, but a great segment o the radio signal propagates over path in the same time zone where the receiver is.
We investigated VLF signal changes and observed that the shape of the signal curve changes dramatically with time. Characteristic examples of SID event that happened dur ing occurrence of strong SF is presented in Figure 2. SF X3.6-class on 9 October 2005 lasted from 09:42-10:08 UT with a peak at 09:59 UT. The dramatical changes of GQD amplitude and phase values at peak time can be observed. The monitoring and investigation of VLF data was carried out simultaneously with the examination of the corresponding incoming solar X-ray fluxes [36,37]. The intensity of X-ray flux Ix is recorded by GOES satellite in the band 0.1-0.8 nm, available from [38] and [36]. SFs can be classified in categories based on the peak emission recorded by GOES satellite in the 0.1-0.8 nm band, namely from A-, B-, C-, M-, and X-classes.
In the presence of SIDs, a standard numerical procedure for the estimation of plasma parameters is based on comparison of the recorded changes of amplitude and phase with the matching values acquired in simulations by the LWPC numerical software package [27,30] as explained in [26,29,39]. Furthermore, for details, see Section 2.

Recorded Variations of VLF Amplitude and Phase Data
In this paper we analyzed around 200 solar X-ray flares events starting from the year 2003. Intensity of those X-ray flares were in range from C1-to X17-class, 1 × 10 −6 Wm −2 to 1.7 × 10 −3 Wm −2 , respectively. For each flare, we measured amplitude and phase perturbations on VLF signal, except when the transmitter was off-air. The amplitude perturbations were estimated as a difference between values of the perturbed amplitude induced by X-ray flux and amplitude in the normal ionospheric condition ∆A = A per − A nor . Phase perturbations were estimated as ∆P = P per − P nor .
Data points of the observed GQD/22.10 kHz are shown in Figure 3-amplitude and phase perturbations as a function of solar X-ray intensity. The range of size in amplitude and phase perturbations varies for different X-ray solar flares. The amplitude perturbations of GQD/22.10 kHz radio signal are distributed between enhancement and attenuation, which depends on signal frequency and the angle of reflection in a function of intensity of solar X-ray flux (see also [40]). On the GQD-BEL path, amplitude decreases during occurrences for C-class solar flares (see Figure 3). The size of amplitude perturbations is proportional to the intensity of solar X-ray flux. A phase change on this path is complicated, displaying an increase or decrease during independent to SF-class. During M-class solar flares, amplitude perturbation is between rising and falling but there are more cases of rising, while during X class SF amplitudes always have positive enhancement. In order to visualize the distribution of the amplitude and phase data, we give Figure  4. In Figure 4, a rug plots of VLF data are displayed as marks along an axis. A small tick on the edge of the rug plot represents each individual observation. The influence of different classes of SF explains the shape of the rug plots. Some grouping of ticks around ΔA = −2, −1, −0.3, 0.8 dB can be noticed, which correspond to the influence of C-class of SF and around 3.4, 6, 7 dB, which is due to M-class and X-class SF, respectively. Similarly, for ΔP, a grouping of ticks around 3 and 5 deg can be noticed, due to C-class SF and also around −20 and 20 deg which correspond to M-class. Around 40, 60 deg, some grouping exists due to X-class SF. In order to visualize the distribution of the amplitude and phase data, we give Figure 4. In Figure 4, a rug plots of VLF data are displayed as marks along an axis. A small tick on the edge of the rug plot represents each individual observation. The influence of different classes of SF explains the shape of the rug plots. Some grouping of ticks around ∆A = −2, −1, −0.3, 0.8 dB can be noticed, which correspond to the influence of C-class of SF and around 3.4, 6, 7 dB, which is due to M-class and X-class SF, respectively. Similarly, for ∆P, a grouping of ticks around 3 and 5 deg can be noticed, due to C-class SF and also around −20 and 20 deg which correspond to M-class. Around 40, 60 deg, some grouping exists due to X-class SF.

Long Lasting Intense Solar Radiation: An Analysis of Strong SF
M3.7-class SF that occurred on 9 September 2017 is an illustrative example of longlasting intense solar radiation (~60 min) which induced SID and caused variation in VLF signal. Furthermore, this example is good for showing the methodology of our analysis.
In Figure 5, we present simultaneous variations of X-ray flux, amplitude, and phase of GQD/22.10 kHz signals against universal time during the occurrence of M3.7-class SF on 9 September 2017. The lower panel shows observed amplitude perturbations on GQD radio signals measured at the Belgrade station with calculated values (red line). Middle panel shows observed and calculated phase perturbations (red line). Time variation of Xray flux measured by the GOES-15 satellite on 9 September 2017 is on the upper panel. Shapes of the VLF signal parameters are observably modified similar to the shape of the X-ray flux.
We simulated propagation of VLF radio signal through the waveguide in the perturbed D-region induced by additional X-ray radiation for resolution of 1 min during the flare duration using a procedure which relies on LWPC code.

Long Lasting Intense Solar Radiation: An Analysis of Strong SF
M3.7-class SF that occurred on 9 September 2017 is an illustrative example of longlasting intense solar radiation (~60 min) which induced SID and caused variation in VLF signal. Furthermore, this example is good for showing the methodology of our analysis.
In The procedure is explained in detail in Section 2. Parameters were obtained when the calculated amplitude and phase perturbations matched as well as possible with observed data (see lower panels of Figure 5). The obtained Wait's parameters βper and H'per were used for our further calculations of ionosphere parameters.
In Figure 6, calculated parameters are shown (obtained by the method explained in Section 2): time dependent effective reflection height H', sharpness β, and electron density We simulated propagation of VLF radio signal through the waveguide in the perturbed D-region induced by additional X-ray radiation for resolution of 1 min during the flare duration using a procedure which relies on LWPC code.
The procedure is explained in detail in Section 2. Parameters were obtained when the calculated amplitude and phase perturbations matched as well as possible with observed data (see lower panels of Figure 5). The obtained Wait's parameters β per and H' per were used for our further calculations of ionosphere parameters.
In Figure 6, calculated parameters are shown (obtained by the method explained in Section 2): time dependent effective reflection height H', sharpness β, and electron density at reference height during occurrences of M3.7-class solar flare on 9 September 2017 from lower to upper panel, respectively. The variation of H'(t) with time has expected behavior, i.e., after starting SF, it drops dropped rapidly to minimum, as expected, and after peak time of X-ray, it slowly keeps rising to referent value (see Figure 6 middle panel). The behavior of β(t) plot is in correlation with changes of X-ray flux. At 11:04 UT, peak of the solar X-ray flux reached Ix = 3.78•10 −5 Wm −2 , and induced perturbations of amplitude and phase ΔA = 5.6 dB and ΔP = −5.6 deg, respectively. For the sharpness, we obtained values β = 0.49 km −1 and, for reflection, height H'= 62.5 km. This means that during this SF the reflection height drops by 11.5 km. We obtained that the electron density increased at reference height h = 74 km from 2.2•10 8 m −3 to 6.7•10 10 m −3 .

Intense Solar Radiation: Consecutive Strong SF
We analyzed data for 10 June 2014 as an example of a day with several strong and successive SFs that induced SIDs and dramatically affected the VLF radio signal. Specifically, on 10 June 2014, the Sun released several consecutive SFs including two major Xclass solar flares that occurred one after the other as shown in Table 1.

Time [UT]
SID VLF Data at Peak Time The variation of H'(t) with time has expected behavior, i.e., after starting SF, it drops dropped rapidly to minimum, as expected, and after peak time of X-ray, it slowly keeps rising to referent value (see Figure 6 middle panel). The behavior of β(t) plot is in correlation with changes of X-ray flux. At 11:04 UT, peak of the solar X-ray flux reached Ix = 3.78 × 10 −5 Wm −2 , and induced perturbations of amplitude and phase ∆A = 5.6 dB and ∆P = −5.6 deg, respectively. For the sharpness, we obtained values β = 0.49 km −1 and, for reflection, height H'= 62.5 km. This means that during this SF the reflection height drops by 11.5 km. We obtained that the electron density increased at reference height h = 74 km from 2.2 × 10 8 m −3 to 6.7 × 10 10 m −3 .

Intense Solar Radiation: Consecutive Strong SF
We analyzed data for 10 June 2014 as an example of a day with several strong and successive SFs that induced SIDs and dramatically affected the VLF radio signal. Specifically, on 10 June 2014, the Sun released several consecutive SFs including two major X-class solar flares that occurred one after the other as shown in Table 1.  Figure 7 presents simultaneous variations of X-ray flux, amplitude, and phase of GQD radio signal against universal time during five successive flares on 10 June 2014. From Figure 7, one can see that VLF and X-ray flux peaks happened simultaneously, i.e., intense radiation instantly disturbed ionosphere, creating SIDs and resulting in perturbed VLF signal. Figure 7 presents simultaneous variations of X-ray flux, amplitude, and phase of GQD radio signal against universal time during five successive flares on 10 June 2014. From Figure 7, one can see that VLF and X-ray flux peaks happened simultaneously, i.e., intense radiation instantly disturbed ionosphere, creating SIDs and resulting in perturbed VLF signal.
Accurately, variation of amplitude and phase on GQD signals during X-class SFs react as a well-defined enhancement that follows the X-ray radiation maximum. SID VLF changes at the time of the maximum of SF of X2.2-and X1.5-classes are ΔA ~ 7.3 dB, ΔP ~ 32° and ΔA ~ 7.4 dB, ΔP ~ 15°, respectively. During C-class SFs, perturbations of amplitude and phase are not well-defined but they exist. These C-class SFs with similar characteristics produce increased ionization in the ionospheric D-region, generating amplitude perturbations of GQD radio signal with a decrease in amplitude up to 2 dB with changes in phase up to ΔP = 10 deg.
For the first X-class flare, it should be emphasized that the development to the maximum flux of X-ray lasted 6 minutes, and after 2 minutes, the additional radiation was completed. The graph of measured phase changes is characterized by a turning point when the phase stops declining and begins to rise, but not so well-defined as for another event. During this X2.2-class flare, the height of reflection dropped from 74 km to 59.1 km following the increase in β from 0.30 km −1 to 0.53 km −1 .
During the time interval of 16 minutes, X-ray radiation of the second X-class flare reached a maximum. During the same interval, the shape of the phase changes is more noticeable with well-defined turning points on the graph of measured phase values as a function of time. Throughout this X1.5-class flare, the height of reflection dropped from 74 km to 60.85 km following the increase in β from 0.30 km −1 to 0.54 km −1   Accurately, variation of amplitude and phase on GQD signals during X-class SFs react as a well-defined enhancement that follows the X-ray radiation maximum. SID VLF changes at the time of the maximum of SF of X2.2-and X1.5-classes are ∆A~7.3 dB, ∆P~32 • and ∆Ã 7.4 dB, ∆P~15 • , respectively. During C-class SFs, perturbations of amplitude and phase are not well-defined but they exist. These C-class SFs with similar characteristics produce increased ionization in the ionospheric D-region, generating amplitude perturbations of GQD radio signal with a decrease in amplitude up to 2 dB with changes in phase up to ∆P = 10 deg.
For the first X-class flare, it should be emphasized that the development to the maximum flux of X-ray lasted 6 minutes, and after 2 minutes, the additional radiation was completed. The graph of measured phase changes is characterized by a turning point when the phase stops declining and begins to rise, but not so well-defined as for another event. During this X2.2-class flare, the height of reflection dropped from 74 km to 59.1 km following the increase in β from 0.30 km −1 to 0.53 km −1 .
During the time interval of 16 minutes, X-ray radiation of the second X-class flare reached a maximum. During the same interval, the shape of the phase changes is more noticeable with well-defined turning points on the graph of measured phase values as a function of time. Throughout this X1.5-class flare, the height of reflection dropped from 74 km to 60.85 km following the increase in β from 0.30 km −1 to 0.54 km −1 Figure 8 shows During occurrences of X2.2-and X1.5-class SFs, the variation of H'(t) with time has expected behavior: after SF starts, it drops rapidly to a minimum as expected, and after Xray flux peaks, it keeps rising to a referent value (see Figure 8 middle panel). Behavior of β(t) is correlated with the shape of recorded increase in amplitude. Practically during this SF of X-class, it rises rapidly to a maximum as expected after the peak time of X-ray flux and drops to a referent value. At 11:42 UT, during peak time of X2.2-class SF, with the solar X-ray flux Ix = 2.22•10 −4 Wm −2 , the reflection height lowers by 14.9 km to a value of H'= 59.1 km and the sharpness values rises to β = 0.53 km −1 . At peak time 12:52 UT, during X1.5-class SF with the solar X-ray flux Ix = 1.55•10 −4 Wm −2 , the reflection height lowers to H'= 60.85 km and the sharpness values rises to a same value as for X2.2. During C-class SFs, the variation of reflection height H' and the sharpness β are in correlation with X-ray flux. It can be noted that it is very difficult to analyze and model the D-region in the case of strong SFs. Therefore, the results related to ionospheric parameters of some authors have some deviations (see Section 3.5). An even more difficult case for modeling is when we have consecutive flares which are continuously perturbing the D-region. During consecutive huge SFs which repeatedly perturb this region for a few hours, the ionosphere lacks time for regeneration and returns to a normal condition. The flares follow each other, literally sticking together and amplifying their influence in the D-region. Therefore, we have the case that a weaker SF, i.e., X1.5-class, induces similar changes in amplitude as the stronger flare of class X2.2 (ΔA ~ 7.4 dB and ΔA ~ 7.3, respectively). It is also necessary to point out that the solar zenith angle and the duration of flare have an impact on the parameters of the D-region.
Parameters H' and β are used to determine the spatial structures of disturbances in the lower ionosphere. The change in distribution of ionization at the D-region due to flare existence can be illustrated by the electron density altitude profile changes. We present the height profile of electron density at peak time of five successive flares on 10 June 2014 in Figure 9. At unperturbed ionospheric conditions (black line), there is a moderate increase in electron density from 2 × 10 7 m −3 at 60 km height, to about 2 × 10 9 m −3 at the upper barrier of the D-region i.e., at 90 km height. These three lines have almost the same slope. During occurrences of X2.2-and X1.5-class SFs, the variation of H'(t) with time has expected behavior: after SF starts, it drops rapidly to a minimum as expected, and after X-ray flux peaks, it keeps rising to a referent value (see Figure 8 middle panel). Behavior of β(t) is correlated with the shape of recorded increase in amplitude. Practically during this SF of X-class, it rises rapidly to a maximum as expected after the peak time of X-ray flux and drops to a referent value. At 11:42 UT, during peak time of X2.2-class SF, with the solar X-ray flux Ix = 2.22 × 10 −4 Wm −2 , the reflection height lowers by 14.9 km to a value of H' = 59.1 km and the sharpness values rises to β = 0.53 km −1 . At peak time 12:52 UT, during X1.5-class SF with the solar X-ray flux Ix = 1.55 × 10 −4 Wm −2 , the reflection height lowers to H' = 60.85 km and the sharpness values rises to a same value as for X2.2. During C-class SFs, the variation of reflection height H' and the sharpness β are in correlation with X-ray flux.
It can be noted that it is very difficult to analyze and model the D-region in the case of strong SFs. Therefore, the results related to ionospheric parameters of some authors have some deviations (see Section 3.5). An even more difficult case for modeling is when we have consecutive flares which are continuously perturbing the D-region. During consecutive huge SFs which repeatedly perturb this region for a few hours, the ionosphere lacks time for regeneration and returns to a normal condition. The flares follow each other, literally sticking together and amplifying their influence in the D-region. Therefore, we have the case that a weaker SF, i.e., X1.5-class, induces similar changes in amplitude as the stronger flare of class X2.2 (∆A~7.4 dB and ∆A~7.3, respectively). It is also necessary to point out that the solar zenith angle and the duration of flare have an impact on the parameters of the D-region.
Parameters H' and β are used to determine the spatial structures of disturbances in the lower ionosphere. The change in distribution of ionization at the D-region due to flare existence can be illustrated by the electron density altitude profile changes. We present the height profile of electron density at peak time of five successive flares on 10 June 2014 in Figure 9. At unperturbed ionospheric conditions (black line), there is a moderate increase in electron density from 2 × 10 7 m −3 at 60 km height, to about 2 × 10 9 m −3 at the upper barrier of the D-region i.e., at 90 km height. These three lines have almost the same slope. Totally different behavior and slope have lines at peak times of X-class SF (dashed orange and pink lines) are already 2 × 10 9 m −3 at 60 km height. These variation in altitude profile of electron density are essential for the VLF propagation characteristics. Totally different behavior and slope have lines at peak times of X-class SF (dashed orange and pink lines) are already 2•10 9 m −3 at 60 km height. These variation in altitude profile of electron density are essential for the VLF propagation characteristics.    To enable the better and more adequate use of data, we give electron density results obtained using simple approximative formula, which is logarithmic and represented by a second-degree polynomial with height-dependent coefficients a1(h), a2(h), a3(h): To enable the better and more adequate use of data, we give electron density results obtained using simple approximative formula, which is logarithmic and represented by a second-degree polynomial with height-dependent coefficients a 1 (h), a 2 (h), a 3 (h):

Approximative Expressions
Here, Ix is solar X-ray flux (Wm −2 ), and h is height (km). Table 2 presents data of the selected fits. The whole table with smaller step in altitudes is available in the Supplementary material (see Table S2.) To obtain a formula that applies to all conditions and input parameters, we analyzed all investigated SF cases from this study and obtained coefficients. We have selected the most appropriate ones ( Table 2 and Table S2) that cover all the examples and can honestly satisfy the accuracy. Moreover, in order to check the validity of the formula, we presented the values of the approximative formula together with the calculations from this work and also of other authors in Figure S4. One can see that it is in good agreement with all relevant data. Surely a more complicated formula could have been chosen, but there was no need for such accuracy because it is known that the values of the D-region parameters can vary literally by an order of magnitude. The idea was to make the formula easy to use as well as to enable rapid calculations and analysis.
During this study, two hundred events of solar X-ray flares were examined. The results were obtained on the basis of VLF and satellite data and static analysis of numerous events. One of the main results is that we can obtain values of electron density at any height, against maximum intensity of X-ray flux by expression (3). Figure 11 shows observed amplitude and phase perturbations on GQD/22.10kHz radio signals, due to solar X-ray flares, measured at the Belgrade station (during 2003-2017), as a function of X-ray flux. Values of electron density at reference height h = 74 km, against maximum intensity of X-ray flux obtained by Equation (3), are presented. The color of the point is proportional to the electron density values. This makes it easier to analyze the D-region during different solar disturbances.
In modern times, one of the largest solar flare measured with instruments occurred on 28 October 2003 (X17.2-class, Ix = 1.72 × 10 −3 Wm −2 ). Coronal mass ejection (CME) was associated with the 28 October 2003 X17.2 solar flare. In Figure S3 in Supplementary Material, VLF and satellite data are presented. The approximative Equation (3) allows evaluation of electron density even for this extremely large SF. Lower panels in Figure S3 present perturbations of amplitude and perturbation of phase on the same day. At the peak time, amplitude changes ∆A = 7.2 dB and phase changes ∆P = 76 deg (it is the farthest point on the right in Figure 11). In Figure S3 height, against maximum intensity of X-ray flux by expression (3). Figure 11 shows observed amplitude and phase perturbations on GQD/22.10kHz radio signals, due to solar X-ray flares, measured at the Belgrade station (during 2003-2017), as a function of X-ray flux. Values of electron density at reference height h = 74 km, against maximum intensity of X-ray flux obtained by Equation (3), are presented. The color of the point is proportional to the electron density values. This makes it easier to analyze the Dregion during different solar disturbances.  Figure S3 in Supplementary Material, VLF and satellite data are presented. The approximative Equation (3) allows evaluation of electron density even for this extremely large SF. Lower panels in Figure S3 present perturbations of amplitude and perturbation of phase on the same day. At the peak time, amplitude changes ΔA = 7.2 dB and phase changes ΔP = 76 deg (it is the farthest point on the right in Figure 11). In Figure S3, the upper panel variations of X-ray irradiance are shown, as measured by GOES-15 satellite, as well as the evaluated electron density at reference height versus universal time UT on 28 October 2003 during occurrence of huge X17.2-class SF.

Comparison of Obtained Ionospheric Parameters
As it was said in Section 2, with the sharpness parameter β, which usually takes values 0.3-0.6 km −1 , and with a reference height H' with values in range 74-50 km, one can describe the electron density height profile, rate coefficients, etc. in the lower ionosphere. We noted that there is a lack of agreement for sharpness values obtained by different models and techniques. Although the values of the β parameter are in the range between 0.3-0.6 km −1 , in some papers β values are higher, i.e., 1.0 km −1 at daytime and~3.0 km −1 at nighttime (see e.g., [35]). Such a deviation illustrates the general deficiency of precision in the ionosphere D-region research and the need for new and more precise methods and techniques for its investigation. Therefore, we presented a comparison of obtained ionospheric parameters with the other authors.
The sharpness parameter β and reflection height H' (red squares) calculated in present work are plotted in Figures 12 and 13 together with results of other authors [41][42][43][44]. In Table S3, one can find the data of amplitude and phase perturbations of the VLF signal and the ionospheric parameters induced by different SF events which were also analyzed in this study. In Figure S4, we also give the electron density at reference height evaluated from the data presented in Figures 12 and 13. We can point out that the results of the sharpness parameter β and reflection height H' presented in this paper are in good agreement with the results of most authors found in the available literature. and the ionospheric parameters induced by different SF events which were also analyzed in this study. In Figure S4, we also give the electron density at reference height evaluated from the data presented in Figures 12 and 13. We can point out that the results of the sharpness parameter β and reflection height H' presented in this paper are in good agreement with the results of most authors found in the available literature.  It can be noticed that, when the intense solar flux increases during SF, values of the reflection height H' decrease. Looking at Figure 12, generally speaking, there is a good agreement in the values of the reflection height from the literature [41][42][43][44]. When the intensity of radiation increases from 10 −6 to 10 −3 Wm −2 (i.e., when we go from normal conditions of the ionosphere to the perturbed ionosphere due to strong X-class SF), the reflections go down from 74 km to 50 km. Figure 13 gives the simple dependence for H' with the well-defined slope Looking at the Figure 13, one can see a totally different situation for the results of Figure 13. Same as in Figure 12 but for the sharpness parameter.
It can be noticed that, when the intense solar flux increases during SF, values of the reflection height H' decrease. Looking at Figure 12, generally speaking, there is a good agreement in the values of the reflection height from the literature [41][42][43][44]. When the intensity of radiation increases from 10 −6 to 10 −3 Wm −2 (i.e., when we go from normal conditions of the ionosphere to the perturbed ionosphere due to strong X-class SF), the reflections go down from 74 km to 50 km. Figure 13 gives the simple dependence for H' with the well-defined slope where k 1 = 36.85032 and k 2 =−6.12059. A larger deviation can be noticed only in the results from Pandey et al., where a very strong flare of X-class with the Ix~10 −3 Wm −2 induce a very small decrease in the reflection height (H'~65 km), like the M3-class gives. Looking at the Figure 13, one can see a totally different situation for the results of sharpness parameter. Generally speaking, sharpness parameter β increases with the increase in X radiation intensity. However, the results are more scattered with the huge disagreement for the cases of strong X-class solar flares (right corner of the Figure) due to results of Gavrilov et al. [44]. These results of Gavrilov et al. deviate a lot from the slope obtained from the results of other authors. Excluding the results from Gavrilov et al., we obtain for the β a linear dependence, as in Equation (4), where now k 1 = 0.89756 and k 2 = 0.09654.
Generally, the results show that electron density can be obtained by processing the observed amplitude and phase perturbations with a specially developed method. More precisely, we obtained parameters H' and β which were further used to obtain the electron density for different classes of SF starting from usual to extreme. The parameters are in good agreement with the most of the existing data from the literature. As expected, the electron density time variation during solar-induced SIDs follows the variation of the registered solar GOES flux. As a conclusion, Ne varies in time and height (few orders of magnitude) during SFs.

Discussion and Conclusions
In this paper we analyzed around 200 solar X-ray flares events during the 2003-2017 period. The perturbed daytime VLF amplitude and phase data are obtained by recording the GQD/22.10 kHz radio signal for a GCP = 1982 km between the transmitter in Anthorn, UK and the receiver in Belgrade, Serbia. We inspect the magnitude of impact of thesẽ 200 SFs on Earth atmosphere and consequences of these explosive events. The VLF data and important Wait's parameters-during the enhancements of X-ray flux due to these SFs-are presented and discussed.
We would like to emphasize that, for the first time, we analyzed the D-region in detail; we obtained ionospheric parameters and electron density profile during consecutive huge SFs which continuously perturbed this region for a few hours. The results show that the impact of subsequent flares can be higher than expected and cannot be analyzed as a stand-alone case without analysis of the state of the ionosphere in the moments before their starts. Furthermore, the results of the ionospheric parameters presented in this paper are compared with the results of other authors in literature and we conclude that they are in good agreement with the findings of most authors. It can be noticed that the intense solar radiation, namely solar extreme events, lead to an increased electron production rate and can increase electron density up to a few orders of magnitude, depending on the flare intensity with distortion of the amplitude and phase VLF signal. In addition, we give a simple approximative formula for altitude electron density profile as a function of X-ray flux which is valid for nonperturbed and also for the perturbed D-region.
The results confirmed the advantageous usage of the presented method for investigation of solar-terrestrial coupling processes as well as for detection and analysis of space weather phenomena such as solar explosive events. The analyzed data can be used in studies of the radiation characteristics before its impact on the atmosphere, modeling of VLF radio propagation in the lower ionosphere, modeling of disturbed ionosphere, etc. We note that the research requires multidisciplinary approach together with application of various models [45,46].

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/app11167194/s1, Figure S1: X-ray flux, amplitude and phase of GQD signals during occurrence of M3.7 class solar flare on 09 Sep 2017, Figure S2: Schematic presentation of a method for atmosphere electron density modeling. Figure S3: X-ray flux, perturbation of amplitudes and perturbation of phases and the corresponding electron density at reference height on 28 Oct 2003 during X17.2 class solar flare, Figure S4: The electron density at reference height calculated for this study together with the results of other authors, Table S1: Data of amplitude and phase perturbations of VLF signal and ionospheric parameters induced by different SF events analyzed in this study, Table S2: VLF  receiving and transmitting sites, Table S3: Data needed for approximative function. Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: https://hesperia.gsfc.nasa.gov/goes/goes_event_listings/, accessed on 10 June 2021, https://satdat.ngdc.noaa.gov/sem/goes/data, accessed on 10 June 2021.