LIR-ACheM: Modelling of the D-Region Response to Solar Flares

1. Introduction

The D-region is the lowest ionosphere layer, situated between 60 km and 90 km. This region is principally composed of $N{O}^{+}$ resulting from the ionisation of nitric oxide ($NO$) by solar Lyman-$\alpha$ radiation (e.g., [2]). A secondary source of ions comes from the ionisation of $N_2$, O and N by solar X-rays (XR, 0.05 - 0.4 nm). During solar flares, the X-ray fluxes may increase by several orders of magnitude in a few minutes. This surge of X-rays becomes the dominant ionisation source, thus increasing the D-region’s electron density by photoionisation. Modelling the D-region presents some challenges ([3], p. 173). First, it is dominated by neutrals, and the variety of species allows diverse ions to form, multiplying the number of possible interactions. Contrary to the upper layers, the relatively high pressure in the D-region also allows three-body reactions to occur, which further increases the chemistry complexity of this region. Distinct ionisation sources are responsible for the ionisation of this region: EUV, Lyman-$\alpha$ and X-ray fluxes, cosmic rays (CR) at low altitudes and electron precipitations at high latitudes (e.g., [1,4,5,6,7]). They operate on different timescales (e.g. minutes to hours for XR during solar flares, while CR and electron precipitations vary on the timescales of hours for magnetic storms or years for the solar cycle), at different altitudes (Lyman-$\alpha$ ionise above 65 km, while CR impact more efficiently at lower altitudes), and on different species (Lyman-$\alpha$ on $NO$, XR on O, N and $N_2$). This adds to the complexity of this region and the difficulties in its accurate modelling.

Multiple chemistry schemes exist, differing in their treatment of the ion species. GPI models [8,9] group ion species into categories such as negative ions and positive cluster ions. Though these models are efficient, the use of ion categories removes the possibility of identifying the role of each ion in the evolution of the ionosphere. In addition, some chemical reactions are neglected. In particular, they do not consider the reaction from $O_4^{+}$ to $O_2^{+}$ that decreases the effective recombination coefficient at 80 km [1], resulting in longer recovery times. On the other hand, more complete models, such as SIC [10] or WACCM-D [11] consider dozens of ion species and hundreds of reactions, allowing for more precise estimates of the D-region chemistry, but requiring longer computational time and implying knowledge of some minor species (e.g. $Cl$, $N{O}_2$). Alternative schemes exist, among which the Mitra-Rowe scheme [1,12,13] that describes the ionosphere with six ion species and the electrons, and was created for solar flare perturbations.

We present here a new model based on the Mitra-Rowe scheme described in Section 2. The analysis of the D-region sluggishness and recovery time after a short duration M1 flare illustrates the capabilities of the code (Section 3).

2. Chemistry Modelling of the D-Region

The D-region is often described by the Wait profile [14,15]:

$$N_e\left(z\right)=1.43\times {10}^{13}exp(-0.15h^{\prime })exp\left((\beta -0.15)(z-h^{\prime })\right),$$

(1)

with $h^{\prime }$ the D-region equivalent height and $\beta$ the density gradients. This two-parameter representation of electron density is easy to integrate into propagation codes (such as LWPC [16] or LMP [17]), improves computational efficiency, and facilitates interpretation. Nonetheless, this simplicity comes at the expense of losing subtle features within the D-region and lacks a detailed representation of the underlying chemical processes. To achieve a more accurate understanding of the ionosphere’s response to the two flares, a more comprehensive modelling of the D-region chemistry is required.

Modelling the D-region chemistry is a three-step process: (i) defining the neutral background, (ii) computing the external forcing from the different sources and (iii) applying the Mitra-Rowe scheme to derive the ion and electron densities in quiescent time and their time evolution. These steps are detailed in the following.

2.1. Neutral Background

