Daedalus Ionospheric Profile Continuation (DIPCont): Monte Carlo studies assessing the quality of in situ measurement extrapolation

In situ satellite exploration of the lower thermosphere–ionosphere system (LTI) as anticipated in the recent Daedalus mission proposal to ESA will be essential to advance the understanding of the interface between the Earth's atmosphere and its space environment. To address physical processes also below perigee, in situ measurements are to be extrapolated using models of the LTI. Motivated by the need for assessing how cost-critical mission elements such as perigee and apogee distances as well as the number of spacecraft affect the accuracy of scientific inference in the LTI, the Daedalus Ionospheric Profile Continuation (DIPCont) project is concerned with the attainable quality of in situ measurement extrapolation for different mission parameters and configurations. This report introduces the methodological framework of the DIPCont approach. Once an LTI model is chosen, ensembles of model parameters are created by means of Monte Carlo simulations using synthetic measurements based on model predictions and relative uncertainties as specified in the Daedalus Report for Assessment. The parameter ensembles give rise to ensembles of model altitude profiles for LTI variables of interest. Extrapolation quality is quantified by statistics derived from the altitude profile ensembles. The vertical extent of meaningful profile continuation is captured by the concept of extrapolation horizons defined as the boundaries of regions where the deviations remain below a prescribed error threshold. To demonstrate the methodology, the initial version of the DIPCont package presented in this paper contains a simplified LTI model with a small number of parameters. As a major source of variability, the pronounced change in temperature across the LTI is captured by self-consistent non-isothermal neutral-density and electron density profiles, constructed from scale height profiles that increase linearly with altitude. The resulting extrapolation horizons are presented for dual-satellite measurements at different inter-spacecraft distances but also for the single-satellite case to compare the two basic mission scenarios under consideration. DIPCont models and procedures are implemented in a collection of Python modules and Jupyter notebooks supplementing this report.

1 Introduction

The lower thermosphere–ionosphere system (LTI) at altitudes between about 100 and 200 km is characterized by transitions of several atmospheric attributes. It is the lower part of the heterosphere where atmospheric constituents are no longer mixed by turbulence and start to follow separate barometric laws (e.g., Picone et al., 2002; Izakov, 2007). As part of the thermosphere, the temperature profile shows a significant increase with altitude throughout the whole LTI (e.g., Chamberlain and Hunten, 1987). As part of the ionosphere, it includes the E-layer peak in electron density and the bottom side of the F layer (e.g., Hargreaves, 1992). With strongly altitude-dependent neutral–ion and neutral–electron collision frequencies, the LTI supports an anisotropic conductivity tensor that gives rise to a complex interplay of electric fields and currents. The conductivity tensor components affecting the directions perpendicular to the ambient magnetic field, namely Pedersen and Hall conductivities, show pronounced maxima in the LTI (e.g., Baumjohann and Treumann, 1996). A key variable quantifying its energetics is the Joule heating rate. Particular rich dynamics can be observed in the auroral region at high latitudes where energy and momentum from the magnetosphere are fed into the ionosphere through currents flowing parallel to the ambient magnetic field lines (e.g., Vogt et al., 1999). A comprehensive review of LTI features, measurement techniques, and models is provided by Palmroth et al. (2021).

Since the early 20th century, the LTI has been studied extensively using ground-based remote sensing facilities such as ionosondes and radars, but in all aspects requiring in situ observations, it remains underexplored territory. Rocket flights (e.g., Sangalli et al., 2009; Pfaff et al., 2022a) can offer only local and temporally confined information. Major technical challenges have so far prevented a satellite mission to the deep, dense part of the LTI, despite scientific interest, community proposals, and feasibility studies by major space agencies (e.g., Grebowsky and Gervin, 2001; Pfaff et al., 2022b). An early conception of the TIMED mission (e.g., Yee et al., 1999) considered dipper options for in situ investigations of the LTI. A recent initiative along this line is the Daedalus mission proposal (Sarris et al., 2020), submitted to ESA in response to the Explorer 10 call under the Earth Observation Program and selected together with two other proposals for a Phase-0 science and technical study. Daedalus aims to perform in situ measurements in the LTI from an elliptical orbit, with a nominal perigee of 150 km and an apogee on the order of 2000 km. Very low altitudes down to 120 km will be sampled by use of propulsion, through a series of short excursions in the form of perigee descent maneuvers. These are planned to be performed at high latitudes (> 65^∘ magnetic latitude), where Pedersen conductivity and Joule heating maximize. The highly elliptical orbit of Daedalus leads to a natural precession of the orbit's semi-major axis, in both magnetic latitude and magnetic local time; this means that Daedalus will perform measurements along its elliptical orbit down to the nominal perigee of 150 km throughout all magnetic latitudes. The geophysical observables sampled by Daedalus will enable a series of derived products to be obtained, as described in Table 1 of the Daedalus Report for Assessment (ESA, 2020), which, among many other geophysical variables, include the calculation of Pedersen conductivity and Hall conductivity.

The range of accessible perigees will be particular critical for any future LTI satellite mission, with severe impact on the propellant budget and other mission performance parameters (Sarris et al., 2020). With nominal perigees not much less than 150 km, peak conductivities and currents controlling E-region electrodynamics typically lie below the orbits. Physically meaningful downward continuation of in situ satellite measurements is desired, ideally using state-of-the-art models of the LTI (e.g., Sarris et al., 2023b). Another critical element of LTI mission conception is the number of spacecraft (Sarris et al., 2023a). The Daedalus Ionospheric Profile Continuation (DIPCont) project is concerned with vertical profiles of LTI variables and their reconstruction from dual-spacecraft and single-spacecraft observations. More specifically, the focus of the project is on the quality of profile continuation towards the lower LTI with its maxima in conductivities and current intensities, as given by the accuracy, the resolution, and the coverage of the reconstructions obtained from in situ measurements. Inspired by early work to extrapolate vertical profiles carried out under the Daedalus Phase-0 Science Study, DIPCont introduces a systematic probabilistic approach to the problem.

The DIPCont procedure to assess the quality of in situ measurement downward continuation is detailed in Sect. 3. In brief, after choosing an LTI model, representative ensembles of altitude profiles are generated by means of Monte Carlo simulations. The altitude profile ensembles give rise to statistical measures of relative deviation, which in turn allow for estimation of extrapolation horizons, effectively capturing the altitude range where deviations remain within given error thresholds. The basic ideas are illustrated in Fig. 1, displaying electron density and Pedersen conductivity extrapolation horizons for a range of relative error thresholds along the orbits of a dual-satellite mission. Horizontal distance corresponds to the latitudinal (north–south) direction, with the origin of the horizontal axis centered at the highest latitude along the satellite orbits. In the LTI model runs leading to Fig. 1, latitudinal inhomogeneity parameters are set to reproduce the two electron density maxima observed by a polar-orbiting satellite when crossing the auroral oval. See Sect. 4 and Appendix D for details. It is important to note that the filled contour representations of electron density and Pedersen conductivity model distributions mainly serve to provide contextual information, while the essential results of the DIPCont modeling procedure are the extrapolation horizons represented as plain contour lines, in response to the satellite orbit configuration (white lines). The extrapolation horizons of the model run shown in Fig. 1 suggest that for a dual-satellite mission as anticipated in the Daedalus Report for Assessment (ESA, 2020), downward continuation yields relative errors in the range of about 10–50 % at altitudes where electron density and Pedersen conductivity maximize. Implications are discussed in more detail further below in Sect. 4 and contrasted with the single-satellite case.

Figure 1 Extrapolation horizons and orbit configuration displayed on top of a two-dimensional section of the modeled LTI. (a) Electron density N_e. (b) Pedersen conductivity σ_P. Horizontal distance corresponds to the latitudinal (north–south) direction, with the origin of the horizontal axis centered at the highest latitude along the satellite orbits. In the LTI model runs leading to this figure, latitudinal inhomogeneity parameters are set to reproduce the two electron density maxima observed by a polar-orbiting satellite when crossing the auroral oval. Synthetic measurements are produced along the two satellite orbits (dashed white lines). The parameters of vertical profiles are estimated using measurements within a window (solid white rectangle) 2Δx in width around the nodes of a horizontal grid (dashed gray lines). Extrapolation horizons (solid and dotted colored lines) for a set of relative error levels are displayed as contours of a relative deviation measure, here the root-mean-square deviation of the ensemble of extrapolated profiles from the synthetic model prediction.

The LTI model used to introduce and demonstrate the DIPCont methodology in this paper is presented in Sect. 2. The parametric model captures the whole LTI temperature range and thus addresses a main source of variability. To limit the number of model parameters and thus also instabilities during model inversion in this initial DIPCont study, LTI variables showing less pronounced changes and ionization source mechanisms are treated in a simplified manner. Furthermore, since the quality of downward continuation is the focus of our study, the LTI model is restricted to E-region physics, with the influence of the F region left for future work.

Further first results are presented in Sect. 4, including a brief comparison between the single-spacecraft and the dual-spacecraft scenario. In Sect. 5, our findings are discussed in the context of important technical parameters and constraints relevant for a low-altitude mission. The body of the paper is concluded in Sect. 6 with prospects for upcoming work. Model derivations and technical details are presented in the Appendices, with particular emphasis on the incorporation of a non-isothermal temperature profile varying linearly with altitude.

2 Parametric models of LTI variables

Probabilistic measures of extrapolation quality produced by the DIPCont procedure detailed in Sect. 3 are based on synthetic in situ observations predicted by a model of the LTI. As emphasized in space physics textbooks and reviews of the LTI (e.g., Pfaff, 2012; Richmond, 1995), the full complexity of LTI variability and dynamics calls for a full multi-species description, taking into account source and loss processes varying in importance and efficiency as functions of magnetic latitude and local time and further factors. In the future, DIPCont functionality is planned to be included in the Daedalus MASE (Mission Assessment through Simulation Exercise) toolset (Sarris et al., 2023b), designed with the purpose to assess and demonstrate the closure of the mission objectives of the proposed Daedalus mission.

The more complex the LTI model of choice, however, the larger the number of parameters that are to be estimated with a downward continuation of in situ satellite measurements, which in turn tend to negatively affect the stability of model inversion. With these implications in mind, the initial version of the DIPCont package contains a simplified LTI description based on a limited set of parameters. Extrapolation quality of a single but important process, namely the formation of Pedersen conductivity σ_P, is supposed to be studied in a self-consistent manner. To this end, only a single particle species is considered, and classical photoionization physics is applied to parametrize ionospheric layer formation. Furthermore, as explained in reviews of ionospheric physics (e.g., Rishbeth, 1997), contributions from electron–neutral collisions to the Pedersen conductivity σ_P peak in the D region but are unimportant at higher altitudes; see also Fig. 4 in Sarris et al. (2023b). We thus arrive at the expression

$$\begin{array}{}\text{(1)}& {\mathit{\sigma }}_{\mathrm{P}}={\displaystyle \frac{N_{\mathrm{e}}{e}^{\mathrm{2}}}{m_{\mathrm{i}}}}{\displaystyle \frac{{\mathit{\nu }}_{\mathrm{in}}}{{\mathit{\nu }}_{\mathrm{in}}^{\mathrm{2}}+{\mathrm{\Omega }}_{\mathrm{i}}^{\mathrm{2}}}}\end{array}$$

(e: elementary charge; m_i: ion particle mass; Ω_i: ion gyrofrequency), suggesting that the altitude variabilities in electron number density N_e and ion–neutral collision frequency ν_in need to be modeled carefully. Less critical is the dependence of ion gyrofrequency Ω_i on magnetic field strength as it does not vary much over the LTI altitude range and is captured with sufficient accuracy by well-established empirical models. Different parametrizations exist for the ion–neutral collision frequency ν_in (e.g., Palmroth et al., 2021; Huba, 2019; Evans et al., 1977). In general the expressions are directly proportional to the number density N_n of neutral particles. As presented in Appendices A and B, also the self-consistent construction of electron density N_e rests prominently on the N_n profile, which in turn is conveniently modeled in terms of the density scale height $H_{\mathrm{n}}^N$. This aspect is chosen as a starting point below in Sect. 2.1 to further explain and motivate the LTI modeling approach.

The LTI model can be summarized in the form $\mathit{m}=\mathit{m}\left(z\right|\mathit{p})$, with a vector m of LTI observables and derived functions and a vector p=p(x) of model parameters, separating the primary (strong) dependence on altitude z from the secondary (weak) dependence on horizontal location x in a numerically efficient manner. Note the model functions are local representations of altitude profiles in the sense that they refer to a flexible reference altitude, z₀, that can be adapted to the locations where measurements are taken. In the DIPCont development phase it was observed that parameters of model functions in local representations typically showed weaker correlations and could be estimated more reliably, in particular as compared to the regional representations, relying on parameters at some fixed altitude, like the peak electron density height.

