© P. Teysseyre et al., Published by EDP Sciences 2026

1 Introduction

The D-region is the lowest layer of the ionosphere, between 60 and 90 km. It regulates much of the High Frequency (HF, 3–30 MHz) wave absorption (Zawdie et al. 2017), a frequency range for civil, broadcasting and military use. The altitude of this region makes continuous satellite or balloon monitoring impossible. Instead, it relies on narrow Very Low Frequency/Low Frequency (VLF/LF, 15–45 kHz) transmitters, emitting waves propagating in the Earth-ionosphere waveguide. The propagating waves are affected by any changes to the waveguide’s boundaries (e.g. Wait & Spies 1965; Ferguson 1992; Thomson 1993; thomsonSPSUNDSCRlowlatitudeSPSUNDSCR2014; Teysseyre et al. 2025). Propagation models like LMP (Gasdia & Marshall 2021) or LWPC (Ferguson 1998) provide modeled VLF amplitude and phase, if a certain ionospheric electron density condition is assumed. Therefore, their comparison with measurements yields estimates of the electron density profiles of the D-region (Thomson 1993; McRae & Thomson 2000; Thomson et al. 2022; Xu et al. 2023). Solar flares produce increases in the Sun’s X-ray flux on timescales from tens of minutes to hours. This energetic radiation increases the ionisation of the neutrals down to the D-region, resulting in an increase of the electron density. HF absorption is stronger when the electron density is higher. Solar flares cause HF absorption to rise on short timescales (ten minutes to an hour). Fast and accurate real-time detection of solar flares is therefore crucial in mitigating their impact on radio operation.

Different flare detection systems already exist. Most issued alerts are based on solar X-ray flux measured by the Geostationary Operational Environmental Satellite (GOES) from the National Oceanic and Atmospheric Administration (NOAA). The NOAA’s Space Weather Prediction Center (SWPC) provides both post-flare reports and alerts when the solar X-ray flux attains an M5 (5 × 10⁻⁵ W.m⁻²) threshold (Rockville 2019). The reliability of GOES data and alerts stems from the satellites’ continuous operation and direct measurement of the Sun’s X-ray flux. It also presents a latency of four to six minutes (George et al. 2019; Hudson 2025), which is significant compared to solar flare onset timescales. Solar flares also cause UV and EUV emissions to rise, along with the Sun’s soft X-ray flux. Therefore, alternative alert systems use different wavelengths (e.g. Bonte et al. 2013; Kraaikamp & Verbeeck 2015; Lu et al. 2024). Those alerts detect 80–88% of M and X-class flares, 5 min before the peak flux on average (Kraaikamp & Verbeeck 2015). While those systems present robust algorithms allowing the detection of most flares, they rely on satellite data. Flare detection is made more complex in particular cases (e.g. the South Atlantic Anomaly; Bonte et al. 2013) and could be influenced by eruptive events like filaments (Lu et al. 2024). Neural network methods (Fernandez Borda et al. 2002; Qu et al. 2003) were also investigated, with improved detection rates and capabilities for real-time processing. However, those attempts similarly relied on satellite data.

Developing an alert system based on ground-based data presents several advantages. First, ground-based instruments are less vulnerable to strong space-weather conditions than satellites, which can suffer from electronic glitches, loss of position and communication. Building a resilient flare detection service thus involves collecting data from various sources. It also requires relying on diverse data analysis methods to act as backups should one fail. VLF in particular show potential for real-time flare alerts. GOES’s direct measurements of solar X-ray flux are invaluable, but they saturate at X17, contrary to VLF (Thomson et al. 2005). VLF measurements also offer a shorter data latency (less than one minute, compared to four to six minutes for satellites). The receiver’s low cost enables the deployment of worldwide networks. Note also that the VLF instruments allow the detection of Transient Luminous Events (e.g. Mika & Haldoupis 2008) and electron precipitations from the magnetosphere (e.g. Potemra & Rosenberg 1973; Clilverd et al. 2009), two other forcing sources of the D-region electron density. The aim is to maximise both the probability of a flare detection and the time between the detection and the flare peak. Two ground-based VLF alert systems were thus suggested. The first one relied on the Global Ionospheric Flare Detection System (GIFDS) network of VLF electric antennas (Wenzel et al. 2016; Banys 2017). This method depends on the capability of determining a quiet-day profile relying on the previous days, which may be an issue in active periods. George et al. (2019) suggested another flare nowcasting approach. An estimate of the Sun’s X-ray flux was provided from a linear regression of VLF phase data and additional parameters such as the solar zenith angle. However, the authors noted that this idea presented two limitations: the need to compute the start of the solar flare, and the reliance on near real-time GOES measurements. To the best of our knowledge, none of these methods are yet implemented for real-time applications. We therefore developed a real-time flare detection system based on VLF measurements alone, with the following aims:

1. Detect the occurrence of solar flares very early after its onset, and more crucially before its peak;

2. Estimate the Sun’s X-ray flux;

3. Estimate the duration of solar flares;