In the D-region, the main neutral species are O, $O_2$ and $N_2$. They are obtained from the NRLMISE-2.1 model [18], which also provides the neutral temperature $T_n$, common to all species. O, $O_2$ and $N_2$ are species being ionised. Other such species include N, which is taken from NRLMISE-00 [19] and an excited state of $O_2$ ($O_2{(}^1{\Delta }_g)$), which remains almost constant in the D-region at ${10}^{10}$ cm⁻³ [2,20].

The dominant source of ions in quiet time is the ionisation of $NO$ by Lyman-$\alpha$ radiation. $NO$ is a highly variable minor species; as such its density remains one of the largest sources of uncertainty in the model. The $\left[NO\right]$ reported from measurements (e.g., [21,22,23,24,25]) or models (e.g., [26,27,28,29,30]) vary from $5\times {10}^5$ cm⁻³ [28] to ${10}^9$ ${\mathrm{cm}}^{-3}$ [22] at 70 km. Another recent example of interest in $\left[NO\right]$ is the incorporation of $\left[NO\right]$ into NRLMSISE-2.1 [18] above 73 km. Here, $NO$ is parameterised from [26]:

$$\left[NO\right]=2\times {10}^{-2}exp(-300/T_n)\left[O_2\right]+5\times {10}^{-7}\left[O\right].$$

(2)

Since the ionisation rate of $NO$ is directly proportional to its concentration, the uncertainty of the $NO$ concentration leads to very different atmospheric compositions, even during quiescent periods. This is solved by a two-step ionosphere initialisation (Section 2.4).

Some neutral species are not directly ionised by external sources, but participate in the reaction rates between ions. A higher $\left[H_{2}O\right]$ thus increases the conversion rate from $O_4^{+}$ to cluster ions. [$H_{2}O$] is computed from [31], assuming a mixing ratio of 3 ppm. $\left[H\right]$ increases the detachment of electrons from $O_2^{-}$, and is provided by the NRLMSISE-2.1 model. Similarly, $O_3$ facilitates the transition from $O_2^{-}$ to other negative ions. It is initialised from [32], p. 178]:

$$\left[O_3\right]\approx {\displaystyle \frac{6\times {10}^{-34}}{1.02\times {10}^{-5}}}{\left(T_n/300\right)}^{-2.3}\left[M\right]\left[O_2\right]\left[O\right],$$

(3)

where $\left[M\right]=\left[O_2\right]+\left[N_2\right]$. More details are provided in Appendix B of the appendix.

Figure 1 displays the altitude profiles of different species that compose the neutral ionosphere. $N_2$ and $O_2$ are the dominant species.

2.2. Forcing Sources and Ionisation Rates

The variety of external ionisation sources is partially responsible for the complexity of the D-region. Since they have diverse energy ranges, they act differently on the neutral species (e.g., [1,4,5,6,7])

• The Lyman-$\alpha$ radiation acts on nitric oxide ($NO$),
• The EUV and UV radiation ionise $O_2{(}^1{\Delta }_g)$,
• The X-ray fluxes impact $O_2$, O and $N_2$ during solar flares,
• The cosmic rays ionise all species below 65 km.