As in the Daedalus Report for Assessment (ESA, 2020), the vertical boundaries of the LTI region are assumed to be at z_B=100 km (base or bottom-side altitude) and z_T=200 km (topside altitude).

2.1 Scale height parameters

As demonstrated in Appendix A, profiles of (neutral gas) pressure P_n and neutral (number) density N_n are conveniently constructed using

$$\begin{array}{}\text{(2)}& {\displaystyle }{P}_{\mathrm{n}}\left(z\right)=P_{\mathrm{n}\mathrm{0}}\exp \left\{-\underset{z_{\mathrm{0}}}{\overset{z}{\int }}{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{̃}}{z}}{H_{\mathrm{n}}^P\left(\stackrel{\mathrm{̃}}{z}\right)}}\right\},\text{(3)}& {\displaystyle }{N}_{\mathrm{n}}\left(z\right)=N_{\mathrm{n}\mathrm{0}}\exp \left\{-\underset{z_{\mathrm{0}}}{\overset{z}{\int }}{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{̃}}{z}}{H_{\mathrm{n}}^N\left(\stackrel{\mathrm{̃}}{z}\right)}}\right\},\end{array}$$

where $H_{\mathrm{n}}^P$ and $H_{\mathrm{n}}^N$ denote the pressure scale height and the density scale height, respectively. Furthermore, it is shown that

$$\begin{array}{}\text{(4)}& H_{\mathrm{n}}^N=H_{\mathrm{n}}^P{\left(\mathrm{1}+{\displaystyle \frac{\mathrm{d}{H}_{\mathrm{n}}^P}{\mathrm{d}z}}\right)}^{-\mathrm{1}}\end{array}$$

if the pressure scale height

$$\begin{array}{}\text{(5)}& H_{\mathrm{n}}^P={\displaystyle \frac{k{T}_{\mathrm{n}}}{m_{\mathrm{n}}g}}={\displaystyle \frac{R_{\mathrm{gas}}{T}_{\mathrm{n}}}{M_{\mathrm{n}}g}}\end{array}$$

(k: Boltzmann constant; T_n: neutral temperature; g: gravitational acceleration; R_gas: universal gas constant; m_n: average particle mass; M_n: average molar mass) changes only with temperature T_n. Equations (4) and (5) further imply that, if temperature T_n varies linearly with altitude z, then $H_{\mathrm{n}}^P$ and $H_{\mathrm{n}}^N$ do too. In Appendix A it is shown that the constant inverse gradients

$$\begin{array}{}\text{(6)}& \mathit{\gamma }={\left({\displaystyle \frac{\mathrm{d}{H}_{\mathrm{n}}^P}{\mathrm{d}z}}\right)}^{-\mathrm{1}}\end{array}$$

and

$$\begin{array}{}\text{(7)}& \mathit{\eta }={\left({\displaystyle \frac{\mathrm{d}{H}_{\mathrm{n}}^N}{\mathrm{d}z}}\right)}^{-\mathrm{1}}\end{array}$$

are related through

$$\begin{array}{}\text{(8)}& \mathit{\eta }=\mathit{\gamma }+\mathrm{1}.\end{array}$$

Variations in gravity g across the LTI are in the range of a few percent and can be neglected in this context. Profiles of T_n, M_n, and $H_{\mathrm{n}}^P$ as predicted by the empirical atmospheric model NRLMSIS 2.0 (Emmert et al., 2021) for different seasons and latitudes are displayed in Fig. S1a–d as part of the Supplement to this paper, indicating that relative variations in average molar mass are indeed significantly smaller than those of neutral temperature. We thus disregard altitude changes in average molar mass M_n as imposed by changes in atmospheric composition and further assume that temperature T_n, pressure scale height $H_{\mathrm{n}}^P$, and density $H_{\mathrm{n}}^N$ vary linearly with altitude in a self-consistent manner as described by Eqs. (4) and (5). According to Eq. (A11), the local density scale height $H_{\mathrm{n}\mathrm{0}}^N$ can be obtained using the inverse scale height gradients and the local pressure scale height $H_{\mathrm{n}\mathrm{0}}^N$ from the expression

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{H}_{\mathrm{n}\mathrm{0}}^N& ={\displaystyle \frac{H_{\mathrm{n}\mathrm{0}}^P}{\mathrm{1}+{\mathit{\gamma }}^{-\mathrm{1}}}}={\displaystyle \frac{\mathit{\gamma }}{\mathit{\gamma }+\mathrm{1}}}{H}_{\mathrm{n}\mathrm{0}}^P={\displaystyle \frac{\mathit{\eta }-\mathrm{1}}{\mathit{\eta }}}{H}_{\mathrm{n}\mathrm{0}}^P\\ \\ & ={\displaystyle \frac{\mathit{\eta }-\mathrm{1}}{\mathit{\eta }}}{\displaystyle \frac{R_{\mathrm{gas}}{T}_{\mathrm{n}\mathrm{0}}}{M_{\mathrm{n}}g}},\end{array}\end{array}$$

where the subscript 0 indicates that the respective variable is taken at the reference (measurement) altitude z₀.

Since in the construction of neutral-density and electron density profiles (Appendix A and B) the scale height parameters are of central importance, their relative change is reflected in the following reference values: at z_B, T_nB∼200 K, and $H_{\mathrm{nB}}^P\sim \mathrm{6}\mathrm{km}$. Figure S1a–d suggest that the pressure scale height varies by a factor of ∼ 5 across the LTI; thus $H_{\mathrm{nT}}^P\sim \mathrm{5}\cdot H_{\mathrm{nB}}^P\sim \mathrm{30}\mathrm{km}$, and T_nT∼1000 K at z_T= 200 km. We obtain $\mathrm{d}{H}_{\mathrm{n}}^P/\mathrm{d}z=\mathrm{0.24}={\mathit{\gamma }}^{-\mathrm{1}}$, γ∼4, $\mathit{\eta }=\mathit{\gamma }+\mathrm{1}\sim \mathrm{5}$, and $H_{\mathrm{nB}}^N=\frac{\mathit{\eta }-\mathrm{1}}{\mathit{\eta }}{H}_{\mathrm{nB}}^P\sim \mathrm{5}\mathrm{km}$ for the density scale height at the base of the LTI.

2.2 Neutral temperature

Neutral temperature T_n is assumed to vary linearly with altitude z:

$$\begin{array}{}\text{(10)}& T_{\mathrm{n}}\left(z\right|z_{\mathrm{0}},T_{\mathrm{n}\mathrm{0}},L_{\mathrm{n}\mathrm{0}})=T_{\mathrm{n}\mathrm{0}}\cdot \left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{L_{\mathrm{n}\mathrm{0}}}}\right).\end{array}$$

The parameters T_n0 and L_n0 are the neutral temperature and the gradient length scale, respectively, at a reference altitude z₀. The constant temperature gradient is given by

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}{T}_{\mathrm{n}}}{\mathrm{d}z}}={\displaystyle \frac{T_{\mathrm{n}\mathrm{0}}}{L_{\mathrm{n}\mathrm{0}}}}.\end{array}$$

Using Eq. (A13), the neutral-temperature profile can be expressed by means of the parameters η and $H_{\mathrm{n}\mathrm{0}}^N$ as follows:

$$\begin{array}{}\text{(12)}& T_{\mathrm{n}}\left(z\right|z_{\mathrm{0}},T_{\mathrm{n}\mathrm{0}},H_{\mathrm{n}\mathrm{0}}^N,\mathit{\eta })=T_{\mathrm{n}\mathrm{0}}\cdot \left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}}\right).\end{array}$$

2.3 Neutral density

The altitude dependence of neutral density N_n for linear scale height profiles is derived in Appendix A, resulting in the following local representation:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& N_{\mathrm{n}}\left(z\right|z_{\mathrm{0}},N_{\mathrm{n}\mathrm{0}},H_{\mathrm{n}\mathrm{0}}^N,\mathit{\eta })\\ & =N_{\mathrm{n}\mathrm{0}}\cdot \exp \left\{-\mathit{\eta }\ln \left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}}\right)\right\}\\ & =N_{\mathrm{n}\mathrm{0}}\cdot {\left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}}\right)}^{-\mathit{\eta }};\end{array}\end{array}$$

see also Eq. (A15). The parameter N_n0=N_n(z₀) is the local neutral density, i.e., its value at the reference altitude z₀.

2.4 Electron density

The altitude dependence of electron density N_e for linear scale height profiles is derived in Appendix B, resulting in the following local representation:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& N_{\mathrm{e}}\left(z\right|z_{\mathrm{0}},N_{\mathrm{e}\mathrm{0}},L_{\mathrm{r}\mathrm{0}}\cos \mathit{\chi },H_{\mathrm{n}\mathrm{0}}^N,\mathit{\eta })\\ & =N_{\mathrm{e}\mathrm{0}}\exp \left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}}\left[-{\mathit{\theta }}_{\mathrm{0}}+{\displaystyle \frac{H_{\mathrm{n}\mathrm{0}}^N}{L_{\mathrm{r}\mathrm{0}}\cos \mathit{\chi }}}\left(\mathrm{1}-e^{-{\mathit{\theta }}_{\mathrm{0}}}\right)\right]\right\},\end{array}\end{array}$$

with ${\mathit{\theta }}_{\mathrm{0}}={\mathit{\theta }}_{\mathrm{0}}\left(z\right)=(\mathit{\eta }-\mathrm{1})\ln \left(\mathrm{1}+\frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}\right)$; see Eq. (B16). The parameter N_e0=N_e(z₀) gives electron density at the chosen reference altitude z₀. Note that L_r0 and χ, the angle of incident radiation with the atmospheric layer normal direction, cannot be estimated separately but only combined as L_r0cos χ. The parameters $H_{\mathrm{n}\mathrm{0}}^N$ and η can be inherited from estimations using neutral-temperature and/or neutral-density data, effectively reducing the number of electron density parameters and thus stabilizing the estimation procedure.

The non-isothermal electron density model can also be expressed in terms of the ionization peak parameters, namely the altitude z_* and the electron density value $N_{{\mathrm{e}}^{*}}=N_{\mathrm{e}}\left(z_{*}\right)$:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& N_{\mathrm{e}}\left(z\right|z_{*},N_{{\mathrm{e}}^{*}},H_{{\mathrm{n}}^{*}}^N,\mathit{\eta })\\ & =N_{{\mathrm{e}}^{*}}\exp \left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}}\left[-{\mathit{\theta }}_{*}+\mathrm{1}-e^{-{\mathit{\theta }}_{*}}\right]\right\},\end{array}\end{array}$$

with ${\mathit{\theta }}_{*}={\mathit{\theta }}_{*}\left(z\right)=(\mathit{\eta }-\mathrm{1})\ln \left(\mathrm{1}+\frac{z-z_{*}}{\mathit{\eta }{H}_{{\mathrm{n}}^{*}}^N}\right)$ and $H_{{\mathrm{n}}^{*}}^N$ denoting the density scale height at $z=z_{*}$. See Appendix B2 for details. Electron density profiles for identical peak parameters but different values of η are displayed in Fig. 2.

Figure 2 Altitude dependence of the non-isothermal electron density model for different values of the inverse neutral-density scale height gradient η. Common electron density peak parameters are $z_{*}=\mathrm{110}\mathrm{km}$, $N_{{\mathrm{e}}^{*}}={\mathrm{10}}^{\mathrm{11}}{\mathrm{m}}^{-\mathrm{3}}$, and $H_{{\mathrm{n}}^{*}}^N=\mathrm{7}\mathrm{km}$. The case η→∞ (in the legend, η=inf) corresponds to the isothermal limit.

The electron density model is designed to describe the ionospheric E layer, assuming that contributions from the F layer are modeled separately and subtracted from the measurements. To account for residuals that may remain after subtraction, the DIPCont package contains a parameter N_eF.

2.5 Ion temperature

Temperature profiles obtained by the International Reference Ionosphere (IRI) 2.0 model (Bilitza et al., 2022) indicate that ion and neutral temperatures are very similar throughout the LTI; see Fig. S2a–d. In analogy with the neutral-temperature case, ion temperature T_i is assumed to vary linearly with altitude z:

$$\begin{array}{}\text{(16)}& T_{\mathrm{i}}\left(z\right|z_{\mathrm{0}},T_{\mathrm{i}\mathrm{0}},L_{\mathrm{i}\mathrm{0}})=T_{\mathrm{i}\mathrm{0}}\cdot \left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{L_{\mathrm{i}\mathrm{0}}}}\right).\end{array}$$