4. Estimate the HF absorption resulting from the flare and the impacted regions.

Here is presented methods and a tool to address the first two points. The latter two will be implemented in a future version of this tool. The preparation of the data (amplitude calibration and phase correction) is explained in Section 2. Section 3 details the flare detection method based on breakpoint detection and its performance. Section 4 addresses the question of solar flux estimation using a probabilistic approach. Section 5 presents the method for real-time estimation of the electron density in the D-region. Finally, Section 6 presents the Python package performing these different steps. In particular, a general workflow is presented, summarising the various steps described previously. Section 7 details some of the possibilities for improvement of the methods’ performances.

2 Instrumentation and calibration

Amplitude calibration and phase correction constitute the initial step in analysis. Amplitude calibration involves converting voltage into magnetic field strength. This study will examine data from an AWESOME receiver (Cohen et al. 2018) situated in Nançay, Sologne, France. The AWESOME instrument consists of two orthogonal air-core loop antennas. In Nançay, those loops are triangles with a 4.47 m base and a 2.28 m height, and made of 11 turns of copper wire. It is thus sensitive to magnetic fields above 357 fT.√Hz, with absolute timing accuracy of 15–20 ns and 96 dB of dynamic range. It records broadband (waveforms) and narrowband measurements which are used here. Narrowband data consist of amplitude and phase data within a 200 Hz bandwidth centred around a transmitter frequency, at a 1 Hz sampling rate. Four transmitters are used for the real-time analysis (Table 1).

The instrument’s response to a known signal is used to calculate the amplitude calibration factors through an in-built amplitude-calibration device. VLF transmitters encode communication signals typically with minimum shift keying (MSK) which alternates between one of two frequencies to communicate bits. The phase is obtained from the real-time demodulation of the signal (Gross et al. 2018). Phase correction involves unwrapping, then detrending the raw phase. Unwrapping compensates for spurious jumps introduced by constraining the phase between −180° and 180°. The phase presents an artificial slope coming from a small drift of the transmitter’s frequency with time. Post-processing analysis compensates for this slope by fitting a straight line to the data and eliminating the trend. Yet, this relies on previous data, which may not be immediately available if the VLF transmission is disrupted. Instead, another method is employed for real-time processing. The total phase is written as:

$$\begin{array}{c}\hfill \varphi ={\varphi }_0+\delta f\frac{2\pi }{360}t\end{array}$$(1)

where ϕ is the unwrapped phase measured and ϕ ₀ is the stable phase. δ f is the small shift in frequency causing the phase drift, and t is the time (in s) since the start of the day. Given δ f, ϕ ₀ is derived from ϕ with no need for past data. δ f depends on the transmitter and varies along the year, following adjustments from the operators (Gross et al. 2018; Cannon & Marshall 2025). It is approximated by unwrapping the phase on consecutive days and minimising the phase jumps between days. The values of δ f for the transmitters monitored at Nançay for the period from November 2023 to December 2024 are presented in Table 1. The users are invited to check those values monthly, as they may vary in time. Once the two-step correction is complete, the phase trend is accurate, though its absolute value is still determined by the first point of each day. As the subsequent analysis relies exclusively on phase variations, absolute phase calibration is not required.

To further mitigate noise, a moving median filter is applied to the VLF phase and amplitude, with a time-resolution t _r . In addition, a flag system is implemented to highlight usable periods. Flags denote intervals during which a transmitter is inactive, identified by a low amplitude signal that varies according to environmental noise levels. In Nançay, this noise is estimated to be 0.5 pT. Flags also eliminate nighttime periods, defined by the solar zenith angle value: nighttime starts when the solar zenith angle is above 85° at either the receiver or the transmitter. This value of 85° corresponds to the minima in the VLF data due to sunrise or sunset. This is a conservative value that satisfies all seasons. Since we consider the two extremes of the path, we are not constrained by its length or obliquity.

3 Flare detection

After applying corrections and flags, solar flares are identified. The phase always shows a peak during a flare, while amplitude can present successions of peaks and dips, depending on the distance between the transmitter and receiver (Briand et al. 2022). Thus, the detection is based on phase data, as it is better correlated with the solar X-ray flux (George et al. 2019). In post-processing, flares are identified through phase peaks and high solar X-ray flux. Our goal is to detect flares before their peak. The solution is to seek changes in the phase trend (breakpoints hereafter). The primary challenge in breakpoint detection lies in optimising the trade-off for detection sensitivity. A comprehensive detection of all breakpoints inherently risks the inclusion of spurious signals arising from instrumental or environmental noise. In contrast, focusing solely on large trend changes misses smaller events. Another limitation arises from the need for fast, low-delay computation. The incremental algorithm meets those criteria (Guralnik & Srivastava 1999) and is described in Section 3.1.