The main source of ions outside solar flare periods comes from Lyman-$\alpha$ radiation. Solar X-ray fluxes only dominate during flares. It should be noted that other sources of ionisation are not considered in this study: they are the particle precipitation [e.g [33] which mainly occurs at high latitudes, as well as the impact of lightning strikes (e.g., [34]) and Lyman-$\alpha$ scattering by the geocorona [e.g. [3], p168], which play a more important role during nighttime only.

The calculation of the ionisation from the various sources depends significantly on the source considered. The next subsections describe these different ionisation terms.

2.2.1. Ionisation by Lyman-$\alpha$ Radiation and EUV Fluxes

The main ionisation source in quiet time is the Lyman-$\alpha$ acting on $NO$. The electron production $Q_{NO}$ is written:

$$Q_{NO}=F(\lambda =121.6\mathrm{nm},z){\sigma }_{NO}^i\left[NO\right]$$

(4)

where ${\sigma }_{NO}^i=2\times {10}^{-18}{\mathrm{cm}}^2$ represents the ionisation cross-section of $NO$ [35]. The main unknown in this case is $\left[NO\right]$, as mentioned in Section 2.1. $F(\lambda ,z)$ is the solar flux at wavelength $\lambda$ and altitude z. Its computation is detailed in Section 2.2.3.

Another major source of ions during quiet times is the ionisation of $O_2{(}^1{\Delta }_g)$ by wavelengths between 102.8-111.8 nm. This has been parameterised as [36]:

$$\begin{array}{cc}\hfill Q_{O_2{(}^1{\Delta }_g)}& =\left[O_2{(}^1{\Delta }_g)\right](5.49\times {10}^{-10}exp(-2.406\times {10}^{-20}\left[O_2\right]HCh(z,\chi ))\hfill \\ \hfill & +2.614\times {10}^{-9}exp(-8.508\times {10}^{-20}\left[O_2\right]HCh(z,\chi )))\hfill \end{array}$$

(5)

The Chapman function $Ch(z,\chi )$ is used to compute this integral. It depends on the altitude and the solar zenith angle $\chi$ (Appendix A.1 in the appendix). H is the atmospheric scale height. This expression considers the absorption of those wavelengths by $O_2$ and, implicitly through the numerical coefficients, the absorption by $C{O}_2$.

The EUV flux depends on solar activity. This can be modelled by multiplying Equation 5 by the following [37]:

$$C=0.764+3.693\times {10}^{-3}\left(F_{10.7}-80\right),$$

(6)

where we replaced the $A10.7$ index, which is an 81-day average on the $F_{10.7}$ index, with the $F_{10.7}$ index itself, to avoid relying on future values of this index.

2.2.2. Ionisation by X-Rays

Electron production by X-rays only dominates during solar flares, by ionisating $N_2$, O and N to form $O_2^{+}$, $O^{+}$ and $N_2^{+}$. In the D-region, only the HXR ionise species below 85 km (see Section 2.2.3). The production of ions from neutral species n is written:

$$Q_n=F(\lambda ,z){\sigma }_n^i\left[n\right]k_{\lambda }.$$

(7)

In this equation, $k_{\lambda }$ represents the number of electron-ion pairs created by one photon of wavelength $\lambda$ [38,39]. Indeed, the electrons created by direct ionisation from X-rays are energetic enough to participate in reactions, including excitation, ionisation and dissociation of neutral species. The $k_{\lambda }$ values have been tabulated to account for the different photon energies and species involved, and range from 67 to 555 [40]. Both $k_{\lambda }$ and ${\sigma }_n^i$ are presented in the appendix (Appendix A.3 and Appendix A.2).

2.2.3. Absorption of Solar Radiation Above the D-Region

Before reaching the D-region, solar fluxes cross the upper ionospheric regions. The neutrals at high altitudes absorb a fraction of this radiation, described by the Beer-Lambert law [e.g [41,42]:

$$F(\lambda ,z)=F(\lambda ,\infty )exp(-\tau ),$$

(8)

with $F(\lambda ,\infty )$ the solar flux at wavelength $\lambda$ before its absorption, $F(\lambda ,z)$ the same radiation at altitude z (z is measured from the ground), and $exp(-\tau )$ the transmission factor. $\tau$ is given by:

$$\tau \left(\lambda \right)=\sum _n{\int }_z^{\infty }\left[n\right]{\sigma }_n^a\left(\lambda \right)ds,$$

(9)

with $\left[n\right]$ the density of neutral species n and ${\sigma }_n^a$ its absorption cross-section. The integration runs along the radiation path, which depends on the solar zenith angle $\chi$. Equation 9 then becomes:

$$\tau \left(\lambda \right)=\sum _n\left[n\right]{\sigma }_n^a\left(\lambda \right)HCh(z,\chi ),$$

(10)

with H the atmosphere’s scale height and $Ch(z,\chi )$ the Chapman function. Further details on the Chapman function, the absorption cross sections of the different species at various wavelengths and the absorption of Lyman-$\alpha$ are given in Appendix A.1 in the appendix.

The transmission factor given by Equation 10 is represented in Figure 2. From this figure, several statements can be drawn. First, the SXR flux does not penetrate below 85 km. This implies that studies of the D-region response to solar flares should avoid classifying flares according to their soft X-ray flux. Then, the GOES HXR band comprises wavelengths of 0.05 nm to 0.4 nm. Using the integrated HXR flux over this entire band is not accurate enough, as the transmission factor for the entire band prevents the HXR from ionising the lower D-region (light-blue curve, Figure 2). This has already been noted in [40].

2.2.4. Ionisation by Cosmic Rays

The cosmic rays ionise $N_2$ and $O_2$ to produce $N_2^{+}$ and $O_2^{+}$ [e.g [4], and are the dominant ionisation term below 65 km. This ionisation is written [9]:

$$Q_{CR}=Q_p{\displaystyle \frac{N}{N_p}}exp(1-N/N_p)$$

(11)

with N the total number of neutral species, $N_p=N\left(15\mathrm{km}\right)$ and $Q_p=10{\mathrm{cm}}^{-3}{\mathrm{s}}^{-1}$.

All different ionisation terms are shown in Figure 3. Ionisation by Lyman-$\alpha$ radiation only concerns $NO$, but it is dominant in quiet times. X-rays ionise N, $N_2$ and $O_2$, and are mostly relevant during solar flares, above 70 km. The EUV ionisation of $O_2{(}^1{\Delta }_g)$ is the second most important source in quiet time, and operates above 75 km. The lowest altitudes of the D-region are maintained by the cosmic ray ionisation of the major neutral species $O_2$ and $N_2$.

2.3. Mitra-Rowe Scheme

To reproduce the time evolution of each species while maintaining a low numerical cost, the Mitra-Rowe chemistry scheme [1,12,13] is adopted. This scheme (Figure 4) includes seven species constituting the D-region: the electrons $e^{-}$, four positive ion species and two negative ion types. Those positive ions are the primary $N{O}^{+}$ and $O_2^{+}$, $O_4^{+}$ (which is important to treat separately due to its back reaction to $O_2^{+}$ [1]), and the positive cluster ions of the form $H^{+}{\left(H_{2}O\right)}_m$, grouped into $Y^{+}$. All negative ions are lumped together, which is the main assumption of the model, except for $O_2^{-}$, the main negative ion being produced [3,4].

Eighteen reactions are considered (Figure 4):

• Recombination between positive and negative species (rates ${\alpha }_1$ to ${\alpha }_4$ for the recombination with electrons, ${\alpha }_5$ for the recombination between $O_4^{+}$ and $X^{-}$ and ${\alpha }_i$ for other ion-ion recombination).
• Attachment of electrons on neutrals (reaction from $e^{-}$ to $O_2^{-}$)
• Detachment of electrons from negative ions (rates ${\gamma }_1$ and ${\gamma }_2$)
• Charge exchange reactions between positive species (reactions between $N{O}^{+}$, $O_2^{+}$ and $O_4^{+}$).
• Conversion of ions into heavy clusters (rates A, B, and reaction from $O_2^{-}$ to $X^{-}$)

Those reactions may also be written as a set of seven coupled differential equations. As an example, the electron density time evolution is described by:

$${\displaystyle \frac{d{N}_e}{dt}}={\gamma }_1\left[O_2^{-}\right]+{\gamma }_2\left[X^{-}\right]-f(O_2,N_2)N_e-\left({\alpha }_1\left[N{O}^{+}\right]+{\alpha }_2\left[O_2^{+}\right]+{\alpha }_3\left[O_4^{+}\right]+{\alpha }_4\left[Y^{+}\right]\right)N_e,$$

(12)

with $f(O_2,N_2)$ and the other coefficients given in Table A3 of the appendix. $f(O_2,N_2)$ represents the attachment rate of electrons on $O_2$ to form $O_2^{-}$.

The ionisation of neutral constituents by solar X-rays produces $O^{+}$, $N^{+}$ and $N_2^{+}$. No $O_2^{+}$ were created, since the ionisation cross-section for SXR and HXR wavelengths is zero [43]. However, the Mitra-Rowe scheme only takes $O_2^{+}$ and $N{O}^{+}$ as primary ions. $N_2^{+}$ ions are instantaneously transformed into $O_2^{+}$ [13,44], so their ionisation rate is directly added to the one of $O_2^{+}$. Both $O^{+}$ and $N^{+}$ are converted into $N{O}^{+}$ and $O_2^{+}$. To consider this conversion, the reaction rates in reactions 1 and 2, Tables 6 for $O^{+}$ and 14 to 17 for $N^{+}$ of [45] are computed. The ratio between these rates approximates the ion fractions converted into $N{O}^{+}$ or $O_2^{+}$.

The Mitra-Rowe scheme was notably validated by [6] during quiet times. They compared the outputs from this scheme to the SIC model and EISCAT-derived profiles, after fitting the solar X-ray flux to match E-region densities and adjusting $\left[NO\right]$, showing good agreement between the different approaches above 80 km. [12] also validated this scheme against observational data from rocket flights, incoherent scatter radar and ionosondes. They show that, with correct ionising fluxes and $NO$ profiles, agreements between the Mitra-Rowe scheme and observational profiles are reached. Note that this scheme solves the coupled differential equations at every altitude separately. It thus does not take into account the transport of species or diffusion.

2.4. Initialisation of the Ionosphere

All ion and electron concentrations are initially set to zero. The ionosphere is then initialised a first time from the actions of the various sources on the neutral background, using the Mitra-Rowe scheme (Section 2.3). The set of differential equations are numerically solved using the explicit Runge-Kutta method of order 4(5) implemented in the SciPy Python package. The electron density is compared to the Faraday International Reference Ionosphere (FIRI-2018, [46]; Figure 6, left panel). This model describes the electron density in the D-region between 60 km and 150 km, for solar zenith angles between 0 and 130° and latitudes below 60°. It has been compiled from various rocket flight measurements. As such, this is the reference data for the D-region. After this first initialisation, the modelled electron density significantly diverges from the FIRI model, especially in the middle layers, due to the low $\left[NO\right]$ profile (${10}^6$ ${\mathrm{cm}}^{-3}$ at 70 km in the example case). Indeed, in quiet periods, the ionisation of $NO$ dominates, especially at 75 km. To correct the $\left[NO\right]$ values and thus the quiet time ionisation of the D-region, $\left[NO\right]$ is scaled by a height-independent factor, defined at 75 km:

$$k={\displaystyle \frac{N_{e_{FIRI}}\left(75\mathrm{km}\right)}{N_{e_{ini}}\left(75\mathrm{km}\right)}},$$

(13)

with $N_{e_{FIRI}}$ and $N_{e_{ini}}$ from the FIRI model and the first initialisation, respectively. The variation of k with time for this first initialisation is in Figure 5. After the $NO$ renormalisation, the densities are then initialised a second time, from the same neutral background and forcing sources. As in the previous step, they are initialised to zero and then computed from the actions of the forcing sources on the neutral atmosphere. The density agrees more closely to the FIRI model at all altitudes (Figure 6, right), with errors between 0.4% (67 km) and 109% (60 km). The simulation retrieves the four areas described by [47]:

Figure 5. Variation of k (Equation 13 with time during the first initialisation. This scaling factor is constant after 3000 s, thus defining the time required to initialise the D-region.)

Figure 5. Variation of k (Equation 13 with time during the first initialisation. This scaling factor is constant after 3000 s, thus defining the time required to initialise the D-region.)

Figure 6. Comparison of the electron density to FIRI with initial $\left[NO\right]$ (left) or corrected $\left[NO\right]$ with scaling factor $k=7.5$ (right).

Figure 6. Comparison of the electron density to FIRI with initial $\left[NO\right]$ (left) or corrected $\left[NO\right]$ with scaling factor $k=7.5$ (right).

• above 85 km, dominated by $N{O}^{+}$ and $O_2^{+}$;
• between 82 km and 85 km, defined by strong gradients;
• between 70 km and 82 km, characterised by a nearly constant electron density and the presence of water clusters (included through $Y^{+}$);
• below 70 km, where negative ions are important.

The results presented in Figure 7 showed good agreement for most species with published results [48]. Once the initial ionosphere is set, time-varying external sources may be applied to model the density time evolution during perturbed times (Figure 7, for example, shows the densities at a peak of an M1 flare). This new chemistry model is named the Lower Ionospheric Region – Absorption and Chemistry Modelling (LIR-ACheM). On a personal laptop, modelling the response to a flare takes about a quarter of an hour, with about 90% due to the ionosphere initialisation.

3. D-Region Response to a Flare

In the following, the D-region response to a short M1 solar flare occurring on 2023/11/05 and peaking at 11:43 UT is analysed. The modelling is run at a latitude of 42° and a longitude of 9° to study the D-region behaviour at middle latitudes. One important aspect of the D-region response to a flare is the time between the X-ray peak and the electron density peak, named slugishness and denoted $\Delta t$ [49,50,51,52,53,54]. Indeed, the ionisation of neutral species is not instantaneous, and occurs with a small-time delay. Appleton [49] showed that

$$\Delta t={\displaystyle \frac{1}{2{\alpha }_{eff}{N}_e}},$$

(14)

with ${\alpha }_{eff}$ the effective recombination rate. $\Delta t$ varies with latitude and season [53], and more importantly, with altitudes (Figure 8). Here, $\Delta t_{SXR}$ (resp. $\Delta t_{HXR}$) denotes the delay between the SXR (resp. HXR) peak and the electron density modelled at each altitude. As presented in [54], $\Delta t_{HXR}$ is always higher than $\Delta t_{SXR}$, as the HXR peaks before the SXR. Figure 9 shows the computed ${\alpha }_{eff}$ variation with altitude, and displays the same variation with altitude observed in Figure 5.b. of [52], though the values are higher. In Equation 14, $N_e$ is evaluated at the flux peak. Figure 9 highlights the importance of defining this peak flux: the reference should be the HXR below 85 km and the SXR above this altitude.

Another important metric for the D-region’s response to a solar flare is its recovery time ${\tau }_r$, defined as the time for the electron density to decrease to half its peak value. It is height-dependent, as seen in Figure 8. Three areas defining both $\Delta t$ and ${\tau }_r$ are defined:

• Above 85 km, with a $\Delta t$ of less than a minute and ${\tau }_r$ of a few minutes, decreasing with height,
• Between 70 km and 85 km, where $\Delta t$ and ${\tau }_r$ increase with altitude,
• Below 70 km, with longer $\Delta t$ and ${\tau }_r$, decreasing with altitudes.

First, at low altitudes, the main source of electrons is not the photoionisation of neutral species, but the detachment of electrons from $X^{-}$ (Figure 10). This detachment operates on longer timescales than the direct ionisation by solar radiation. The creation of negative ions by attachment and subsequent release of electrons by detachment is thus perceived as a longer recovery, because the electron production by detachment is an additional source of ionisation and reaches its peak after the XR flux. This also causes a rise in $\Delta t$.

Above 80 km, the SXR flux is not completely absorbed. This is important, as the ionisation cross-sections for SXR wavelengths are much higher than those for the HXR fluxes (Appendix A in the appendix), and the SXR flux is approximately one order of magnitude above the HXR flux. A part of the ionisation at those altitudes is thus due to the SXR flux. A consequence of this is a rise in ${\tau }_r$ above 80 km, as the SXR peak occurs after the HXR (e.g. [54] and Figure 11). At 85 km, the ionisation rises with the (early) HXR flux, and falls with the (late) SXR flux. Thus, the perceived ${\tau }_r$ is longer, and $\Delta t$ increases.

Above 85 km, ${\tau }_r$ drops with altitude. At those altitudes, the main source of electrons is from the SXR flux (Figure 11). The peak of ionisation occurs later as the ionisation profile starts to follow the SXR flux. ${\tau }_r$ thus decreases, as the electron density responds to SXR only. Another factor to explain the altitude variation of ${\tau }_r$ is the high $\left[NO\right]$ at high altitude. Indeed, at 85 km and above, one of the main electron sinks is the recombination with $N{O}^{+}$. A higher $\left[NO\right]$ thus entails lower $\Delta t$ and ${\tau }_r$ as recombination occurs faster.

4. Discussion and Conclusions

A new chemistry model is presented here, LIR-ACheM, based on the Mitra-Rowe scheme. The LIR-ACheM model’s only required inputs are the geographic latitude and longitude, the time, and the forcing (given by the $F_{10.7}$ index for the EUV fluxes, the SXR flux in the band 0.1 - 0.8 nm and the HXR flux in the 0.05 - 0.4 nm band, as is given in the GOES data files). Those parameters are used to define the background neutral profiles and the ionisation rates. The outputs contain the variations of the different electron and ion densities for the duration of the modelling, as well as the scaling parameter k for $\left[NO\right]$.

Through this model, two important characteristics of the D-region response to solar flares were analysed: the peak delay $\Delta t$ and the recovery time ${\tau }_r$. Both were found to follow similar altitude profiles. At low altitudes, the main source of electrons is due to a photodetachment. This implies longer $\Delta t$ and apparent ${\tau }_r$, as this reaction operates on longer timescales. At higher altitudes, however, the ionisation by SXR becomes relevant. At the top of the D-region, it results in a longer apparent ${\tau }_r$ at 85 km, which decreases with altitude.

Further improvements in the model include adding more forcing sources (especially the scattering of Lyman-$\alpha$ radiation by the Earth’s geocorona, which is dominant in nighttime, and the forcing from electron precipitations for high latitudes). This model is mainly limited in altitude by the production rate of $O_2^{+}$ from EUV, parametrised by Paulsen et al. [36]. Another description of this rate could therefore extend the validity of this model to higher altitudes.

The user is free to provide their own neutral background. This model version (v1.3.1) indeed initialises the background neutral atmosphere from NRLMSISE-00 by default. [55] showed that using a different neutral background causes deviations in the computed electron density of the order of 34%. In the present study, the initialisation of the background was done, for example, through the NRLMSISE-2.1 model and the various references presented in Section 2.1.

Appendix A.1. Chapman Function

The Chapman function is approximated by [56]. If $\chi \le 90$°, then:

$$\begin{array}{cc}\hfill Ch(X,\chi )& =\sqrt{0.5\pi X}f\left(Y\right)\hfill \\ \hfill X& ={\displaystyle \frac{R_E+z}{H}}\hfill \\ \hfill Y& =|cos\chi |\sqrt{0.5X}\hfill \\ \hfill f\left(Y\right)& =exp\left(Y^2\right)erfcY\hfill \\ \hfill erfc& =1-{\displaystyle \frac{2}{\sqrt{\pi }}}{\int }_0^{Y}exp(-t^2)dt\hfill \end{array}$$

with $R_E$ the Earth’s radius. For $\chi >90$° (the Sun is below the horizon), $Ch$ becomes [56]:

$$Ch(X,\chi )=\sqrt{2\pi X}\left(\sqrt{sin\chi }expX(1-sin\chi )-0.5f\left(Y\right)\right)$$

Appendix A.2. Absorption and Ionisation Cross Sections

The absorption of the different wavelengths is due to the major species O, $O_2$ and $N_2$ (e.g., [43,45]). The absorption cross-sections are represented in Table A1.

Table A1. Absorption cross-sections (${\sigma }_n^a$, ${10}^{-18}$ cm²) of the main neutral species $N_2$, $O_2$ and O [43].

	SXR (0.1-0.8 nm)	HXR (0.05-0.4 nm)
$N_2$	0.0201	0.0025
$O_2$	0.0340	0.0045
O	0.0170	0.0023

Table A1. Absorption cross-sections (${\sigma }_n^a$, ${10}^{-18}$ cm²) of the main neutral species $N_2$, $O_2$ and O [43].

	SXR (0.1-0.8 nm)	HXR (0.05-0.4 nm)
$N_2$	0.0201	0.0025
$O_2$	0.0340	0.0045
O	0.0170	0.0023

The Lyman-$\alpha$ radiation, however, is mainly absorbed by $O_2$ [57]. The mean absorption cross-section is given by:

$$\begin{array}{cc}\hfill {\sigma }_{O_2}^a& =7.65\times {10}^{-21}\left(1+0.35exp(-5\times {10}^{-21}N\left(O_2\right))\right)\hfill \\ \hfill & \times \left(1.1+0.1tanh(3\times {10}^{-21}N\left(O_2\right)max((0.8-cos\chi ),0)-2.4)\right)\hfill \\ \hfill & \times \left(1.16-0.0021T_n+6\times {10}^{-6}{T}_n^2\right),\hfill \end{array}$$

(A1)

where the neutral temperature $T_n$ is in K and $N\left(O_2\right)=\left[O_2\right]HCh(X,\chi ){\mathrm{cm}}^{-2}$.

If the HXR range is discretised in 0.05 nm bins, then the cross-sections need to be adapted. They are presented in Table A2.

Table A2. Absorption cross-sections (${10}^{-18}$ ${\mathrm{cm}}^2$) of the main neutral species $N_2$, $O_2$ and O for different wavelengths in the HXR range [40].

Wavelength (nm)	${\mathit{N}}_2$	${\mathit{O}}_2$	O
0.05-0.10	$3.38$	$6.21$	$3.1$
0.10-0.15	$1.49$	$2.73$	$1.4$
0.15-0.20	$4.55$	$8.32$	$4.2$
0.20-0.25	$1.03$	$1.87$	$9.3$
0.25-0.30	$2.04$	$3.66$	$1.8$
0.30-0.40	$3.91$	$6.94$	$3.5$

Table A2. Absorption cross-sections (${10}^{-18}$ ${\mathrm{cm}}^2$) of the main neutral species $N_2$, $O_2$ and O for different wavelengths in the HXR range [40].

Wavelength (nm)	${\mathit{N}}_2$	${\mathit{O}}_2$	O
0.05-0.10	$3.38$	$6.21$	$3.1$
0.10-0.15	$1.49$	$2.73$	$1.4$
0.15-0.20	$4.55$	$8.32$	$4.2$
0.20-0.25	$1.03$	$1.87$	$9.3$
0.25-0.30	$2.04$	$3.66$	$1.8$
0.30-0.40	$3.91$	$6.94$	$3.5$

In Equation 7, ${\sigma }_n^i$ represents the ionisation cross-section of species n. This ionisation comes from the direct interaction of solar radiation with the neutral species. This interaction can result in ionisation, dissociation or dissociative ionisation, and the ratio $\beta$ between each of those possibilities is called the branching ratio. Specifically, ${\beta }_{ionisation}$ is the branching ratio for ionisation. Those values are tabulated in [40,43], so that:

$${\sigma }_n^i={\sigma }_n^a{\beta }_{ionisation}.$$

(A2)

Appendix A.3. Computation of k λ

In [43], the value of the ratio of photoelectron impact ionisation to ionisation by external sources $pe/pi$ is given for the three major neutral species, for different wavelengths. Here, $k_{\lambda }=pe/pi+1$. The $pe/pi$ values are later reevaluated by [40], for narrower HXR bins, and range between 66 and 554. Most of the D-region ionisation during flares thus comes from electron impact ionisation.