The parameters T_i0 and L_i0 are the ion temperature and the gradient length scale, respectively, at the chosen reference altitude z₀.

2.6 Ion–neutral collision frequency

In quantitative terms, collision processes in the partially ionized LTI medium remain inadequately described and are major sources of uncertainties in empirical models (e.g., Palmroth et al., 2021; Heelis and Maute, 2020; Sarris, 2019). At this stage, the DIPCont project is less concerned with optimizing the quantitative description of the LTI, but rather with the quality of parameter estimation extrapolation. While the choice of the best LTI model is certainly important for recovering the real values of targeted observables, further work will be needed by parametric studies, comparison with previous work, and data analysis when a low-perigee mission such as Daedalus (Sarris et al., 2020) provides in situ measurements in the LTI. For our goal here, the chosen variant among the models for ion–neutral collision frequency ν_in should not matter too much as long as the underlying variability associated with erroneous measurements is captured. To this end, we follow the description of Huba (2019) and write

$$\begin{array}{}\text{(17)}& {\mathit{\nu }}_{\mathrm{in}}={\mathit{\sigma }}_{\mathrm{in}}{N}_{\mathrm{n}}\sqrt{{\displaystyle \frac{k{T}_{\mathrm{i}}}{m_{\mathrm{i}}}}},\end{array}$$

with the collision cross-section ${\mathit{\sigma }}_{\mathrm{in}}\sim \mathrm{5}\cdot {\mathrm{10}}^{-\mathrm{15}}{\mathrm{cm}}^{\mathrm{2}}$. An even simpler expression could neglect the variation with ion temperature T_i so that ν_in becomes directly proportional to the neutral density N_n.

2.7 Pedersen conductivity

Using the approximations explained at the beginning of Sect. 2, Pedersen conductivity is given by

$$\begin{array}{}\text{(18)}& {\mathit{\sigma }}_{\mathrm{P}}={\displaystyle \frac{N_{\mathrm{e}}{e}^{\mathrm{2}}}{m_{\mathrm{i}}}}{\displaystyle \frac{{\mathit{\nu }}_{\mathrm{in}}}{{\mathit{\nu }}_{\mathrm{in}}^{\mathrm{2}}+{\mathrm{\Omega }}_{\mathrm{i}}^{\mathrm{2}}}}\end{array}$$

for a quasi-neutral two-component plasma when the contribution from electron–neutral collisions is neglected; see also Eq. (1), reproduced here for convenience. Compared to other variables and parameters of the LTI models presented here, the dependence of ion gyrofrequency ${\mathrm{\Omega }}_{\mathrm{i}}=q_{\mathrm{i}}B/m_{\mathrm{i}}$ (q_i: ion charge; m_i: ion mass) on magnetic field strength B can be determined from measurements or models of the magnetic field with very good accuracy; hence the associated variability should not much affect our results. Furthermore, in the logic of the LTI model constructed for the initial version of the DIPCont package, changes in atmospheric composition and thus average ion mass are disregarded. Inspection of Fig. S2a–d indicates that in the lower part of the LTI (altitudes below about 150 km), being the focus of downward continuation quality in the current study, variations in average ion mass with altitude are relatively small. Hence, altitude variations in ion gyrofrequency are neglected. In the same way as for other LTI model variables, namely through the dependence of the parameters in the vector p=p(x) (see Sect. 2.8 below) on the coordinate x, horizontal variations in magnetic field strength B and thus ion gyrofrequency Ω_i can be modeled and are planned to be considered in future work.

2.8 LTI model in compact form

Parameters of model functions in local representation are listed in Table 1.

The description of the DIPCont modeling procedure in Sect. 3 benefits from summarizing the LTI model in compact form as $\mathit{m}=\mathit{m}\left(z\right|\mathit{p})$, with parameters $T_{\mathrm{n}\mathrm{0}},H_{\mathrm{n}\mathrm{0}}^N,\mathit{\eta }$, etc. entering the vector p. The parametric functions T_n(z), N_n(z), N_e(z), T_i(z), ${\mathit{\nu }}_{\mathrm{in}}\left(z\right)={\mathit{\nu }}_{\mathrm{in}}\left(N_{\mathrm{n}}\right(z),T_{\mathrm{i}}(z\left)\right)$, and ${\mathit{\sigma }}_{\mathrm{P}}\left(z\right)={\mathit{\sigma }}_{\mathrm{P}}\left(N_{\mathrm{e}}\right(z),{\mathit{\nu }}_{\mathrm{in}}(z\left)\right)$ constitute the components of the vectorial function m.

3 DIPCont modeling procedure

The DIPCont modeling procedure is as follows.

• Synthetic noise-free measurements ${\mathit{m}}_j=\mathit{m}\left(z_j\right|\mathit{p})$ are created along anticipated Daedalus satellite orbit sections around perigee at altitudes z_j=z(t_j) and horizontal distances x_j=x(t_j). The chosen model parameters are defined by vectors p=p(x_#) on a grid of horizontal distances x_#. The integration and approximation methods employed for constructing the satellite orbits are described in Sect. 3.1 and in Appendix C.

• Using the multiplicative noise model presented in Sect. 3.2, synthetic measurements are contaminated by random errors in accordance with relative uncertainties specified in the Daedalus Report for Assessment (ESA, 2020), yielding ensembles $\mathit{\left\{}{\stackrel{\mathrm{̃}}{\mathit{m}}}_j^k\mathit{\right\}}$ of noisy synthetic data sets.

• For a point x_# on the horizontal grid, synthetic data with horizontal distances x_j in $[x_{\mathit{\#}}-\mathrm{\Delta }x,x_{\mathit{\#}}+\mathrm{\Delta }x]$ are considered to produce a least squares estimate ${\widehat{\mathit{p}}}^k\left(x_{\mathit{\#}}\right)$ of the parameter vector p(x_#). Repeating the estimation procedure for all members k of the ensemble $\mathit{\left\{}{\stackrel{\mathrm{̃}}{\mathit{m}}}_j^k\mathit{\right\}}$ of synthetic data sets yields ensembles of model parameters $\mathit{\left\{}{\widehat{\mathit{p}}}^k\right(x_{\mathit{\#}}\left)\mathit{\right\}}$ for all horizontal grid points x_#. Specifics of the estimation procedure are discussed in Sect. 3.3.

• With parameter vectors $\mathit{p}\in \mathit{\left\{}{\widehat{\mathit{p}}}^k\right(x_{\mathit{\#}}\left)\mathit{\right\}}$, the parametric model function $\mathit{m}=\mathit{m}\left(z\right|\mathit{p})$ can be evaluated to obtain ensembles $\mathit{\left\{}{\widehat{\mathit{m}}}^k\right(z,x_{\mathit{\#}}\left)\mathit{\right\}}=\mathit{\left\{}\mathit{m}\right(z\left|{\widehat{\mathit{p}}}^k\right(x_{\mathit{\#}}\left)\right)\mathit{\}}$, representing altitude profiles of LTI observables and derived variables such as ν_in and σ_P over the entire range of LTI altitudes and for all horizontal grid points x_#. The resulting altitude profiles form a representative ensemble in the sense that their statistics are compatible with the model functions and the set of given relative errors. Relative deviation measures of observables and derived variables as functions of altitude are constructed. Finally, the concept of extrapolation horizons, introduced in Sect. 3.4, captures the altitude range where errors are tolerable according to predefined thresholds.

Table 1 Parameters of model functions in local representation. The list is partially redundant, e.g., $L_{\mathrm{n}\mathrm{0}}=\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N$. The parameters L_r0 and χ cannot be estimated independently but only in combination (L_r0cos χ). Boundary data (neutral and ion temperatures at the base of the LTI) are used to constrain the parameters L_n0, η, $H_{\mathrm{n}\mathrm{0}}^N$, and L_i0; see Sect. 3.3.

3.1 Satellite orbits around perigee

The DIPCont model offers two options for computing altitudes and horizontal distances along the orbits of satellites around perigee, namely numerical integration by means of the Störmer–Verlet method (e.g., Hairer et al., 2003) and the polynomial approximation

$$\begin{array}{}\text{(19)}& {\displaystyle }z\left(t\right)=z_{\mathrm{per}}+{\displaystyle \frac{a_{\mathrm{per}}}{\mathrm{2}}}{t}^{\mathrm{2}},\text{(20)}& {\displaystyle }x\left(t\right)={\displaystyle \frac{R_{\mathrm{E}}{V}_{\mathrm{per}}t}{R_{\mathrm{per}}}}\left(\mathrm{1}-{\displaystyle \frac{a_{\mathrm{per}}}{\mathrm{3}{R}_{\mathrm{per}}}}{t}^{\mathrm{2}}\right),\end{array}$$

with the acceleration a_per at perigee given by

$$\begin{array}{}\text{(21)}& a_{\mathrm{per}}={\displaystyle \frac{G{M}_{\mathrm{E}}}{R_{\mathrm{per}}^{\mathrm{2}}}}{\displaystyle \frac{R_{\mathrm{apo}}-R_{\mathrm{per}}}{R_{\mathrm{apo}}+R_{\mathrm{per}}}}=g_{\mathrm{per}}\mathit{\epsilon };\end{array}$$

see Appendix C. Here R_E is the Earth's radius; M_E is the Earth's mass; G is the gravitational constant; z_per, R_per, and V_per are the altitude, geocentric distance, and satellite velocity at perigee; $g_{\mathrm{per}}=\frac{G{M}_{\mathrm{E}}}{R_{\mathrm{per}}^{\mathrm{2}}}$ is the Earth's gravitational acceleration at geocentric distance R_per; R_apo is the geocentric distance at apogee; and $\mathit{\epsilon }=\frac{R_{\mathrm{apo}}-R_{\mathrm{per}}}{R_{\mathrm{apo}}+R_{\mathrm{per}}}$ is the orbital eccentricity. For the parameter range considered in this study, the deviation of the polynomial approximation from the more precise orbit integration is on the order of a few hundred meters; see Fig. S4b of this report.

The observation time T_obs spent by Daedalus in the LTI during a perigee pass controls the amount of data that can be gathered for statistical investigations. Using the quadratic orbital approximation around perigee, T_obs is twice the time needed to move from z=z_per to the upper boundary at z=z_T; thus $z_{\mathrm{T}}-z_{\mathrm{per}}=\frac{a_{\mathrm{per}}}{\mathrm{2}}(T_{\mathrm{obs}}/\mathrm{2}{)}^{\mathrm{2}}$, $T_{\mathrm{obs}}^{\mathrm{2}}=\frac{\mathrm{8}(z_{\mathrm{T}}-z_{\mathrm{per}})}{a_{\mathrm{per}}}$, and

$$\begin{array}{}\text{(22)}& T_{\mathrm{obs}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{8}(z_{\mathrm{T}}-z_{\mathrm{per}})R_{\mathrm{per}}^{\mathrm{2}}}{G{M}_{\mathrm{E}}}}{\displaystyle \frac{R_{\mathrm{apo}}+R_{\mathrm{per}}}{R_{\mathrm{apo}}-R_{\mathrm{per}}}}={\displaystyle \frac{\mathrm{8}(z_{\mathrm{T}}-z_{\mathrm{per}})}{g_{\mathrm{per}}\mathit{\epsilon }}}.\end{array}$$

The variations in T_obs with apogee altitude in the range of $\mathrm{1500}\mathrm{km}\le z_{\mathrm{apo}}\le \mathrm{3000}\mathrm{km}$ for the three perigee altitudes z_per=115, 130, and 150 km are displayed in Fig. 3. Raising the perigee from 115 to 130 km yields a small reduction in observation time by about 10 %. Within the range of orbital parameters considered here, the overall amount of data gathered during a perigee pass turns out to depend only moderately on apogee altitude z_apo, with a relative difference of not more than about 20 % for changes in z_apo between 2000 and 3000 km.

Figure 3 Observation time T_obs in the LTI versus apogee altitude z_apo for three different values of perigee altitude z_per. The topside of the LTI is assumed to be at z_T=200 km.

When dual-satellite missions to the LTI are considered, the question arises as to how synchronous the measurements are with respect to ground horizontal distance x, assuming the two spacecraft share the same orbital plane, have identical semi-major axes and thus orbital periods, and pass through their perigees at the same time. Figure S4a illustrates how visit times of ground horizontal distances are expected to differ for two satellites with perigee altitudes of 130 and 150 km. Differences in satellite visit times turn out to be on the order of seconds.

3.2 Synthetic measurements and positivity constraints

Synthetic measurements $\mathit{\{}{\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathrm{1}},{\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathrm{2}},{\stackrel{\mathrm{̃}}{\mathit{\mu }}}_{\mathrm{3}},\mathrm{\dots }\mathit{\}}$ of an observable at altitudes $\mathit{\{}{z}_{\mathrm{1}},z_{\mathrm{2}},z_{\mathrm{3}},\mathrm{\dots }\mathit{\}}$ are constructed from a parametric model function $\mathit{\mu }=\mathit{\mu }\left(z\right|\mathit{p})$ producing predictions that are contaminated by random errors $\mathit{\{}{\mathit{\sigma }}_{\mathrm{1}},{\mathit{\sigma }}_{\mathrm{2}},{\mathit{\sigma }}_{\mathrm{3}},\mathrm{\dots }\mathit{\}}$ from a suitable probability distribution. The model parameter vector p is estimated through minimization of a cost function. Following the standard least squares approach, the cost function is chosen to be the error-scaled square deviation

$$\begin{array}{}\text{(23)}& {\mathit{\chi }}^{\mathrm{2}}\left(\mathit{p}\right)=\sum _j{\left({\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\mu }}}_j-\mathit{\mu }\left(z_j\right|\mathit{p})}{{\mathit{\sigma }}_j}}\right)}^{\mathrm{2}}.\end{array}$$

