Open Access
Issue
J. Space Weather Space Clim.
Volume 14, 2024
Article Number 36
Number of page(s) 20
DOI https://doi.org/10.1051/swsc/2024016
Published online 28 November 2024

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

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 (RE) 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 Lpp 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 Lpp. 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 ≤ Lpp ≤ 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−3 density thresholds for the plasmasphere outer edge (L100). The L100 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 Lpp-function is embedded in a 3D modelling algorithm allowing to fit of the function parameters directly from observed plasmapause torus position vectors Rpp$ {\vec{R}}_{{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 Neo$ {N}_e^o$ are generally replaced by smoothed measurements via a 3-point moving average approach

Ne(lk) :=13i=-11Neo(lk+i),$$ {N}_e\left({l}_k\right)\enspace:=\frac{1}{3}{\sum }_{i=-1}^1{N}_e^o\left({l}_{k+i}\right), $$(2.1)

where lk runs through the profile on which Lpp 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 Ne to an average of 50 values per profile in total in order to apply the same plasmapause detection algorithm.

The point

LargminlllogNe(l)$$ L\mathbf{Error:02254} \underset{l}{\mathrm{argmin}}{\partial }_l\mathrm{log}{N}_e(l) $$(2.2)

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

llogNe(lk):=logNe(lk+1)- logNe(lk-1)lk+1-lk-1.$$ {\partial }_l\mathrm{log}{N}_e\left({l}_k\right):=\frac{\mathrm{log}{N}_e\left({l}_{k+1}\right)-\enspace \mathrm{log}{N}_e\left({l}_{k-1}\right)}{{l}_{k+1}-{l}_{k-1}}. $$(2.3)

thumbnail 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 ̅$ \overline{\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 Lpp. 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.

llogNe(L)< llogNe̅-2σ(llogNe),$$ {\partial }_l\mathrm{log}{N}_e(L) < \enspace \overline{{\partial }_l\mathrm{log}{N}_e}-2\sigma \left({\partial }_l\mathrm{log}{N}_e\right), $$(2.4)

where ̅$ \overline{\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 Ne is expected to be small in relation to the other changes, so we assume

ΔlogNe(L)< ΔlogNe̅-2σ(ΔlogNe)$$ \Delta \mathrm{log}{N}_e(L) < \enspace \overline{\Delta \mathrm{log}{N}_e}-2\sigma (\Delta \mathrm{log}{N}_e) $$(2.5)

with ΔlogNe(lk) ≔ logNe(lk) − logNe(lk−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 Lpp. 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.,

llogNe(l)  llogNe̅-1.5 σ(llogNe).$$ {\partial }_l\mathrm{log}{N}_e(l)\enspace \ge \enspace \overline{{\partial }_l\mathrm{log}{N}_e}-1.5\enspace {\sigma }\left({\partial }_l\mathrm{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 Lpp is given by a single sharp gradient and change in Ne 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, Lpp 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 Lpp 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, Lppmax$ {L}_{{{pp}}_{{max}}}$, typically points to the evening-afternoon hours and the Earth-nearest point, Lppmin$ {L}_{{{pp}}_{{min}}}$, is located in the opposite direction. The ellipse radii, among them Lppmax$ {L}_{{{pp}}_{{max}}}$ and Lppmin$ {L}_{{{pp}}_{{min}}}$, are equal to the distances from the focal point coinciding with the Earth to the ellipse periphery, in L-units.

thumbnail Figure 2

Meaning of the orientation angle ϖ, the semi-major axis a, Lppmax$ {L}_{{{pp}}_{{max}}}$ and Lppmin$ {L}_{{{pp}}_{{min}}}$.

3.1 Modelling of Lpp 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 Lpp, 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 Lpp-ellipse is modelled with the formula describing the orbital radius of a satellite (e.g., Escobal, 1965)

Lpp=a·(1-e2)1+e·cosυ,$$ {L}_{{pp}}=\frac{a\middot \left(1-{e}^2\right)}{1+e\middot \mathrm{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, Lpp the true anomaly is expressed as

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

where

x … magnetic local time expressed in radians, i.e. x=π12 h ·MLT$ x=\frac{\pi }{12\enspace \mathrm{h}\enspace }\middot {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 Lpp 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):

Lpp=a·(1-e2)1+e·cos(x+ϖ).$$ {L}_{{pp}}=\frac{a\middot \left(1-{e}^2\right)}{1+e\middot \mathrm{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

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

where

MMLT 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:

Table 1

Dst levels and their ranges used for ellipse parameters fitting.

Table 2

Empirical functions f(Dst) fitted to represent the three ellipse parameters: semi-major axis a, eccentricity e, and M denoting the MLT of the plasmapause bulge as obtained by reformulating equation (3.4).

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 Lpp-parameters with a 3D modelling approach

The IMAGE and Van Allen Probes recorded 3D plasmapause positions in the SM system. Of these, Lpp-values were derived, indicating corresponding plasmapause locations in the geomagnetic equatorial plane. Both 3D plasmapause positions and derived Lpp-values are listed in the IMAGE and Van Allen Probes data records. It is not fully clear how the Lpp were projected from higher latitudes to the equatorial plane, probably by applying the dipole assumption Lpp = Rpp/cos2φ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 Lpp-function, equation (3.3), and allows to fit its parameters directly from observed plasmapause torus position vectors Rpp$ {\vec{R}}_{{pp}}$. This allows to account for a non-dipole geomagnetic field and avoids a projection of nonzero latitudinal measured radii Rpp into the geomagnetic equatorial plane according to the Lpp = Rpp/cos2φ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 |B|$ \left|\vec{B}\right|$,

  2. the invariant I=QQ'P|| ds$ I={\oint }_Q^{{Q\prime}{P}_{||}\enspace \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, |B|$ \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, (|B|,I)=const$ \left(\left|\vec{B}\right|,I\right)={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(|B|,I)$ 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

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

It is obtained by solving the relation Lpp = Rpp/cos2φ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 Lpp values are determined during plasmapause detection, cf. Section 2, and a cos2φm term is computed now with a given B$ \vec{B}$ vector (dipole or non-dipole). In this way, both approaches, invariant latitude and modifying the cos2φm term, are somehow equivalent ways to take the non-dipole aspect into account.

In the 3D modelling approach, the full 3D position vectors Rpp$ {\vec{R}}_{{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

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

where

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

  • LL-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 R$ \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 ε=R|R|$ \overrightarrow{\epsilon }=\raisebox{1ex}{$\vec{R}$}\!\left/ \!\raisebox{-1ex}{$\left|\vec{R}\right|$}\right.$ pointing to the position in space, and the other two are scalar quantities. When describing the plasmapause position Rpp$ {\vec{R}}_{{pp}}$, the first factor (1) is substituted by Lpp according to equation (3.3). In principle, any Lpp 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

φm=arctan(-0.75·Bz+0.5625·Bz2+0.5·(Bx2+By2)Bx2+By2)withB={BxByBz},$$ \begin{array}{ccc}{\phi }_m=\mathrm{arctan}\left(\frac{-0.75\middot {B}_z+\sqrt{0.5625\middot {{B}_z}^2+0.5\middot \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 B$ \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 Rpp$ {\vec{R}}_{{pp}}$ and related epoch. However, the B$ \vec{B}$ vector components need to be converted into the SM system before entering into equation (3.7) since the Lpp-function (3.3) is referred to SM (see Annex B).

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

Rx+vx=Lppa|0·a·cos2φm·εx+Lppe|0·e·cos2φm·εx+Lppϖ|0·ϖ ·cos2φm·εx,Ry+vy=Lppa|0·a·cos2φm·εy+Lppe|0·e·cos2φm·εy+Lppϖ|0·ϖ ·cos2φm·εy,Rz+vz=Lppa|0·a·cos2φm·εz+Lppe|0·e·cos2φm·εz+Lppϖ|0·ϖ ·cos2φm·εz .$$ \begin{array}{ccc}{\Delta R}_x+{v}_x& =& \frac{\partial {L}_{{pp}}}{{\partial a}}{|}_0\middot \Delta a\middot {\mathrm{cos}}^2{\phi }_m\middot {\epsilon }_x+\frac{\partial {L}_{{pp}}}{{\partial e}}{|}_0\middot \Delta e\middot {\mathrm{cos}}^2{\phi }_m\middot {\epsilon }_x+\frac{\partial {L}_{{pp}}}{\partial \varpi }{|}_0\middot \Delta \varpi \enspace \middot {\mathrm{cos}}^2{\phi }_m\middot {\epsilon }_x,\\ \Delta {R}_y+{v}_y& =& \frac{\partial {L}_{{pp}}}{{\partial a}}{|}_0\middot \Delta a\middot {\mathrm{cos}}^2{\phi }_m\middot {\epsilon }_y+\frac{\partial {L}_{{pp}}}{{\partial e}}{|}_0\middot \Delta e\middot {\mathrm{cos}}^2{\phi }_m\middot {\epsilon }_y+\frac{\partial {L}_{{pp}}}{\partial \varpi }{|}_0\middot \Delta \varpi \enspace \middot {\mathrm{cos}}^2{\phi }_m\middot {\epsilon }_y,\\ \Delta {R}_z+{v}_z& =& \frac{\partial {L}_{{pp}}}{{\partial a}}{|}_0\middot \Delta a\middot {\mathrm{cos}}^2{\phi }_m\middot {\epsilon }_z+\frac{\partial {L}_{{pp}}}{{\partial e}}{|}_0\middot \Delta e\middot {\mathrm{cos}}^2{\phi }_m\middot {\epsilon }_z+\frac{\partial {L}_{{pp}}}{\partial \varpi }{|}_0\middot \Delta \varpi \enspace \middot {\mathrm{cos}}^2{\phi }_m\middot {\epsilon }_z\enspace.\end{array} $$(3.8)

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

[Rx1Ry1Rz1RxnRynRzn]+[vx1vy1vz1vxnvynvzn]=[Lpp1a|0·cos2φ1·εx1Lpp1a|0·cos2φ1·εy1Lpp1a|0·cos2φ1·εz1Lpp1e|0·cos2φ1·εx1Lpp1e|0·cos2φ1·εy1Lpp1e|0·cos2φ1·εz1Lpp1ϖ|0·cos2φ1·εx1Lpp1ϖ|0·cos2φ1·εy1Lpp1ϖ|0·cos2φ1·εz1Lppna|0·cos2φn·εxnLppna|0·cos2φn·εynLppna|0·cos2φn·εznLppne|0·cos2φn·εxnLppne|0·cos2φn·εynLppne|0·cos2φn·εznLppnϖ|0·cos2φn·εxnLppnϖ|0·cos2φn·εynLppnϖ|0·cos2φn·εzn]*[aeϖ],$$ \begin{array}{ccc}\left[\begin{array}{c}\begin{array}{c}{\Delta R}_{{x}_1}\\ {\Delta R}_{{y}_1}\\ {\Delta R}_{{z}_1}\end{array}\\ \vdots \\ \begin{array}{c}{\Delta R}_{{x}_n}\\ {\Delta R}_{{y}_n}\\ {\Delta 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}\\ \vdots \\ \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}_{{{pp}}_1}}{{\partial a}}{|}_0\middot {\mathrm{cos}}^2{\phi }_1\middot {\epsilon }_{{x}_1}\\ \frac{\partial {L}_{{{pp}}_1}}{{\partial a}}{|}_0\middot {\mathrm{cos}}^2{\phi }_1\middot {\epsilon }_{{y}_1}\\ \frac{\partial {L}_{{{pp}}_1}}{{\partial a}}{|}_0\middot {\mathrm{cos}}^2{\phi }_1\middot {\epsilon }_{{z}_1}\end{array}& \begin{array}{c}\frac{\partial {L}_{{{pp}}_1}}{{\partial e}}{|}_0\middot {\mathrm{cos}}^2{\phi }_1\middot {\epsilon }_{{x}_1}\\ \frac{\partial {L}_{{{pp}}_1}}{{\partial e}}{|}_0\middot {\mathrm{cos}}^2{\phi }_1\middot {\epsilon }_{{y}_1}\\ \frac{\partial {L}_{{{pp}}_1}}{{\partial e}}{|}_0\middot {\mathrm{cos}}^2{\phi }_1\middot {\epsilon }_{{z}_1}\end{array}& \begin{array}{c}\frac{\partial {L}_{{{pp}}_1}}{\partial \varpi }{|}_0\middot {\mathrm{cos}}^2{\phi }_1\middot {\epsilon }_{{x}_1}\\ \frac{\partial {L}_{{{pp}}_1}}{\partial \varpi }{|}_0\middot {\mathrm{cos}}^2{\phi }_1\middot {\epsilon }_{{y}_1}\\ \frac{\partial {L}_{{{pp}}_1}}{\partial \varpi }{|}_0\middot {\mathrm{cos}}^2{\phi }_1\middot {\epsilon }_{{z}_1}\end{array}\\ \vdots & \vdots & \vdots \\ \begin{array}{c}\frac{\partial {L}_{{{pp}}_n}}{{\partial a}}{|}_0\middot {\mathrm{cos}}^2{\phi }_n\middot {\epsilon }_{{x}_n}\\ \frac{\partial {L}_{{{pp}}_n}}{{\partial a}}{|}_0\middot {\mathrm{cos}}^2{\phi }_n\middot {\epsilon }_{{y}_n}\\ \frac{\partial {L}_{{{pp}}_n}}{{\partial a}}{|}_0\middot {\mathrm{cos}}^2{\phi }_n\middot {\epsilon }_{{z}_n}\end{array}& \begin{array}{c}\frac{\partial {L}_{{{pp}}_n}}{{\partial e}}{|}_0\middot {\mathrm{cos}}^2{\phi }_n\middot {\epsilon }_{{x}_n}\\ \frac{\partial {L}_{{{pp}}_n}}{{\partial e}}{|}_0\middot {\mathrm{cos}}^2{\phi }_n\middot {\epsilon }_{{y}_n}\\ \frac{\partial {L}_{{{pp}}_n}}{{\partial e}}{|}_0\middot {\mathrm{cos}}^2{\phi }_n\middot {\epsilon }_{{z}_n}\end{array}& \begin{array}{c}\frac{\partial {L}_{{{pp}}_n}}{\partial \varpi }{|}_0\middot {\mathrm{cos}}^2{\phi }_n\middot {\epsilon }_{{x}_n}\\ \frac{\partial {L}_{{{pp}}_n}}{\partial \varpi }{|}_0\middot {\mathrm{cos}}^2{\phi }_n\middot {\epsilon }_{{y}_n}\\ \frac{\partial {L}_{{{pp}}_n}}{\partial \varpi }{|}_0\middot {\mathrm{cos}}^2{\phi }_n\middot {\epsilon }_{{z}_n}\end{array}\end{array}\right]& *& \left[\begin{array}{c}\Delta a\\ \Delta e\\ \Delta \varpi \end{array}\right]\end{array}\end{array}\end{array}, $$(3.9)

where for the non-linear fit

  • Rpp=RppObs-RppMod=RppObs-LppMod·cos2φm·εpp$ \Delta {\vec{R}}_{{pp}}={\vec{R}}_{{ppObs}}-{\vec{R}}_{{ppMod}}={\vec{R}}_{{ppObs}}-{L}_{{ppMod}}\middot {\mathrm{cos}}^2{\phi }_m\middot {\overrightarrow{\epsilon }}_{{pp}}$,

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

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

  • cos2φm … computed from B$ \vec{B}$ vector, equation (3.7),

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

  • v … fitted correction to the observable,

  • Lppa|0$ \frac{\partial {L}_{{pp}}}{{\partial a}}{|}_0$ … partial of Lpp w.r.t. semi-major axis a, evaluated with initial values,

  • a … estimated correction to initial value of a,

  • Lppe|0$ \frac{\partial {L}_{{pp}}}{{\partial e}}{|}_0$ … partial of Lpp w.r.t. eccentricity e, evaluated with initial values,

  • e … estimated correction to initial value of e,

  • Lppϖ|0$ \frac{\partial {L}_{{pp}}}{\partial \varpi }{|}_0$ … partial of Lpp 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 Rpp$ {\vec{R}}_{{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 εg={cosφ·cos λcosφ·sinλsinφ}$ {\overrightarrow{\epsilon }}_g=\left\{\begin{array}{c}\mathrm{cos}\phi \middot \mathrm{cos}\enspace \lambda \\ \mathrm{cos}\phi \middot \mathrm{sin}\lambda \\ \mathrm{sin}\phi \end{array}\right\}$.

  3. Transform εg$ {\overrightarrow{\epsilon }}_g$ from geographic into SM ➧ unit vector εpp$ {\overrightarrow{\epsilon }}_{{pp}}$.

  4. Compute MLT = Λ·12 hπ+12 h$ {MLT}\enspace =\enspace \Lambda \middot \frac{12\enspace \mathrm{h}}{\pi }+12\enspace \mathrm{h}$, with Λ=arctan(εppyεppx)$ \Lambda =\mathrm{arctan}(\frac{{\epsilon }_{{{pp}}_y}}{{\epsilon }_{{{pp}}_x}})$.

  5. Compute cos2φm term for dipole from the components of εpp$ {\overrightarrow{\epsilon }}_{{pp}}$ via

cos 2 φ m = ε pp x 2 + ε pp y 2 ε pp x 2 + ε pp y 2 + ε pp z 2 . $$ {\mathrm{cos}}^2{\phi }_m=\frac{{{\epsilon }_{{{pp}}_x}}^2+{{\epsilon }_{{{pp}}_y}}^2}{{{\epsilon }_{{{pp}}_x}}^2+{{\epsilon }_{{{pp}}_y}}^2+{{\epsilon }_{{{pp}}_z}}^2}. $$

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

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

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

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

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

  6. Transform non-dipole Rpp$ {\vec{R}}_{{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 Lpp-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(Dst)=A·exp(DstmedDscal)+C,$$ f\left({Dst}\right)=A\middot \mathrm{exp}\left(\frac{D{stmed}}{{Dscal}}\right)+C, $$(3.10)

where

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

  • DscalDst 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, Ca = 1.7 RE, Ce = 0.27, CM = 12 h, Table 2, also boundary values for Dst = 0 had to be defined, which were 5 RE, 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.

thumbnail 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 Lppmin1.4$ {L}_{{{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 Dstmin ≈ −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 Lppmin=1.3$ {L}_{{{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 Lpp by the relations

Lppmin=a·(1-e)$$ {L}_{{{pp}}_{{min}}}=a\middot (1-e) $$(3.11)

and

Lppmax=a·(1+e).$$ {L}_{{{pp}}_{{max}}}=a\middot \left(1+e\right). $$(3.12)

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

Lppmax-Lppmin=2a·e.$$ {L}_{{{pp}}_{{max}}}-{L}_{{{pp}}_{{min}}}=2a\middot {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 Lppmax$ {L}_{{{pp}}_{{max}}}$ and Lppmin$ {L}_{{{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 Lppmax$ {L}_{{{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 CM = 12 h. As with the other boundary values Ca and Ce 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.

thumbnail 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 Lpp 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.

thumbnail 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 Lpp values from IMAGE and Van Allen Probes as well as the modelled Lpp values for NEPPM (left) and GCPM (right), in the equatorial plane (top) as well as by Lpp over MLT (middle), and its corresponding RMS (bottom). The Lpp 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 RE, GCPM shows a higher mean RMS of 0.9 RE. 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 Lpp values per Dst level reveals more information. Figure 6 shows for each Dst level, the extracted Lpp values next to the resulting mean Lpp computed with both the NEPPM (black) and the GCPM (grey). At all Dst levels, GCPM reveals a consistent overestimation of the extracted Lpp values. This overestimation expands even further with strongly negative Dst levels. The bulge in GCPM is generally higher than the extracted Lpp 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.

thumbnail 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.

thumbnail 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 RE around the bulge at 19 MLT and the lowest RMS values with 0.4 RE 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 Lpp data sources for the RMS calculations. So, the pre-selection of our Lpp 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 Lppdatabase 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 (Lpp 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 Lpp, detecting a sharp gradient of the electron density is essential. Naturally, some measurements do not fulfill the requirements for reliable Lpp 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 Rpp$ {\vec{R}}_{{pp}}$ to a Lpp-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 Lpp = 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 Lpp 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 RE. 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 Rpp$ {\vec{R}}_{{pp}}$ that obey the representation equation (3.6). When computing them with a constant Lpp value in an equidistant grid, a plot of these Rpp$ {\vec{R}}_{{pp}}$ would display a torus-like picture, blue in Figure A1 below, corresponding to the dipole assumption inherent in the definition of L-shells.

thumbnail Figure A1

The vector fields Rpp$ {\vec{R}}_{{pp}}$ (blue) and Tpp$ {\vec{T}}_{{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

cos2φm=Rx2+Ry2Rx2+Ry2+Rz2withRpp={RxRyRz},$$ \begin{array}{ccc}{\mathrm{cos}}^2{\phi }_m=\frac{{{R}_x}^2+{{R}_y}^2}{{{R}_x}^2+{{R}_y}^2+{{R}_z}^2}& \mathrm{with}& {\vec{R}}_{{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 B$ \vec{B}$ if, like in the IGRF, a non-dipole field is the base. For this, B$ \vec{B}$ has to be computed with the selected geomagnetic field model at the position indicated by Rpp$ {\vec{R}}_{{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 cos2φm from B$ \vec{B}$ are once more derived from Rpp$ {\vec{R}}_{{pp}}$, i.e. for a dipole. Nevertheless, as will be shown below, the components from a non-dipole B$ \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 B$ \vec{B}$ vector at that point in the magnetic field.

First, reformulate equation (3.6) to

Rpp=Lpp·{cos3φm·cosΛcos3φm·sinΛcos2φm·sinφm}$$ {\vec{R}}_{{pp}}={L}_{{pp}}\middot \left\{\begin{array}{c}{\mathrm{cos}}^3{\phi }_m\middot \mathrm{cos}\Lambda \\ {\mathrm{cos}}^3{\phi }_m\middot \mathrm{sin}\Lambda \\ {\mathrm{cos}}^2{\phi }_m\middot \mathrm{sin}{\phi }_m\end{array}\right\} $$(A.2)

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

Tpp=Rppφm=Lpp·{-3·cos2φm·sinφm·cosΛ-3·cos2φm·sinφm·sinΛ-2·cosφm·sin2φm+cos3φm},$$ {\vec{T}}_{{pp}}=\frac{\partial {\vec{R}}_{{pp}}}{\partial {\phi }_m}={L}_{{pp}}\middot \left\{\begin{array}{c}-3\middot {\mathrm{cos}}^2{\phi }_m\middot \mathrm{sin}{\phi }_m\middot \mathrm{cos}\Lambda \\ -3\middot {\mathrm{cos}}^2{\phi }_m\middot \mathrm{sin}{\phi }_m\middot \mathrm{sin}\Lambda \\ -2\middot \mathrm{cos}{\phi }_m\middot {\mathrm{sin}}^2{\phi }_m+{\mathrm{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 B$ \vec{B}$. The geometry behind Rpp$ {\vec{R}}_{{pp}}$ and Tpp$ {\vec{T}}_{{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 TppB$ {\vec{T}}_{{pp}}\uparrow \uparrow \vec{B}$, analogously to spherical latitude computation from Cartesian coordinates, one can build the ratio

TzTx2+Ty2=BzBx2+By2withTpp={TxTyTz}, B={BxByBz}.$$ \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}}_{{pp}}=\left\{\begin{array}{c}{T}_x\\ {T}_y\\ {T}_z\end{array}\right\},\enspace \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 cos2Λ + sin2Λ = 1, see e.g. Sigl (1977), the denominator of the left-hand side in (A.4) can be condensed to

Tx2+Ty2=3·cos2φm·sinφm.$$ \sqrt{{{T}_x}^2+{{T}_y}^2}=3\middot {\mathrm{cos}}^2{\phi }_m\middot \mathrm{sin}{\phi }_m. $$(A.5)

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

BzBx2+By2=-2·cosφm·sin2φm3·cos2φm·sinφm+cos3φm3·cos2φm·sinφm=-23·tanφm+13·1tanφm,$$ \frac{{B}_z}{\sqrt{{{B}_x}^2+{{B}_y}^2}}=-\frac{2\middot \mathrm{cos}{\phi }_m\middot {\mathrm{sin}}^2{\phi }_m}{3\middot {\mathrm{cos}}^2{\phi }_m\middot \mathrm{sin}{\phi }_m}+\frac{{\mathrm{cos}}^3{\phi }_m}{3\middot {\mathrm{cos}}^2{\phi }_m\middot \mathrm{sin}{\phi }_m}=-\frac{2}{3}\middot \mathrm{tan}{\phi }_m+\frac{1}{3}\middot \frac{1}{\mathrm{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

(tanφm)2+32·BzBx2+By2p·(tanφm)-12q=0.$$ {\left(\mathrm{tan}{\phi }_m\right)}^2+\underset{p}{\underbrace{\frac{3}{2}\middot \frac{{B}_z}{\sqrt{{{B}_x}^2+{{B}_y}^2}}}}\middot \left(\mathrm{tan}{\phi }_m\right)\underset{q}{\underbrace{-\frac{1}{2}}}=0. $$(A.7)

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

x1,2=-p2±(p2)2-q.$$ {x}_{\mathrm{1,2}}=-\frac{p}{2}\pm \sqrt{{\left(\frac{p}{2}\right)}^2-{q}.} $$(A.8)

Leading to

(tanφm)1,2=-34·BzBx2+By2±916·Bz2Bx2+By2+12.$$ {\left(\mathrm{tan}{\phi }_m\right)}_{\mathrm{1,2}}=-\frac{3}{4}\middot \frac{{B}_z}{\sqrt{{{B}_x}^2+{{B}_y}^2}}\pm \sqrt{\frac{9}{16}\middot \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 Bz = 0 the dipole loop reaches its maximal/minimal elevation above/below the geomagnetic equator, i.e. (tanφm)1,2=±12$ {\left(\mathrm{tan}{\phi }_m\right)}_{\mathrm{1,2}}=\pm \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\sqrt{2}$}\right.$, 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 Tpp$ {\vec{T}}_{{pp}}$ have to be expressed as a function of the components {Rx Ry Rz} of the plasmapause position vector Rpp$ {\vec{R}}_{{pp}}$. Therefore in equation (A.3) the cosine and sine terms are formulated as functions of the plasmapause position vector components

cosφm=Rx2+Ry2Rx2+Ry2+Rz2, sinφm=RzRx2+Ry2+Rz2, cosΛ=RxRx2+Ry2, sinΛ=RyRx2+Ry2.$$ \mathrm{cos}{\phi }_m=\frac{\sqrt{{{R}_x}^2+{{R}_y}^2}}{\sqrt{{{R}_x}^2+{{R}_y}^2+{{R}_z}^2}},\enspace \mathrm{sin}{\phi }_m=\frac{{R}_z}{\sqrt{{{R}_x}^2+{{R}_y}^2+{{R}_z}^2}},\enspace \mathrm{cos}\Lambda =\frac{{R}_x}{\sqrt{{{R}_x}^2+{{R}_y}^2}},\enspace \mathrm{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

Tpp=Lpp·Rx2+Ry2(Rx2+Ry2+Rz2)3·{-3·Rz·Rx-3·Rz·Ry-2·Rz2+Rx2+Ry2}.$$ {\vec{T}}_{{pp}}={L}_{{pp}}\middot \frac{\sqrt{{{R}_x}^2+{{R}_y}^2}}{{\left(\sqrt{{{R}_x}^2+{{R}_y}^2+{{R}_z}^2}\right)}^3}\middot \left\{\begin{array}{c}-3\middot {R}_z\middot {R}_x\\ -3\middot {R}_z\middot {R}_y\\ -2\middot {{R}_z}^2+{{R}_x}^2+{{R}_y}^2\end{array}\right\}. $$(A.11)

Whereby for describing purely the direction of Tpp$ {\vec{T}}_{{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 Rpp$ {\vec{R}}_{{pp}}$ is computed for a given geomagnetic latitude φm and solar longitude Λ according to equation (A.2). Reversely, φm can be computed from the Rpp$ {\vec{R}}_{{pp}}$ Cartesian components in a backwards transformation

tanφm=RzRx2+Ry2.$$ \mathrm{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 B$ \vec{B}$ vector components can be replaced by the tangent vector Tpp$ {\vec{T}}_{{pp}}$ components. This approach will, after some algebra, allow to express equation (A.9) in terms of the Rpp$ {\vec{R}}_{{pp}}$ vector Cartesian components and finally lead to tan φm formulae giving both φm solutions of the quadratic equation as a function of {Rx Ry Rz}. Thus equation (A.9) is now expressed in terms of the Tpp$ {\vec{T}}_{{pp}}$ components {Tx Ty Tz}

(tanφm)1,2=-34·TzTx2+Ty2±916·Tz2Tx2+Ty2+12.$$ {\left(\mathrm{tan}{\phi }_m\right)}_{\mathrm{1,2}}=-\frac{3}{4}\middot \frac{{T}_z}{\sqrt{{{T}_x}^2+{{T}_y}^2}}\pm \sqrt{\frac{9}{16}\middot \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

±916·Tz2Tx2+Ty2+12=±916·Tz2Tx2+Ty2+816·Tx2+Ty2Tx2+Ty2=±9·Tz2+8·Tx2+8·Ty24·Tx2+Ty2.$$ \pm \sqrt{\frac{9}{16}\middot \frac{{{T}_z}^2}{{{T}_x}^2+{{T}_y}^2}+\frac{1}{2}}=\pm \sqrt{\frac{9}{16}\middot \frac{{{T}_z}^2}{{{T}_x}^2+{{T}_y}^2}+\frac{8}{16}\middot \frac{{{T}_x}^2+{{T}_y}^2}{{{T}_x}^2+{{T}_y}^2}}=\pm \frac{\sqrt{9\middot {{T}_z}^2+8\middot {{T}_x}^2+8{{\middot T}_y}^2}}{4\middot \sqrt{{{T}_x}^2+{{T}_y}^2}}. $$(A.14)

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

(tanφm)1,2=-34·TzTx2+Ty2±9·Tz2+8·Tx2+8·Ty24·Tx2+Ty2=-3·Tz±9·Tz2+8·Tx2+8·Ty24·Tx2+Ty2.$$ {\left(\mathrm{tan}{\phi }_m\right)}_{\mathrm{1,2}}=-\frac{3}{4}\middot \frac{{T}_z}{\sqrt{{{T}_x}^2+{{T}_y}^2}}\pm \frac{\sqrt{9\middot {{T}_z}^2+8\middot {{T}_x}^2+8{{\middot T}_y}^2}}{4\middot \sqrt{{{T}_x}^2+{{T}_y}^2}}=\frac{-3\middot {T}_z\pm \sqrt{9\middot {{T}_z}^2+8\middot {{T}_x}^2+8{{\middot T}_y}^2}}{4\middot \sqrt{{{T}_x}^2+{{T}_y}^2}}. $$(A.15)

Now, the expressions of {Tx Ty Tz} in terms of {Rx Ry Rz}, equation (A.11), must be substituted into equation (A.15). Thereby the term Lpp·Rx2+Ry2(Rx2+Ry2+Rz2)3$ \raisebox{1ex}{${L}_{{pp}}\middot \sqrt{{{R}_x}^2+{{R}_y}^2}$}\!\left/ \!\raisebox{-1ex}{${\left(\sqrt{{{R}_x}^2+{{R}_y}^2+{{R}_z}^2}\right)}^3$}\right.$ of equation (A.11) is an overall scaling factor which cancels when establishing a ratio between Tpp$ {\vec{T}}_{{pp}}$ components (similar to equation (A.4)), as given in equation (A.15). Therefore, it is sufficient to write here

Tx=-3·Rz·Rx,Ty=-3·Rz·Ry,Tz=-2·Rz2+Rx2+Ry2.$$ \begin{array}{c}{T}_x=-3\middot {R}_z\middot {R}_x,\\ {T}_y=-3\middot {R}_z\middot {R}_y,\\ {T}_z=-2\middot {{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

Tx2=9·Rz2·Rx2,Ty2=9·Rz2·Ry2,Tz2=(-2·Rz2+Rx2+Ry2)2=4·Rz4+Rx4+Ry4-4·Rz2·Rx2-4·Rz2·Ry2+2·Rx2·Ry2.$$ \begin{array}{c}{{T}_x}^2=9\middot {{R}_z}^2\middot {{R}_x}^2,\\ {{T}_y}^2=9\middot {{R}_z}^2\middot {{R}_y}^2,\\ {{T}_z}^2={\left(-2\middot {{R}_z}^2+{{R}_x}^2+{{R}_y}^2\right)}^2=4\middot {{R}_z}^4+{{R}_x}^4+{{R}_y}^4-4\middot {{R}_z}^2\middot {{R}_x}^2-4\middot {{R}_z}^2\middot {{R}_y}^2+2\middot {{R}_x}^2\middot {{R}_y}^2.\end{array} $$(A.17)

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

9·Tz2+8·Tx2+8·Ty2=36·Rz4+9·Rx4+9·Ry4-36·Rz2·Rx2-36·Rz2·Ry2+18·Rx2·Ry2+72·Rz2·Rx2+72·Rz2·Ry2 =9·(+2·Rz2+Rx2+Ry2)2.$$ 9\middot {{T}_z}^2+8\middot {{T}_x}^2+8{{\middot T}_y}^2=36\middot {{R}_z}^4+{9\middot {R}_x}^4+{{9\middot R}_y}^4-36\middot {{R}_z}^2\middot {{R}_x}^2-36\middot {{R}_z}^2\middot {{R}_y}^2+18\middot {{R}_x}^2\middot {{R}_y}^2+72\middot {{R}_z}^2\middot {{R}_x}^2+72\middot {{R}_z}^2\middot {{R}_y}^2\enspace =9\middot {\left(+2\middot {{R}_z}^2+{{R}_x}^2+{{R}_y}^2\right)}^2. $$

Compare with equation (A.17) for Tz2$ {{T}_z}^2$ and

9·Tz2+8·Tx2+8·Ty2=3·(2·Rz2+Rx2+Ry2).$$ \sqrt{9\middot {{T}_z}^2+8\middot {{T}_x}^2+8{{\middot T}_y}^2}=3\middot \left(2\middot {{R}_z}^2+{{R}_x}^2+{{R}_y}^2\right). $$(A.18)

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

tanφm1,2=-3·Tz±9·Tz2+8·Tx2+8·Ty24·Tx2+Ty2=-3·(-2·Rz2+Rx2+Ry2)±3·(2·Rz2+Rx2+Ry2)4·9·Rz2·Rx2+9·Rz2·Ry2 =6·Rz2-3·Rx2-3·Ry2±6·Rz2±3·Rx2±3·Ry212·Rz·Rx2+Ry2.$$ \mathrm{tan}{{\phi }_m}_{\mathrm{1,2}}=\frac{-3\middot {T}_z\pm \sqrt{9\middot {{T}_z}^2+8\middot {{T}_x}^2+8{{\middot T}_y}^2}}{4\middot \sqrt{{{T}_x}^2+{{T}_y}^2}}=\frac{-3\middot \left(-2\middot {{R}_z}^2+{{R}_x}^2+{{R}_y}^2\right)\pm 3\middot \left(2\middot {{R}_z}^2+{{R}_x}^2+{{R}_y}^2\right)}{4\middot \sqrt{9\middot {{R}_z}^2\middot {{R}_x}^2+9\middot {{R}_z}^2\middot {{R}_y}^2}}\enspace =\frac{6\middot {{R}_z}^2-3\middot {{R}_x}^2-3\middot {{R}_y}^2\pm 6\middot {{R}_z}^2\pm 3\middot {{R}_x}^2\pm 3\middot {{R}_y}^2}{12\middot {R}_z\middot \sqrt{{{R}_x}^2+{{R}_y}^2}}. $$(A.19)

Now the two φm solutions of the quadratic equation (A.9) in tan φm, as function of {Rx Ry Rz}, can finally be computed:

Positive square root

tanφm1=6·Rz2-3·Rx2-3·Ry2+6·Rz2+3·Rx2+3·Ry212·Rz·Rx2+Ry2=RzRx2+Ry2,$$ \mathrm{tan}{{\phi }_m}_1=\frac{6\middot {{R}_z}^2-3\middot {{R}_x}^2-3\middot {{R}_y}^2+6\middot {{R}_z}^2+3\middot {{R}_x}^2+3\middot {{R}_y}^2}{12\middot {R}_z\middot \sqrt{{{R}_x}^2+{{R}_y}^2}}=\frac{{R}_z}{\sqrt{{{R}_x}^2+{{R}_y}^2}}, $$(A.20a)

Negative square root

tanφm2=6·Rz2-3·Rx2-3·Ry2-6·Rz2-3·Rx2-3·Ry212·Rz·Rx2+Ry2=-12·Rx2+Ry2Rz=-12·1tanφm1.$$ \mathrm{tan}{{\phi }_m}_2=\frac{6\middot {{R}_z}^2-3\middot {{R}_x}^2-3\middot {{R}_y}^2-6\middot {{R}_z}^2-3\middot {{R}_x}^2-3\middot {{R}_y}^2}{12\middot {R}_z\middot \sqrt{{{R}_x}^2+{{R}_y}^2}}=\frac{-\frac{1}{2}\middot \sqrt{{{R}_x}^2+{{R}_y}^2}}{{R}_z}=-\frac{1}{2}\middot \frac{1}{\mathrm{tan}{{\phi }_m}_1}. $$(A.20b)

When comparing equations (A.12) and (A.20a), one can see that the positive square root provides the sought-after solution.

And obviously between the two solutions the following relation holds

tanφm1·tanφm2=-12.$$ \mathrm{tan}{{\phi }_m}_1\middot \mathrm{tan}{{\phi }_m}_2=-\frac{1}{2}. $$(A.20c)

A2 Meaning of the negative square root solution

For completeness, also the plasmapause position vector Rpp2$ {\vec{R}}_{{{pp}}_2}$, corresponding to the negative square root of equation (A.9), shall be explained. It can be obtained by putting φm2$ {{\phi }_m}_2$, equation (A.20b), into equation (A.2). However, it can also be expressed in terms of the plasmapause position vector Rpp$ {\vec{R}}_{{pp}}$. Taking tanφm2$ \mathrm{tan}{\phi }_{{m}_2}$ from equation (A.20b) and regarding the following trigonometric relations (e.g. Sigl, 1977)

sinφm2=tanφm21+ tan2φm2, cosφm2=11+tan2φm2,$$ \mathrm{sin}{\phi }_{{m}_2}=\frac{\mathrm{tan}{\phi }_{{m}_2}}{\sqrt{1+\enspace {\mathrm{tan}}^2{\phi }_{{m}_2}}},\enspace \mathrm{cos}{\phi }_{{m}_2}=\frac{1}{\sqrt{1+{\mathrm{tan}}^2{\phi }_{{m}_2}}}, $$(A.21)

one gets by substituting the expression (A.20b) for tanφm2$ {{tan}{\phi }_m}_2$

sinφm2=-12·Rx2+Ry2Rz2+14·(Rx2+Ry2), cosφm2=RzRz2+14·(Rx2+Ry2).$$ \mathrm{sin}{\phi }_{{m}_2}=\frac{-\frac{1}{2}\middot \sqrt{{{R}_x}^2+{{R}_y}^2}}{\sqrt{{{R}_z}^2+\frac{1}{4}\middot \left({{R}_x}^2+{{R}_y}^2\right)}},\enspace \mathrm{cos}{\phi }_{{m}_2}=\frac{{R}_z}{\sqrt{{{R}_z}^2+\frac{1}{4}\middot \left({{R}_x}^2+{{R}_y}^2\right)}}. $$(A.22)

And since Rpp$ {\vec{R}}_{{pp}}$ and Rpp2$ {\vec{R}}_{{{pp}}_2}$ are lying on the same dipole loop, they share the same meridian, having the same solar longitude Λ, i.e. the cosΛ and sinΛ expressions provided by equations (A.10) hold also for Rpp2$ {\vec{R}}_{{{pp}}_2}$. Thus, by substituting equations (A.22) and (A.10) into equation (3.6) one obtains, after some algebra, Rpp2$ {\vec{R}}_{{{pp}}_2}$ as a function of the Rpp$ {\vec{R}}_{{pp}}$ components

Rpp2=Lpp·Rz2Rz2+14·(Rx2+Ry2)(1)·1Rz2+14·(Rx2+Ry2)·1Rx2+Ry2·{Rz·RxRz·Ry-12·(Rx2+Ry2)}(2).$$ {\vec{R}}_{{{pp}}_2}=\underset{(1)}{\underbrace{{L}_{{pp}}\middot \frac{{{R}_z}^2}{{{R}_z}^2+\frac{1}{4}\middot \left({{R}_x}^2+{{R}_y}^2\right)}}}\middot \underset{(2)}{\underbrace{\frac{1}{\sqrt{{{R}_z}^2+\frac{1}{4}\middot \left({{R}_x}^2+{{R}_y}^2\right)}}\middot \frac{1}{\sqrt{{{R}_x}^2+{{R}_y}^2}}\middot \left\{\begin{array}{c}{R}_z\middot {R}_x\\ {R}_z\middot {R}_y\\ -\frac{1}{2}\middot \left({{R}_x}^2+{{R}_y}^2\right)\end{array}\right\}}}. $$(A.23)

Thereby term (1) in equation (A.23) gives the Rpp2$ {\vec{R}}_{{{pp}}_2}$ magnitude, (1)=Lpp·cos2φm2$ {(1)=L}_{{pp}}\middot {\mathrm{cos}}^2{\phi }_{{m}_2}$, compare with equation (A.22), bottom, and term (2) is the unit vector of Rpp2$ {\vec{R}}_{{{pp}}_2}$, (2)=εpp2$ (2)={\overrightarrow{\epsilon }}_{{{pp}}_2}$. When computing the magnitude of the rightmost braced vector term {⋮} in (2) one recognizes that this magnitude is identical to the product of the denominators of the two preceding terms in (2), i.e. the product of all three terms in (2) results in the unit vector εpp2$ {\overrightarrow{\epsilon }}_{{{pp}}_2}$.

By the relation Lpp=Rpp/cos2φm1$ {L}_{{pp}}={R}_{{pp}}/{\mathrm{cos}}^2{\phi }_{{m}_1}$ and equation (A.1), also Lpp can be expressed as a function of the Rpp$ {\vec{R}}_{{pp}}$ components

Lpp=(Rx2+Ry2+Rz2)3Rx2+Ry2.$$ {L}_{{pp}}=\frac{{\left(\sqrt{{{R}_x}^2+{{R}_y}^2+{{R}_z}^2}\right)}^3}{{{R}_x}^2+{{R}_y}^2}. $$(A.24)

Thus with equation (A.24) Rpp2$ {\vec{R}}_{{{pp}}_2}$ can be fully formulated as a function of the Rpp$ {\vec{R}}_{{pp}}$ components. Starting from equation (A.23) by merging the two reciprocal square roots of the term (2) with the term (1), some algebra leads to

Rpp2=Rz2·(Rx2+Ry2+Rz2(Rx2+Ry2)·(Rz2+14·(Rx2+Ry2)))32·{Rz·RxRz·Ry-12·(Rx2+Ry2)}.$$ {\vec{R}}_{{{pp}}_2}={{R}_z}^2\middot {\left(\frac{{{R}_x}^2+{{R}_y}^2+{{R}_z}^2}{\left({{R}_x}^2+{{R}_y}^2\right)\middot \left({{R}_z}^2+\frac{1}{4}\middot \left({{R}_x}^2+{{R}_y}^2\right)\right)}\right)}^{\frac{3}{2}}\middot \left\{\begin{array}{c}{R}_z\middot {R}_x\\ {R}_z\middot {R}_y\\ -\frac{1}{2}\middot \left({{R}_x}^2+{{R}_y}^2\right)\end{array}\right\}. $$(A.25)

And in an analogous exercise also a formula for the angle enclosed by Rpp$ {\vec{R}}_{{pp}}$ and Rpp2$ {\vec{R}}_{{{pp}}_2}$ can be established, by using the unit vectors εpp$ {\overrightarrow{\epsilon }}_{{pp}}$ and εpp2$ {\overrightarrow{\epsilon }}_{{{pp}}_2}$,

cosϑ=(εpp·εpp2)=12·Rz·Rx2+Ry2(Rx2+Ry2+Rz2)·(Rz2+14·(Rx2+Ry2)).$$ \mathrm{cos}\vartheta =\left({\overrightarrow{\epsilon }}_{{pp}}\middot {\overrightarrow{\epsilon }}_{{{pp}}_2}\right)=\frac{1}{2}\middot {R}_z\middot \sqrt{\frac{{{R}_x}^2+{{R}_y}^2}{\left({{R}_x}^2+{{R}_y}^2+{{R}_z}^2\right)\middot \left({{R}_z}^2+\frac{1}{4}\middot \left({{R}_x}^2+{{R}_y}^2\right)\right)}}. $$(A.26)

By an appropriate exchange of signs of the Rpp$ {\vec{R}}_{{pp}}$ and Rpp2$ {\vec{R}}_{{{pp}}_2}$ components, these two vectors can also be mirrored symmetrically to other points on their magnetic loop as well as on the opposite – and perpendicular loops, as far as the magnetic field is considered to be perfect dipole.

The Rpp$ {\vec{R}}_{{pp}}$ vs. Rpp2$ {\vec{R}}_{{{pp}}_2}$ behaviour has been investigated in some detail. However, a complete presentation of the outcome would be out of the scope of this paper. Nevertheless, the most essential results, found so far, can be summarized in short as follows:

  1. The angle ϑ enclosed by Rpp$ {\vec{R}}_{{pp}}$ and Rpp2$ {\vec{R}}_{{{pp}}_2}$ ranges between 70°.5 at φm=±35°.3=arctan(±12)$ {\phi }_m=\pm 35{}^{\circ}.3=\mathrm{arctan}\left(\raisebox{1ex}{$\pm 1$}\!\left/ \!\raisebox{-1ex}{$\sqrt{2}$}\right.\right)$ and 90° at φm = 0°, ±90°.

  2. φm = ±90° is a singularity point. Here the dipole curve passes through the geocentre and the plasmapause position vectors vanish.

  3. Rpp2=0$ {\vec{R}}_{{{pp}}_2}=\vec{0}$ when Rpp$ {\vec{R}}_{{pp}}$ is pointing to the equator (Rz = 0). Thus φm = 0° has been skipped in Rpp2$ {\vec{R}}_{{{pp}}_2}$ computation.

  4. If Rpp$ {\vec{R}}_{{pp}}$ has a negative z-component Rz, the resulting Rpp2$ {\vec{R}}_{{{pp}}_2}$, computed with equation (A.25), is pointing into the opposite dipole loop. To get in this case Rpp2$ {\vec{R}}_{{{pp}}_2}$ into the correct dipole loop, it has simply to be turned in sign, i.e. this problem can be overcome by multiplying with sign(Rz).

  5. Also cosϑ, computed with equation (A.26), has to be turned in sign when Rz is negative, i.e. by multiplying with sign(Rz), due to the same reason as for Point 4.

  6. When equation (A.20b) is solved with the atan function instead of atan2, the obtained Rpp2$ {\vec{R}}_{{{pp}}_2}$ latitude is always correct, also due to the same reason as for Point 4.

Annex B

This annex describes an efficient transformation between SM – and geographic coordinates, working without the need to explicitly set up rotation matrices and also without the transformation into geomagnetic coordinates as an intermediate step. The method is based on the three unit vectors indicating the coordinate axes of the SM system. These vectors must be known in the geographic system as well as in the SM system. To compute them, the following input parameters are needed:

  1. Geographic latitude φM and longitude λM of the northern geomagnetic pole (e.g., v. Biel, 1990).

  2. Geographic direction to the Sun, Annex B1.

Then the required unit vectors are established (whereby in the SM system, the geomagnetic dipole axis coincides with the SM z-axis, Laundal & Richmond, 2016):

Geographic system

εMgg={cosφM·cosλMcosφM·sinλMsinφM}, εYgg=εMgg×εSgg|εMgg×εSgg|, εXgg=εYgg×εMgg,$$ {\overrightarrow{\epsilon }}_{{M}_{{gg}}}=\left\{\begin{array}{c}\mathrm{cos}{\phi }_M\middot \mathrm{cos}{\lambda }_M\\ \mathrm{cos}{\phi }_M\middot \mathrm{sin}{\lambda }_M\\ \mathrm{sin}{\phi }_M\end{array}\right\},\enspace {\overrightarrow{\epsilon }}_{{Y}_{{gg}}}=\frac{{\overrightarrow{\epsilon }}_{{M}_{{gg}}}\times {\overrightarrow{\epsilon }}_{{S}_{{gg}}}}{\left|{\overrightarrow{\epsilon }}_{{M}_{{gg}}}\times {\overrightarrow{\epsilon }}_{{S}_{{gg}}}\right|},\enspace {\overrightarrow{\epsilon }}_{{X}_{{gg}}}={\overrightarrow{\epsilon }}_{{Y}_{{gg}}}\times {\overrightarrow{\epsilon }}_{{M}_{{gg}}}, $$(B.1a)

where

  • εMgg$ {\overrightarrow{\epsilon }}_{{M}_{{gg}}}$ … geocentric unit vector pointing to the geographic position of the northern geomagnetic pole,

  • εSgg$ {\overrightarrow{\epsilon }}_{{S}_{{gg}}}$ … geocentric unit vector pointing into the geographic direction of the Sun,

SM system

εMSM={001}, εYSM={010}, εXSM={100}.$$ {\overrightarrow{\epsilon }}_{{M}_{{SM}}}=\left\{\begin{array}{c}0\\ 0\\ 1\end{array}\right\},\enspace {\overrightarrow{\epsilon }}_{{Y}_{{SM}}}=\left\{\begin{array}{c}0\\ 1\\ 0\end{array}\right\},\enspace {\overrightarrow{\epsilon }}_{{X}_{{SM}}}=\left\{\begin{array}{c}1\\ 0\\ 0\end{array}\right\}. $$(B.1b)

In the SM system, the position vector of the Sun is by the way

εsSM={cosφmS0sinφmS}$ {\overrightarrow{\epsilon }}_{{s}_{{SM}}}=\left\{\begin{array}{c}\mathrm{cos}{\phi }_{{m}_S}\\ 0\\ \mathrm{sin}{\phi }_{{m}_S}\end{array}\right\}$ as per SM definition the Sun’s SM longitude is always zero

where

φmS$ {\phi }_{{m}_S}$ … geomagnetic latitude of the Sun, which can be computed from the geographic unit vectors as

sinφmS=(εSgg·εMgg),cosφmS=+1-sin2φmS.$$ \begin{array}{ccc}\mathrm{sin}{\phi }_{{m}_S}=\left({\overrightarrow{\epsilon }}_{{S}_{{gg}}}\middot {\overrightarrow{\epsilon }}_{{M}_{{gg}}}\right)&,& \mathrm{cos}{\phi }_{{m}_S}=+\end{array}\sqrt{1-{\mathrm{sin}}^2{\phi }_{{m}_S}.} $$

To transform a given position vector R$ \vec{R}$ from one system into the other, its projections onto the unit vectors are computed in the system in which it is given, and then the unit vectors in the other system are scaled by these projections to obtain the components of R$ \vec{R}$ in the other system:

Geographic ➙ SM

RSM=(εXgg·Rgg)·εXSM+(εYgg·Rgg)·εYSM+(εMgg·Rgg)·εMSM.$$ {\vec{R}}_{{SM}}=\left({\overrightarrow{\epsilon }}_{{X}_{{gg}}}\middot {\vec{R}}_{{gg}}\right)\middot {\overrightarrow{\epsilon }}_{{X}_{{SM}}}+\left({\overrightarrow{\epsilon }}_{{Y}_{{gg}}}\middot {\vec{R}}_{{gg}}\right)\middot {\overrightarrow{\epsilon }}_{{Y}_{{SM}}}+\left({\overrightarrow{\epsilon }}_{{M}_{{gg}}}\middot {\vec{R}}_{{gg}}\right)\middot {\overrightarrow{\epsilon }}_{{M}_{{SM}}}. $$(B.2a)

SM ➙ geographic

Rgg=(εXSM·RSM)·εXgg+(εYSM·RSM)·εYgg+(εMSM·RSM)·εMgg.$$ {\vec{R}}_{{gg}}=\left({\overrightarrow{\epsilon }}_{{X}_{{SM}}}\middot {\vec{R}}_{{SM}}\right)\middot {\overrightarrow{\epsilon }}_{{X}_{{gg}}}+\left({\overrightarrow{\epsilon }}_{{Y}_{{SM}}}\middot {\vec{R}}_{{SM}}\right)\middot {\overrightarrow{\epsilon }}_{{Y}_{{gg}}}+\left({\overrightarrow{\epsilon }}_{{M}_{{SM}}}\middot {\vec{R}}_{{SM}}\right)\middot {\overrightarrow{\epsilon }}_{{M}_{{gg}}}. $$(B.2b)

Regarding that the unit vectors in the SM system consist only of Ones and Zeros, equations (B.2a) and (B.2b)) can be simplified as follows:

Geographic ➙ SM

RSM=(εXgg·Rgg)·{100}+(εYgg·Rgg)·{010}+(εMgg·Rgg)·{001}={(εXgg·Rgg)(εYgg·Rgg)(εMgg·Rgg)},$$ {\vec{R}}_{{SM}}=\left({\overrightarrow{\epsilon }}_{{X}_{{gg}}}\middot {\vec{R}}_{{gg}}\right)\middot \left\{\begin{array}{c}1\\ 0\\ 0\end{array}\right\}+\left({\overrightarrow{\epsilon }}_{{Y}_{{gg}}}\middot {\vec{R}}_{{gg}}\right)\middot \left\{\begin{array}{c}0\\ 1\\ 0\end{array}\right\}+\left({\overrightarrow{\epsilon }}_{{M}_{{gg}}}\middot {\vec{R}}_{{gg}}\right)\middot \left\{\begin{array}{c}0\\ 0\\ 1\end{array}\right\}=\left\{\begin{array}{c}\left({\overrightarrow{\epsilon }}_{{X}_{{gg}}}\middot {\vec{R}}_{{gg}}\right)\\ \left({\overrightarrow{\epsilon }}_{{Y}_{{gg}}}\middot {\vec{R}}_{{gg}}\right)\\ \left({\overrightarrow{\epsilon }}_{{M}_{{gg}}}\middot {\vec{R}}_{{gg}}\right)\end{array}\right\}, $$(B.3a)

SM  geographic

Rgg=({100}·{RSMxRSMyRSMz})·εXgg+({010}·{RSMxRSMyRSMz})·εYgg+({001}·{RSMxRSMyRSMz})·εMgg=RSMx·εXgg+RSMy·εYgg+RSMz·εMgg.$$ {\vec{R}}_{{gg}}=\left(\left\{\begin{array}{c}1\\ 0\\ 0\end{array}\right\}\middot \left\{\begin{array}{c}{R}_{{{SM}}_x}\\ {R}_{{{SM}}_y}\\ {R}_{{{SM}}_z}\end{array}\right\}\right)\middot {\overrightarrow{\epsilon }}_{{X}_{{gg}}}+\left(\left\{\begin{array}{c}0\\ 1\\ 0\end{array}\right\}\middot \left\{\begin{array}{c}{R}_{{{SM}}_x}\\ {R}_{{{SM}}_y}\\ {R}_{{{SM}}_z}\end{array}\right\}\right)\middot {\overrightarrow{\epsilon }}_{{Y}_{{gg}}}+\left(\left\{\begin{array}{c}0\\ 0\\ 1\end{array}\right\}\middot \left\{\begin{array}{c}{R}_{{{SM}}_x}\\ {R}_{{{SM}}_y}\\ {R}_{{{SM}}_z}\end{array}\right\}\right)\middot {\overrightarrow{\epsilon }}_{{M}_{{gg}}}={R}_{{{SM}}_x}\middot {\overrightarrow{\epsilon }}_{{X}_{{gg}}}+{R}_{{{SM}}_y}\middot {\overrightarrow{\epsilon }}_{{Y}_{{gg}}}+{R}_{{{SM}}_z}\middot {\overrightarrow{\epsilon }}_{{M}_{{gg}}}. $$(B.3b)

B1 Simplified computation of geographic Sun unit vector

In Annex B1, a very simple way of geographic Sun direction computation is worked out, which is considered accurate enough when applying the NEPPM.

Similar to the unit vector of the northern geomagnetic pole in the geographic system, equation (B.1a), left, the geographic unit vector pointing into the direction of the Sun can be computed if the Sun’s geographic latitude φS and longitude λS are known. To conduct the four computation steps below, the user has to provide the following parameters as input:

Computation steps:

  1. Compute local time, and from that solar time, from the user’s longitude λ and MJD (MJD always refers to Universal Time, i.e., to the local time of the zero meridian; thus, the user’s geographic longitude λ has to be added onto MJD to obtain the user’s local time)

LT   =   mod ( MJD , 1 ) · 24   h + λ 15 °   [ h ours ] . $$ \begin{array}{cc}{LT}\enspace =\enspace {mod}\left({MJD},1\right)\middot 24\enspace \mathrm{h}+\raisebox{1ex}{$\lambda $}\!\left/ \!\raisebox{-1ex}{$15{}^{\circ} $}\right.\enspace & \left[\mathrm{h}\mathrm{ours}\right]\end{array}. $$(B.4)

where

mod(MJD,1)… modulo operation, in this case, returning the day-fractional of MJD.

The local time obtained by equation (B.4) is the mean solar time, referring to the uniform motion of a fictitious Sun. What is needed here is the location of the true Sun, i.e., true solar time. The difference between both solar times can, depending on the season, be up to 17 min and is caused by 1) the Earth’s elliptic motion around the Sun and 2) due to the Earth’s rotation axis tilt w.r.t. the ecliptic plane. This periodically varying offset between true solar time and mean solar time is described by the equation of time (e.g., https://adsabs.harvard.edu/full/1989MNRAS.238.1529H). Davies (1990) provides the following simple formula giving the equation of time (ET) with an accuracy of 2 min

ET=7.75·sin{0.9856·(doy-3.0)}-9.94·sin{1.971·(doy-80.7)} [minutes],$$ \begin{array}{cc}{ET}=7.75\middot \mathrm{sin}\left\{0.9856\middot \left({doy}-3.0\right)\right\}-9.94\middot \mathrm{sin}\left\{1.971\middot \left({doy}-80.7\right)\right\}\enspace & \left[\mathrm{minutes}\right]\end{array}, $$(B.5)

where

  • (360°)/365.25 = 0.9856 … conversion factor from [days] to [degrees],

  • 1.971 = 2·0.9856,

  • doy … day-of-the-year [days], can be computed from MJD$ {MJD}$, should be doy-fractional,

  • 3.0 … doy of 3 January, i.e. perihelion transit,

  • 80.7 … doy of 21/22 March, i.e. vernal equinox.

Tests revealed that the result of equation (B.5) has to be subtracted from the LT value computed with equation (B.4) to get the true solar time ST

ST = LT-ET60 [hours],$$ \begin{array}{cc}{ST}\enspace =\enspace {LT}-\raisebox{1ex}{${ET}$}\!\left/ \!\raisebox{-1ex}{$60$}\right.\enspace & \left[\mathrm{h}\mathrm{ours}\right],\end{array} $$(B.6)

where ET, as provided by equation (B.5), has to be converted from [minutes] to [hours] in equation (B.6).

  1. Compute the Sun’s geographic longitude

λ S   =   λ - ( ST - 12   h ) · 180 ° 12   h   [ degree s ] . $$ \begin{array}{cc}{\lambda }_S\enspace =\enspace \lambda -\left({ST}-12\enspace \mathrm{h}\right)\middot \frac{180{}^{\circ} }{12\enspace \mathrm{h}}\enspace & \left[\mathrm{degree}\mathrm{s}\right]\end{array}. $$(B.7)

  1. In addition, Davies (1990) provides a low-precision formula for the solar declination, by which, for the application here, the Sun’s geographic latitude can be approximated as

δ S = 23 ° . 44 · sin { 0.9856 · ( doy - 80.7 ) }   [ degree s ] , $$ \begin{array}{cc}{\delta }_S=23{}^{\circ}.44\middot \mathrm{sin}\left\{0.9856\middot \left({doy}-80.7\right)\right\}\enspace & \left[\mathrm{degree}\mathrm{s}\right]\end{array}, $$(B.8)

where

  • 23°.44 … obliquity of the ecliptic (low-precision value),

  • (360°)/365.25 = 0.9856 … conversion factor from [days] to [degrees],

  • doy … day-of-the-year [days], can be computed from MJD, should be doy-fractional,

  • 80.7 … doy of 21/22 March, i.e. vernal equinox.

  1. Here δS can be set equal to the Sun’s geographic latitude φS, and the Sun’s low-precision geographic longitude λS was computed with equation (B.7). Therefore, the geographic direction to the Sun can, for the application here, be described by the following unit vector with an accuracy of about 1 degree

ε S gg = { cos δ S · cos λ S cos δ S · sin λ S sin δ S }   with   φ S δ S . $$ {\overrightarrow{\epsilon }}_{{S}_{{gg}}}=\left\{\begin{array}{c}\mathrm{cos}{\delta }_S\middot \mathrm{cos}{\lambda }_S\\ \mathrm{cos}{\delta }_S\middot \mathrm{sin}{\lambda }_S\\ \mathrm{sin}{\delta }_S\end{array}\right\}\enspace \mathrm{with}\enspace {\phi }_S\approx {\delta }_S. $$(B.9)

References

Cite this article as: Banyś D, Feltens J, Jakowski N, Berdermann J, Hoque MM, et al. 2024. A new model for plasmapause locations derived from IMAGE RPI and Van Allen Probes data. J. Space Weather Space Clim. 14, 36. https://doi.org/10.1051/swsc/2024016.

All Tables

Table 1

Dst levels and their ranges used for ellipse parameters fitting.

Table 2

Empirical functions f(Dst) fitted to represent the three ellipse parameters: semi-major axis a, eccentricity e, and M denoting the MLT of the plasmapause bulge as obtained by reformulating equation (3.4).

All Figures

thumbnail 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 ̅$ \overline{\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).

In the text
thumbnail Figure 2

Meaning of the orientation angle ϖ, the semi-major axis a, Lppmax$ {L}_{{{pp}}_{{max}}}$ and Lppmin$ {L}_{{{pp}}_{{min}}}$.

In the text
thumbnail 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.

In the text
thumbnail 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).

In the text
thumbnail 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).

In the text
thumbnail 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).

In the text
thumbnail Figure 7

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

In the text
thumbnail Figure A1

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

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.