To test this method, a full year of data was processed, covering the period from 2023/11/01 to 2024/10/31. This period includes highly active months like May 2024 (which included numerous X-flares), as well as quieter periods such as November and December 2023. Additionally, this period includes instances when the transmitters were inactive, likely due to maintenance. To give more specific illustrations, the example of three consecutive flares on 2024/07/28, is taken hereafter. The NRK phase and the X-ray flux defining those flares are displayed in Figure 1. Hereafter, GOES X-ray flux will always refer to the soft X-ray band, even though hard X-ray are more crucial for strong flares (Briand et al. 2022).

Figure 1. Phase (Φ, in °, black curve) and soft X-ray flux (in W/m2, green curve) on 2024/07/28. Three flares are presented and taken as illustration for the reminder of the paper. They are denoted by their soft X-ray flux peak time: 10.71, 12.01 and 12.84 UT.

3.1 The incremental algorithm

Breakpoints are detected through the incremental algorithm. This method minimises the distance of the data points to a fit. We assume that a straight line correctly fits the phase, as it is constant during quiet periods. This approach ensures the fast computation required for real-time treatment, compared to using a more complex fit. This distance to minimise is measured by the likelihood criteria L:

$$\begin{array}{c}\hfill L=mlog\frac{R}{m},\end{array}$$(2)

with m the number of data points, and R the sum of residues of data points to the linear fit. It is assumed that different segments present different variances (i.e. heteroscedastic hypothesis). This was checked on data from 2024/07/28, where the variance of different data segments varies by about 30%.

The first data point is automatically declared to be a breakpoint. Afterwards, each new data point is a potential breakpoint. For each one, two hypotheses are examined:

1. It follows the previous phase trend. The likelihood criterion is then computed and denoted L _no change, representing the distance to the fit over the entire dataset.

2. This data point introduces a break in the phase trend. It is further assumed that this break occurs five data points prior, as five points are needed to accurately fit a straight line. This assumption ensures fast detection while guaranteeing that the new trend is accurate. Moreover, it limits the computation time arising from checking each data point. Last, the new global likelihood criterion is calculated by adding the local likelihood criteria from both segments (separated by the breakpoint) and is denoted L _change.

δ is a user-defined threshold; the new breakpoint is kept when it satisfies

$$\begin{array}{c}\hfill \frac{(L_{\mathrm{change}}-L_{\mathrm{no}\mathrm{change}})}{L_{\mathrm{no}\mathrm{change}}}>\delta \end{array}$$(3)

Then, the phase slope associated with the breakpoint is named s and is used to determine the state of the flare evolution (see Sect. 3.3).

3.2 Parameter tuning for the detection method

The detection performance is controlled by two parameters: δ and t _r (Eq. (3) and Sect. 2 respectively). The goal is to find the optimal values to detect the maximum of flares and minimize the false detections. The optimal values of δ and t_r were estimated by analysing data from the NRK transmitter (Table 1) over one-year period of reference. During this period, 390 M and 17 X-flares (407 total) occurred when the transmitter was on and the path remained entirely within daylight. Each flare above C3 was initially found using GOES X-ray flux measurements. Periods when NRK was off were identified during the flagging step (Sect. 2). Of the 407 flares, 271 have a noticeable response in the data. The others were too close to nighttime or to another more powerful flare. All X-flares were observed in the VLF phase measurement. Figure 2 shows an example of breakpoint detection from phase measurement of NRK performed on 2024/07/28. For each flare, a visual inspection confirmed the identification. In parallel, the phase maximum time is determined. A detection threshold is also defined as the 75th quantile of quiet-time slopes. A flare is ‘detected’ when the breakpoint’s slope exceeds this detection threshold.

Figure 2. NRK phase (black line, also shown in Fig. 1) and detected breakpoints (red crosses) for 2024/07/28. The right vertical axis gives the scale of the slope. Four periods are indicated: flare detection (red vertical lines), flare onset (orange area), flare decay (light purple) and quiet (white). The breakpoint classification is detailed in Section 3.3.

The forecast time is defined as the time difference between the flare detection and phase maximum. Optimal parameter values maximize flare detection and forecast time, while minimizing computation time. The forecast time is flare-dependent: the mean forecast time is presented in Figure 3. Detecting more flares was prioritized over a longer forecast time. δ = 0.1 ensures the maximum number of detected flare: this parameter value improves the number of detected flare by about 20. t_r was then determined to ensure a balance between this number and the forecast time. t_r is thus set to 60 s, a rate that optimizes the code running-time and makes it more manageable. With these parameters, 82.7% M and X flares visible in the training set were detected. All X-flares were detected. 24 of the 47 undetected flares occurred during the decay of a nearby flare. The other undetected flares were 8 faint flares (the phase increase Δ Φ is below 10°), one short-duration flare (less than 3 min), and 14 flares presenting a positive slope below the detection threshold. The optimal parameters values were determined from NRK data. However, they depend on both environmental noise and the difference in slope between flares and quiet times. Thus, they vary for each transmitter-receiver path. Nevertheless, as long as δ <  0.2 and t_r <  120 s, the detection performance are comparable (Fig. 3). These parameter values should thus hold for any transmitter-receiver path without high noise levels or weak phase slopes during flares.