The observables of interest T_n, N_n, N_e, and T_i are all positive; hence a straightforward additive noise model would not be appropriate as it may produce negative synthetic data. Furthermore, instrumental uncertainties as provided in the Daedalus Report for Assessment (ESA, 2020) are typically specified as relative (multiplicative) errors. Both issues are addressed by considering as model predictions ${\mathit{\mu }}_j=\mathit{\mu }\left(z_j\right|\mathit{p})$ and data ${\stackrel{\mathrm{̃}}{\mathit{\mu }}}_j$ not the positive observables as such but their (natural) logarithms and relative uncertainties for the random errors $\mathit{\{}{\mathit{\sigma }}_{\mathrm{1}},{\mathit{\sigma }}_{\mathrm{2}},{\mathit{\sigma }}_{\mathrm{3}},\mathrm{\dots }\mathit{\}}$. In the case of a (neutral or electron) density N, one obtains

$$\begin{array}{}\text{(24)}& \ln {\stackrel{\mathrm{̃}}{N}}_j=\ln N\left(z_j\right)+{\mathit{\sigma }}_j{r}_j,\end{array}$$

where $r_j\sim \mathcal{N}(\mathrm{0},\mathrm{1})$ represents Gaussian noise (normally distributed random numbers with zero mean and unit variance), and N=N(z) refers to the (positive) density model. Then

$$\begin{array}{}\text{(25)}& {\stackrel{\mathrm{̃}}{N}}_j=e^{{\mathit{\sigma }}_j{r}_j}\cdot N\left(z_j\right)\end{array}$$

so that positivity is guaranteed. Furthermore,

$$\begin{array}{}\text{(26)}& e^{{\mathit{\sigma }}_j{r}_j}\approx \mathrm{1}+{\mathit{\sigma }}_j{r}_j,\end{array}$$

showing that the parameters σ_j correspond to relative error levels. Table 2 summarizes the values used in this report.

Table 2 Relative error levels used in this study, according to Table 2 of the Daedalus Report for Assessment (ESA, 2020).

In general, the parameters enter the logarithms of model functions nonlinearly, and an iterative estimation procedure is required.

3.3 Parameter estimation strategies

The model parameters listed in Table 1 are estimated from observations of neutral temperature T_n, neutral density N_n, electron density N_e, and ion temperature T_i as follows.

• For a given horizontal grid location x_#, data within the interval $[x_{\mathit{\#}}-\mathrm{\Delta }x,x_{\mathit{\#}}+\mathrm{\Delta }x]$ are considered. The effective window width is 2Δx; see the solid white rectangles in Figs. 1 and 4.

• From T_n data and constraining the neutral-temperature profile at the LTI lower boundary z_B as explained below, infer T_n0, $H_{\mathrm{n}\mathrm{0}}^N$, and η. See Eq. (10) and Sect. 2.1.

• Using $H_{\mathrm{n}\mathrm{0}}^N$ and η, estimate N_n0 from N_n data. See Eq. (13).

• Using $H_{\mathrm{n}\mathrm{0}}^N$ and η, estimate N_e0 and L_r0cos χ from N_e data. See Eq. (14).

• From T_i data and a suitable constraint at the lower LTI boundary z_B in analogy to the neutral-temperature case, infer T_i0 and L_i0. See Eq. (16).

Altitude profiles of these observables allow for construction of the height dependence of derived variables such as ion–neutral collision frequency ν_in and the Pedersen conductivity σ_P; see Eqs. (17) and (18), respectively.

3.3.1 Lower LTI boundary constraints

As explained in Appendix A, Eqs. (A10) and (A13), the linear density scale height profile can be parametrized using $H_{\mathrm{n}\mathrm{0}}^N$ and η in the form

$$\begin{array}{}\text{(27)}& H_{\mathrm{n}}^N\left(z\right|z_{\mathrm{0}},H_{\mathrm{n}\mathrm{0}}^N,\mathit{\eta })=H_{\mathrm{n}\mathrm{0}}^N\cdot \left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}}\right).\end{array}$$

It is important to note that the $H_{\mathrm{n}}^N$ profile takes center stage in the LTI models of the observables T_n, N_n, and N_e. While the local temperature amplitude T_n0 is essentially an average of local temperature data around an altitude z₀, and the same applies to the local pressure scale height $H_{\mathrm{n}\mathrm{0}}^P$ obtained from T_n0 by simple multiplication, the inverse density scale height gradient η and thus also the local density scale height parameter $H_{\mathrm{n}\mathrm{0}}^N=\frac{\mathit{\eta }-\mathrm{1}}{\mathit{\eta }}{H}_{\mathrm{n}\mathrm{0}}^P$ are very challenging to estimate from purely local data with little variance in altitude, as suggested already by the standard error of the slope in linear regression analysis. Fortunately, neutral temperature at the base of the LTI is known with reasonable tolerances from atmospheric models (e.g., Picone et al., 2002; Emmert et al., 2021). This remote data point constitutes a valuable constraint for estimating the density scale height profile. To incorporate model uncertainties and expected deviations from actual values, boundary data at the base of the LTI are contaminated by random errors according to the approach described in Sect. 3.2.

To be specific, the pressure scale height gradient, constant under the assumptions discussed in Sect. 2.1 and Appendix A, can be obtained from its values $H_{\mathrm{n}\mathrm{0}}^P$ and $H_{\mathrm{nB}}^P$ at z₀ and the LTI base altitude z_B, respectively, as follows:

$$\begin{array}{}\text{(28)}& {\displaystyle \frac{\mathrm{d}{H}_{\mathrm{n}}^P}{\mathrm{d}z}}={\displaystyle \frac{H_{\mathrm{n}\mathrm{0}}^P-H_{\mathrm{nB}}^P}{z_{\mathrm{0}}-z_{\mathrm{B}}}}.\end{array}$$

The inverse gradients γ and η of pressure scale height and density scale height, respectively, are related by Eq. (8) through $\mathit{\eta }=\mathit{\gamma }+\mathrm{1}$; thus the parameter η is given by

$$\begin{array}{}\text{(29)}& \mathit{\eta }={\displaystyle \frac{z_{\mathrm{0}}-z_{\mathrm{B}}}{H_{\mathrm{n}\mathrm{0}}^P-H_{\mathrm{nB}}^P}}+\mathrm{1}={\displaystyle \frac{M_{\mathrm{n}}g}{R_{\mathrm{gas}}}}{\displaystyle \frac{z_{\mathrm{0}}-z_{\mathrm{B}}}{T_{\mathrm{n}\mathrm{0}}-T_{\mathrm{nB}}}}+\mathrm{1},\end{array}$$

where T_nB denotes the neutral temperature at z_B. The local density scale height $H_{\mathrm{n}\mathrm{0}}^N$ can now be obtained from Eq. (9) as

$$\begin{array}{}\text{(30)}& H_{\mathrm{n}\mathrm{0}}^N={\displaystyle \frac{H_{\mathrm{n}\mathrm{0}}^P}{\mathrm{1}+{\mathit{\gamma }}^{-\mathrm{1}}}}={\displaystyle \frac{\mathit{\eta }-\mathrm{1}}{\mathit{\eta }}}{\displaystyle \frac{R_{\mathrm{gas}}{T}_{\mathrm{n}\mathrm{0}}}{M_{\mathrm{n}}g}}.\end{array}$$

Figure 4 Model distributions of electron density N_e (a) and Pedersen conductivity σ_P (b) in the LTI. Synthetic measurements are produced along the two satellite orbits (dashed white lines). The parameters of vertical profiles are estimated using measurements within a window (solid white rectangle) around two locations in the horizontal direction (dashed blue and green lines).

3.3.2 Linear estimation of electron density parameters

The logarithm of the electron density model considered here,

$$\begin{array}{}\text{(31)}& \begin{array}{rl}\ln N_{\mathrm{e}}\left(z\right)& =\ln N_{\mathrm{e}\mathrm{0}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}}\left[-{\mathit{\theta }}_{\mathrm{0}}+{\displaystyle \frac{H_{\mathrm{n}\mathrm{0}}^N}{L_{\mathrm{r}\mathrm{0}}\cos \mathit{\chi }}}\left(\mathrm{1}-e^{-{\mathit{\theta }}_{\mathrm{0}}}\right)\right],\end{array}\end{array}$$

can be combined with the logarithm of the neutral-density model,

$$\begin{array}{}\text{(32)}& \ln N_{\mathrm{n}}\left(z\right)=\ln N_{\mathrm{n}\mathrm{0}}-{\displaystyle \frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}}{\mathit{\theta }}_{\mathrm{0}},\end{array}$$

to find

$$\begin{array}{}\text{(33)}& \begin{array}{rl}& \ln N_{\mathrm{e}}\left(z\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\ln N_{\mathrm{n}}\left(z\right)\\ & =\ln N_{\mathrm{e}\mathrm{0}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\ln N_{\mathrm{n}\mathrm{0}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}}{\displaystyle \frac{H_{\mathrm{n}\mathrm{0}}^N}{L_{\mathrm{r}\mathrm{0}}\cos \mathit{\chi }}}\left(\mathrm{1}-e^{-{\mathit{\theta }}_{\mathrm{0}}}\right)\\ & =a+b\left(\mathrm{1}-e^{-{\mathit{\theta }}_{\mathrm{0}}}\right),\end{array}\end{array}$$

showing that $a=\ln N_{\mathrm{e}\mathrm{0}}-\frac{\mathrm{1}}{\mathrm{2}}\ln N_{\mathrm{n}\mathrm{0}}$ and $b=\frac{\mathrm{1}}{\mathrm{2}}\frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}\frac{H_{\mathrm{n}\mathrm{0}}^N}{L_{\mathrm{r}\mathrm{0}}\cos \mathit{\chi }}$ can be obtained from linear regression of $\ln N_{\mathrm{e}}-\frac{\mathrm{1}}{\mathrm{2}}\ln N_{\mathrm{n}}$ versus $\mathrm{1}-e^{-{\mathit{\theta }}_{\mathrm{0}}}$ with ${\mathit{\theta }}_{\mathrm{0}}={\mathit{\theta }}_{\mathrm{0}}\left(z\right)=(\mathit{\eta }-\mathrm{1})\ln \left(\mathrm{1}+\frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}\right)$. Since the parameters η and $H_{\mathrm{n}\mathrm{0}}^N$ are available as estimates from T_n modeling, and N_n0 is known from N_n modeling, N_e0 and L_r0cos χ can be computed from the linear coefficients a and b; hence this special case does not necessitate an iterative parameter estimation approach.

3.4 Error profiles and extrapolation horizons

