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 |
Research Article
A new model for plasmapause locations derived from IMAGE RPI and Van Allen Probes data
1
German Aerospace Center (DLR), Institute for Solar-Terrestrial Physics, Kalkhorstweg 53, 17235 Neustrelitz, Germany
2
Telespazio Germany GmbH c/o European Space Operations Centre (ESA/ESOC), Robert-Bosch-Str. 5, 64293 Darmstadt, Germany (Retired)
3
European Space Operations Centre (ESA/ESOC), Robert-Bosch-Str. 5, 64293 Darmstadt, Germany
* Corresponding author: Daniela.Banys@dlr.de
Received:
20
August
2023
Accepted:
18
May
2024
The outer boundary of the plasmasphere, the plasmapause, is characterized by a sharp electron density gradient that changes under varying space weather conditions. We developed a new model, called the Neustrelitz/ESOC PlasmaPause Model (NEPPM), for providing plasmapause location in terms of L-shell utilizing electron density measurements from the Van Allen Probes from 2012 to 2018 and the IMAGE satellite data from 2001 to 2005. Both datasets were preprocessed, and algorithms were developed for the automatic detection of plasmapause location Lpp where L denotes the McIlwain parameter. The suggested model provides a simple ellipse-based approach determined by the semi-major axis, the eccentricity, and the orientation angle of the semi-major axis. The modelled Lpp varies as a function of the Dst index and magnetic local time MLT. The NEPPM results are compared with the Global Core Plasma Model (GCPM). The plasmapause bulge in the evening hours follows the level of geomagnetic activity. The NEPPM will complete the NPSM (Neustrelitz PlasmaSphere Model), which was derived from dual-frequency GPS measurements onboard the CHAMP satellite mission.
Key words: Plasmapause / Plasmasphere / Ionosphere / Electron density / Empirical model
© D. Banyś et al., Published by EDP Sciences 2024
This 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 . 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 are generally replaced by smoothed measurements via a 3-point moving average approach
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
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
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 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.
where 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
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.,
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, , typically points to the evening-afternoon hours and the Earth-nearest point, , is located in the opposite direction. The ellipse radii, among them and , are equal to the distances from the focal point coinciding with the Earth to the ellipse periphery, in L-units.
Figure 2 Meaning of the orientation angle ϖ, the semi-major axis a, and . |
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)
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
where
x … magnetic local time expressed in radians, i.e. ,
ϖ … 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):
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
where
M… MLT of bulge.
At low to moderate solar activities, the bulge is typically located in the evening hours while moving more and more into the early afternoon hours when solar activity increases. However, the bulge will never appear in the morning and forenoon hours, i.e. M > 12 h, and thus ϖ is always negative, with decreasing magnitude at increasing solar activity.
Different Dst levels were defined, and then for each Dst level, a set of ellipse parameters (a,e,ϖ) was fitted to the plasmapause values obtained from the plasmapause detection, Section 2.3, whereby ϖ=ϖ(MLT). The selected Dst levels were, in Table 1:
Dst levels and their ranges used for ellipse parameters fitting.
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 . 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:
the magnitude of the geomagnetic field vector ,
the invariant .
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, 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, 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 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
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 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 , 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
where
-
… 3D Cartesian vector indicating a geocentric plasmaspheric position in the SM,
-
L … L-shell value at plasmaspheric position,
-
φm … geomagnetic latitude,
-
Λ … solar longitude, i.e., longitude in the plane of the geomagnetic equator relative to the longitude of the Sun.
When plotting 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 pointing to the position in space, and the other two are scalar quantities. When describing the plasmapause position , 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
an ansatz that is developed in detail in Annex A. Principally, a 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 and related epoch. However, the 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 , three component-wise observation equations are set up from equation (3.6)
Expressed in matrix form, a set of 3n observation equations can be established for n observed vectors based on equations (3.8),
where for the non-linear fit
-
,
-
… plasmapause position vector detected from the IMAGE & Van Allen Probes data, Sect. 2.3,
-
LppMod … Lpp-value computed with equation (3.3) using some initial values for a, e, ϖ,
-
cos2φm … computed from vector, equation (3.7),
-
… unit vector of ,
-
v … fitted correction to the observable,
-
… partial of Lpp w.r.t. semi-major axis a, evaluated with initial values,
-
∆a … estimated correction to initial value of a,
-
… partial of Lpp w.r.t. eccentricity e, evaluated with initial values,
-
∆e … estimated correction to initial value of e,
-
… 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 is requested. Then the computation is conducted in the following steps (see also Annex B):
Set transformation parameters geographic ➧ SM.
Compute geographic position unit vector .
Transform from geographic into SM ➧ unit vector .
Compute , with .
Compute cos2φm term for dipole from the components of via
-
The fitted Lpp function is evaluated with equation (3.3).
-
The dipole position is computed with equation (3.6), where term (3) is here .
-
With the IGRF (or another non-dipole model) a vector is computed at the position indicated by dipole position vector and for the requested epoch T, transform into the SM system.
-
Compute from the non-dipole vector components an improved cos2φm term, equation (3.7).
-
Compute non-dipole vector, again by using equation (3.6), but now with the non-dipole cos2φm term obtained from Step 9).
-
Transform non-dipole 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
where
-
Dstmed … median Dst of considered Dst level represented by a certain Dst range, Table 1,
-
Dscal … Dst scaling factor,
-
A … amplitude to be fitted,
-
C … boundary value fixed for extreme storms by assuming Dstmed → −∞.
It turned out, that exponential functions are best suited to represent the ellipse parameters in dependency of Dst. This approach has the advantage that besides the fitting using the available data set, extreme storms that have never been observed can be considered in a plausible way. It is evident that the boundary values fixed for all three fits are not approved by observations. We have defined these values following plausible indications of observations and physical arguments. Consequently, the given boundary values might be slightly modified using an extended database in future studies. Because such a modification has a negligible impact on the model approach for weak and moderate storms (i.e., Dst > −100 nT) the model is applicable in common operational services. Besides Dscal and strength A obtained by least squares fitting, also the estimated boundary values C are listed in Table 2. For the semi-major axis, the eccentricity, and the orientation angle of the ellipse.
While fitting the strength A of the f(Dst) function, the input Dst values had to be scaled by an empirical Dscal quotient. This procedure converted the range of Dst values into a reasonable argument range for the exponential function. Similar to the boundary values for extreme storms, 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.
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.
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 and a ≈ 1.9, respectively.
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.
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 , 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
and
Thus, the difference between both extreme plasmapause positions is given by
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 and 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 . 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.
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.
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.
Figure 6 NEPPM (black) and GCPM (grey) for each Dst level, cf. Table 1. The extracted Lpp values are marked by dots in red (IMAGE), and blue (RBSP). |
Figure 7 displays the specific deviations of the two models from the IMAGE and Van Allen Probes data with respect to Dst and MLT. NEPPM displays the highest RMS in the afternoon for Dst < −50 nT. GCPM shows two higher RMS values, one in the afternoon for very low Dst values and another one in the evening hours for Dst > −40 nT. Furthermore, the overall dayside shows enhanced RMS values for GCPM, which are not seen for NEPPM. Since we used IMAGE and Van Allen Probes data to build the NEPPM, our model performs well. Certainly, further data is necessary to validate our approach. As we can see, GCPM does not represent the recent data properly, though it was a good and, most importantly, consistent global approach back then.
Figure 7 RMS (top) and relative RMS (bottom) for NEPPM (left) and GCPM (right). |
Meanwhile, much more detailed data has become available for improved accuracy and higher resolution. Various plasmapause models were proposed, and especially neural networks (NN) are becoming of higher interest as proposed by Guo et al. (2021). A comparison of the RMS with this NN model shows similar results. Both NEPPM and NN models show the highest RMS values of about 0.9 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 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 that obey the representation equation (3.6). When computing them with a constant Lpp value in an equidistant grid, a plot of these would display a torus-like picture, blue in Figure A1 below, corresponding to the dipole assumption inherent in the definition of L-shells.
Figure A1 The vector fields (blue) and (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
whereas it must be linked to the geomagnetic field vector if, like in the IGRF, a non-dipole field is the base. For this, has to be computed with the selected geomagnetic field model at the position indicated by 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 are once more derived from , i.e. for a dipole. Nevertheless, as will be shown below, the components from a non-dipole 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 vector at that point in the magnetic field.
First, reformulate equation (3.6) to
and determine the first derivative of w.r.t. the geomagnetic latitude. By this, one has the tangent vector
which always points tangentially northwards along the dipole field line, i.e., into the same direction as . The geometry behind and 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 , analogously to spherical latitude computation from Cartesian coordinates, one can build the ratio
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
Thus, linking the dipole and φm, we get
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
This can be solved according to the p-q-formula with
Leading to
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. , 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 have to be expressed as a function of the components {Rx Ry Rz} of the plasmapause position vector . Therefore in equation (A.3) the cosine and sine terms are formulated as functions of the plasmapause position vector components
After having substituted equations (A.10) into equation (A.3), one obtains after some simple algebra
Whereby for describing purely the direction of 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 is computed for a given geomagnetic latitude φm and solar longitude Λ according to equation (A.2). Reversely, φm can be computed from the Cartesian components in a backwards transformation
The proof starts with equation (A.9), in which, due to equation (A.4), the vector components can be replaced by the tangent vector components. This approach will, after some algebra, allow to express equation (A.9) in terms of the 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 components {Tx Ty Tz}
In an intermediate step, the square root term of equation (A.13) is reformulated
Substituting equation (A.14) into (A.13)
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 of equation (A.11) is an overall scaling factor which cancels when establishing a ratio between components (similar to equation (A.4)), as given in equation (A.15). Therefore, it is sufficient to write here
Equation (A.15) needs the squares of the equation (A.16) terms
Thus is (numerator square root term of equation (A.15)):
Compare with equation (A.17) for and
Substituting Tz from equation (A.16) and equations (A.17) and (A.18) into equation (A.15)
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
Negative square root
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
A2 Meaning of the negative square root solution
For completeness, also the plasmapause position vector , corresponding to the negative square root of equation (A.9), shall be explained. It can be obtained by putting , equation (A.20b), into equation (A.2). However, it can also be expressed in terms of the plasmapause position vector . Taking from equation (A.20b) and regarding the following trigonometric relations (e.g. Sigl, 1977)
one gets by substituting the expression (A.20b) for
And since and 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 . Thus, by substituting equations (A.22) and (A.10) into equation (3.6) one obtains, after some algebra, as a function of the components
Thereby term (1) in equation (A.23) gives the magnitude, , compare with equation (A.22), bottom, and term (2) is the unit vector of , . 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 .
By the relation and equation (A.1), also Lpp can be expressed as a function of the components
Thus with equation (A.24) can be fully formulated as a function of the 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
And in an analogous exercise also a formula for the angle enclosed by and can be established, by using the unit vectors and ,
By an appropriate exchange of signs of the and 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 vs. 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:
The angle ϑ enclosed by and ranges between 70°.5 at and 90° at φm = 0°, ±90°.
φm = ±90° is a singularity point. Here the dipole curve passes through the geocentre and the plasmapause position vectors vanish.
when is pointing to the equator (Rz = 0). Thus φm = 0° has been skipped in computation.
If has a negative z-component Rz, the resulting , computed with equation (A.25), is pointing into the opposite dipole loop. To get in this case 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).
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.
When equation (A.20b) is solved with the atan function instead of atan2, the obtained 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:
Geographic latitude φM and longitude λM of the northern geomagnetic pole (e.g., v. Biel, 1990).
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
where
-
… geocentric unit vector pointing to the geographic position of the northern geomagnetic pole,
-
… geocentric unit vector pointing into the geographic direction of the Sun,
SM system
In the SM system, the position vector of the Sun is by the way
as per SM definition the Sun’s SM longitude is always zero
where
… geomagnetic latitude of the Sun, which can be computed from the geographic unit vectors as
To transform a given position vector 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 in the other system:
Geographic ➙ SM
SM ➙ geographic
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
SM geographic
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:
His/her geographic longitude λ
Modified Julian Date (MJD) epoch, incl. day-fractional, for which the Sun direction is requested (e.g., http://www.csgnetwork.com/julianmodifdateconv.html)
Computation steps:
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)
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
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 , 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
where ET, as provided by equation (B.5), has to be converted from [minutes] to [hours] in equation (B.6).
Compute the Sun’s geographic longitude
-
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
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.
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
References
- v. Biel HA. 1990. The geomagnetic time and position of a terrestrial station. J Atm Terr Phys 52(9): 687–694. https://doi.org/10.1016/0021-9169(90)90001-4. [CrossRef] [Google Scholar]
- Botek E, Pierrard V, Darrouzet F. 2021. Assessment of the Earth’s cold plasmatrough modeling by using Van Allen Probes/EMFISIS and Arase/PWE electron density data. J Geophys Res Space Phys 126: e2021JA029737. https://doi.org/10.1029/2021JA029737. [CrossRef] [Google Scholar]
- Carpenter DL. 1963. Whistler evidence of a ‘knee’ in the magnetospheric ionization density profile. J Geophys Res 68(6): 1675–1682. https://doi.org/10.1029/JZ068i006p01675. [CrossRef] [Google Scholar]
- Carpenter DL, Anderson RR. 1992. An ISEE/whistler model of equatorial electron density in the magnetosphere. J Geophys Res 97: 1097–1108. https://doi.org/10.1029/91JA01548. [CrossRef] [Google Scholar]
- Carpenter DL, Smith RL. 1964. Whistler measurements of electron density in the magnetosphere. Rev Geophys 2(3): 415–441. https://doi.org/10.1029/RG002i003p00415. [CrossRef] [Google Scholar]
- Cliver EW, Svalgaard L. 2004. The 1859 solar-terrestrial disturbance and the current limits of extreme space weather activity. Sol Phys 224: 407–422. https://doi.org/10.1007/s11207-005-4980-z. [CrossRef] [Google Scholar]
- Darrouzet F, De Keyser J, Décréau PME, El Lemdani-Mazouz F, Vallières X. 2008. Statistical analysis of plasmaspheric plumes with Cluster/WHISPER observations. Ann Geophys 26(8): 2403–2417. https://doi.org/10.5194/angeo-26-2403-2008. [CrossRef] [Google Scholar]
- Darrouzet F, De Keyser J, Pierrard V. 2009. The Earth’s plasmasphere: a CLUSTER and IMAGE perspective. Springer. p. 296. ISBN: 978-1-4419-1322-7. [Google Scholar]
- Davies K. 1990. Ionospheric radio. Peter Peregrinus Ltd., London, United Kingdom. ISBN: 978-0-86341-186-1. [CrossRef] [Google Scholar]
- Denton RE, Wang Y, Webb PA, Tengdin PM, Goldstein J, Redfern JA, Reinisch BW. 2012. Magnetospheric electron density long-term (>1 day) refilling rates inferred from passive radio emissions measured by IMAGE RPI during geomagnetically quiet times. J Geophys Res 117: A3. http://dx.doi.org/10.1029/2011JA017274. [Google Scholar]
- Escobal PR. 1965. Methods of orbit determination. In: Robert E. Krieger Publishing Company, Malabar, Florida, 2nd Edition, John Wiley & Sons, Inc., Reprinted by Arrangement, ISBN-10: 0882753193/ISBN-13: 978-0882753195. [Google Scholar]
- Gallagher DL, Craven PD, Comfort RH. 1988. An empirical model of the Earth’s plasmasphere. Adv Space Res 8(8): 15–24. https://doi.org/10.1016/0273-1177(88)90258-X. [CrossRef] [Google Scholar]
- Gallagher DL, Craven PD, Comfort RH. 2000. Global core plasma model. J Geophys Res 105(A8): 18819–18833. https://doi.org/10.1029/1999JA000241. [CrossRef] [Google Scholar]
- Galkin IA, Reinisch BW, Grinstein G, Khmyrov G, Kozlov A, Huang X, Fung SF. 2004. Automated exploration of the radio plasma imager data. J Geophys Res 109: A12. http://dx.doi.org:10.1029/2004JA010439. [Google Scholar]
- Gerzen T, Feltens J, Jakowski N, Galkin I, Denton R, Reinisch BW, Zandbergen R. 2014. Validation of plasmasphere electron density reconstructions derived from data on board CHAMP by IMAGE/RPI data. Adv Space Res 55(1): 170–183. http://dx.doi.org/10.1016/j.asr.2014.08.005. [Google Scholar]
- Goldstein J, Spasojević M, Reiff PH, Sandel BR, Forrester WT, Gallagher DL, Reinisch BW. 2003. Identifying the plasmapause in IMAGE EUV data using IMAGE RPI in situ steep density gradients. J Geophys Res 108(A4): 1147. http://dx.doi.org/10.1029/2002JA009475. [Google Scholar]
- Guo D, Fu S, Xiang Z, Ni B, Guo Y, Feng M, et al. 2021. Prediction of dynamic plasmapause location using a neural network. Space Weather 19(5): 1–13. https://doi.org/10.1029/2020SW002622. [Google Scholar]
- Hayakawa H, Ebihara Y, Willis DM, Toriumi S, Iju T, Hattori K, et al. 2019. Temporal and spatial evolutions of a large sunspot group and great auroral storms around the Carrington event in 1859. Space Weather 17(11): 1553–1569. https://doi.org/10.1029/2019SW002269. [CrossRef] [Google Scholar]
- He F, Zhang XX, Lin RL, Fok MC, Katus RM, Liemohn MW, et al. 2017. A new solar wind-driven global dynamic plasmapause model: 2. Model and validation. J Geophys Res Space Phys 122(7): 7172–7187. https://doi.org/10.1002/2017JA023913. [CrossRef] [Google Scholar]
- Heilig B, Lühr H. 2013. New plasmapause model derived from CHAMP field-aligned current signatures. Ann Geophys 31(3): 529–539. https://doi.org/10.5194/angeo-31-529-2013. [CrossRef] [Google Scholar]
- Hoque MM, Jakowski N, Prol FS. 2022. A new climatological electron density model for supporting space weather services. J Space Weather Space Clim 12: 1. https://doi.org/10.1051/swsc/2021044. [CrossRef] [EDP Sciences] [Google Scholar]
- Huba JD, Krall J. 2013. Modeling the plasmasphere with SAMI3. Geophys Res Lett 40: 6–10. https://doi.org/10.1029/2012GL054300. [CrossRef] [Google Scholar]
- Jakowski N, Hoque MM. 2018. A new electron density model of the plasmasphere for operational applications and services. J Space Weather Space Clim 8: A16. https://doi.org/10.1051/swsc/2018002. [CrossRef] [EDP Sciences] [Google Scholar]
- Kimball DS. 1960. A study of the aurora of 1859, Sci. Rpt. 6, UAG-R109, Univ. of Alaska, Fairbanks, Alaska. [Google Scholar]
- Kletzing CA, Kurth WS, Acuna M, MacDowall RJ, Torbert RB, et al. 2013. The electric and magnetic fields instrument suite and integrated science (EMFISIS) on RBSP. Space Sci Rev 179: 127–181. https://doi.org/10.1007/s11214-013-9993-6. [CrossRef] [Google Scholar]
- Laundal KM, Richmond AD. 2016. Magnetic coordinate systems. Space Sci Rev 206: 27–59. https://doi.org/10.1007/s11214-016-0275-y. [Google Scholar]
- Lemaire JF, Gringauz KI. 1998. The earth’s plasmasphere. Cambridge University Press. ISBN 769 0521430917. https://doi.org/10.1017/CBO9780511600098. [CrossRef] [Google Scholar]
- Liu X, Liu W, Cao JB, Fu HS, Yu J, Li X. 2015. Dynamic plasmapause model based on THEMIS measurements. J Geophys Res Space Phys 120(12): 10543–10556. https://doi.org/10.1002/2015JA021801. [Google Scholar]
- Lunt N, Kersley L, Bishop GJ, Mazzella AJ Jr. 1999. The contribution of the protonosphere to GPS total electron content: Experimental measurements. Radio Sci 34(5): 1273–1280. https://doi.org/10.1029/1999RS900016. [CrossRef] [Google Scholar]
- Matzka J, Bronkalla O, Tornow K, Elger K, Stolle C. 2021. Geomagnetic Kp index. V. 1.0. GFZ Data Services. https://doi.org/10.5880/Kp.0001. [Google Scholar]
- McIlwain CE. 1961. Coordinates for mapping the distribution of magnetically trapped particles. J Geophys Res 66(11): 3681–3691. https://doi.org/10.1029/JZ066i011p03681. [CrossRef] [Google Scholar]
- McIlwain CE. 1965. Magnetic coordinates. Space Sci Rev 5: 585–598. https://doi.org/10.1007/BF00167327. [Google Scholar]
- Moldwin MB, Downward L, Rassoul HK, Amin R, Anderson RR. 2002. A new model of the location of the plasmapause: CRRES results. J Geophys Res 107(A11): 1339. https://doi.org/10.1029/2001JA009211. [Google Scholar]
- Nose M, Iyemori T, Sugiura M, Kamei T. 2015. Geomagnetic Dst index, World Data Center for Geomagnetism, Kyoto. https://doi.org/10.17593/14515-74000. [Google Scholar]
- Obana Y, Maruyama N, Shinbori A, Hashimoto KK, Fedrizzi M, Nosé M, et al. 2019. Response of the ionosphere-plasmasphere coupling to the September 2017 storm: what erodes the plasmasphere so severely? Space Weather 17: 861–876. https://doi.org/10.1029/2019SW002168. [CrossRef] [Google Scholar]
- O’Brien BJ, Laughlin CD, Allen JAV, Frank LA. 1962. Measurements of the intensity and spectrum of electrons at 1000-kilometer altitude and high latitudes. J Geophys Res 67(4): 1209–1225. https://doi.org/10.1029/JZ067i004p01209. [CrossRef] [Google Scholar]
- O’Brien TP, Moldwin MB. 2003. Empirical plasmapause models from magnetic indices. Geophys Res Lett 30(4): 1152. https://doi.org/10.1029/2002GL016007. [Google Scholar]
- Park W, Lee J, Kim KC, Lee JK, Park K, Miyashita Y, et al. 2021. Operational Dst index prediction model based on combination of artificial neural network and empirical model. J Space Weather Space Clim 11: 38. https://doi.org/10.1051/swsc/2021021. [CrossRef] [EDP Sciences] [Google Scholar]
- Pierrard V, Goldstein J, André N, Jordanova VK, Kotova GA, Lemaire JF, Liemohn MW, Matsui H. 2009. Recent progress in physics-based models of the plasmasphere. Space Sci Rev 145(1–2): 193–229. https://doi.org/10.1007/s11214-008-9480-7. [CrossRef] [Google Scholar]
- Reinisch BW, Huang X, Haines DM, Galkin IA, Green JL, Benson RF, Fung SF, Taylor WWL, Reiff PH, Gallagher DL, Bougeret JL, Manning R, Carpenter DL, Boardsen SA. 2001a. First results from the radio plasma imager on IMAGE. Geophys Res Lett 28(6): 1167–1170. https://doi.org/10.1029/2000GL012398. [CrossRef] [Google Scholar]
- Reinisch BW, Huang X, Song P, Sales GS, Fung SF, Green JL, Gallagher DL, Vasyliunas VM. 2001b. Plasma density distribution along the magnetospheric field: RPI observations from IMAGE. Geophys Res Lett 28(24): 4521–4524. https://doi.org/10.1029/2001GL013684. [CrossRef] [Google Scholar]
- Reinisch BW, Mark W, Moldwin B, Denton RE, Gallagher DL, Matsui H, Pierrard V, Tu J. 2009. Augmented empirical models of plasmaspheric density and electric field using IMAGE and CLUSTER data. Space Sci Rev 145: 231–261. http://dx.doi.org/10.1007/s11214-008-9481-6. [CrossRef] [Google Scholar]
- Ripoll JF, Thaller SA, Hartley DP, Cunningham GS, Pierrard V, Kurth WS, et al. 2022. Statistics and empirical models of the plasmasphere boundaries from the Van Allen Probes for radiation belt physics. Geophys Res Lett 49(21), e2022GL101402. https://doi.org/10.1029/2022GL101402. [CrossRef] [Google Scholar]
- Sigl R. 1977. Ebene und Sphärische Trigonometrie. Herbert Wichmann Verlag, Karlsruhe, ISBN 10: 3400000620/ISBN-13: 9783400000627. [Google Scholar]
- Tsurutani BT, Gonzalez WD, Lakhina GS, Alex S. 2003. The extreme magnetic storm of 1–2 September 1859. J Geophys Res 108: A7. https://doi.org/10.1029/2002JA009504. [Google Scholar]
- Yizengaw E, Moldwin MB, Galvan D, Iijima BA, Komjathy A, Mannucci AJ. 2008. Global plasmaspheric TEC and its relative contribution to GPS TEC. J Atmos Sol-Terr Phys 70(11–12): 1541–1548. https://doi.org/10.1016/j.jastp.2008.04.022. [CrossRef] [Google Scholar]
- Zhelavskaya IS, Spasojevic M, Shprits YY, Kurth WS. 2016. Automated determination of electron density from electric field measurements on the Van Allen Probes spacecraft. J Geophys Res Space Phys 121(5): 4611–4625. http://doi.org/10.1002/2015JA022132. [CrossRef] [Google Scholar]
- Zhelavskaya IS, Shprits YY, Spasojevic M. 2017. Empirical modeling of the plasmasphere dynamics using neural networks. J Geophys Res Space Phys 122(11): 227–244. https://doi.org/10.1002/2017JA024406. [CrossRef] [Google Scholar]
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
All Figures
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 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 |
Figure 2 Meaning of the orientation angle ϖ, the semi-major axis a, and . |
|
In the text |
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 |
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 |
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 |
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 |
Figure 7 RMS (top) and relative RMS (bottom) for NEPPM (left) and GCPM (right). |
|
In the text |
Figure A1 The vector fields (blue) and (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.