Figure 3. Number of flares for which a positive phase slope was detected (left) and median time delay (right) of the detection method applied to one year of phase data from NRK for different combinations of δ and t r .

The goal is to reduce the detection time, that is the time between the flare start (when the VLF-phase slopes become positive) and the flare detection. The detection time depends heavily on the rise time (the time between the flare start and the phase maximum). A ratio of detection time to rise time greater than one indicates that the flare was detected after the phase reaches its peaks (Fig. 4). Most of the flares are characterised by a rise time shorter than 6.5 min. If the flares are constrained to isolated flares above M3, the ratio presents a median of 0.22 and a mean of 0.24 (red crosses in Fig. 4). Strong isolated flares are thus detected at one fourth of their rise time on average. Also, as anticipated, a longer rise time leads to an earlier detection.

Figure 4. Ratio between the detection time and the phase rise time (r) vs. phase rise time (in minutes).

To quantify the quality of the detection, we check the delay between the soft X-ray flux increase and the detection. The ionosphere has a slight lag in its response, as reflected in the time-delay from amplitude measurements (Appleton 1953; Basak & Chakrabarti 2013; Briand et al. 2022). The VLF-flare detection is on average 7.1 min after the start of an X-ray flux increase. However, the average time from VLF phase increase to flare detection is 2.6 min. For comparison, GOES alerts are triggered by M5 threshold with a latency of 4–6 min (George et al. 2019; Hudson 2025). Our method relies on the detection of the flare start for any flare strength.

3.3 Breakpoint analysis

At this stage, a breakpoint has been detected, but not classified. Therefore, each new breakpoint must be analysed to determine whether it is quiet time, flare onset, flare maximum, or flare decay.

The first step of this classification is based on the computation ΔΦ, the variations in phase during a solar flare. It only depends on the previous measured breakpoints. ΔΦ is calculated through:

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathrm{\Phi }=\mathrm{\Phi }-({\mathrm{\Phi }}_q+s\times (t-t_q))\end{array}$$(4)

Φ _q , s and t _q are respectively the phase measured at the last breakpoint, the slope, and time of the last breakpoint. Φ and t are the current phase and time. The breakpoint is then classified using the algorithm presented in Figure 5.

Figure 5. Breakpoint identification scheme. It starts from a “New breakpoint detection” (left side), and continue with the slope s sign, which will determine whether the breakpoint is due to a flare onset, flare maximum, flare decay or is still in quiescent periods.

After building a catalogue of all breakpoints from the one-year period described previously, they are sorted by ‘flare’ and ‘quiet period’. Comparison with GOES-X catalogue shows that 94% of our identification was correct. Let’s call ‘positive’ and ‘negative’ times categorised as flare or quiet from our identification. ‘True’ and ‘false’ refer to the correctness of our identification. ‘False’ means that the identification is wrong. Thus, a ‘true positive’ denotes a flare that has been correctly detected, while a ‘false negative’ represents a flare identified as a quiet period. Any undetected C flare was considered as quiet time, as those faint flares will not greatly impact the ionosphere; but detected C flares were considered ‘true positives’ since they induced an actual change in the phase data. Our identification includes 51% for true positives and 43% of true negatives. The incorrect identifications consist of 2% of false negatives and 4% of false positives. Some false negatives were due to flares producing a weak phase slope, and false positives were occasionally present during quiet periods.

Figure 2 presents the different detected breakpoints and their classifications on the 2024/07/28. The flare peaking at 10.71 UT (resp. 12.84 UT) is detected 12.6 min (resp. 4.1 min) before the phase peak, and 8.4 min (resp. 2.1 min) before the X-ray flux peak. The flare at 12.01 UT remains undetected, as the flare at 10.71 UT did not recovered its quiescent level when it occurred. Each period (onset, quiet and decay) is delimited by a breakpoint. Thus, the onset period of the flare is overestimated for both the M7 and M2 flare.

4 Solar flux estimation

4.1 Method

Once a breakpoint has been classified as a flare the solar X-ray flux is estimated from ΔΦ. Daily variability in the ionosphere depends on physical parameters (such as the X-ray spectrum, the solar zenith angle or the season). Thus, there is no bijection between ΔΦ and the X-ray flux. Instead, there is a dispersion of ΔΦ regarding the solar X-ray flux (Fig. 6, top panel and Thomson et al. 2005). This implies that a probabilistic approach must be employed instead of a direct derivation of the flux from the phase variation.

Figure 6. Probability density functions (PDF) for NRK. Top panel: ΔΦ vs. GOES X-ray flux. Bottom left panel: PDFs for discretized phase. Bottom right: PDF fitted by a normalized Gaussian distribution.

The probability P of the flux value being F is computed from Bayesian statistics, given ΔΦ₁ and ΔΦ₂ the phase variation for two transmitters (as measured with a single receiver), and background information I:

$$\begin{array}{cc}\hfill P(F|\mathrm{\Delta }{\mathrm{\Phi }}_1,\mathrm{\Delta }{\mathrm{\Phi }}_2,I)& \propto P(\mathrm{\Delta }{\mathrm{\Phi }}_2,\mathrm{\Delta }{\mathrm{\Phi }}_1|F,I)\times P(F|I)\hfill \\ \hfill & \propto P(\mathrm{\Delta }{\mathrm{\Phi }}_1|F,I)\times P(\mathrm{\Delta }{\mathrm{\Phi }}_2|F,I)\times P(F|I)\hfill \\ \hfill & \propto \frac{P(F|\mathrm{\Delta }{\mathrm{\Phi }}_1,I)}{P(F|I)}\frac{P(F|\mathrm{\Delta }{\mathrm{\Phi }}_2,I)}{P(F|I)}\times P(F|I)\hfill \end{array}$$(5)

P(F|ΔΦ₁, ΔΦ₂, I) is called the posterior density function, and represents the probability of the X-ray flux being F, knowing ΔΦ₁ and ΔΦ₂. Statistical independence between ΔΦ₁ and ΔΦ₂ is assumed to simplify this computation, and probabilities independent of F are constants. This is justified, since the variations of ΔΦ depend on the interactions of the different propagating modes in the Earth-Ionosphere waveguide, which is notably parametrised by the transmitter frequency and ground conductivity. The phase variations also depend on the solar zenith angle. In most cases, the transmitter-receiver paths are far enough that the solar zenith angles and ground conductivities are distinct (except for the GQD and GVT transmitters). This simplification can be checked with the ΔΦ data for various transmitters. However, the limited number of flares in the data set, particularly X-class ones, would make this calculation unreliable at the moment. This option remains available for later execution, pending further data acquisition.

Equation (5) can be generalised to any number of transmitters. Therefore, the X-ray flux is estimated from the prior P(F|I) and the successive conditional probability density functions associated with each ΔΦ _i . Since the prior represents the probability density of the flux without knowledge of the ΔΦ values, it is taken as uniform from 10⁻⁶ to 10⁻³ W.m⁻².

4.2 PDFs from a single transmitter

To compute the PDF in equation (5), it is necessary to discretise the values of ΔΦ and the X-ray flux. ΔΦ is computed from fully resolved M and X flares following equation (4), and discretized into five bins. If A is the maximum of ΔΦ, then the edges of these bins are 0, A/5, A/2, 2× A/3, 4× A/5 and A. A depends on each transmitter-receiver path. As A depends on each propagation path, the ΔΦ binning will depend on each transmitter and receiver. The X-ray flux values are similarly discretized, with bin edges M, M3, M6, M9, X2, X5 and X10. The choice of the bin edges is motivated by trying to limit the impact of the overabundance of weak M flares compared to the stronger X flares, by having larger bins for low ΔΦ values. However, it is also important to distinguish between M1 and M5 or M9 flares, especially since NOAA’s Space Weather scales for Radio Blackout change at M5 and X1.

The probability density functions thus computed are then fitted by normalised Gaussian distributions to smooth them. The normalisation first ensures that the density functions never go to zero, leaving a (small) probability that the flux differs from the best value computed from a single transmitter. It also facilitates the storage of these probability functions, since only the average and standard error are needed for each ΔΦ bin.

Figure 6 shows the application of the method to 183 isolated flares observed from NRK. Phase variations between 0 and 97° can hardly distinguish between M-flare strength. Furthermore, since the training only contained 13 X-class flares, the PDF are not well-prepared for very strong flare. We will discuss these limitations in Section 8.

PDF are computed for every transmitter monitored in Nançay (Fig. 7). For most of them, PDFs behave as expected: larger ΔΦ corresponds to higher X-flux. However, GQD phase changes do not strongly depend on the solar flare X-ray flux: the probability densities are almost similar for all ΔΦ values. Hence, transmitters behaving like GQD should not be included in the X-flux estimate. The phase behaviour is not an intrinsic property of the transmitter and must be verified for each propagation path. In Nançay, the flux estimation no longer relies on GQD data.

Figure 7. Probability density functions for the four transmitters monitored in Nançay. The PDFs shown in the bottom left panel for NRK were already present in Figure 6.