With $\mathit{\{}{\stackrel{\mathrm{̃}}{\mathit{m}}}_j^k{\mathit{\}}}_{j\in \left[\mathit{\#}\right]}^{k=\mathrm{1}}$ being a single set (k=1) of synthetic measurements and j∈[#] indicating that horizontal distances are selected to be within ±Δx around a predefined grid point x_#, the estimation procedure yields a specific estimate ${\widehat{\mathit{p}}}^k$ of the parameter vector p(x_#). In a Monte Carlo setup, different instances of random errors are applied to the model predictions to produce data sets $\mathit{\{}{\stackrel{\mathrm{̃}}{\mathit{m}}}_j^{\mathrm{1}},{\stackrel{\mathrm{̃}}{\mathit{m}}}_j^{\mathrm{2}},{\stackrel{\mathrm{̃}}{\mathit{m}}}_j^{\mathrm{3}},\mathrm{\dots }{\mathit{\}}}_{j\in \left[\mathit{\#}\right]}$. The ensemble of data sets gives rise to an ensemble of parameter vectors $\mathit{\left\{}{\widehat{\mathit{p}}}^k\mathit{\right\}}=\mathit{\{}{\widehat{\mathit{p}}}^{\mathrm{1}},{\widehat{\mathit{p}}}^{\mathrm{2}},{\widehat{\mathit{p}}}^{\mathrm{3}},\mathrm{\dots }\mathit{\}}$, which in turn, when entered in $\mathit{m}=\mathit{m}\left(z\right|\mathit{p})$, yields an ensemble of profiles $\mathit{\left\{}{\widehat{\mathit{m}}}^k\right(z,x_{\mathit{\#}}\left)\mathit{\right\}}=\mathit{\left\{}{\widehat{\mathit{m}}}^{\mathrm{1}}\right(z,x_{\mathit{\#}}),{\widehat{\mathit{m}}}^{\mathrm{2}}(z,x_{\mathit{\#}}),{\widehat{\mathit{m}}}^{\mathrm{3}}(z,x_{\mathit{\#}}),\mathrm{\dots }\mathit{\}}$ for the entire range of altitudes z and at each point x_# of the horizontal coordinate grid.

The procedure is illustrated in Figs. 4 and 5. Figure 4 shows the model functions and the satellite orbits used for computing the predictions that enter the Monte Carlo simulation. The ensemble of altitude profiles generated from the Monte Carlo distributions of model parameters is visualized in Fig. 5 by means of selected quantiles evaluated at the vertical grid of LTI altitudes.

Figure 5 Visualization of the ensemble of altitude profiles generated from the Monte Carlo distributions of model parameters. Shown are selected quantiles evaluated at the vertical grid of LTI altitudes. (a, c) Electron density N_e. (b, d) Pedersen conductivity σ_P. (a, b) Center position (dashed blue line) in Fig. 4. (c, d) Right position (dashed green line) in Fig. 4.

The ensemble of altitude profiles forms the basis for quantifying extrapolation quality through measures of relative deviation from a model prediction. Suppressing altitude and horizontal grid dependencies and considering only a single model variable μ with ensemble members ${\widehat{\mathit{\mu }}}^{\mathrm{1}},{\widehat{\mathit{\mu }}}^{\mathrm{2}},{\widehat{\mathit{\mu }}}^{\mathrm{3}},\mathrm{\dots },{\widehat{\mathit{\mu }}}^K$, the root-mean-square deviation is given by

$$\begin{array}{}\text{(34)}& \mathit{\delta }\mathit{\mu }=\sqrt{〈(\widehat{\mathit{\mu }}-\mathit{\mu }{)}^{\mathrm{2}}〉}=\sqrt{{\displaystyle \frac{\mathrm{1}}{K}}\sum _{k=\mathrm{1}}^K({\widehat{\mathit{\mu }}}^k-\mathit{\mu }{)}^{\mathrm{2}}}.\end{array}$$

Figure 6 shows the altitude profiles of relative root-mean-square deviation $\mathit{\delta }\mathit{\mu }/\mathit{\mu }=\sqrt{〈(\widehat{\mathit{\mu }}-\mathit{\mu }{)}^{\mathrm{2}}〉}/\mathit{\mu }$ for the variables and horizontal locations as in Figs. 4 and 5.

Figure 6 Solid lines (blue and green) give the relative root-mean-square deviations of Monte Carlo altitude profiles from the respective input model profiles at two horizontal locations. Vertical dotted and solid lines represent a set of chosen error levels, ranging from 0.5 % and 1 % (yellow) to 32 % and 64 % (magenta). The corresponding horizontal lines show the extrapolation horizons indicating at which altitude the relative deviation equals the respective error level. (a, c,) Electron density N_e. (b, d) Pedersen conductivity σ_P. (a, b) Center position (dashed blue line) in Fig. 4. (c, d) Right position (dashed green line) in Fig. 4.

Figure S3 provides additional information on this DIPCont model run; visualizing model distributions; ensembles of altitude profiles; and extrapolation horizons also for neutral temperature T_n, neutral density N_n, ion temperature T_i, and ion–neutral collision frequency ν_in.

Alternative relative deviation measures considered in the DIPCont package are based on the empirical distribution of absolute deviations $|\widehat{\mathit{\mu }}-\mathit{\mu }|$, e.g., the average absolute deviation from the model prediction μ,

$$\begin{array}{}\text{(35)}& (\mathit{\delta }\mathit{\mu }{)}_{\mathrm{abs}}=〈|\widehat{\mathit{\mu }}-\mathit{\mu }|〉={\displaystyle \frac{\mathrm{1}}{K}}\sum _{k=\mathrm{1}}^K|{\widehat{\mathit{\mu }}}^k-\mathit{\mu }|,\end{array}$$

or selected quantiles of the distribution.

3.5 Implementation

The DIPCont model is implemented as a bundle of Python instructions and functions collected in three modules.

In the module DIPContBas.py, the basic setup of the DIPCont framework is defined, e.g., LTI region boundaries and boundary values, satellite orbit parameters, horizontal grid locations, and auxiliary plot parameters. Furthermore, it also provides configurational variables that are exchanged between DIPCont functions and modules, e.g., parameters shared by different parametric models.

The module DIPContMod.py provides parametric model functions of LTI variables and plot routines.

The module DIPContEst.py is concerned with Monte Carlo parameter estimation and profile continuation. Estimation of parameters that enter the model functions nonlinearly is accomplished by the function curve_fit() from the module scipy.optimize, whereas linear parameter estimation is performed using the function linregress from the module scipy.stats. Monte Carlo ensembles of parameters and altitude profiles are stored in pandas DataFrames.

The three DIPCont modules are provided in the Supplement to this report, together with Jupyter notebooks to explain and illustrate their usage.

4 First results

The major ingredients of the DIPCont processing chain, namely generation of synthetic in situ measurements along satellite orbits, Monte Carlo simulations of vertical profiles, and construction of extrapolation horizons, are summarized in Figs. 4–6, displaying electron density N_e and Pedersen conductivity σ_P as two variables of key importance for the structure and the dynamics of the LTI. As indicated by Eq. (1) and the respective profiles in Fig. 5, electron density makes the main contribution to the peaked height variation in Pedersen conductivity, with secondary contributions of neutral density and possibly ion temperature through the parametric form chosen for the ion–neutral collision frequency; see Sect. 2.6 and also Fig. S3. Furthermore, Pedersen conductivity controls the height variation in Joule heating, whose characterization is one of the main scientific targets of the proposed Daedalus mission (ESA, 2020). In the neutral wind reference frame, Joule heating is ${\mathit{j}}_{⟂}\cdot {\mathit{E}}_{⟂}={\mathit{\sigma }}_{\mathrm{P}}{\left|{\mathit{E}}_{⟂}\right|}^{\mathrm{2}}$, where the subscript ⟂ indicates a vectorial component perpendicular to the ambient magnetic field direction $\widehat{\mathit{B}}$. Height variations in E_⟂ are negligible according to the following rationale; see, e.g., Rishbeth (1997). Due to high parallel conductivity, the electric field component $E_{\Vert }=E_s$ parallel to $\widehat{\mathit{B}}$ vanishes, i.e., $\mathrm{0}=E_s=-\frac{\partial \mathrm{\Phi }}{\partial s}$, where s is the magnetic field line coordinate, and Φ denotes the electric potential. The electric field component E_q in a direction perpendicular to $\widehat{\mathit{B}}$ captured by a coordinate q then satisfies $\frac{\partial E_q}{\partial s}=-\frac{\partial }{\partial s}\frac{\partial \mathrm{\Phi }}{\partial q}=-\frac{\partial }{\partial q}\frac{\partial \mathrm{\Phi }}{\partial s}=\frac{\partial E_s}{\partial q}=\mathrm{0}$.

When instead of two selected horizontal locations as in Fig. 6 an equidistant grid of horizontal coordinates is defined for DIPCont simulations and the construction of extrapolation horizons, the results can be displayed together with the underlying model distributions and satellite orbits as in Fig. 1. In the following examples, such displays are used to visualize DIPCont results for different spacecraft configurations. Section 4.1 offers a first qualitative assessment of extrapolation quality in terms of varying inter-spacecraft distance. Section 4.2 contrasts the performance of the dual-spacecraft configuration considered so far with the results of the single-spacecraft case.

Note that the horizontal axis corresponds to the latitudinal (north–south) direction. In the simulations that led to Figs. 4–6, horizontal variations were disregarded for better comparability. In Fig. 1 and in the following, latitudinal inhomogeneity of electron density is meant to reproduce the two maxima observed by a polar-orbiting satellite when crossing the auroral oval. The highest latitude corresponds to the origin of the horizontal axis. Since the physics of energetic particle precipitation is not incorporated in this initial version of the DIPCont package, the horizontal variation in electron density expected for an auroral oval crossing is prescribed through ad hoc choices of horizontal electron density peak parameter profiles; see the option LTIModelType='NeAuroralZoneCrossing' in the DIPCont code as part of the Supplement to this report. The functional forms of horizontal electron density peak parameters are given in Appendix D.

4.1 Varying inter-spacecraft distance

Extrapolation of two-point measurements is expected to perform best if the spatial separation matches the relevant physical length scale. In the LTI this should be the (local) density scale height, in the range of 10–20 km for altitudes above 130 km, as in our example of a dual-spacecraft setup with perigee altitudes of 130 and 150 km; see Fig. 1. The inter-spacecraft distance remains close to 20 km throughout the whole orbit section and thus also to the density scale height as the relevant physical scale. Note that in all dual-satellite DIPCont model runs presented in this paper, apogee distances of the second satellite have been adjusted such that the sum of perigee and apogee distances are identical for both satellites and thus also the semi-major axes and the orbital periods.

Figure 7 displays extrapolation horizons for the same simulation setup except that the perigee altitude of the second satellite is reduced to 135 km, producing an inter-spacecraft distance at perigee of only 5 km. The separation is now smaller than the local density scale height with values of about 15 km at altitudes around 150 km. Compared to Fig. 1, the errors are increased and the extrapolation horizons reduced. The changes are not dramatic but enough to show that inter-spacecraft distance is a parameter to be considered when extrapolation quality is supposed to be optimized.

Figure 7 Same as Fig. 1 but for an inter-spacecraft separation of 5 km at perigee.

4.2 Single-satellite case

To check how much a second satellite improves extrapolation quality, the Monte Carlo simulations summarized in Fig. 1 are repeated for the single-spacecraft case, with all other parameters left unchanged. The resulting extrapolation horizons are shown in Fig. 8. Compared to ionospheric profile continuation from dual-spacecraft observations, the single-spacecraft case yields significantly worse results, with extrapolation horizons collapsing into the orbit near the perigee due to lacking variability in altitudes. Away from the perigee, the orbital motion of the satellite during the time corresponding to the horizontal window width 2Δx yields some height range that allows for profile reconstruction but with significant errors. The peaks in electron density and Pedersen conductivity are clearly outside the largest considered error level of 64 %, while Fig. 1 shows that in the dual-spacecraft case the peaks are between the 16 % and 32 % error levels.

Figure 8 Same as Fig. 1 but for the single-satellite case.

5 Discussion

Our first results suggest that altitude profiles of key LTI variables can be reconstructed with sufficient accuracy from in situ measurements if the effective altitude range covers relevant physical scales such as the local density scale height $H_{\mathrm{n}\mathrm{0}}^N$. This is the case for a dual-spacecraft configuration with an inter-spacecraft separation of 20 km at perigee; see Figs. 1 and 4–6. By two-point sampling, one can retrieve the vertical profiles of electron density and Pedersen conductivity essentially down to the bottom of the LTI region, a few scale heights under the lower satellite and including the peak altitudes. For Pedersen conductivity, errors are expected in the range of about 10–50 %, with the peak altitude and most of the conductivity within the 32 % extrapolation horizon in the chosen example, consistent with rocket observations (Sangalli et al., 2009).

Given the current knowledge of key LTI variables, error levels of about 10–50 % may well improve the situation. An important motivation behind the Daedalus proposal was the large error margin in Joule heating estimates, with a major contribution by errors in conductance (height-integrated conductivity). Thus, Sarris et al. (2020) pointed out that for a substorm event investigated by Palmroth et al. (2005), there were differences of up to 500 % between three proxies for the Joule heating rates integrated over the Northern Hemisphere. Even if this setup cannot be directly compared to our virtual environment, the order-of-magnitude difference between the two error margins looks encouraging for follow-up work on ionospheric profile continuation.

The DIPCont framework allows economical and technical questions regarding the impact of different LTI mission cost factors to be addressed. On the one hand, a dual-spacecraft mission seems to automatically imply higher costs because a second satellite needs to be built. On the other hand, a major cost driver of any deep LTI mission is the necessary amount of propellant that is required in order to maintain a spacecraft in orbit due to enhanced atmospheric drag at very low perigee altitudes. Since in a dual-spacecraft setup the role of the lower perigee satellite can be shared, each of the two probes would have to carry half of the total amount of propellant required to maintain the same total observation time required by a single-satellite mission. Moreover, the necessary amount of thermal shielding depends as well on perigee altitude, and each of the two probes would have to withstand the maximum thermal stress at perigee less often. Our findings show that the two-point setup allows for a more effective extrapolation to lower altitudes, which in turn means that a higher perigee may well be a meaningful option.

Data processing would also benefit from raising the perigee. As shown by simulations carried out for the technical assessment of Daedalus (ESA, 2020), a hydrodynamic shock develops in front of the spacecraft at altitudes under ∼120–130 km, complicating the retrieval of unperturbed data from the observed ones. Another LTI mission parameter considered in this paper is the apogee altitude controlling the proximity to the Van Allen belts and thus the necessary amount of radiation shielding but also affecting the available LTI observation time near perigee. The analysis presented in Sect. 3.1 shows that the amount of data gathered for statistical studies depends only moderately on apogee altitude.

The current version of the DIPCont framework concentrates on the E layer, assuming that contributions from the F layer can be disregarded or subtracted before processing, e.g., using the NeQuick approach to model topside ionospheric sounding data (Pignalberi et al., 2020). The DIPCont package contains a parameter N_eF to study the effect of F-layer residuals on extrapolation quality in future work.

The first results presented here are planned to be validated and extended in more extensive studies. Besides varying orbital parameters such as perigee altitude and inter-spacecraft distance, the impact of numerical parameters such as the horizontal selection window, 2Δx, needs further investigation. As already commented in Sect. 2.6, alternative functional forms for modeling ion–neutral collision frequencies or other variables may also be considered.

6 Conclusions and outlook

The DIPCont methodology introduced in this paper is designed to assess the quality of downward continuation of LTI variables using in situ satellite measurements and parametric models. While first results have been obtained with a simplified LTI description based on a single particle species, the Monte Carlo simulation machinery in DIPCont is not constrained to a particular model setup. By quantifying the quality of extrapolated in situ measurements, DIPCont can help to assess the science return of specific configurations and thus to optimize the parameters of upcoming LTI missions.

First DIPCont tests, performed on electron density and Pedersen conductivity, show promising results, to be consolidated by further parametric studies. Application of DIPCont to a modeled event, like the geomagnetic storm event of March 2015 addressed in the Daedalus Report for Assessment (ESA, 2020), is an upcoming target. This could be performed using the capabilities of the Daedalus MASE toolset (Sarris et al., 2023b). Future studies are planned to include Joule heating, which was a major driver of the Daedalus mission proposal. To investigate auroral processes and the electrodynamics of magnetosphere–ionosphere coupling, ionization through energetic particle precipitation needs to be incorporated. The Hall current nature of auroral electrojets calls for including electron–neutral collisional interaction as a major contributor to Hall conductivity formation.

Coordination between an LTI mission, like Daedalus, and a topside mission, like Swarm or DMSP (Defense Meteorological Satellite Program), would enhance the return of both missions. As an example, reconstruction of vertical profiles of ionospheric conductivity based on LTI observations could help to calibrate topside estimates of the conductance, while topside electron density could provide upper continuation and constrain the height-integrated total electron content (TEC) inferred from LTI data. Combination with ground-based observatory data such as ionosondes would offer further valuable constraints to DIPCont and thus enable more comprehensive modeling of the LTI.

Appendix A: Neutral-density profile for linear variations in scale height

Consider an atmospheric layer dominated by possibly several neutral constituents with an average or representative particle mass m_n, total pressure P_n, mass density ϱ_n, effective neutral number density $N_{\mathrm{n}}={\mathit{\varrho }}_{\mathrm{n}}/m_{\mathrm{n}}$, and temperature T_n. Under hydrostatic conditions, $\mathrm{d}{P}_{\mathrm{n}}=-{\mathit{\varrho }}_{\mathrm{n}}g\mathrm{d}z$, where z is altitude, and g is gravity (gravitational acceleration), assumed to vary so little within the layer that it can be safely considered constant. Using the ideal gas law P_n = N_nkT_n, where k denotes the Boltzmann constant, one obtains $\mathrm{d}{P}_{\mathrm{n}}=-P_{\mathrm{n}}\frac{m_{\mathrm{n}}g}{k{T}_{\mathrm{n}}}\mathrm{d}z=-P_{\mathrm{n}}\frac{\mathrm{d}z}{H_{\mathrm{n}}^P}$ with the pressure scale height

$$\begin{array}{}\text{(A1)}& H_{\mathrm{n}}^P={\displaystyle \frac{k{T}_{\mathrm{n}}}{m_{\mathrm{n}}g}}.\end{array}$$

Rearranging $-\frac{\mathrm{d}z}{H_{\mathrm{n}}^P}=\frac{\mathrm{d}{P}_{\mathrm{n}}}{P_{\mathrm{n}}}=\mathrm{d}\ln P_{\mathrm{n}}$ and integrating leads to

$$\begin{array}{}\text{(A2)}& P_{\mathrm{n}}\left(z\right)=P_{\mathrm{n}\mathrm{0}}\exp \left\{-\underset{z_{\mathrm{0}}}{\overset{z}{\int }}{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{̃}}{z}}{H_{\mathrm{n}}^P\left(\stackrel{\mathrm{̃}}{z}\right)}}\right\},\end{array}$$

where the altitude dependence of $H_{\mathrm{n}}^P$ directly reflects the change in temperature T_n with z.

Analogous differential and integral expressions for the neutral density, namely $\mathrm{d}\ln N_{\mathrm{n}}=-\frac{\mathrm{d}z}{H_{\mathrm{n}}^N}$ and

$$\begin{array}{}\text{(A3)}& N_{\mathrm{n}}\left(z\right)=N_{\mathrm{n}\mathrm{0}}\exp \left\{-\underset{z_{\mathrm{0}}}{\overset{z}{\int }}{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{̃}}{z}}{H_{\mathrm{n}}^N\left(\stackrel{\mathrm{̃}}{z}\right)}}\right\},\end{array}$$

