© D. Banyś et al., Published by EDP Sciences 2024

1 Introduction

The plasmasphere and its upper boundary region, the plasmapause, are important parts of the near-Earth space environment (cf. Darrouzet et al., 2009). The co-rotating plasma extends from the ionosphere up to the plasmapause characterized by a sharp decrease of the electron density towards the non-co-rotating magnetosphere (e.g. Carpenter, 1963). Mainly driven by the dynamics of the solar wind and the geomagnetic field conditions, the plasmapause position and shape are highly variable. The plasmasphere is refilled after geomagnetic storm-induced erosion by ionospheric ions, primarily protons, up to the plasmapause at about 5–7 Earth radii (R_E) in the geomagnetic equatorial plane (cf. Lemaire & Gringauz, 1998). Under perturbed conditions, the enhanced magnetospheric convection electric field causes a strong contraction of the plasmasphere measurable as an inward motion of the plasmapause position L_pp down to about L = 2 (cf. Obana et al., 2019 and references therein).

The plasmasphere is not only of interest for studying solar-terrestrial relationships but also for space-based applications in communication, navigation and remote sensing. A better understanding of its behaviour will enhance the accuracy and reliability of these applications, e.g., by estimating total electron content (TEC) induced range errors in applications of Global Navigation Satellite Systems (GNSS) as reported e.g. by Lunt et al. (1999), Yizengaw et al. (2008) and Jakowski & Hoque (2018).

Early plasmapause models by Carpenter & Anderson (1992) and Gallagher et al. (1988, 2000) addressed this behaviour by referring to the geomagnetic Kp index (Matzka et al., 2021). The Carpenter & Anderson (1992) plasmapause model, based on whistler measurements (Carpenter & Smith, 1964), has been largely used for radiation belt studies over many years (cf. Ripoll et al., 2022). Gallagher et al. (2000) established a unified model, the so-called Global Core Plasma Model (GCPM), which is based on comprehensive previous studies and considers the space weather impact by relying on the geomagnetic Kp index. Different geomagnetic indices such as Kp, Dst and AE have been used by O’Brien & Moldwin (2003) to present an empirical model based on 900 Combined Release and Radiation Effects Satellite (CRRES) electron density profiles. Unfortunately, the methodology to define reliable plasmapause profiles remains unclear. It is worth mentioning that Dst has been shown to correlate slightly better than Kp.

Considerable progress has been made in theoretical and empirical modelling of the plasmasphere and plasmapause position in the recent two decades (cf. Ripoll et al., 2022 and references therein). Reinisch et al. (2009) subjected the augmentation of empirical plasmasphere models by the inclusion of data from IMAGE (Imager for Magnetopause-to-Aurora Global Exploration) and Cluster. Pierrard et al. (2009) investigated the performance and progress of physics-based models in terms of modelling approaches, behaviour at geomagnetic storms and couplings of the plasmasphere to ionosphere and magnetosphere. Heilig & Lühr (2013) analyzed more than 20,000 plasmapause crossings obtained from the CHAllenging Minisatellite Payload (CHAMP). Their model depends on magnetic local time MLT and Kp. Liu et al. (2015) developed a dynamic plasmapause model based on the Time History of Events and Macroscale Interactions during Substorms (THEMIS) data from 2009 to 2013, including 5878 crossing events over different MLT sectors. The model utilizes 5 different parameters having a high correlation with L_pp. Huba & Krall (2013) simulated for the first time the 3D plasmasphere with the first-principles physics-based model SAMI3, with a special focus on plasmaspheric behaviour during storms investigating the evolvement of plumes. Zhelavskaya et al. (2017) have developed a Plasma density in the Inner magnetosphere Neural network-based Empirical (PINE) model for reconstructing the plasmasphere density distribution and its dynamics for 2 ≤ L_pp ≤ 6 at all local times. The neural network utilizes a wide variety of input parameters such as geomagnetic indices, solar wind data and solar wind coupling functions. A New Solar Wind driven Global Dynamic Plasmapause (NSW-GDP) model has been developed by He et al. (2017) based on a large database utilizing multiple sources from 1977 to 2015 covering four solar cycles. The model refers to solar wind parameters as well as geomagnetic indices. Guo et al. (2021) have developed a neural network model of the plasmapause location using Van Allen Probes data obtained during the period from 2012 to 2017. The model achieves good results when using only AE or Kp indices as input parameters. Botek et al. (2021) applied the Space Weather Integrated Forecasting Framework (SWIFF) Plasmasphere Model (SPM), a 3D kinetic plasmasphere model, to Van Allen Probes data in order to improve the model’s plasma trough equations and compared the results, among others, with Arase data. The plasmapause model developed by Ripoll et al. (2022) utilizes NASA Van Allen Probes data and Integrated Science measurements for extracting the plasmasphere boundaries during 2012–2019. The gradient method for locating the plasmapause is equivalent to the 100 cm⁻³ density thresholds for the plasmasphere outer edge (L₁₀₀). The L₁₀₀ boundary starts varying with MLT for Kp > 2.

The Neustrelitz/ESOC PlasmaPause Model (NEPPM) presented here has been developed to serve as a robust and easy-to-use background model in operational space weather services. For detecting the plasmapause position, an effective automatic procedure is proposed by formulating clear conditions for electron density gradient shapes. The newly applied criteria for automated filtering and subsequent model fitting shall be illustrated in the paper by using IMAGE and Van Allen Probes data. The database is explained in Section 2. The therein-introduced conditions for automatically selecting reliable profile data are described in detail. Applying this procedure is expected to obtain a conveniently restricted set of electron density profiles as input for building the plasmapause model. Section 3 describes the ellipse-based modelling approach assumed for the NEPPM. The ellipse-based L_pp-function is embedded in a 3D modelling algorithm allowing to fit of the function parameters directly from observed plasmapause torus position vectors ${\vec{R}}_{\mathrm{pp}}$. In addition, this way of 3D modelling allows accounting for a non-dipole Earth magnetic field. Section 4 summarizes the results obtained with the NEPPM and compares them with other modelling approaches, in particular the GCPM, before concluding the paper in Section 5. Attached are two annexes: Annex A details formulae developments of algorithms needed for the 3D plasmapause modelling, and Annex B provides an efficient coordinate transformation between geographic and solar-magnetic (SM) coordinates.

2 Database and preparation