The quality of the PDF is limited by the dataset used for the training: the small number of flares above X5 limits their detectability. Indeed, the PDFs built on NRK data (bottom left panel, Fig. 7 show considerable variation if ΔΦ is above 116°: the 95% confidence interval is thus [M1.7, X3.6], due to there being only 13 X-class flares in the training period. Similarly, this confidence interval is [X1.0, X4.2] (resp. [M2.6, X9.3]) for NSY (resp. NAA), relying on 15 (resp. 5) X-class flares. This can be improved by generating synthetic data obtained by paring a chemistry scheme and a VLF propagation model. We will discuss this point in more details in Section 8.

4.3 PDFs from multiple transmitters

The PDF were built on data from November 2023 to October 2024, and the validation sample consists of data from November and December 2024. During those two months, two M2 and one M9 flares were detected by NAA, NSY and NRK. For each flare, the flux is evaluated using NRK only, NRK + NAA or the three transmitters (Fig. 8).

Figure 8. Estimation of the X-flux using one, two or three transmitters for: top panel, M2.9 flare (14.32 UT on 2024/11/05); middle panel, M9.5 flare (12.11 UT on 2024/11/10); bottom panel, M2.4 flare (13.14 UT on 2024/12/07). The posterior probability density functions obtained with equation (5) are indicated, combining information from one (black), two (brown) or three (red) transmitters.

As the transmitter number increases, the confidence interval reduces. In all but one case, the true X-flux is included within the 2σ confidence interval (shaded area) deduced from the VLF measurements. However, the estimate does not improve significantly with the combination of transmitters: for the M9 flare, the best flux estimate is even found using a single transmitter. This is partly due to the limited number of samples of intense flares available for training, which affects the phase grid discretization and the training process itself. Yet, maintaining several transmitters in the evaluation reduces the impact of a transmitter’s failure, reduces the number of false positives (Sect. 8), and provides a comprehensive overview of electron density.

5 Electron density estimate

Computing the D-region electron density is necessary to estimate the solar flare impact on HF absorption, and is achieved through the Longwave Mode Propagator (LMP, Gasdia & Marshall 2021). The electron density profile is described by the Wait profile (Wait & Spies 1964; Thomson 1993):

$$\begin{array}{c}\hfill \mathit{Ne}=1.43\times {10}^{13}exp(-0.15h^{\prime }+(\beta -0.15)(z-h^{\prime }))\end{array}$$(6)

h′ is the effective ionosphere height, β is the inverse of the typical gradient length and z is the altitude in the ionosphere, with z = 0 on the ground.

For each transmitter, amplitude and phase lookup tables are generated for the range of (h′, β) expected during quiet times to intense flares (Teysseyre et al. 2025, Fig. 1). Given a measured amplitude and phase, the values of (h′, β) minimising the difference between modelled and measured amplitude and phase are sought, according to the method in Teysseyre et al. (2025). Electron densities are thus deduced for each time step of an observed profile, from which maps of electron density along each path are generated (Fig. 9). These maps show each propagation path colour-coded for its electron density at 70 km in altitude. Note that in practice, the colour-code adapts to high electron densities without saturating at red as the right panel of Figure 9.

Figure 9. Maps representing the electron density at 70 km in quiet times (left) or at the peak of the M7 flare on 2024/07/28 (right).

6 The vlf4ions Python package

Based on the methods presented in previous sections, we designed a Python package vlf4ions. While designed for magnetic loop antennas, it is also suitable for VLF electric-antenna types. Open-source access is provided for the source code (version 1.3.6 at the time of this paper) and its documentation.

Figure 10 depicts the step-by-step functionality of this library’s methods. VLF data analysis using the methods from earlier sections yields solar flare detection, solar X-ray flux and D-region electron density estimates. The ‘Data correction’ and ‘Flag the data’ blocks correspond to Section 2, ‘Breakpoint detection’ is described in Section 3 and the ‘Nowcast time’ branch is detailed in Sections 4 and 5. Different branches are called at given times (specified by the user) to do separate tasks (see Fig. 10 from ‘Start’):

Figure 10. Diagram of the main steps of the real-time data correction, flare detection and electron density computation from the vlf4ions package. The different methods and results are detailed in the next sections.

• Newsletter time: A short email is sent to the VLF receiver maintainer as proof that the script is still running. A newsletter summarising daily flare detections (start time, end time and estimated peak flux) is also sent.

• Transmitter reading time: Each transmitter’s data is read to detect flares (with a frequency given by t _r ).

• Nowcast time: Evaluations of the solar X-ray flux and electron density profiles along each monitored propagation path are produced at the nowcast time (with a frequency given by t _r ).

Every necessary parameter (δ, t _r , the threshold for the flags and the thresholds for the flare detection) are given by the user in a separate script. Alerts are sent at each flare detection to various recipients. Each alert includes the solar X-ray flux estimation and the electron density profile along monitored paths. vlf4ions users can customise these alerts, as detailed in the documentation. Maintenance warnings are sent if the receiver is off so that the instrument maintainer knows to restart the instrument, and newsletters recapitulating each flare detection are crafted at the end of each day.

7 Improving the detection lead-time

The detection time is improved by decreasing the time between two data reading, δt: instead of 60 s, data can be read every 15 s which maximises the chance of detecting a flare early. For illustration, Figure 11 shows the difference in detection times when δt is 1 min or 15 s, using the data from three transmitters recorded in Nançay. Lowering δt to 15 s improves the detection time by 51 s, and the detection of the M5 threshold by 31 s.

Figure 11. Detection of the X5 flare occuring on 2025/11/11. The soft and hard X-rays fluxes are represented in black (full and dashed lines respectively) and the phase from NRK in green. The red vertical lines represent the first detections of this flare, and the blue vertical lines the moment at which the flux estimation reaches an M5 threshold. Full lines show the time detection with δ t = 60 s and dashed lines those with δ t = 15 s. The blue dotted line is the time of the issue of the SWPC alert for the M5 threshold.

Similarly, lowering the probability threshold at which the flare alerts are sent (Fig. 10) increases the chance to detect all flares and improves the lead-time. However, there is a trade-off with the false positive times, defined here as times when the flux estimation was above M5 (with a given probability threshold, P _t ). As an illustration, Table 2 below presents the percentage of time with false positives, for the months of November and December 2024. The only days considered for this study were those when all three transmitters were working, and without flares above M5, which represents 19 days in total. Neither changing δt nor P_t significantly impact the false positive time, though some variations are observed with the parameters. For the envisioned usage of this package, it is thus advised to choose δt = 15 s and P_t = 0.1.

8 Discussion and conclusion

vlf4ions is a Python package for VLF receiver owners that analyses VLF data in real-time. An incremental algorithm detects solar flares in real-time by seeking trend breaks in VLF phase data. 82.7% of M and X flares are detected within one fourth of their rise times. Solar X-ray flux estimate is performed using combined probability density functions derived from the phase variations of several transmitters. The computation is accurate within a factor of two with the current PDF built, and as assessed on our validation dataset. Finally, the D-region’s electron density profile is estimated by comparing VLF changes to those predicted by models like LMP.

Note that the method presented here detects changes in the D-region ionisation. Those were attributed to solar flares in the context of this work. However, any change in D-region ionisation will be detected if it is intense enough to affect the VLF phase data. This also implies that flares that do not greatly impact the D-region will be missed. This includes both faint flares and flares occurring close to sunrise or sunset, which do not affect the D-region very much probably due to their very low solar zenith angle.

This approaches complements the flare detection systems already in place. Compared to the previous suggestions of a VLF ground-based flare detection systems, the main improvement is the reliance solely on the instrument’s short term data. Thus, access to GOES X-ray flux or to past days of data is not required, which increases the resilience of these methods. We also provide estimates of the D-region electron density and the Sun’s X-ray flux. This method is adaptable to other VLF receivers; the steps to achieve this are presented in the package’s documentation.

The most recent versions of this package (v1.3.4 to v1.3.6) have been running in Nançay since 2025/07/29. Previous versions were tested since April 2025. A computationally efficient method is imperative to work in real-time on the instruments’ computers, without disturbing the data recording. The data analysis for a single transmitter is completed in less than one second. Both flux estimate and electron density computations require a maximum of two seconds. Consequently, the method demonstrates suitability for real-time application.

Various improvements to the flare detection and X-flux estimate are still possible. The analysis still relies on past data for the flux estimate. An alternative solution for deriving the different probability density functions is to employ synthetic data provided by a chemistry model. Indeed, electron and ion profiles can be derived from various chemistry models, for example WACCM-D (Verronen et al. 2016), GPI (Glukhov et al. 1992; Lehtinen & Inan 2007) or other chemistry schemes (e.g. Mitra & Rowe 1972). These models are able to reproduce the ions species and electrons for different forcing sources. The densities can then be inputted into a propagation code for different forcing conditions to estimate the VLF wave time profile both in amplitude and phase. This would have a second crucial improvement: it would allow a finer discretization in ΔΦ. It would also provide an extension of the probability density functions towards C-flares, which were not included here in the pdf as measuring and validating their ΔΦ values is very time-consuming. It thus would be unrealistic to perform this task every time this method is applied to a new antenna. In addition, it would extend the flux estimation to strong X-flares, which were not included in the pdf due to the lack of data. Finally, this flux estimate running in real-time depends only on the values of ΔΦ for each transmitter. Thus, this method works similarly in post-processing.

Phase measurement was preferred in the present version of the code, as it presents a simpler response to flares compared to amplitude, which exhibits more diverse responses (e.g. Briand et al. 2022). Yet, inverted or through profiles also carry information on the flare strength. Seeking negative slopes in the amplitude data may offer a different solution for flare detection and flux estimation. Moreover, relying on ΔA as well as ΔΦ for the flux estimation would reduce the errors and produce estimates with better confidence. There is significant work to achieve this, since the amplitude response to solar flares is more diverse than the phase’s. This is thus planned for a future version of the tool.

Phase measurements also require careful correction. An alternative to them would be the use of polarisation parameters (Gross et al. 2018). Indeed, three out of the four polarisation parameters are unaffected by the phase drifts. Using the tilt angle as well as major and minor axis would provide three quantities on which to base the flare detection and subsequent steps.

In this implementation of the electron density computation, the ionosphere was assumed to follow the Wait profile, as the most commonly-used model for the D-region. It also guarantees a fast computation-time and an easy comparison with the existing literature (e.g. Thomson et al. 2011, 2022). This two-parameter model however misses some fine features of the actual electron density profile during flares (e.g. Gasdia & Marshall 2023). This could be improved in later versions of the vlf4ions package by relying on pre-computed profiles from GPI-4 or GPI-5 instead.

Since the slopes in the phase data are detected, a short-term forecast of a few minutes may also be produced by extrapolation of the phase slope on short timescales. Thus, alerts could be sent sooner, although this may result in more false positives. This will be studied and implemented in a future version of the package.

The performance of this system is quantified using the alerts issued by the SWPC if the X-ray flux reaches M5¹. It should be emphasised that these alerts are not sent at the beginning of the flares and only concern flares above M5. On average, there are 11.3 min between the flare alerts obtained from VLF data and the alerts issued by the SWPC (the median is 5.2 min, for 39 flares). However, at this point, alerts from VLF data do not distinguish between weak M flares and stronger ones. Though there is not enough data yet to systematically compare the alerts at the M5 level with our system and from GOES data, it is presented for a specific example. For the flare in Figure 11, the M5 alert was issued first with VLF data with δ t = 15 s. Then, 23 s later (resp. 36 s), it was sent by SWPC (resp. our system, with δ t = 60 s). The detection of this M5 flare is thus equivalent between the different systems. Note that the ionosphere reacts to the flare with a small delay (e.g. Appleton 1953; Basak & Chakrabarti 2013; Briand et al. 2022). However, this delay decreases with flare strength. Indeed, VLF data is sensitive to hard X-rays, as soft X-rays do not penetrate in the D-region (Ratcliffe 1960, p. 186, Fig. 25). The stronger the flare, the earlier the hard X-ray peak occurs compared to soft X-rays. The VLF phase will thus react faster for strong flares, providing earlier alerts.

The main improvement to this system will be the computation of the HF absorption caused by the flares. Since the electron density is already estimated, this could be achieved through a different chemistry model or by implementing tabulated ion profiles and collision frequencies. This would allow a single receiver to produce HF absorption estimates over an entire continent or ocean, without relying on satellite data.

Acknowledgments

The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

Funding

The CNRS supported the present work through its PNST/ATST program, as well as the CSAA. The CNES also supported this work through its support to Space Weather programs.

Conflicts of interest

The authors declare no Conflict of Interest.

Data availability statement

The open-source code and documentation for the vlf4ions library are found at https://codeberg.org/pteysseyre/vlf4ions_library. At the time of the paper, the library’s version is 1.3.6. The VLF data from the AWESOME receiver in Nançay is accessible by emailing the authors.

References

All Tables

All Figures

	Figure 1. Phase (Φ, in °, black curve) and soft X-ray flux (in W/m2, green curve) on 2024/07/28. Three flares are presented and taken as illustration for the reminder of the paper. They are denoted by their soft X-ray flux peak time: 10.71, 12.01 and 12.84 UT.
In the text

	Figure 2. NRK phase (black line, also shown in Fig. 1) and detected breakpoints (red crosses) for 2024/07/28. The right vertical axis gives the scale of the slope. Four periods are indicated: flare detection (red vertical lines), flare onset (orange area), flare decay (light purple) and quiet (white). The breakpoint classification is detailed in Section 3.3.
In the text

	Figure 3. Number of flares for which a positive phase slope was detected (left) and median time delay (right) of the detection method applied to one year of phase data from NRK for different combinations of δ and t r .
In the text

	Figure 4. Ratio between the detection time and the phase rise time (r) vs. phase rise time (in minutes).
In the text

	Figure 5. Breakpoint identification scheme. It starts from a “New breakpoint detection” (left side), and continue with the slope s sign, which will determine whether the breakpoint is due to a flare onset, flare maximum, flare decay or is still in quiescent periods.
In the text

	Figure 6. Probability density functions (PDF) for NRK. Top panel: ΔΦ vs. GOES X-ray flux. Bottom left panel: PDFs for discretized phase. Bottom right: PDF fitted by a normalized Gaussian distribution.
In the text

	Figure 7. Probability density functions for the four transmitters monitored in Nançay. The PDFs shown in the bottom left panel for NRK were already present in Figure 6.
In the text

	Figure 8. Estimation of the X-flux using one, two or three transmitters for: top panel, M2.9 flare (14.32 UT on 2024/11/05); middle panel, M9.5 flare (12.11 UT on 2024/11/10); bottom panel, M2.4 flare (13.14 UT on 2024/12/07). The posterior probability density functions obtained with equation (5) are indicated, combining information from one (black), two (brown) or three (red) transmitters.
In the text

	Figure 9. Maps representing the electron density at 70 km in quiet times (left) or at the peak of the M7 flare on 2024/07/28 (right).
In the text

	Figure 10. Diagram of the main steps of the real-time data correction, flare detection and electron density computation from the vlf4ions package. The different methods and results are detailed in the next sections.
In the text

Figure 11. Detection of the X5 flare occuring on 2025/11/11. The soft and hard X-rays fluxes are represented in black (full and dashed lines respectively) and the phase from NRK in green. The red vertical lines represent the first detections of this flare, and the blue vertical lines the moment at which the flux estimation reaches an M5 threshold. Full lines show the time detection with δ t = 60 s and dashed lines those with δ t = 15 s. The blue dotted line is the time of the issue of the SWPC alert for the M5 threshold.

In the text