are derived as follows. Combining the differential of the ideal gas law $\mathrm{d}{P}_{\mathrm{n}}=N_{\mathrm{n}}k\mathrm{d}{T}_{\mathrm{n}}+k{T}_{\mathrm{n}}\mathrm{d}{N}_{\mathrm{n}}$ with the hydrostatic condition yields $-N_{\mathrm{n}}{m}_{\mathrm{n}}g\mathrm{d}z=N_{\mathrm{n}}k\mathrm{d}{T}_{\mathrm{n}}+k{T}_{\mathrm{n}}\mathrm{d}{N}_{\mathrm{n}}$ and thus $-\frac{m_{\mathrm{n}}g}{k{T}_{\mathrm{n}}}\mathrm{d}z-\frac{\mathrm{1}}{T_{\mathrm{n}}}\mathrm{d}{T}_{\mathrm{n}}=\frac{\mathrm{1}}{N_{\mathrm{n}}}\mathrm{d}{N}_{\mathrm{n}}=\mathrm{d}\ln N_{\mathrm{n}}$. Since $\frac{\mathrm{d}{T}_{\mathrm{n}}}{T_{\mathrm{n}}}=\mathrm{d}\ln T_{\mathrm{n}}=\mathrm{d}\ln H_{\mathrm{n}}^P=\frac{\mathrm{d}{H}_{\mathrm{n}}^P}{H_{\mathrm{n}}^P}$, one obtains

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\mathrm{d}\ln N_{\mathrm{n}}}{\mathrm{d}z}}=-{\displaystyle \frac{\mathrm{1}}{H_{\mathrm{n}}^P}}\left(\mathrm{1}+{\displaystyle \frac{\mathrm{d}{H}_{\mathrm{n}}^P}{\mathrm{d}z}}\right).\end{array}$$

Therefore, the density scale height $H_{\mathrm{n}}^N$ in the expression $\mathrm{d}\ln N_{\mathrm{n}}=-\frac{\mathrm{d}z}{H_{\mathrm{n}}^N}$ is given by

$$\begin{array}{}\text{(A5)}& H_{\mathrm{n}}^N=H_{\mathrm{n}}^P{\left(\mathrm{1}+{\displaystyle \frac{\mathrm{d}{H}_{\mathrm{n}}^P}{\mathrm{d}z}}\right)}^{-\mathrm{1}}.\end{array}$$

To be more specific, we suppose the neutral temperature T_n varies linearly with altitude z,

$$\begin{array}{}\text{(A6)}& T_{\mathrm{n}}\left(z\right)=T_{\mathrm{n}\mathrm{0}}\cdot \left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{L_{\mathrm{n}\mathrm{0}}}}\right),\end{array}$$

where T_n0 is the temperature at a reference altitude z₀, and $L_{\mathrm{n}\mathrm{0}}=\frac{T_{\mathrm{n}\mathrm{0}}}{\mathrm{d}{T}_{\mathrm{n}}/\mathrm{d}z}$ denotes the local gradient length. Then

$$\begin{array}{}\text{(A7)}& H_{\mathrm{n}}^P\left(z\right)=H_{\mathrm{n}\mathrm{0}}^P\cdot \left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{L_{\mathrm{n}\mathrm{0}}}}\right),\end{array}$$

with

$$\begin{array}{}\text{(A8)}& H_{\mathrm{n}\mathrm{0}}^P={\displaystyle \frac{k{T}_{\mathrm{n}\mathrm{0}}}{mg}},\end{array}$$

so that the pressure scale height gradient $\frac{\mathrm{d}{H}_{\mathrm{n}}^P}{\mathrm{d}z}=\frac{H_{\mathrm{n}\mathrm{0}}^P}{L_{\mathrm{n}\mathrm{0}}}$ is constant and thus also the gradient of density scale height:

$$\begin{array}{}\text{(A9)}& {\displaystyle \frac{\mathrm{d}{H}_{\mathrm{n}}^N}{\mathrm{d}z}}={\displaystyle \frac{\mathrm{d}{H}_{\mathrm{n}}^P}{\mathrm{d}z}}\cdot {\left(\mathrm{1}+{\displaystyle \frac{\mathrm{d}{H}_{\mathrm{n}}^P}{\mathrm{d}z}}\right)}^{-\mathrm{1}}.\end{array}$$

The linear profile of density scale height is given by

$$\begin{array}{}\text{(A10)}& H_{\mathrm{n}}^N\left(z\right)=H_{\mathrm{n}\mathrm{0}}^N\cdot \left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{L_{\mathrm{n}\mathrm{0}}}}\right),\end{array}$$

with

$$\begin{array}{}\text{(A11)}& H_{\mathrm{n}\mathrm{0}}^N={\displaystyle \frac{H_{\mathrm{n}\mathrm{0}}^P}{\mathrm{1}+H_{\mathrm{n}\mathrm{0}}^P/L_{\mathrm{n}\mathrm{0}}}}={\displaystyle \frac{H_{\mathrm{n}\mathrm{0}}^P}{\mathrm{1}+{\mathit{\gamma }}^{-\mathrm{1}}}}.\end{array}$$

Here

$$\begin{array}{}\text{(A12)}& \mathit{\gamma }={\left({\displaystyle \frac{\mathrm{d}{H}_{\mathrm{n}}^P}{\mathrm{d}z}}\right)}^{-\mathrm{1}}\end{array}$$

denotes the inverse gradient of pressure scale height. The inverse gradient of density scale height

$$\begin{array}{}\text{(A13)}& \mathit{\eta }={\left({\displaystyle \frac{\mathrm{d}{H}_{\mathrm{n}}^N}{\mathrm{d}z}}\right)}^{-\mathrm{1}}={\displaystyle \frac{L_{\mathrm{n}\mathrm{0}}}{H_{\mathrm{n}\mathrm{0}}^N}}\end{array}$$

is related to γ through $\mathit{\eta }=\mathit{\gamma }+\mathrm{1}$.

In the non-isothermal case L_n0<∞, integrating $\mathrm{1}/H_{\mathrm{n}}^N$ gives the expression