For developing the NEPPM, over 3000 plasmapause crossings have been utilized. The selected electron density data have been recorded during plasmapause crossings onboard the IMAGE satellite (cf. https://image.gsfc.nasa.gov/) and the Van Allen Probes (https://www.nasa.gov/van-allen-probes).

2.1 IMAGE RPI data

On 25 March 2000, NASA launched the IMAGE (Imager for Magnetopause-to-Aurora Global Explorer) satellite, which operated successfully until 18 December 2005. The Radio-Plasma Imager (RPI) on board took passive in-situ plasma wave measurements suitable for reconstructing electron densities and deriving plasmapause positions. Detailed information on the instrument is given in Reinisch et al. (2001a, 2001b, 2009) and Goldstein et al. (2003). The passive observations can be displayed as a function of frequency over time, i.e., a dynamic spectrum (cf. Galkin et al., 2004). Considering the upper-hybrid band, the continuum edge, and band emissions, a semi-automatic fitting technique in the dynamic spectra was then applied to derive electron density values from the dynamic RPI spectrogram (Gerzen et al., 2014). Data in the Level Zero Telemetry L0 format is limited to the beginning of 2001 until the end of 2005. This summarizes over 200,000 individual electron density values together with their respective observation time, location, and L-shell distance. For further information, we refer to Denton et al. (2012), Gerzen et al. (2014), and the website of the Space Science Lab of the University of Massachusetts Lowell (UML, http://ulcar.uml.edu/rpi.html). Our NEPPM model is based on these electron densities derived from the passive RPI data, their respective geographic position vectors and times, and furthermore, the equivalent Van Allen Probes data, detailed in the next section.

2.2 Van Allen Probes reconstructed electron density data

The Van Allen Probes, Radiation Belt Storm Probes RBSP-A and B, were launched on 20 August 2012 and operated for over 7 years. Onboard, the EMFISIS instrument suite (Electric and Magnetic Fields Instrument Suite and Integrated Science; Kletzing et al., 2013) provided electric field measurements in the frequency range of 10–487 kHz, so that the upper hybrid resonance band could be identified allowing precise electron density estimations. The NEPPM uses the electron densities, position vectors, and times that are provided by the NURD (Neural-network-based Upper Hybrid Resonance Determination) algorithm published in Zhelavskaya et al. (2016). The feedforward neural network NURD was trained on the Van Allen Probes data to derive the upper hybrid resonance frequency and finally the electron number density. The database consists of over 24,000,000 single electron density values for Van Allen Probes RBSP-A and over 16,000,000 for RBSP-B. To ensure a reliable database, we only use orbits with a quality flag of 1, representing good quality, while disregarding orbits with questionable or missing data.

2.3 Plasmapause detection

Due to the highly eccentric orbits of IMAGE and the Van Allen Probes, one can relate electron density measurements to the altitude. For retrieving distinguished electron density profiles, in particular, for IMAGE, the measurements are split, first, at gaps larger than 5 min, and second, at their global maxima and minima in case these are not at the beginning or end of the respective time series. Each section obtained this way is considered physically plausible and is treated as an electron density profile. Nevertheless, too short profiles may produce unwanted artefacts. We hence only use profiles with at least 10 measurement points, considering a common IMAGE profile of about 50 measurements.

Many publications suggest looking for density drops of factor 5 within an interval of 0.5L (Moldwin et al., 2002; Liu et al., 2015; Guo et al., 2021). However, this is quite a conservative criterion excluding rapid density fluctuations during active conditions (Moldwin et al., 2002). Also, this bound will not detect smaller but well-shaped plasmapause gradients. For including these potentially valuable profiles, we use a refined set of criteria based on averaging to a comparable step size and relative bounds.

Predefinitions

To avoid outliers, raw electron density measurements $N_e^o$ are generally replaced by smoothed measurements via a 3-point moving average approach

$$N_e\left(l_k\right) :=\frac{1}{3}{\sum }_{i=-1}^1{N}_e^o\left(l_{k+i}\right),$$(2.1)

where l_k runs through the profile on which L_pp is considered. Since the provided electron density profiles of the Van Allen Probes are of higher resolution than the IMAGE ones, we are furthermore reducing the resolution of N_e to an average of 50 values per profile in total in order to apply the same plasmapause detection algorithm.

The point

$$L≔\underset{l}{\mathrm{argmin}}{\partial }_l\log N_e\left(l\right)$$(2.2)

shall represent the position of the steepest slope and possible L_pp (cf. red dot in Fig. 1a–f), where the derivative is defined via central differences

$${\partial }_l\log N_e\left(l_k\right):=\frac{\log N_e\left(l_{k+1}\right)- \log N_e\left(l_{k-1}\right)}{l_{k+1}-l_{k-1}}.$$(2.3)

Figure 1 Sample profiles obtained from IMAGE demonstrating the use of the different conditions for Lpp detection (red dots): Each case a-f shows the original profile in blue and its average in black (upper left), condition 1 (upper right), condition 2 (lower left), and condition 3 (lower right). The red dotted lines represent the mean values $\bar{\cdot }$ for each condition, the red dashed lines the standard deviation 2σ(⋅) of conditions 1 and 2, and the orange dashed lines indicate the standard deviation 1.5 σ(⋅) of condition 3. The orange dots are the next local minima. Condition 3 requires that they are above the orange dashed line if they are outside of ℒ (indicated green).

Now, all of the following conditions must be satisfied to regard L as L_pp. Measurements that are in conflict with one or more of them are excluded. Typical profiles for each condition are shown in Figure 1.

Condition 1

At the L value of the steepest decrease, the derivative must be much smaller than for the remainder, i.e.

$${\partial }_l\log N_e\left(L\right)< \overline{{\partial }_l\log N_e}-2\sigma \left({\partial }_l\log N_e\right),$$(2.4)

where $\bar{\cdot }$ stands for the arithmetic mean and σ(⋅) for the standard deviation (cf. the small gradients in Fig. 1b and c).

Condition 2

The absolute change of N_e is expected to be small in relation to the other changes, so we assume

$$\Delta \log N_e\left(L\right)< \overline{\Delta \log N_e}-2\sigma (\Delta \log N_e)$$(2.5)

with ΔlogN_e(l_k) ≔ logN_e(l_k) − logN_e(l_k−1) being the difference of the neighboring values (see Fig. 1b and d). This covers more global aspects than the very local Condition 1 and is attributed mainly to the fact that we are stuck with an inhomogeneous grid.

Condition 3

Finally, we ensure that there is only one distinctive peak in negative slopes, i.e. a unique L_pp. Still, a certain amount of noise is normal and acceptable. To this end, ℒ shall be the largest interval containing L so that Condition 1 is valid for all points inside.

Then for all local minima l outside of ℒ, a sharpened reversal of Condition 1 shall be satisfied, i.e.,

$${\partial }_l\log N_e\left(l\right) \ge \overline{{\partial }_l\log N_e}-1.5 \sigma \left({\partial }_l\log N_e\right).$$(2.6)

If there are two separate peaks or the data is too noisy, Condition 3 is not fulfilled (see Fig. 1e and f). However, a plateau or two very close peaks—so that the second minimum is contained in ℒ — will be counted as one (deformed) peak.

Thus, a good profile showing a plasmapause position L_pp is given by a single sharp gradient and change in N_e by meeting all 3 conditions, such as presented in Figure 1a. By this method, over 500 plausible plasmapause crossings for IMAGE, over 1500 for Van Allen Probes RBSP-A and over 1100 for Van Allen Probes RBSP-B are extracted and used as the basis for our model.

3 Modelling approach

As mentioned before, the plasmapause location depends on a broad spectrum of geophysical and external space weather conditions and is, therefore, highly dynamic. As O’Brien & Moldwin (2003) have shown, L_pp correlates quite well with several geomagnetic indices such as Kp, AE or Dst. The best correlation was obtained with Dst. Perhaps this result might be explained by the fact that Dst is a prominent measure of the ring current intensity. The ring current intensity is closely related to the driving force, the convection electric field, which essentially forms the shape of the plasmasphere, including the plasmapause location. Consequently, we use Dst as the main driver for the model. With regard to the influence of the Earth’s magnetic field on the plasmapause torus, we refer to solar-magnetic (SM) coordinates (e.g. Laundal & Richmond, 2016) in the following. Given the general shape of the plasmapause position in the magnetic equatorial plane as shown in many papers (e.g., Carpenter & Anderson, 1992; Heilig & Lühr, 2013; O’Brien & Moldwin, 2003; Ripoll et al., 2022) we define an ellipse in the equatorial plane with the Earth residing in one of the two focal points. The ellipse parameters, i.e. semi-major axis, eccentricity, and orientation of the ellipse in the MLT plane, are described as functions of the Dst index. Then, this approach enables an easy formulation of L_pp as a function of MLT by estimating the semi-major axis, the eccentricity, and the alignment of the ellipse with respect to the MLT axis along the midnight-noon line from the measurements (Fig. 2). Thereby, the diurnal maximum of the plasmapause position, $L_{{\mathrm{pp}}_{\max }}$, typically points to the evening-afternoon hours and the Earth-nearest point, $L_{{\mathrm{pp}}_{\min }}$, is located in the opposite direction. The ellipse radii, among them $L_{{\mathrm{pp}}_{\max }}$ and $L_{{\mathrm{pp}}_{\min }}$, are equal to the distances from the focal point coinciding with the Earth to the ellipse periphery, in L-units.

Figure 2 Meaning of the orientation angle ϖ, the semi-major axis a, $L_{{\mathrm{pp}}_{\max }}$ and $L_{{\mathrm{pp}}_{\min }}$.

3.1 Modelling of L_pp by an ellipse approach

As stated above, in our approach, an ellipse is employed to describe in the geomagnetic equatorial plane the principal plasmapause shape in terms of L_pp, whereby one ellipse focal point coincides with the centre of the Earth and the ellipse maximum radius points into the direction of the plasmapause bulge. Size (semi-major axis), shape (eccentricity), and orientation (orientation angle of the line of apsides) of the ellipse is, in turn, described as functions of Dst and MLT. For the anticipated task, the L_pp-ellipse is modelled with the formula describing the orbital radius of a satellite (e.g., Escobal, 1965)

$$L_{\mathrm{pp}}=\frac{a·\left(1-e^2\right)}{1+e·\cos \upsilon },$$(3.1)

where

a … semi-major axis,

e … eccentricity,

υ … true anomaly.

In the above formulation, the bulge location coincides then with the orbital apogee, i.e. υ = 180°. For the application, L_pp the true anomaly is expressed as

$$\upsilon =x+\varpi ,$$(3.2)

where

x … magnetic local time expressed in radians, i.e. $x=\frac{\pi }{12 \mathrm{h} }·\mathrm{MLT}$,

ϖ … orientation angle to align the ellipse line of apsides (υ = 0°, 180°) to the bulge direction in the MLT system, i.e. w.r.t. the direction of the Sun (midnight-noon line), Figure 2.

Due to the eccentricity, the progress of true anomaly along the ellipse periphery is not linear, while MLT is. However, with regard to the low eccentricities to be dealt with in L_pp modelling, this effect is neglected here, and the sum described by equation (3.2) is considered sufficient. Substitution of equation (3.2) into equation (3.1) leads to equation (3.3):

$$L_{\mathrm{pp}}=\frac{a·\left(1-e^2\right)}{1+e·\cos \left(x+\varpi \right)}.$$(3.3)

The orientation angle ϖ can then be interpreted as follows, Figure 2: When viewing from the North Pole, MLT is counted counter-clockwise from midnight via noon (direction to the Sun) back to midnight. In that system, the bulge peak appears typically around MLT = 18 h (in Fig. 2 around MLT = 16 h, enhanced solar activity). On the other hand, in the ellipse system, the apogee, corresponding to the bulge location, is achieved at υ = 180°. And the task of the orientation angle ϖ is to get υ = 180° coinciding with a given bulge location, in Figure 2 at M = 16 h. Thus, the relation between ϖ and M is

$$\varpi =\frac{\pi }{12 \mathrm{h}}·\left(12 \mathrm{h}-M\right),$$(3.4)

where

M… MLT of bulge.

At low to moderate solar activities, the bulge is typically located in the evening hours while moving more and more into the early afternoon hours when solar activity increases. However, the bulge will never appear in the morning and forenoon hours, i.e. M > 12 h, and thus ϖ is always negative, with decreasing magnitude at increasing solar activity.

Different Dst levels were defined, and then for each Dst level, a set of ellipse parameters (a,e,ϖ) was fitted to the plasmapause values obtained from the plasmapause detection, Section 2.3, whereby ϖ=ϖ(MLT). The selected Dst levels were, in Table 1:

For Dst > −10 nT, the shape of the plasmapause in the equatorial plane is almost circular (as will be shown in Section 4). Therefore, the eccentricity and direction of the ellipse can no longer be determined properly. Hence, this range is not considered in the fitting process of a(Dst), e(Dst), and M(Dst) in Section 3.3, but is nonetheless presented in the overall results of Section 4. Due to the lack of strong geomagnetic storms (here indicated with −500 nT < Dst < −100 nT) and extreme events (given by Dst < −500 nT), which would prevent proper fitting, Dst level 6 has been given an overlap with the richer level 5. Hence, level 6 should be regarded with caution.

3.2 Fitting L_pp-parameters with a 3D modelling approach

The IMAGE and Van Allen Probes recorded 3D plasmapause positions in the SM system. Of these, L_pp-values were derived, indicating corresponding plasmapause locations in the geomagnetic equatorial plane. Both 3D plasmapause positions and derived L_pp-values are listed in the IMAGE and Van Allen Probes data records. It is not fully clear how the L_pp were projected from higher latitudes to the equatorial plane, probably by applying the dipole assumption L_pp = R_pp/cos²φ_m. Here φ_m is the geomagnetic latitude, counted from the geomagnetic equator, positively along the local geomagnetic meridian towards the northern geomagnetic pole, and negatively towards the southern geomagnetic pole. Thus, a 3D approach, that will be described in the following, has been introduced, which embeds the L_pp-function, equation (3.3), and allows to fit its parameters directly from observed plasmapause torus position vectors ${\vec{R}}_{\mathrm{pp}}$. This allows to account for a non-dipole geomagnetic field and avoids a projection of nonzero latitudinal measured radii R_pp into the geomagnetic equatorial plane according to the L_pp = R_pp/cos²φ_m rule, i.e., assuming a dipole, where φ_m is the geomagnetic latitude. In this way, relations of the physics in the background can be captured. The attempt presented in the following is initially derived for the dipole case but can very easily also be applied to a non-dipole field, as will be shown later.

McIlwain (1961, 1965) proposes the following set of 2D coordinates for displaying phenomena in the geomagnetic field:

1. the magnitude of the geomagnetic field vector $\left|\vec{B}\right|$,

2. the invariant $I={\oint }_Q^{\mathrm{Q\text{'}}}{P}_{\left|\right|} \mathrm{d}s$.

I is an invariant along a magnetic field line, and states that charged particles spiral forth and back, i.e., they oscillate along the field line between two mirror points Q and Q'. The quantity P_|| is a momentum which, if no force acts from the exterior to the system, will not change. Thus, also the invariant I remains constant between the two mirror points. At the geomagnetic equator, $\left|\vec{B}\right|$ is minimal along a field line, and the mirror points coincide. At the same time, typically, Q and Q' are located symmetrically and equidistantly on both sides of the equatorial plane. Therefore, $\left(\left|\vec{B}\right|,I\right)=\mathrm{const}$ establishes rings along latitude circles, of which mirrored rings share the same value pairs, delimiting the oscillation along the magnetic field lines. These lines are lying on shells, represented by the parameter L.

Since I is too cumbersome for practical use, a new coordinate $L=F\left(\left|\vec{B}\right|,I\right)$ was introduced by McIlwain (1961) instead. For the computation of L, McIlwain (1961, 1965) used polynomial representations based on coefficient tables listing required numbers according to the best available geomagnetic field model at the time of his publications. For practical applications, O’Brien et al. (1962) and McIlwain (1965) defined the so-called invariant latitude

$${\phi }_m=\arccos \left(\sqrt{\frac{1}{L}}\right).$$(3.5)

It is obtained by solving the relation L_pp = R_pp/cos²φ_m for the geomagnetic latitude φ_m with a unit radius R = 1, and can therefore be considered as a geomagnetic latitude, too. Laundal & Richmond (2016) employ similar definitions in some of their coordinate systems. While with the invariant latitude formula, a latitude value is computed for a given L-value, in the application here L_pp values are determined during plasmapause detection, cf. Section 2, and a cos²φ_m term is computed now with a given $\vec{B}$ vector (dipole or non-dipole). In this way, both approaches, invariant latitude and modifying the cos²φ_m term, are somehow equivalent ways to take the non-dipole aspect into account.

In the 3D modelling approach, the full 3D position vectors ${\vec{R}}_{\mathrm{pp}}$, obtained from the plasmapause detection are taken as observables. Regarding the above considerations, the vectorised dependency between latitude-determined and equatorial radii can be written with spherical coordinates (φ_m, Λ) induced by the SM system and radius based on the McIlwain (1961) L-parameter as

$$\vec{R}=\underset{\left(1\right)}{\underbrace{L}}·\underset{\left(2\right)}{\underbrace{{\cos }^2{\phi }_m}}·\underset{\left(3\right)}{\underbrace{\left\{\begin{array}{c}\cos {\phi }_m·\cos \Lambda \\ \cos {\phi }_m·\sin \Lambda \\ \sin {\phi }_m\end{array}\right\}}},$$(3.6)

where

• $\vec{R}$ … 3D Cartesian vector indicating a geocentric plasmaspheric position in the SM,

• L … L-shell value at plasmaspheric position,

• φ_m … geomagnetic latitude,

• Λ … solar longitude, i.e., longitude in the plane of the geomagnetic equator relative to the longitude of the Sun.

When plotting $\vec{R}$ vectors with constant L over an equidistant grid of SM latitudes and longitudes (φ_m, Λ), one obtains a torus-like object (see Fig. A1 in Annex A).

Equation (3.6) consists of three factors (1)–(3). The last, braced term (3) is the unit vector $\vec{\epsilon }=\vec{R}\left|\vec{R}\right|$ pointing to the position in space, and the other two are scalar quantities. When describing the plasmapause position ${\vec{R}}_{\mathrm{pp}}$, the first factor (1) is substituted by L_pp according to equation (3.3). In principle, any L_pp function could be put in here. The same holds true for the cosine term in the middle (2), but we here opted to compute it from the geomagnetic field vector’s components via

$$\begin{array}{ccc}{\phi }_m=\arctan \left(\frac{-0.75·B_z+\sqrt{0.5625·{B_z}^2+0.5·\left({B_x}^2+{B_y}^2\right)}}{\sqrt{{B_x}^2+{B_y}^2}}\right)& \mathrm{with}& \vec{B}=\left\{\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right\}\end{array},$$(3.7)

an ansatz that is developed in detail in Annex A. Principally, a $\vec{B}$ vector from any dipole and non-dipole field could be used here, for instance, computed with the International Geomagnetic Reference Field (IGRF, https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html) at the position indicated by ${\vec{R}}_{\mathrm{pp}}$ and related epoch. However, the $\vec{B}$ vector components need to be converted into the SM system before entering into equation (3.7) since the L_pp-function (3.3) is referred to SM (see Annex B).

In order to fit the ellipse parameters of the L_pp-function, equation (3.3), from plasmapause position vectors, per observed ${\vec{R}}_{\mathrm{pp}}$, three component-wise observation equations are set up from equation (3.6)

$$\begin{array}{ccc}{∆R}_x+v_x& =& \frac{\partial L_{\mathrm{pp}}}{\partial a}{|}_0·∆a·{\cos }^2{\phi }_m·{\epsilon }_x+\frac{\partial L_{\mathrm{pp}}}{\partial e}{|}_0·∆e·{\cos }^2{\phi }_m·{\epsilon }_x+\frac{\partial L_{\mathrm{pp}}}{\partial \varpi }{|}_0·∆\varpi ·{\cos }^2{\phi }_m·{\epsilon }_x,\\ ∆R_y+v_y& =& \frac{\partial L_{\mathrm{pp}}}{\partial a}{|}_0·∆a·{\cos }^2{\phi }_m·{\epsilon }_y+\frac{\partial L_{\mathrm{pp}}}{\partial e}{|}_0·∆e·{\cos }^2{\phi }_m·{\epsilon }_y+\frac{\partial L_{\mathrm{pp}}}{\partial \varpi }{|}_0·∆\varpi ·{\cos }^2{\phi }_m·{\epsilon }_y,\\ ∆R_z+v_z& =& \frac{\partial L_{\mathrm{pp}}}{\partial a}{|}_0·∆a·{\cos }^2{\phi }_m·{\epsilon }_z+\frac{\partial L_{\mathrm{pp}}}{\partial e}{|}_0·∆e·{\cos }^2{\phi }_m·{\epsilon }_z+\frac{\partial L_{\mathrm{pp}}}{\partial \varpi }{|}_0·∆\varpi ·{\cos }^2{\phi }_m·{\epsilon }_z .\end{array}$$(3.8)

Expressed in matrix form, a set of 3n observation equations can be established for n observed ${\vec{R}}_{\mathrm{pp}}$ vectors based on equations (3.8),

$$\begin{array}{ccc}\left[\begin{array}{c}\begin{array}{c}{∆R}_{x_1}\\ {∆R}_{y_1}\\ {∆R}_{z_1}\end{array}\\ ⋮\\ \begin{array}{c}{∆R}_{x_n}\\ {∆R}_{y_n}\\ {∆R}_{z_n}\end{array}\end{array}\right]& +& \begin{array}{ccc}\left[\begin{array}{c}\begin{array}{c}{v}_{x_1}\\ v_{y_1}\\ v_{z_1}\end{array}\\ ⋮\\ \begin{array}{c}{v}_{x_n}\\ v_{y_n}\\ v_{z_n}\end{array}\end{array}\right]& =& \begin{array}{ccc}\left[\begin{array}{ccc}\begin{array}{c}\frac{\partial L_{{\mathrm{pp}}_1}}{\partial a}{|}_0·{\cos }^2{\phi }_1·{\epsilon }_{x_1}\\ \frac{\partial L_{{\mathrm{pp}}_1}}{\partial a}{|}_0·{\cos }^2{\phi }_1·{\epsilon }_{y_1}\\ \frac{\partial L_{{\mathrm{pp}}_1}}{\partial a}{|}_0·{\cos }^2{\phi }_1·{\epsilon }_{z_1}\end{array}& \begin{array}{c}\frac{\partial L_{{\mathrm{pp}}_1}}{\partial e}{|}_0·{\cos }^2{\phi }_1·{\epsilon }_{x_1}\\ \frac{\partial L_{{\mathrm{pp}}_1}}{\partial e}{|}_0·{\cos }^2{\phi }_1·{\epsilon }_{y_1}\\ \frac{\partial L_{{\mathrm{pp}}_1}}{\partial e}{|}_0·{\cos }^2{\phi }_1·{\epsilon }_{z_1}\end{array}& \begin{array}{c}\frac{\partial L_{{\mathrm{pp}}_1}}{\partial \varpi }{|}_0·{\cos }^2{\phi }_1·{\epsilon }_{x_1}\\ \frac{\partial L_{{\mathrm{pp}}_1}}{\partial \varpi }{|}_0·{\cos }^2{\phi }_1·{\epsilon }_{y_1}\\ \frac{\partial L_{{\mathrm{pp}}_1}}{\partial \varpi }{|}_0·{\cos }^2{\phi }_1·{\epsilon }_{z_1}\end{array}\\ ⋮& ⋮& ⋮\\ \begin{array}{c}\frac{\partial L_{{\mathrm{pp}}_n}}{\partial a}{|}_0·{\cos }^2{\phi }_n·{\epsilon }_{x_n}\\ \frac{\partial L_{{\mathrm{pp}}_n}}{\partial a}{|}_0·{\cos }^2{\phi }_n·{\epsilon }_{y_n}\\ \frac{\partial L_{{\mathrm{pp}}_n}}{\partial a}{|}_0·{\cos }^2{\phi }_n·{\epsilon }_{z_n}\end{array}& \begin{array}{c}\frac{\partial L_{{\mathrm{pp}}_n}}{\partial e}{|}_0·{\cos }^2{\phi }_n·{\epsilon }_{x_n}\\ \frac{\partial L_{{\mathrm{pp}}_n}}{\partial e}{|}_0·{\cos }^2{\phi }_n·{\epsilon }_{y_n}\\ \frac{\partial L_{{\mathrm{pp}}_n}}{\partial e}{|}_0·{\cos }^2{\phi }_n·{\epsilon }_{z_n}\end{array}& \begin{array}{c}\frac{\partial L_{{\mathrm{pp}}_n}}{\partial \varpi }{|}_0·{\cos }^2{\phi }_n·{\epsilon }_{x_n}\\ \frac{\partial L_{{\mathrm{pp}}_n}}{\partial \varpi }{|}_0·{\cos }^2{\phi }_n·{\epsilon }_{y_n}\\ \frac{\partial L_{{\mathrm{pp}}_n}}{\partial \varpi }{|}_0·{\cos }^2{\phi }_n·{\epsilon }_{z_n}\end{array}\end{array}\right]& *& \left[\begin{array}{c}∆a\\ ∆e\\ ∆\varpi \end{array}\right]\end{array}\end{array}\end{array},$$(3.9)

where for the non-linear fit

• $∆{\vec{R}}_{\mathrm{pp}}={\vec{R}}_{\mathrm{ppObs}}-{\vec{R}}_{\mathrm{ppMod}}={\vec{R}}_{\mathrm{ppObs}}-L_{\mathrm{ppMod}}·{\cos }^2{\phi }_m·{\vec{\epsilon }}_{\mathrm{pp}}$,

• ${\vec{R}}_{\mathrm{ppObs}}$ … plasmapause position vector detected from the IMAGE & Van Allen Probes data, Sect. 2.3,

• L_ppMod … L_pp-value computed with equation (3.3) using some initial values for a, e, ϖ,

• cos²φ_m … computed from $\vec{B}$ vector, equation (3.7),

• ${\vec{\epsilon }}_{\mathrm{pp}}$ … unit vector of ${\vec{R}}_{\mathrm{ppObs}}$,

• v … fitted correction to the observable,

• $\frac{\partial L_{\mathrm{pp}}}{\partial a}{|}_0$ … partial of L_pp w.r.t. semi-major axis a, evaluated with initial values,

• ∆a … estimated correction to initial value of a,

• $\frac{\partial L_{\mathrm{pp}}}{\partial e}{|}_0$ … partial of L_pp w.r.t. eccentricity e, evaluated with initial values,

• ∆e … estimated correction to initial value of e,

• $\frac{\partial L_{\mathrm{pp}}}{\partial \varpi }{|}_0$ … partial of L_pp w.r.t. orientation angle ϖ, evaluated with initial values,

• ∆ϖ … estimated correction to initial value of ϖ.

A 3D plasmapause position computation with the NEPPM can be summarized as follows. The user has to provide the geographic latitude φ and longitude λ, epoch T and Dst, for which a 3D plasmapause position ${\vec{R}}_{\mathrm{pp}}$ is requested. Then the computation is conducted in the following steps (see also Annex B):

1. Set transformation parameters geographic ➧ SM.

2. Compute geographic position unit vector ${\vec{\epsilon }}_g=\left\{\begin{array}{c}\cos \phi ·\cos \lambda \\ \cos \phi ·\sin \lambda \\ \sin \phi \end{array}\right\}$.

3. Transform ${\vec{\epsilon }}_g$ from geographic into SM ➧ unit vector ${\vec{\epsilon }}_{\mathrm{pp}}$.

4. Compute $\mathrm{MLT} = \Lambda ·\frac{12 \mathrm{h}}{\pi }+12 \mathrm{h}$, with $\Lambda =\arctan \left(\frac{{\epsilon }_{{\mathrm{pp}}_y}}{{\epsilon }_{{\mathrm{pp}}_x}}\right)$.

5. Compute cos²φ_m term for dipole from the components of ${\vec{\epsilon }}_{\mathrm{pp}}$ via

$${\cos }^2{\phi }_m=\frac{{{\epsilon }_{{\mathrm{pp}}_x}}^2+{{\epsilon }_{{\mathrm{pp}}_y}}^2}{{{\epsilon }_{{\mathrm{pp}}_x}}^2+{{\epsilon }_{{\mathrm{pp}}_y}}^2+{{\epsilon }_{{\mathrm{pp}}_z}}^2}.$$

1. The fitted L_pp function is evaluated with equation (3.3).

2. The dipole position ${\vec{R}}_{\mathrm{pp}}$ is computed with equation (3.6), where term (3) is here ${\vec{\epsilon }}_{\mathrm{pp}}$.

3. With the IGRF (or another non-dipole model) a $\vec{B}$ vector is computed at the position indicated by dipole position vector ${\vec{R}}_{\mathrm{pp}}$ and for the requested epoch T, transform $\vec{B}$ into the SM system.

4. Compute from the non-dipole $\vec{B}$ vector components an improved cos²φ_m term, equation (3.7).

5. Compute non-dipole ${\vec{R}}_{\mathrm{pp}}$ vector, again by using equation (3.6), but now with the non-dipole cos²φ_m term obtained from Step 9).

6. Transform non-dipole ${\vec{R}}_{\mathrm{pp}}$ vector from SM ➧ geographic.

3.3 Fitting of empirical functions to ellipse parameters

Per Dst level listed in Table 1, a least squares fit was run with the L_pp-function, equation (3.3), in combination with the 3D modelling scheme, equation (3.9), to get one set of ellipse parameters (a,e,ϖ) for that Dst level. In order to obtain representations of each ellipse parameter in dependency of Dst, empirical exponential functions f(Dst) were fitted to the ellipse parameters

$$f\left(\mathrm{Dst}\right)=A·\exp \left(\frac{D\mathrm{stmed}}{\mathrm{Dscal}}\right)+C,$$(3.10)

where

• Dstmed … median Dst of considered Dst level represented by a certain Dst range, Table 1,

• Dscal … Dst scaling factor,

• A … amplitude to be fitted,

• C … boundary value fixed for extreme storms by assuming Dstmed → −∞.

It turned out, that exponential functions are best suited to represent the ellipse parameters in dependency of Dst. This approach has the advantage that besides the fitting using the available data set, extreme storms that have never been observed can be considered in a plausible way. It is evident that the boundary values fixed for all three fits are not approved by observations. We have defined these values following plausible indications of observations and physical arguments. Consequently, the given boundary values might be slightly modified using an extended database in future studies. Because such a modification has a negligible impact on the model approach for weak and moderate storms (i.e., Dst > −100 nT) the model is applicable in common operational services. Besides Dscal and strength A obtained by least squares fitting, also the estimated boundary values C are listed in Table 2. For the semi-major axis, the eccentricity, and the orientation angle of the ellipse.

While fitting the strength A of the f(Dst) function, the input Dst values had to be scaled by an empirical Dscal quotient. This procedure converted the range of Dst values into a reasonable argument range for the exponential function. Similar to the boundary values for extreme storms, C_a = 1.7 R_E, C_e = 0.27, C_M = 12 h, Table 2, also boundary values for Dst = 0 had to be defined, which were 5 R_E, 0, 24 h for semi-major axis, eccentricity and MLT of bulge, respectively. Then, with the aid of Dscal, the original Dst scale had to be adapted to an argument scale, such that the exponential function was able to cover the required bandwidths appropriately. Only the blue points in Figure 3 could enter into the fit, and their distribution is so linear, that Dscal had to be fixed manually. An exponential fitting of Dscal together with the amplitude turned out not feasible. Thus, only A could be determined by fitting.

Figure 3 Semi-major axis (a), eccentricity (b), and orientation angle (c) of the Lpp ellipse as a function of Dst (blue lines). The blue dots denote the fitted ellipse parameters for each Dst level. Additional aurora observations for (a) are shown in magenta.

Figure 3a shows the NEPPM representation for the semi-major axis a. The zoom indicates the range where data are available from the data set analyzed in this study. Unfortunately, only a few data for strong and extreme storms (i.e., Dst < −100 nT) are available, therefore, such individual cases are not statistically significant. Moreover, the aurora observations (1)–(3) cited below suffer from the fact that it is often not clear, under which aspect angle these auroras have been seen. So, the points displayed in magenta in Figure 3a for the semi-major axis a are best effort values, also considering the ellipse eccentricity estimates (further explanations see below). These three reference points (1)–(3) indicated in Figure 3a are described below to illustrate the estimation procedure and our exponential fit.

1. Cliver & Svalgaard (2004) have considered largest geomagnetic storms based on Dst between 1932 and 2002. The peak value of the Dst-index observed at the extreme storm on 13/14 March 1989 is given by Dst = −584 nT. Polar lights have been observed at 29°N geomagnetic latitude that would correlate with $L_{{\mathrm{pp}}_{\min }}\approx 1.4$ and a ≈ 1.9, respectively.

2. For the great space weather event in February 1730 Hayakawa et al. (2019) estimated a Dst value in the order of Dst_min ≈ −1200 nT or even less. Auroras have been observed down to 27° magnetic latitude in the Northern sky meaning that the location of the aurora occurred at ionospheric heights in the range of 29°–32°N magnetic latitude, associated with a ≈ 1.85.

3. Tsurutani et al. (2003) have estimated the Dst peak value as Dst ≈ −1760 nT for the Carrington Event 1859 that is related to a plasmapause position at around $L_{{\mathrm{pp}}_{\min }}=1.3$, based on aurora locations analyzed by Kimball (1960). In relation to the eccentricity of the ellipse approach at great storms (see Fig. 3b) we estimate a value of about 1.7 for the semi-major axis, which is considered the boundary value for describing extreme storms by NEPPM. Although there remain uncertainties, it is believed that NEPPM represents aurora observations quite well when using this boundary value.

Besides the semi-major axis a also the eccentricity e is a key parameter of the ellipse. It is connected with L_pp by the relations

$$L_{{\mathrm{pp}}_{\min }}=a·(1-e)$$(3.11)

and

$$L_{{\mathrm{pp}}_{\max }}=a·\left(1+e\right).$$(3.12)

Thus, the difference between both extreme plasmapause positions is given by

$$L_{{\mathrm{pp}}_{\max }}-L_{{\mathrm{pp}}_{\min }}=2a·e.$$(3.13)

Considering the available data set, the boundary value for the eccentricity was fixed at e = 0.27 for the fitting procedure. Higher values would considerably increase the difference between $L_{{\mathrm{pp}}_{\max }}$ and $L_{{\mathrm{pp}}_{\min }}$ towards unrealistic values. The estimate of the semi-major axis, as shown in Figure 3a, is closely related to equation (3.11). To complete the model approach, an expression for the direction of the model ellipse is still needed (see Fig. 3c). It is well-known that a plasma bulge is created in the evening hours due to the superposition of the co-rotation electric field and the solar wind-controlled convection electric field. The enhanced plasma density leads to an outward motion of the plasmapause characterized by $L_{{\mathrm{pp}}_{\max }}$. Considering the poor database for extreme storm events, the definition of a boundary value remains speculative to a certain extent. It is evident that the boundary value must be bigger than MLT = 12 h. Here, the boundary value is set as C_M = 12 h. As with the other boundary values C_a and C_e too, the fixed values can be modified if it is required by results obtained from an extended database.

4 Results and comparisons

The NEPPM clearly demonstrates the torus-shaped plasmapause, as can be seen in Figure 4, which divides the dense co-rotating plasma driven by the electric field and the outer tenuous magnetospheric plasma. Figure 4 illustrates the overall compression of the plasmapause during storm activities, indicated by strongly negative Dst values (red) and an expansion during quiet periods (blue). Moreover, we can clearly see the plasmapause bulge during storm conditions appearing at noon (in the direction of the positive x-axis), indicated by the increased distance from the Earth for the inner red curve. During quiet conditions (outer blue curve), the torus looks equally balanced in distance from the Earth.

Figure 4 Results of the modelled Lpp values obtained by NEPPM in SM coordinates depending on the Dst index varying from quiet (blue) to disturbed (red) conditions: (a) for hourly MLT between 12 h and 24 h and (b) in the noon-midnight meridional plane (i.e. y = 0).

For direct comparison, we use the Global Core Plasma Model (GCPM) by Gallagher et al. (2000). The GCPM is a composition of different region-specific models for the core plasma density with smooth transitions in value and derivative. Intended to match the overall physical appearance, certain specifics of some regions may not be included in high detail. Allowing the estimation of electron densities in a larger region throughout the inner magnetosphere, it considers the space weather impact by relying on the geomagnetic Kp index. In contrast, the NEPPM depends on the Dst index, which has been shown to have better correlations with L_pp than Kp (cf. O’Brien & Moldwin, 2003). Nevertheless, for direct comparisons, the subfigures of Figure 5 are drawn on Dst scale, as we are simply running the GCPM for the same dataset with the associated Kp values but plotting the corresponding Dst values of the given data by IMAGE, and Van Allen Probes.

Figure 5 Results of the NEPPM (left) and GCPM (right) shown for all Dst levels (lines, coloured by the six corresponding Dst levels: [−10, −20], [−20, −30], [−30, −40], [−40, −50], [−50, −70], [−50, −500]), with the extracted Lpp values coloured by Dst (dots) in the equatorial plane (top) and via Lpp over MLT (middle), and its RMS over MLT (bottom).

Figure 5 shows all extracted L_pp values from IMAGE and Van Allen Probes as well as the modelled L_pp values for NEPPM (left) and GCPM (right), in the equatorial plane (top) as well as by L_pp over MLT (middle), and its corresponding RMS (bottom). The L_pp values are generally located between L = 3 and L = 6, and are evenly distributed in MLT. Both models show a plasmapause bulge in the afternoon sector. During high solar activity, the bulge moves towards noon (i.e. westward towards the Sun), and during quiet conditions towards midnight.

However, the size of the bulge is visibly different in the two models. The significant bulge in GCPM is striking, which is considerably less pronounced in NEPPM. As we can see in the bottom panels of Figure 5, the deviation of NEPPM is generally lower than the one of GCPM. While NEPPM shows a mean RMS of about 0.5 R_E, GCPM shows a higher mean RMS of 0.9 R_E. The large overestimation of the bulge by the GCPM is again shown in the afternoon-evening hours (cf. Fig. 5, bottom right).

An additional illustration of extracted and modelled L_pp values per Dst level reveals more information. Figure 6 shows for each Dst level, the extracted L_pp values next to the resulting mean L_pp computed with both the NEPPM (black) and the GCPM (grey). At all Dst levels, GCPM reveals a consistent overestimation of the extracted L_pp values. This overestimation expands even further with strongly negative Dst levels. The bulge in GCPM is generally higher than the extracted L_pp values indicate. Overall, the bulge, which is given by the ellipse in NEPPM, is less pronounced and gets smaller with increasing Dst index. Moreover, for Dst > −10 nT, the ellipse in NEPPM becomes almost circular in shape. Unfortunately, during quiet periods, there is a wide scattering of data at all MLT values making a good fitting of the model difficult, i.e. for Dst > −30 nT.

Figure 6 NEPPM (black) and GCPM (grey) for each Dst level, cf. Table 1. The extracted Lpp values are marked by dots in red (IMAGE), and blue (RBSP).

Figure 7 displays the specific deviations of the two models from the IMAGE and Van Allen Probes data with respect to Dst and MLT. NEPPM displays the highest RMS in the afternoon for Dst < −50 nT. GCPM shows two higher RMS values, one in the afternoon for very low Dst values and another one in the evening hours for Dst > −40 nT. Furthermore, the overall dayside shows enhanced RMS values for GCPM, which are not seen for NEPPM. Since we used IMAGE and Van Allen Probes data to build the NEPPM, our model performs well. Certainly, further data is necessary to validate our approach. As we can see, GCPM does not represent the recent data properly, though it was a good and, most importantly, consistent global approach back then.

Figure 7 RMS (top) and relative RMS (bottom) for NEPPM (left) and GCPM (right).

Meanwhile, much more detailed data has become available for improved accuracy and higher resolution. Various plasmapause models were proposed, and especially neural networks (NN) are becoming of higher interest as proposed by Guo et al. (2021). A comparison of the RMS with this NN model shows similar results. Both NEPPM and NN models show the highest RMS values of about 0.9 R_E around the bulge at 19 MLT and the lowest RMS values with 0.4 R_E between 2 and 3 MLT. However, in addition to small fluctuations, Guo et al. (2021), Figure 7a, exhibit a peak in RMS at 24 MLT. This behaviour is similar to other models (Moldwin et al., 2002; O’Brien & Moldwin, 2003) also shown there. In comparison to this, the NEPPM shows a broader but consistent deviation in RMS from 16 to 20 MLT with only minor fluctuations. Furthermore, the small RMS at 24 MLT there unveils a better suitability of the NEPPM compared to others.

Also note that we use different L_pp data sources for the RMS calculations. So, the pre-selection of our L_pp values is crucial for a good model with a small RMS. Due to the different conditions for plasmapause detection introduced in this paper, we are using a reliable L_ppdatabase without any multiple plasmapause crossings within a profile.

In some cases, still under investigation, a plasma plume may separate from the co-rotating plasma, thus forming a separate plasma structure at even greater radial distances. The current version of the NEPPM does not cover such structures. The NEPPM describes the average plasmapause behaviour related to the fundamental bulge formation process as well as possible in a robust way. This shortcoming is acceptable because the NEPPM has been developed to support, in particular, the current plasmasphere model Neustrelitz PlasmaSphere Model (NPSM) developed at DLR for operational space weather services (Jakowski & Hoque, 2018). Plasma dynamics creating plasma plume structures make it difficult to extract a clear orientation of the ellipse from the data, as can be seen in Figure 6 (L_pp outliers between about 5 and 6 in the afternoon sector) and Figure 5 (enhance RMS in the afternoon sector). Nevertheless, it is evident that with increasing storm intensity, the major ellipse axis turns towards noon in the equatorial plane. This indicates the growing impact of the convection electric field that competes with the co-rotating electric field.

5 Summary and conclusions

In this paper, we present a new approach for modelling the plasmapause position. The database is formed by two comprehensive datasets, the electron densities derived from the passive IMAGE RPI data from 2001 to 2005 and the electron densities provided by the NURD data of the Van Allen Probes between 2012 and 2018. To determine the plasmapause position L_pp, detecting a sharp gradient of the electron density is essential. Naturally, some measurements do not fulfill the requirements for reliable L_pp estimations e.g. if a profile is too short or too noisy. Moreover, different data sources offer different resolutions (in time and space), so uneven grids have to be compared reasonably. Due to the amount of data, an automated detection algorithm has been established considering these uncertainties. Unsuitable parts of the data are filtered out in accordance with well-founded constraints that guarantee the identification of a clear and pronounced plasmapause position. By this, over 3000 plasmapause crossings have been found, which formed the basis for creating our own model.

The presented NEPPM is an ellipse-based approach to model the plasmapause location as a function of the geomagnetic Dst index and the Magnetic Local Time MLT. In this approach, one focal point coincides with the location of the Earth, and the semi-major axis, eccentricity, and orientation of the ellipse are least-squares fitted to the obtained plasmapause locations. By relating an observed 3D plasmapause position vector ${\vec{R}}_{\mathrm{pp}}$ to a L_pp-function, the possibility to account for a non-dipole geomagnetic field can be established in a very easy way. In addition, NEPPM allows the user to directly compute 3D positions on the plasmapause torus for a given latitude, longitude, epoch, and Dst. Due to its short computation times, the NEPPM is an ideal background model in operational space weather services. In the first step, the NEPPM will be integrated into the NPSM so that their currently fixed plasmapause position of L_pp = 5 will be replaced with our dynamic and robust empirical approach. In a second step, the NEPPM will be combined with the three-dimensional Neustrelitz Electron Density Model (NEDM, Hoque et al., 2022), covering the entire ionosphere/plasmasphere systems on a global scale under all space weather conditions. The operational use of such an empirical model requires a synthesis of simplicity and robustness in handling, sufficient accuracy, and fast computation. Considering this, the plasmapause model NEPPM does not claim to provide a comprehensive description of physical processes contributing to the plasmapause. However, the model should agree with current physical knowledge of the plasmapause behaviour and provide an accuracy that is basically comparable with those obtained by state-of-the-art models. The NEPPM is fed by the geomagnetic index Dst and is easily accessible for operational applications via https://wdc.kugi.kyoto-u.ac.jp/dstdir/ (Nose et al., 2015). Furthermore, Dst can be predicted up to 24 h in advance, as described by Park et al. (2021). The predictability and the hourly resolution of Dst fit quite well with current needs in operational services.

Nevertheless, we still lack data for strong and extreme storm events. Therefore, further L_pp observations are required to validate the proposed empirical model, especially during extreme storm conditions. A direct comparison with the GCPM illustrates the closer representation in our model. Moreover, comparing the RMS with the NN model shows similar errors, both ranging between 0.4 and 0.9 R_E. So, our simple approach can compete with such machine learning models and is also based on physical constraints such as limits for the dimension via the semi-major axis a and the bulge location via M (cf. Table 2). Finally, it should be noted that all models presented suffer the largest errors in the afternoon/evening sector between 16 and 20 MLT. This can be attributed to “detached” plasma plumes, commonly observed in this MLT sector (Darrouzet et al., 2008), a phenomenon that still needs further investigation. An extended analysis with additional data is foreseen to capture the cause and development of such structures. Regardless, it is not critical for radio system applications because plasma plume densities do not contribute noticeably to TEC. Hence, the NEPPM serves the aim of supporting operational space weather services.

Acknowledgments

The authors thank R.E. Denton for providing the passive electron density data of IMAGE RPI. We thank Yuri Shprits and Matyas Szabo-Roberts for their support and helpful discussions on the NURD data of the Van Allen Probes. Furthermore, we thank NASA for making Van Allen Probes data available. We would like to thank Kateryna Lubyk for her valuable remarks during her revision. Finally, we are grateful to the reviewers and editors of JSWSC for their valuable comments and their substantial support to improve the paper. The editor thanks Fabien Darrouzet and Scott Thaller for their assistance in evaluating this paper.

Annex A

Plasmapause positions are described in the solar-magnetic (SM) system in terms of Cartesian plasmapause position vectors ${\vec{R}}_{\mathrm{pp}}$ that obey the representation equation (3.6). When computing them with a constant L_pp value in an equidistant grid, a plot of these ${\vec{R}}_{\mathrm{pp}}$ would display a torus-like picture, blue in Figure A1 below, corresponding to the dipole assumption inherent in the definition of L-shells.

Figure A1 The vector fields ${\vec{R}}_{\mathrm{pp}}$ (blue) and ${\vec{T}}_{\mathrm{pp}}$ (red) in a 2D meridional cut (left) and as a 3D illustration (right).

In case the plasmapause modelling shall be conducted in a dipole geomagnetic field, the projection term can be directly computed from the position vector components as

$$\begin{array}{ccc}{\cos }^2{\phi }_m=\frac{{R_x}^2+{R_y}^2}{{R_x}^2+{R_y}^2+{R_z}^2}& \mathrm{with}& {\vec{R}}_{\mathrm{pp}}=\left\{\begin{array}{c}{R}_x\\ R_y\\ R_z\end{array}\right\}\end{array},$$(A.1)

whereas it must be linked to the geomagnetic field vector $\vec{B}$ if, like in the IGRF, a non-dipole field is the base. For this, $\vec{B}$ has to be computed with the selected geomagnetic field model at the position indicated by ${\vec{R}}_{\mathrm{pp}}$ and the related epoch, and it must be transformed into the SM system to be consistent with the plasmapause modelling. This annex aims to give the foundation for the resulting equation previously presented in equation (3.7).

At first sight, the formulae for calculating cos²φ_m from $\vec{B}$ are once more derived from ${\vec{R}}_{\mathrm{pp}}$, i.e. for a dipole. Nevertheless, as will be shown below, the components from a non-dipole $\vec{B}$ vector can enter these, too. The procedure may somehow be comparable with the principle of osculating Keplerian elements in orbit modelling, where at each orbital point another slightly varying set of Keplerian elements is valid. Here, this might be translated into a kind of osculating dipole, always being represented by a dipole magnetic field line put through the direction of the actually used non-dipole $\vec{B}$ vector at that point in the magnetic field.

First, reformulate equation (3.6) to

$${\vec{R}}_{\mathrm{pp}}=L_{\mathrm{pp}}·\left\{\begin{array}{c}{\cos }^3{\phi }_m·\cos \Lambda \\ {\cos }^3{\phi }_m·\sin \Lambda \\ {\cos }^2{\phi }_m·\sin {\phi }_m\end{array}\right\}$$(A.2)

and determine the first derivative of ${\vec{R}}_{\mathrm{pp}}$ w.r.t. the geomagnetic latitude. By this, one has the tangent vector

$${\vec{T}}_{\mathrm{pp}}=\frac{\partial {\vec{R}}_{\mathrm{pp}}}{\partial {\phi }_m}=L_{\mathrm{pp}}·\left\{\begin{array}{c}-3·{\cos }^2{\phi }_m·\sin {\phi }_m·\cos \Lambda \\ -3·{\cos }^2{\phi }_m·\sin {\phi }_m·\sin \Lambda \\ -2·\cos {\phi }_m·{\sin }^2{\phi }_m+{\cos }^3{\phi }_m\end{array}\right\},$$(A.3)

which always points tangentially northwards along the dipole field line, i.e., into the same direction as $\vec{B}$. The geometry behind ${\vec{R}}_{\mathrm{pp}}$ and ${\vec{T}}_{\mathrm{pp}}$ is visualized in Figure A1. The left image shows a cut through the dipole torus (or one dipole loop) with both vector types in their correct lengths. In the right-hand plot, the torus in space is shown, but with tangents drawn as unit vectors, as otherwise the picture would appear too distorted.

Since always ${\vec{T}}_{\mathrm{pp}}\uparrow \uparrow \vec{B}$, analogously to spherical latitude computation from Cartesian coordinates, one can build the ratio

$$\begin{array}{ccc}\frac{T_z}{\sqrt{{T_x}^2+{T_y}^2}}=\frac{B_z}{\sqrt{{B_x}^2+{B_y}^2}}& \mathrm{with}& {\vec{T}}_{\mathrm{pp}}=\left\{\begin{array}{c}{T}_x\\ T_y\\ T_z\end{array}\right\}, \vec{B}=\left\{\begin{array}{c}{B}_x\\ B_y\\ B_z\end{array}\right\}.\end{array}$$(A.4)

In this way, the different magnitudes of both vectors are neutralized. Next, the dependence on longitude Λ will be eliminated. Substituting the components as given in (A.3) and applying the trigonometric Pythagoras law cos²Λ + sin²Λ = 1, see e.g. Sigl (1977), the denominator of the left-hand side in (A.4) can be condensed to

$$\sqrt{{T_x}^2+{T_y}^2}=3·{\cos }^2{\phi }_m·\sin {\phi }_m.$$(A.5)

Thus, linking the dipole $\vec{B}$ and φ_m, we get

$$\frac{B_z}{\sqrt{{B_x}^2+{B_y}^2}}=-\frac{2·\cos {\phi }_m·{\sin }^2{\phi }_m}{3·{\cos }^2{\phi }_m·\sin {\phi }_m}+\frac{{\cos }^3{\phi }_m}{3·{\cos }^2{\phi }_m·\sin {\phi }_m}=-\frac{2}{3}·\tan {\phi }_m+\frac{1}{3}·\frac{1}{\tan {\phi }_m},$$(A.6)

where the definition tanφ_m = sinφ_m/cosφ_m, again Sigl (1977), supports further simplification.

Multiplying the outermost identity of equation (A.6) with tanφ_m and bringing all terms on one side gives a quadratic equation in x = tanφ_m

$${\left(\tan {\phi }_m\right)}^2+\underset{p}{\underbrace{\frac{3}{2}·\frac{B_z}{\sqrt{{B_x}^2+{B_y}^2}}}}·\left(\tan {\phi }_m\right)\underset{q}{\underbrace{-\frac{1}{2}}}=0.$$(A.7)

This can be solved according to the p-q-formula with

$$x_{\mathrm{1,2}}=-\frac{p}{2}\pm \sqrt{{\left(\frac{p}{2}\right)}^2-q.}$$(A.8)

Leading to

$${\left(\tan {\phi }_m\right)}_{\mathrm{1,2}}=-\frac{3}{4}·\frac{B_z}{\sqrt{{B_x}^2+{B_y}^2}}\pm \sqrt{\frac{9}{16}·\frac{{B_z}^2}{{B_x}^2+{B_y}^2}+\frac{1}{2}}.$$(A.9)

As will be demonstrated in Annex A1 below, the positive solution of equation (A.9) is always the correct one. Thus equation (A.9) gives rise to equation (3.7), which is optimized for being processed with the atan2(x, y) function available in many programming languages.

For B_z = 0 the dipole loop reaches its maximal/minimal elevation above/below the geomagnetic equator, i.e. ${\left(\tan {\phi }_m\right)}_{\mathrm{1,2}}=\pm 1\sqrt{2}$, corresponding to φ_m,max/min = ±35°. 264389682754654.

In analogy to equations (A.6), (A.7), a quadratic equation can also be established in cotφ_m.

A1 Rigorous proof that the positive square root is the one sought-after

For this proof, the components of the tangent vector ${\vec{T}}_{\mathrm{pp}}$ have to be expressed as a function of the components {R_x R_y R_z} of the plasmapause position vector ${\vec{R}}_{\mathrm{pp}}$. Therefore in equation (A.3) the cosine and sine terms are formulated as functions of the plasmapause position vector components

$$\cos {\phi }_m=\frac{\sqrt{{R_x}^2+{R_y}^2}}{\sqrt{{R_x}^2+{R_y}^2+{R_z}^2}}, \sin {\phi }_m=\frac{R_z}{\sqrt{{R_x}^2+{R_y}^2+{R_z}^2}}, \cos \Lambda =\frac{R_x}{\sqrt{{R_x}^2+{R_y}^2}}, \sin \Lambda =\frac{R_y}{\sqrt{{R_x}^2+{R_y}^2}}.$$(A.10)

After having substituted equations (A.10) into equation (A.3), one obtains after some simple algebra

$${\vec{T}}_{\mathrm{pp}}=L_{\mathrm{pp}}·\frac{\sqrt{{R_x}^2+{R_y}^2}}{{\left(\sqrt{{R_x}^2+{R_y}^2+{R_z}^2}\right)}^3}·\left\{\begin{array}{c}-3·R_z·R_x\\ -3·R_z·R_y\\ -2·{R_z}^2+{R_x}^2+{R_y}^2\end{array}\right\}.$$(A.11)

Whereby for describing purely the direction of ${\vec{T}}_{\mathrm{pp}}$ the braced term {⋮} of equation (A.11) is sufficient.

With the expression of the tangent vector components in terms of Cartesian components of the plasmapause position vector, equation (A.11), it is now possible to conduct in a rigorous mathematical proof, which of the two square roots of equation (A.9) provides the sought-after tan φ_m. A plasmapause position vector ${\vec{R}}_{\mathrm{pp}}$ is computed for a given geomagnetic latitude φ_m and solar longitude Λ according to equation (A.2). Reversely, φ_m can be computed from the ${\vec{R}}_{\mathrm{pp}}$ Cartesian components in a backwards transformation

$$\tan {\phi }_m=\frac{R_z}{\sqrt{{R_x}^2+{R_y}^2}}.$$(A.12)

The proof starts with equation (A.9), in which, due to equation (A.4), the $\vec{B}$ vector components can be replaced by the tangent vector ${\vec{T}}_{\mathrm{pp}}$ components. This approach will, after some algebra, allow to express equation (A.9) in terms of the ${\vec{R}}_{\mathrm{pp}}$ vector Cartesian components and finally lead to tan φ_m formulae giving both φ_m solutions of the quadratic equation as a function of {R_x R_y R_z}. Thus equation (A.9) is now expressed in terms of the ${\vec{T}}_{\mathrm{pp}}$ components {T_x T_y T_z}

$${\left(\tan {\phi }_m\right)}_{\mathrm{1,2}}=-\frac{3}{4}·\frac{T_z}{\sqrt{{T_x}^2+{T_y}^2}}\pm \sqrt{\frac{9}{16}·\frac{{T_z}^2}{{T_x}^2+{T_y}^2}+\frac{1}{2}}.$$(A.13)

In an intermediate step, the square root term of equation (A.13) is reformulated

$$\pm \sqrt{\frac{9}{16}·\frac{{T_z}^2}{{T_x}^2+{T_y}^2}+\frac{1}{2}}=\pm \sqrt{\frac{9}{16}·\frac{{T_z}^2}{{T_x}^2+{T_y}^2}+\frac{8}{16}·\frac{{T_x}^2+{T_y}^2}{{T_x}^2+{T_y}^2}}=\pm \frac{\sqrt{9·{T_z}^2+8·{T_x}^2+8{{·T}_y}^2}}{4·\sqrt{{T_x}^2+{T_y}^2}}.$$(A.14)

Substituting equation (A.14) into (A.13)

$${\left(\tan {\phi }_m\right)}_{\mathrm{1,2}}=-\frac{3}{4}·\frac{T_z}{\sqrt{{T_x}^2+{T_y}^2}}\pm \frac{\sqrt{9·{T_z}^2+8·{T_x}^2+8{{·T}_y}^2}}{4·\sqrt{{T_x}^2+{T_y}^2}}=\frac{-3·T_z\pm \sqrt{9·{T_z}^2+8·{T_x}^2+8{{·T}_y}^2}}{4·\sqrt{{T_x}^2+{T_y}^2}}.$$(A.15)

Now, the expressions of {T_x T_y T_z} in terms of {R_x R_y R_z}, equation (A.11), must be substituted into equation (A.15). Thereby the term $L_{\mathrm{pp}}·\sqrt{{R_x}^2+{R_y}^2}{\left(\sqrt{{R_x}^2+{R_y}^2+{R_z}^2}\right)}^3$ of equation (A.11) is an overall scaling factor which cancels when establishing a ratio between ${\vec{T}}_{\mathrm{pp}}$ components (similar to equation (A.4)), as given in equation (A.15). Therefore, it is sufficient to write here

$$\begin{array}{c}{T}_x=-3·R_z·R_x,\\ T_y=-3·R_z·R_y,\\ T_z=-2·{R_z}^2+{R_x}^2+{R_y}^2.\end{array}$$(A.16)

Equation (A.15) needs the squares of the equation (A.16) terms

$$\begin{array}{c}{T_x}^2=9·{R_z}^2·{R_x}^2,\\ {T_y}^2=9·{R_z}^2·{R_y}^2,\\ {T_z}^2={\left(-2·{R_z}^2+{R_x}^2+{R_y}^2\right)}^2=4·{R_z}^4+{R_x}^4+{R_y}^4-4·{R_z}^2·{R_x}^2-4·{R_z}^2·{R_y}^2+2·{R_x}^2·{R_y}^2.\end{array}$$(A.17)

Thus is (numerator square root term of equation (A.15)):

$$9·{T_z}^2+8·{T_x}^2+8{{·T}_y}^2=36·{R_z}^4+{9·R_x}^4+{{9·R}_y}^4-36·{R_z}^2·{R_x}^2-36·{R_z}^2·{R_y}^2+18·{R_x}^2·{R_y}^2+72·{R_z}^2·{R_x}^2+72·{R_z}^2·{R_y}^2 =9·{\left(+2·{R_z}^2+{R_x}^2+{R_y}^2\right)}^2.$$

Compare with equation (A.17) for ${T_z}^2$ and

$$\sqrt{9·{T_z}^2+8·{T_x}^2+8{{·T}_y}^2}=3·\left(2·{R_z}^2+{R_x}^2+{R_y}^2\right).$$(A.18)

Substituting T_z from equation (A.16) and equations (A.17) and (A.18) into equation (A.15)

$$\tan {{\phi }_m}_{\mathrm{1,2}}=\frac{-3·T_z\pm \sqrt{9·{T_z}^2+8·{T_x}^2+8{{·T}_y}^2}}{4·\sqrt{{T_x}^2+{T_y}^2}}=\frac{-3·\left(-2·{R_z}^2+{R_x}^2+{R_y}^2\right)\pm 3·\left(2·{R_z}^2+{R_x}^2+{R_y}^2\right)}{4·\sqrt{9·{R_z}^2·{R_x}^2+9·{R_z}^2·{R_y}^2}} =\frac{6·{R_z}^2-3·{R_x}^2-3·{R_y}^2\pm 6·{R_z}^2\pm 3·{R_x}^2\pm 3·{R_y}^2}{12·R_z·\sqrt{{R_x}^2+{R_y}^2}}.$$(A.19)