$$\begin{array}{}\text{(A14)}& \begin{array}{rl}{\mathit{\zeta }}_{\mathrm{0}}& =\underset{z_{\mathrm{0}}}{\overset{z}{\int }}{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{̃}}{z}}{H_{\mathrm{n}}^N\left(\stackrel{\mathrm{̃}}{z}\right)}}=\mathit{\eta }\ln \left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{L_{\mathrm{n}\mathrm{0}}}}\right)\\ & =-\ln {\left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}}\right)}^{-\mathit{\eta }}.\end{array}\end{array}$$

Hence, the altitude profile of number density (Eq. A3) is given by

$$\begin{array}{}\text{(A15)}& N_{\mathrm{n}}\left(z\right)=N_{\mathrm{n}\mathrm{0}}\cdot e^{-{\mathit{\zeta }}_{\mathrm{0}}}=N_{\mathrm{n}\mathrm{0}}\cdot {\left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}}\right)}^{-\mathit{\eta }}.\end{array}$$

In the isothermal limit, η→∞, $\ln \left(\mathrm{1}+\frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}\right)\to \frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}$, and thus ${\mathit{\zeta }}_{\mathrm{0}}\to \frac{z-z_{\mathrm{0}}}{H_{\mathrm{n}\mathrm{0}}^N}$, and $H_{\mathrm{n}\mathrm{0}}^N\to H_{\mathrm{n}\mathrm{0}}^P$ through Eq. (A11).

Appendix B: Electron density profile for linear variations in scale height

Following the approach first presented by Chapman (1931), the ionization rate per unit volume q is expressed in terms of the intensity I of ionizing radiation, the ionization efficiency κ, the angle χ of incident radiation with the atmospheric layer normal vector, the radiation absorption cross-section σ_r, and the neutral density N_n as $q=\mathit{\kappa }\cos \mathit{\chi }\frac{\mathrm{d}I}{\mathrm{d}z}$. Here z is altitude, and the z axis is pointing upwards as before. The function q=q(z) is also called the production function. Although originally proposed for photoionization, the Chapman approach may also be applied to ionization by precipitation of energetic particles as in the auroral region if model variables and coefficients are properly interpreted.

The intensity I satisfies the differential equation

$$\begin{array}{}\text{(B1)}& \mathrm{d}I={\mathit{\sigma }}_{\mathrm{r}}{N}_{\mathrm{n}}I{\displaystyle \frac{\mathrm{d}z}{\cos \mathit{\chi }}}\end{array}$$

with the solution

$$\begin{array}{}\text{(B2)}& I\left(z\right)=I_{\mathrm{\infty }}\exp \left\{{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{r}}}{\cos \mathit{\chi }}}\underset{z_{\mathrm{\infty }}}{\overset{z}{\int }}N\left(\stackrel{\mathrm{̃}}{z}\right)\mathrm{d}\stackrel{\mathrm{̃}}{z}\right\},\end{array}$$

where z_∞ and I_∞ refer to an upper boundary sufficiently remote from the atmospheric layer.

Using $\mathrm{d}I={\mathit{\sigma }}_{\mathrm{r}}{N}_{\mathrm{n}}I\frac{\mathrm{d}z}{\cos \mathit{\chi }}$, the production function q can be rewritten as q = κ σ_rN_nI and thus

$$\begin{array}{}\text{(B3)}& q\left(z\right)=\mathit{\kappa }{\mathit{\sigma }}_{\mathrm{r}}{N}_{\mathrm{n}}\left(z\right)I_{\mathrm{\infty }}\exp \left\{{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{r}}}{\cos \mathit{\chi }}}\underset{z_{\mathrm{\infty }}}{\overset{z}{\int }}N\left(\stackrel{\mathrm{̃}}{z}\right)\mathrm{d}\stackrel{\mathrm{̃}}{z}\right\}.\end{array}$$

The ionization peak altitude z_* is obtained from the condition

$$\begin{array}{}\text{(B4)}& \mathrm{0}={\left.{\displaystyle \frac{\mathrm{d}\ln q}{\mathrm{d}z}}\right|}_{z=z_{*}}={\displaystyle \frac{N_{\mathrm{n}}^{\prime }\left(z_{*}\right)}{N_{\mathrm{n}}\left(z_{*}\right)}}+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{r}}{N}_{\mathrm{n}}\left(z_{*}\right)}{\cos \mathit{\chi }}},\end{array}$$

where the prime denotes differentiation with respect to altitude z. Considering Eqs. (A3) and (A14) gives rise to $N_{\mathrm{n}}\left(z\right)=N_{\mathrm{n}\mathrm{0}}{e}^{-{\mathit{\zeta }}_{\mathrm{0}}}$, and ${\mathit{\zeta }}_{\mathrm{0}}^{\prime }=\mathrm{1}/H_{\mathrm{n}}^N$, and defining the radiation absorption length L_r=L_r(z) by

$$\begin{array}{}\text{(B5)}& L_{\mathrm{r}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathrm{r}}{N}_{\mathrm{n}}}},\end{array}$$

the general ionization peak condition is conveniently expressed as

$$\begin{array}{}\text{(B6)}& H_{\mathrm{n}}^N\left(z_{*}\right)=L_{\mathrm{r}}\left(z_{*}\right)\cos \mathit{\chi }.\end{array}$$

B1 Local representation of electron density

Assuming the neutral temperature T_n varies linearly with altitude z, the altitude dependence of electron density was modeled by Gledhill and Szendrei (1950). Since their formulation does not fit well with the DIPCont nomenclature used in the current report, an independent and extended derivation is presented now. Using $T_{\mathrm{n}}\left(z\right)=T_{\mathrm{n}\mathrm{0}}\cdot \left(\mathrm{1}+\frac{z-z_{\mathrm{0}}}{L_{\mathrm{n}\mathrm{0}}}\right)=T_{\mathrm{n}\mathrm{0}}\cdot \left(\mathrm{1}+\frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}\right)$ and η<∞, the altitude profile of neutral number density can be written in the form

$$\begin{array}{}\text{(B7)}& N_{\mathrm{n}}\left(z\right)=N_{\mathrm{n}\mathrm{0}}{\left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}}\right)}^{-\mathit{\eta }};\end{array}$$

see Appendix A and Eq. (A15). Integration gives

$$\begin{array}{}\text{(B8)}& \underset{z_{\mathrm{\infty }}}{\overset{z}{\int }}{N}_{\mathrm{n}}\left(\stackrel{\mathrm{̃}}{z}\right)\mathrm{d}\stackrel{\mathrm{̃}}{z}=-N_{\mathrm{n}\mathrm{0}}{\displaystyle \frac{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}{\mathit{\eta }-\mathrm{1}}}{\left[{\left(\mathrm{1}+{\displaystyle \frac{\stackrel{\mathrm{̃}}{z}-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}}\right)}^{-(\mathit{\eta }-\mathrm{1})}\right]}_{\stackrel{\mathrm{̃}}{z}=z_{\mathrm{\infty }}}^{\stackrel{\mathrm{̃}}{z}=z}.\end{array}$$

In this LTI modeling context it is safe to assume that the regional temperature increase with altitude is moderate enough to ensure $H_{\mathrm{n}\mathrm{0}}^N<L_{\mathrm{n}\mathrm{0}}$, then η>1. Furthermore, the altitude z_∞ is chosen to be large enough for the contribution from the value at $\stackrel{\mathrm{̃}}{z}=z_{\mathrm{\infty }}$ to be negligible. We obtain

$$\begin{array}{}\text{(B9)}& \underset{z_{\mathrm{\infty }}}{\overset{z}{\int }}{N}_{\mathrm{n}}\left(\stackrel{\mathrm{̃}}{z}\right)\mathrm{d}\stackrel{\mathrm{̃}}{z}=-N_{\mathrm{n}\mathrm{0}}{\displaystyle \frac{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}{\mathit{\eta }-\mathrm{1}}}{\left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}}\right)}^{-(\mathit{\eta }-\mathrm{1})}\end{array}$$

by using Eq. (A14). Defining

$$\begin{array}{}\text{(B10)}& {\mathit{\theta }}_{\mathrm{0}}={\displaystyle \frac{\mathit{\eta }-\mathrm{1}}{\mathit{\eta }}}{\mathit{\zeta }}_{\mathrm{0}}=(\mathit{\eta }-\mathrm{1})\ln \left(\mathrm{1}+{\displaystyle \frac{z-z_{\mathrm{0}}}{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}}\right),\end{array}$$

the radiation intensity profile assumes the form

$$\begin{array}{}\text{(B11)}& {\displaystyle }I\left(z\right)=I_{\mathrm{\infty }}\exp \left\{-{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{r}}{N}_{\mathrm{n}\mathrm{0}}}{\cos \mathit{\chi }}}{\displaystyle \frac{\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N}{\mathit{\eta }-\mathrm{1}}}{e}^{-{\mathit{\theta }}_{\mathrm{0}}}\right\}\text{(B12)}& {\displaystyle }=I_{\mathrm{\infty }}\exp \left\{-{\displaystyle \frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}}{\displaystyle \frac{H_{\mathrm{n}\mathrm{0}}^N}{L_{\mathrm{r}\mathrm{0}}\cos \mathit{\chi }}}{e}^{-{\mathit{\theta }}_{\mathrm{0}}}\right\},\end{array}$$

where L_r0=L_r(z₀). The neutral density (Eq. A15) is rewritten as

$$\begin{array}{}\text{(B13)}& N_{\mathrm{n}}\left(z\right)=N_{\mathrm{n}\mathrm{0}}\exp \left\{-{\displaystyle \frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}}{\mathit{\theta }}_{\mathrm{0}}\right\}\end{array}$$

so that the production function (Eq. B3) assumes the form

$$\begin{array}{}\text{(B14)}& q\left(z\right)={\displaystyle \frac{\mathit{\kappa }{I}_{\mathrm{\infty }}}{L_{\mathrm{r}\mathrm{0}}}}\exp \left\{{\displaystyle \frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}}\left[-{\mathit{\theta }}_{\mathrm{0}}-{\displaystyle \frac{H_{\mathrm{n}\mathrm{0}}^N}{L_{\mathrm{r}\mathrm{0}}\cos \mathit{\chi }}}{e}^{-{\mathit{\theta }}_{\mathrm{0}}}\right]\right\}.\end{array}$$

In the isothermal limit, η→∞, $\frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}\to \mathrm{1}$, ${\mathit{\theta }}_{\mathrm{0}}\to \frac{z-z_{\mathrm{0}}}{H_{\mathrm{n}\mathrm{0}}^N}$, and the isothermal Chapman production function (Chapman, 1931) is recovered.

In static equilibrium of photoionization and quadratic recombination, $q=\mathit{\alpha }{N}_{\mathrm{e}}^{\mathrm{2}}$ with the recombination coefficient α; thus $N_{\mathrm{e}}=\sqrt{q/\mathit{\alpha }}$. Using

$$\begin{array}{}\text{(B15)}& N_{\mathrm{e}\mathrm{0}}=N_{\mathrm{e}}\left(z_{\mathrm{0}}\right)=\sqrt{{\displaystyle \frac{\mathit{\kappa }{I}_{\mathrm{\infty }}}{\mathit{\alpha }{L}_{\mathrm{r}\mathrm{0}}}}}\exp \left\{-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}}{\displaystyle \frac{H_{\mathrm{n}\mathrm{0}}^N}{L_{\mathrm{r}\mathrm{0}}\cos \mathit{\chi }}}\right\},\end{array}$$

we obtain

$$\begin{array}{}\text{(B16)}& \begin{array}{rl}& N_{\mathrm{e}}\left(z\right)=\\ & N_{\mathrm{e}\mathrm{0}}\exp \left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}}\left[-{\mathit{\theta }}_{\mathrm{0}}+{\displaystyle \frac{H_{\mathrm{n}\mathrm{0}}^N}{L_{\mathrm{r}\mathrm{0}}\cos \mathit{\chi }}}\left(\mathrm{1}-e^{-{\mathit{\theta }}_{\mathrm{0}}}\right)\right]\right\}.\end{array}\end{array}$$

B2 Representation of electron density in terms of ionization peak parameters

A meaningful regional representation of the electron density can be constructed by means of the ionization peak parameters. For a given incident radiation angle χ, the altitude z_* of the electron density maximum can be expressed in local parameters as follows:

$$\begin{array}{}\text{(B17)}& z_{*}=z_{\mathrm{0}}+\mathit{\eta }{H}_{\mathrm{n}\mathrm{0}}^N\left[{\mathrm{\Gamma }}^{\mathrm{1}/(\mathit{\eta }-\mathrm{1})}-\mathrm{1}\right],\end{array}$$

where

$$\begin{array}{}\text{(B18)}& \mathrm{\Gamma }={\displaystyle \frac{H_{\mathrm{n}\mathrm{0}}^N}{L_{\mathrm{r}\mathrm{0}}\cos \mathit{\chi }}}.\end{array}$$

The electron density peak value $N_{{\mathrm{e}}^{*}}=N_{\mathrm{e}}(z=z_{*})$ is

$$\begin{array}{}\text{(B19)}& N_{{\mathrm{e}}^{*}}=N_{\mathrm{e}\mathrm{0}}\exp \left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}}\left[-\ln \mathrm{\Gamma }+\mathrm{\Gamma }-\mathrm{1}\right]\right\}.\end{array}$$

With z_* as the reference altitude, $z_{\mathrm{0}}=z_{*}$, we can take advantage of the condition (Eq. B6) $H_{\mathrm{n}\mathrm{0}}^N=H_{{\mathrm{n}}^{*}}^N=L_{r*}\cos \mathit{\chi }=L_{\mathrm{r}\mathrm{0}}\cos \mathit{\chi }$; thus

$$\begin{array}{}\text{(B20)}& N_{\mathrm{e}}\left(z\right)=N_{{\mathrm{e}}^{*}}\exp \left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathit{\eta }}{\mathit{\eta }-\mathrm{1}}}\left[-{\mathit{\theta }}_{*}+\mathrm{1}-e^{-{\mathit{\theta }}_{*}}\right]\right\},\end{array}$$

where ${\mathit{\theta }}_{*}={\mathit{\theta }}_{*}\left(z\right)=(\mathit{\eta }-\mathrm{1})\ln \left(\mathrm{1}+\frac{z-z_{*}}{\mathit{\eta }{H}_{{\mathrm{n}}^{*}}^N}\right)$, and $H_{{\mathrm{n}}^{*}}^N$ denotes the density scale height at $z=z_{*}$. This representation shows that χ is only an implicit parameter of the electron density model and cannot be inferred from knowledge of the peak parameters.

Appendix C: Orbit approximation around perigee

Consider a Kepler orbit with radial distance r=r(t) and azimuth ϕ=ϕ(t), where t denotes time. Distance and velocity at perigee are R_per and V_per, respectively. The corresponding variables at apogee are R_apo and V_apo, the gravitational constant is G, and the planetary mass is M. Combining the conservation laws for angular momentum,

$$\begin{array}{}\text{(C1)}& r^{\mathrm{2}}\dot{\mathit{\varphi }}=R_{\mathrm{apo}}{V}_{\mathrm{apo}}=R_{\mathrm{per}}{V}_{\mathrm{per}},\end{array}$$

and total energy E (here normalized by the test mass m),

$$\begin{array}{}\text{(C2)}& {\displaystyle \frac{E}{m}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\dot{r}}^{\mathrm{2}}+r^{\mathrm{2}}{\dot{\mathit{\varphi }}}^{\mathrm{2}}\right)-{\displaystyle \frac{GM}{r}}\text{(C3)}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_{\mathrm{per}}^{\mathrm{2}}-{\displaystyle \frac{GM}{R_{\mathrm{per}}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{V}_{\mathrm{apo}}^{\mathrm{2}}-{\displaystyle \frac{GM}{R_{\mathrm{apo}}}},\end{array}$$

yields the following expression for the perigee velocity in terms of perigee and apogee distances:

$$\begin{array}{}\text{(C4)}& V_{\mathrm{per}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{2}GM{R}_{\mathrm{apo}}}{R_{\mathrm{per}}(R_{\mathrm{apo}}+R_{\mathrm{per}})}}={\displaystyle \frac{\mathrm{2}{g}_{\mathrm{per}}{R}_{\mathrm{per}}{R}_{\mathrm{apo}}}{(R_{\mathrm{apo}}+R_{\mathrm{per}})}},\end{array}$$

where $g_{\mathrm{per}}=\frac{GM}{R_{\mathrm{per}}^{\mathrm{2}}}$ is the value of Earth's gravitational acceleration at geocentric distance R_per. The radial velocity $\dot{r}$ satisfies

$$\begin{array}{}\text{(C5)}& {\displaystyle }{\dot{r}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{2}E}{m}}+{\displaystyle \frac{\mathrm{2}GM}{r}}-{\displaystyle \frac{(r^{\mathrm{2}}\dot{\mathit{\varphi }}{)}^{\mathrm{2}}}{r^{\mathrm{2}}}}\text{(C6)}& {\displaystyle }={\displaystyle \frac{\mathrm{2}E}{m}}+{\displaystyle \frac{\mathrm{2}GM}{r}}-{\displaystyle \frac{R_{\mathrm{per}}^{\mathrm{2}}{V}_{\mathrm{per}}^{\mathrm{2}}}{r^{\mathrm{2}}}}.\end{array}$$

Differentiating this expression and dividing by $\mathrm{2}\dot{r}$ yields

$$\begin{array}{}\text{(C7)}& \ddot{r}=-{\displaystyle \frac{GM}{r^{\mathrm{2}}}}+{\displaystyle \frac{R_{\mathrm{per}}^{\mathrm{2}}{V}_{\mathrm{per}}^{\mathrm{2}}}{r^{\mathrm{3}}}}.\end{array}$$

Evaluation at perigee r=R_per gives

$$\begin{array}{}\text{(C8)}& {\left.\ddot{r}\right|}_{r=R_{\mathrm{per}}}=-{\displaystyle \frac{GM}{R_{\mathrm{per}}^{\mathrm{2}}}}+{\displaystyle \frac{R_{\mathrm{per}}^{\mathrm{2}}{V}_{\mathrm{per}}^{\mathrm{2}}}{R_{\mathrm{per}}^{\mathrm{3}}}}=-{\displaystyle \frac{GM}{R_{\mathrm{per}}^{\mathrm{2}}}}+{\displaystyle \frac{V_{\mathrm{per}}^{\mathrm{2}}}{R_{\mathrm{per}}}}.\end{array}$$

Inserting the expression for $V_{\mathrm{per}}^{\mathrm{2}}$ yields

$$\begin{array}{}\text{(C9)}& {\left.\ddot{r}\right|}_{r=R_{\mathrm{per}}}={\displaystyle \frac{GM}{R_{\mathrm{per}}^{\mathrm{2}}}}{\displaystyle \frac{R_{\mathrm{apo}}-R_{\mathrm{per}}}{R_{\mathrm{apo}}+R_{\mathrm{per}}}}=g_{\mathrm{per}}\mathit{\epsilon },\end{array}$$

where $\mathit{\epsilon }=\frac{R_{\mathrm{apo}}-R_{\mathrm{per}}}{R_{\mathrm{apo}}+R_{\mathrm{per}}}$ is the orbital eccentricity. The altitude z is related to radial distance r and the Earth's planetary radius R_E through $z=r-R_{\mathrm{E}}$. At perigee, t=0, and z=z_per. The parameter $a_{\mathrm{per}}=\ddot{z}(t=\mathrm{0})$ coincides with the radial acceleration at perigee ${\left.\ddot{r}\right|}_{r=R_{\mathrm{per}}}$. Hence, orbital altitudes around perigee are approximately given by the quadratic function

$$\begin{array}{}\text{(C10)}& z\left(t\right)\simeq z_{\mathrm{per}}+{\displaystyle \frac{a_{\mathrm{per}}}{\mathrm{2}}}{t}^{\mathrm{2}}.\end{array}$$

To the same approximation order, the angular momentum conservation condition $r^{\mathrm{2}}\dot{\mathit{\varphi }}=R_{\mathrm{per}}{V}_{\mathrm{per}}$ can be integrated to yield approximate azimuths ϕ=ϕ(t). In $\dot{\mathit{\varphi }}=R_{\mathrm{per}}{V}_{\mathrm{per}}/r^{\mathrm{2}}$ insert $r=r\left(t\right)=R_{\mathrm{per}}+\frac{a_{\mathrm{per}}}{\mathrm{2}}{t}^{\mathrm{2}}$, then expand

$$\begin{array}{}\text{(C11)}& {\left(R_{\mathrm{per}}+{\displaystyle \frac{a_{\mathrm{per}}}{\mathrm{2}}}{t}^{\mathrm{2}}\right)}^{-\mathrm{2}}\simeq R_{\mathrm{per}}^{-\mathrm{2}}\left(\mathrm{1}-{\displaystyle \frac{a_{\mathrm{per}}}{R_{\mathrm{per}}}}{t}^{\mathrm{2}}\right)\end{array}$$

and integrate $\mathrm{d}\mathit{\varphi }=R_{\mathrm{per}}{V}_{\mathrm{per}}{r}^{-\mathrm{2}}\mathrm{d}t$ to obtain

$$\begin{array}{}\text{(C12)}& {\displaystyle }\mathit{\varphi }\left(t\right)\simeq {\displaystyle \frac{V_{\mathrm{per}}}{R_{\mathrm{per}}}}\underset{\mathrm{0}}{\overset{t}{\int }}\left(\mathrm{1}-{\displaystyle \frac{a_{\mathrm{per}}}{R_{\mathrm{per}}}}{\stackrel{\mathrm{̃}}{t}}^{\mathrm{2}}\right)\mathrm{d}\stackrel{\mathrm{̃}}{t}\text{(C13)}& {\displaystyle }={\displaystyle \frac{V_{\mathrm{per}}}{R_{\mathrm{per}}}}\cdot \left(t-{\displaystyle \frac{a_{\mathrm{per}}}{\mathrm{3}{R}_{\mathrm{per}}}}{t}^{\mathrm{3}}\right).\end{array}$$

The corresponding horizontal distances at the Earth's surface are then given by $x=x\left(t\right)=R_{\mathrm{E}}\mathit{\varphi }\left(t\right)$. By using Eq. (C9), this can be further processed to yield

$$\begin{array}{}\text{(C14)}& x\left(t\right)\simeq {\displaystyle \frac{R_{\mathrm{E}}{V}_{\mathrm{per}}t}{R_{\mathrm{per}}}}\cdot \left(\mathrm{1}-{\displaystyle \frac{\mathit{\epsilon }}{\mathrm{3}}}{\displaystyle \frac{g_{\mathrm{per}}{t}^{\mathrm{2}}}{R_{\mathrm{per}}}}\right).\end{array}$$

The leading term is ground distance for a circular orbit. The correction produced by the second term is proportional to eccentricity.

Appendix D: Parametrization of horizontal electron density variations

In the initial version of the DIPCont package, the horizontal variability in electron density profiles is controlled by the keyword argument LTIModelType. Setting LTIModelType='NeAuroralZoneCrossing' produces two electron density maxima along the horizontal (latitudinal) axis as observed by a polar-orbiting satellite when crossing the auroral oval; see Figs. 1, 7, and 8. More specifically, the horizontal (x) variations in peak altitude $z_{*}=z_{*}\left(x\right)$ and peak electron density $N_{{\mathrm{e}}^{*}}=N_{{\mathrm{e}}^{*}}\left(x\right)$ in Eq. (15) are prescribed by the ad hoc parametrizations

$$\begin{array}{}\text{(D1)}& {\displaystyle }{z}_{*}\left(x\right)=z_{*,\min }+\mathrm{\Delta }{z}_{*}\cdot f\left(x\right),\text{(D2)}& {\displaystyle }{N}_{{\mathrm{e}}^{*}}\left(x\right)=N_{e*,\max }-\mathrm{\Delta }{N}_{{\mathrm{e}}^{*}}\cdot f\left(x\right),\end{array}$$

with

$$\begin{array}{}\text{(D3)}& f\left(x\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left\{\mathrm{1}+\cos \left({\displaystyle \frac{\mathrm{4}\mathit{\pi }x}{x_{\mathrm{R}}-x_{\mathrm{L}}}}\right),\right\}\end{array}$$

so that f=f(x) varies between zero and one. The parameters x_L and x_R are the horizontal boundaries of the modeling domain, here chosen to be x_L= −2000 km and x_R= 2000 km. The values of the electron density peak parameters used in the model runs leading to Figs. 1, 7, and 8 are as follows: $z_{*,\min }=$ 110 km, $\mathrm{\Delta }{z}_{*}=$ 10 km, $N_{e*,\max }=\mathrm{1.5}\cdot {\mathrm{10}}^{\mathrm{11}}{\mathrm{m}}^{-\mathrm{3}}$, $\mathrm{\Delta }{N}_{{\mathrm{e}}^{*}}=\mathrm{0.5}\cdot {\mathrm{10}}^{\mathrm{11}}{\mathrm{m}}^{-\mathrm{3}}$.

All LTI model parameters for the simulation runs of the current report, including the horizontal electron density profile parameters, are provided in the configuration file DIPContBas.py as part of the Supplement.