Width of plasmaspheric plumes related to the level of geomagnetic storm intensity

The plume is a plasma region in the magnetosphere that is detached from the main plasmasphere. It significantly contributes to the dynamic processes in both the inner and outer magnetosphere. In this paper, using Van Allen Probe A (VAP-A), the correlation between plume width and the level of geomagnetic storm intensity is studied. First, through the statistical analysis of all potential plume events, we find that there is almost no correlation between plume width and the level of geomagnetic storm intensity. However, for the plumes in the recovery phase after improved sifting, it seems that there is a negative correlation between the plume width and the absolute value of minimum Dst during a storm. Utilizing test particle simulations, we study the dynamic evolution patterns of plumes during two geomagnetic storms. The simulated structures of the two plasmaspheric plumes are roughly consistent with the structures observed by the Van Allen Probe A. This result suggests that the plasmaspheric particles escape quickly during intense geomagnetic storms, causing the width of the plume to be relatively narrow during the recovery phase of intense geomagnetic storms. These results are helpful for understanding the dynamic evolution of the plasmasphere and plume during geomagnetic storms.

1 Introduction

The plasmasphere is a region of high-density cold particles (at several electron volts) in the inner magnetosphere. The motions of the outer plasmaspheric particles are periodically driven by geomagnetic activity. During geomagnetic storms, the interplanetary magnetic field (IMF) moves southward and leads to geomagnetic reconnection, which subsequently drives the convection electric field (Dungey, 1961). Then, plasmaspheric particles move along the E×B-drift paths in the electric field of the inner magnetosphere and escape from the plasmasphere. The process is known as plasmaspheric erosion. It will force the plasma to extend sunward and produce plasmaspheric plumes that rotate around the Earth during geomagnetic storm intervals (Lakhina et al., 2000).

Previous studies have indicated that the drift paths of plume plasma are not restricted to the innermost magnetosphere (Spasojevic et al., 2005; Spasojevic and Ivan, 2010). This means that the plasmaspheric plume is an important channel for the exchange of mass and energy between the inner magnetosphere and outer magnetosphere (Lakhina et al., 2000). Furthermore, although electromagnetic ion cyclotron (EMIC) waves are not preferentially observed in high-density plumes (Usanova et al., 2013; Grison et al., 2018), the plume may be correlated with the excitation of EMIC waves (Grison et al., 2018; Yu et al., 2016; Yuan et al., 2010), and whistler-mode hiss emissions often exist in plasmaspheric plumes (Su et al., 2018; Kim and Shprits, 2019; Ma et al., 2021). Therefore, understanding the evolution of plumes is essential. When the levels of geomagnetic storm intensity increases, plasmaspheric erosion becomes stronger with the enhancement of the convection electric field (Chen and Grebowsky, 1974; Grebowsky, 1970). However, relatively less research has been conducted regarding the shapes of plasmaspheric plumes. Borovsky and Denton (2008) statistically calculated the linear relationship between the width of the plasmaspheric plume and the levels of geomagnetic storm intensity. Borovsky and Denton (2008) suggested that the linear correlation coefficient between them was almost 0. Since the plasmasphere can be eroded by the enhanced convection electric field during geomagnetic storms (Krall et al., 2017), the enhanced convection field causes low energy plasma drainage to the magnetopause (Denton et al., 2005). We consider that the level of storm intensity may affect the width of the plume in some conditions.

In this paper, we utilized the data recorded by Van Allen Probe A (VAP-A; from 2013 to 2018) to identify plasmaspheric plumes. The correlation coefficient between the width of plasmaspheric plumes and the levels of geomagnetic storm intensity was calculated under different standards. Furthermore, we ran group test particle simulations to support the statistical results.

2 Data and statistical widths of plasmaspheric plumes

The Van Allen Probe A (VAP-A) spacecraft was in a highly elliptical (1.1×5.8 R_E), low-inclination (10^∘) orbit (Mauk et al., 2012) and collected data from August 2012 to October 2019. We use the electron density data provided by the VAP-A (listed in Level-4 data sets), which were calculated from the upper-hybrid resonance frequency (Kurth et al., 2015). In this study, by analyzing the electron density measurements of VAP-A, the structure of plasmaspheric plumes is determined by the two following criteria: (1) the location of the plasmapause is considered to be the position where the electron density decreases to less than 0.2 times within the 0.5 R_E. (2) According to the method commonly used in previous studies, if VAP-A is located outside the location of the plasmapause and the detected density exceeds the result of the plasmaspheric density model given by Sheeley et al. (2001) (the formula is shown below) for more than 10 min, we consider it a plasmaspheric plume (Moldwin et al., 2002, 2004; Zhang et al., 2019):

$$\begin{array}{}\text{(1)}& n_{\mathrm{e}}=\mathrm{1390}{\left({\displaystyle \frac{\mathrm{3}}{L}}\right)}^{\mathrm{4.83}}-\mathrm{240}{\left({\displaystyle \frac{\mathrm{3}}{L}}\right)}^{\mathrm{3.60}}.\end{array}$$

A typical example of a plume is exhibited in Fig. 1. The blue curve indicates the observed electron density from 08:00 to 16:36 UT on 19 October 2016. The red curve indicates the empirical electron density obtained from the plasmaspheric model published in Sheeley et al. (2001). The plasmapause positions are marked by vertical black lines based on the criterion above. In addition, a typical plasmaspheric plume is identified and marked by a gray shadow in Fig. 1.

Figure 1 A typical example of a plume on 19 October 2016. The blue curve represents the detected plasma density, and the red curve displays the density calculated by Sheeley et al. (2001). The position of the plasmaspause is marked by black vertical lines. The gray shadow indicates the plasmaspheric plume detected by VAP-A.

As Li et al. (2022) suggested, the number of plume events in the initial and main phases of geomagnetic storms is very small, and most plumes mainly form in the recovery phase. Partial plume events can still be observed when the geomagnetic activity recovers to quiet conditions after geomagnetic storms; this is because a relatively long time is required for the plasmasphere to recover to normal levels. Besides, referring to Halford et al. (2010) and Wang et al. (2016), the onset time of storm is defined as the time when the Dst index slope becomes negative and remains relatively negative till the minimum of the Dst index in our study. In this paper, we first focus on the relationship between the width of the plume in the interval of 10 d after the storm minimum Dst and the corresponding level of geomagnetic storm intensity (represented by the minimum Dst). According to the plume determination criteria above, we find 423 orbits with plumes (within 10 d after the minimum Dst of storm) by searching 4030 VAP-A orbits from 2013 to 2018. Since several plumes may be identified in some orbits, a total number of 586 plume events are found.

The concept of the observed plume width only has a “fuzzy” definition in the literature. To ensure the accuracy of the analysis, two judgments are adopted to represent the detected plume width. First, the width of the detected plume is defined as the Cartesian distance (ΔP) between the two ends (the entrance and exit) of the VAP-A orbit where the plume is detected:

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }P=\surd \left[(X_{\mathrm{0}}-X_{\mathrm{1}}{)}^{\mathrm{2}}+(Y_{\mathrm{0}}-Y_{\mathrm{1}}{)}^{\mathrm{2}}+(Z_{\mathrm{0}}-Z_{\mathrm{1}}{)}^{\mathrm{2}}\right],\end{array}$$

where X₀, Y₀, and Z₀ represent the Cartesian position of plume entry edge, and X₁, Y₁, and Z₁ represent the Cartesian position of plume exit edge.

Second, the width is considered the difference between the magnetic local times (MLTs) of the two plume ends (ΔMLT). In this study, the corresponding geomagnetic index of the plasmaspheric plume is considered to be the minimum Dst value in the geomagnetic storm. In addition, a complication in measuring the plume width is that the plume is still rotating when VAP-A passes through the plasmasphere, which will lead to more or fewer satellite orbits in the plasmasphere. For statistical significance, the average width can reduce this influence because of the similar behavior between the entry and exit edges of the plume, and it also reminds the reader that the width measurement of the individual plume may be some error (Borovsky and Denton, 2008). To clearly reflect the statistical variation in plume width associated with geomagnetic activity, we use the averaged width of the detected plume in steps of 5 nT (Dst index) to represent the plume width in the corresponding Dst range in the study.

Figure 2 shows the widths of the 586 plumes detected above as a function of the corresponding Dst. The Cartesian distance (ΔP) and ΔMLTs are denoted by the solid black points in Fig. 2a and b, respectively. The red curves connect the averaged widths of plasmaspheric plumes in each step of the 5 nT range (plotted by red asterisks). Then, we fit the red asterisk dots through a linear function, which is drawn by a blue line. It is obvious that the blue curves are almost parallel to the Dst indices for both Fig. 2a and b. This means that the plume width is independent of the level of geomagnetic storm intensity. For Fig. 2a, the calculated Pearson correlation coefficient (marked as R), which indicates the linear relevance between the averaged Cartesian distance of plasmaspheric plumes (indicated by red asterisks) and the corresponding Dst value, is only −0.017. The P value (marked as P), which is generally adopted to express the reliability of their linear correlation, reaches 0.935. For Fig. 2b, the calculated Pearson correlation coefficient between the averaged ΔMLT of plasmaspheric plumes and the corresponding Dst value is only −0.053. The P value is 0.803. The low R and high P values in Fig. 2a and b suggest that there is almost no correlation between the widths of the above plumes and the corresponding Dst values, and their relationship may be very complicated.

Figure 2 (a) The widths of plumes as a function of Dst values are represented by solid black dots, where the width expresses the Cartesian distance between the two ends (the entrance and exit) of the VAP-A orbit where the plume is detected. The linear fitting of the plasmaspheric plume's averaged widths in each step of the 5 nT range (red asterisk points) is indicated by the blue line. (b) The format is similar to (a); however, the ΔMLT of the plasmaspheric plume is adopted to represent the width of the plume.

The formula of the Pearson correlation coefficient is as follows:

$$\begin{array}{}\text{(3)}& R={\displaystyle \frac{N{\sum }_i{x}_i{y}_i-{\sum }_i{x}_i{\sum }_i{y}_i}{\sqrt{N{\sum }_i{x}_i^{\mathrm{2}}-({\sum }_i{x}_i{)}^{\mathrm{2}}}\sqrt{N{\sum }_i{y}_i^{\mathrm{2}}-({\sum }_i{y}_i{)}^{\mathrm{2}}}}},\end{array}$$

where x_i and y_i are two sets of data with the same number (N). As described above, we use the averaged width of the detected plume in steps of 5 nT (Dst index) to represent the plume width in the corresponding Dst range; N is equal to 15 if the geomagnetic activity levels from Dst ∼ −15 nT to Dst ∼ −90 nT are considered in the study.

An explanation of the poor relevance is that some processes may influence the structure of the plasmasphere and plume. For example, the advent of quiet conditions after geomagnetic storms can contribute to the refilling of the plasmasphere because ionospheric particles are drawn upward from low altitudes along magnetic field lines. In addition, the number of samples with geomagnetic storms that are too strong or too weak may be too few to have statistical significance. Therefore, to better understand the influence of storm intensity on the plume width, we further sifted the plume events.

First, considering that the plasmasphere can be obviously refilled after the time of geomagnetic disturbance, we only retain the events during the recovery phase of the geomagnetic storm. This operation ensures that the main factor affecting the structure of the plume is the erosion process of the geomagnetic storm, which is the main topic in the study. Similar to the standards described in Engebretson et al. (2008), Halford et al. (2010), and Bortnik et al. (2008), we define the main phase as the period from onset of the storm until the Dst reaches its minimum value, and the end of the recovery phase is defined as the fifth day after the main phase finishes. The numbers of plume events corresponding to different time intervals (here, 1 d represents the interval of 0–24 h) after the minimum Dst are shown in Table 1. The total number of plume events in the recovery phase (the interval time is 1–5 d) is 377.

Table 1 The number of events corresponding to different intervals after the minimum Dst.

Among the 377 plume events, the minimum Dst value of the most intense geomagnetic storm reaches −209 nT, but 68 % of plume events (256 events) correspond to Dst values ranging from −70 to −15 nT, and 88 % of plume events (333 events) correspond to the minimum Dst values of the most intense geomagnetic storm ranging from −90 to −15 nT. For the accuracy of statistical research, we exclude extremely intense storms in the study, and we only statistically analyze plasmaspheric plume events from −70 to −15 nT (and from −90 to −15 nT) during the recovery phase.

Moreover, in addition to the density exceeding the Sheeley model for more than 10 min outside the plasmapause, we further improve the standard of plume judgment. Referencing the method of Darrouzet et al. (2008), the ΔL of the structure (the difference in the L shell between the entrance and exit of the plume orbit) should be large enough (0.2) to be considered a plume. On the other hand, the events with excessive linear width (>3.5 R_E) are also considered not plumes because they are more likely to be the cross section of the plasmasphere rather than the plume. Based on this standard, we exclude 111 plume events in the range from −70 to −15 nT and 168 plume events in the range from −90 to −15 nT.

The process of deleting events that do not match the standard is shown in Fig. 3. The solid orange dots indicate the events that are not in the recovery phase, which are beyond 5 d after the minimum Dst value. The solid purple dots indicate plume events with corresponding Dst values less than −70 nT or greater than −15 nT. The solid blue and red dots represent the tracks with ΔL less than 0.2 and Cartesian distances greater than 3.5 R_E, respectively. The solid gray and black dots indicate the retained plume events with corresponding Dst ranges from −90 to −70 nT and −70 to −15 nT, respectively. Ultimately, there are 145 retained plume events in the Dst range from −70 to −15 nT, and 165 retained plume events in the Dst range from −90 to −15 nT. In addition, the spatial distribution of 586 events is shown in Fig. 4. The curves represent the VAP-A orbits while the events are observed. The colors of the curves represent the filtering process discussed above, and the corresponding event color is consistent with Fig. 3. It also shows that the retained plumes (indicated by the black curves) are mainly observed on the duskside (MLT ∼ 15:00 to ∼ 21:00).

Figure 3 The category of observed events: the solid orange dots indicate plume events that are not in the recovery phase. The solid purple dots display the events with a corresponding Dst index less than −70 nT or greater than −15 nT. The solid blue and red dots represent the events with ΔL less than 0.2 and Cartesian distances greater than 3.5 R_E, respectively. The events in the Dst range of −70 to −15 nT and −90 to −70 nT eventually retained are indicated by black and gray dots, respectively.

Figure 4 The spatial distribution of 586 plasmaspheric plumes is shown in the MLT–L plane. The color codes are the same as those in Fig. 3.

Figure 5a and b show the correlation analysis between the retained plume width and the level of geomagnetic storm intensity with a minimum Dst of −70 to −15 nT. The formats are similar to Fig. 2a and b. Completely different from the results before the sifting plume events (as shown in Fig. 2), as indicated by the blue lines in Fig. 5, there is a considerable negative correlation between the plume width and level of corresponding geomagnetic storm intensity. This implies that as the minimum Dst value becomes lower, the width of the plasmaspheric plume tends to become narrower. As presented in Fig. 5a, the Pearson correlation coefficient between the averaged Cartesian distance of plumes and the Dst value reaches 0.619. The value of the P value is 0.042, and the low value of the P value means that it is feasible to express the relationship as a linear one. As presented in Fig. 5b, the Pearson correlation coefficient between the averaged ΔMLT and the Dst value reaches 0.546, and the P value is 0.068. The interpretations from Fig. 5a and b both imply that there is a roughly negative correlation between the width of the plume in the recovery phase and the level of geomagnetic storm intensity.

Figure 5 The format is the same as in Fig. 2. However, only the retained plume events that meet more stringent sifting conditions are analyzed.

Similarly, Fig. 5c and d show the correlation analysis between the retained plume width and storm intensity with a minimum Dst from −90 to −15 nT. As exhibited in Fig. 5c, the Pearson correlation coefficient between the averaged Cartesian distance of plumes and the Dst value is 0.580. The value of the P value is 0.021. This also implies that there is a roughly negative correlation between the plume width and levels of geomagnetic storm intensity, with a minimum Dst from −90 to −15 nT. As presented in Fig. 5d, the Pearson correlation coefficient between the averaged ΔMLT and the Dst value is 0.370, and the P value is 0.170. From the perspective of the ΔMLT analysis, it seems that the negative correlation in the Dst range from −90 to −15 nT is weaker.

For the accuracy of statistical research, we exclude extremely intense storms in the study; only the plume events that correspond to intervals of minimum Dst indices from −70 to −15 and −90 to −15 nT are statistically analyzed.

It seems that the negative correlation for the events with minimum Dst from −90 to −15 nT decreases slightly compared to those with minimum Dst from −70 to −15 nT.

In this study, the Cartesian distance and ΔMLT are adopted to represent the detected plume width. Both methods imply that there is a negative correlation between the width of the plume in the recovery phase and the level of corresponding geomagnetic storm intensity. To compare the similarities and differences between the two standards, Fig. 6 exhibits the relationship between $\mathrm{\Delta }L/\langle L\rangle$ and the ΔMLT of plumes, where ΔL indicates the difference between the L shells of the entrance and exit of the plume detected by the VAP-A, and 〈L〉 indicates the average value of L on the plume orbit. There is a positive Pearson correlation coefficient (0.471 for the Dst range from −70 to −15 nT, 0.396 for the Dst range from −90 to −15 nT) between them, which means that when $\mathrm{\Delta }L/\langle L\rangle$ increases, the tendency of ΔMLT also enhances, although this positive correlation is not too strong. The minimal P value (∼ $\mathrm{2.890}\times {\mathrm{10}}^{-\mathrm{10}}$ and ∼ $\mathrm{4.300}\times {\mathrm{10}}^{-\mathrm{7}}$, respectively) shows that the linear relationship between them is very reliability.

Figure 6 The relationship between $\mathrm{\Delta }L/\langle L\rangle$ and ΔMLT.

3 Simulation of plume evolution

To clearly exhibit the effect of the storm intensity on the width of the plasmaspheric plume, the evolutions of the plumes during two geomagnetic storms (with minimum Dst values equal to −39 and −74 nT) are simulated through test particle simulations. In this study, this process differs from the plasmapause test particle (PTP) simulation which only provides the evolution of the plasmapause and plasmaspheric plume boundaries in Goldstein et al. (2004, 2014a, b); the particle simulations in this study also calculate the evolution of density in the plasmasphere and plasmaspheric plume.

3.1 Model inputs

This simulation assumes that all particles move in the dipole magnetic field model. Considering not only that the plasma motion will be driven by the combined action of the corotating electric field and convection electric field during the geomagnetic storm but also that the subauroral polarization stream (SAPS) will play a significant modification of the convection electric field, the magnetospheric electric model is assumed to consist of three parts as follows:

1. The corotation electric potential is Φ_rot with the following formula:

   $$\begin{array}{}\text{(4)}& {\mathrm{\Phi }}_{\mathrm{rot}}=-C{\displaystyle \frac{R_{\mathrm{E}}}{r}},\end{array}$$

   where C indicates a constant of 92, which is provided by Völk and Haerendel (1970), and r indicates the geocentric latitude.

2. The convection electric potential is expressed as follows (Maynard and Chen, 1975):

   $$\begin{array}{}\text{(5)}& {\mathrm{\Phi }}_{\mathrm{VS}}=-E_{\mathrm{IM}}{r}^{\mathrm{2}}\sin \mathit{\phi }(\mathrm{6.6}{R}_{\mathrm{E}}{)}^{-\mathrm{1}},\end{array}$$

   where E_IM is the inner magnetospheric electric field. While the southward IMF turns southward, $E_{\mathrm{IM}}=\mathrm{0.12}\cdot \mathrm{|}{E}_{\mathrm{SW}}\mathrm{|}$, and when the IMF is reverse, it is equal to 0.12⋅0.25 mV m⁻¹. Here, the E_SW is the solar wind electric field. The azimuth angle is indicated by φ, where MLT = $\mathrm{12}(\mathit{\phi }/\mathit{\pi }+\mathrm{1})$.

3. The SAPS electric potential is an intense, radially narrow, westward flow channel that is mainly located in the dusk-to-midnight MLT area (Burke et al., 1998; Foster and Burke, 2002) and is considered to significantly modify convection. As Goldstein et al. (2005a, b, 2014b) suggest, the effects of the SAPS in the equatorial magnetosphere are driven by the Kp index, and the electric potential is described as follows:

   $$\begin{array}{}\text{(6)}& {\mathit{\phi }}_{\mathrm{s}}\left(r\mathit{\phi }t\right)=-F\left(r\mathit{\phi }\right)G\left(\mathit{\phi }\right)V_{\mathrm{S}}\left(t\right),\end{array}$$

   where F(rφ), G(φ), and V_S(t) are functions parameterized by the magnetic latitude, MLT(φ), time (t) and Kp index.

The detail formulas of functions are as follows:

1. The function F(rφ) treats the SAPS flow channel as a potential drop centered at radius R_s:

   $$\begin{array}{}\text{(7)}& F\left(r\mathit{\phi }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}{\tan }^{-\mathrm{1}}\left[{\displaystyle \frac{\mathrm{2}}{\mathit{\alpha }}}\left\{r-R_{\mathrm{s}}\left(\mathit{\phi }\right)\right\}\right],\end{array}$$

   where R_s, α are represented in the following way:

   $$\begin{array}{}\text{(8)}& R_{\mathrm{s}}\left(\mathit{\phi }\right)=R_{\mathrm{0}}{\left[{\displaystyle \frac{\mathrm{1}+\mathit{\beta }}{\mathrm{1}+\mathit{\beta }\cos \left(\mathit{\phi }-\mathit{\pi }\right)}}\right]}^{\mathit{\kappa }},\end{array}$$

   with β=0.97 and κ=0.14; and R₀ is expressed as follows:

   $$\begin{array}{}\text{(9)}& R_{\mathrm{0}}/R_{\mathrm{E}}=\mathrm{4.4}-\mathrm{0.6}\left(\mathrm{Kp}-\mathrm{5}\right),\end{array}$$

   where R_E is the radius of the Earth, 6380 km.

   α is expressed as

   $$\begin{array}{}\text{(10)}& \mathit{\alpha }=\mathrm{0.15}+\left(\mathrm{2.55}-\mathrm{0.27}\mathrm{Kp}\right)\left[\mathrm{1}+\cos \left(\mathit{\phi }-{\displaystyle \frac{\mathrm{7}\mathit{\pi }}{\mathrm{12}}}\right)\right].\end{array}$$

2. Based on the magnitude of the SAPS potential drop decreases eastward across the nightside, azimuthal modulation of SAPS magnitude G(φ) is set to

   $$\begin{array}{}\text{(11)}& \begin{array}{rl}G\left(\mathit{\phi }\right)& =\sum _{m=\mathrm{0}}^{\mathrm{2}}\mathit{\{}{A}_{\mathrm{m}}\cos \left[m\left(\mathit{\phi }-{\mathit{\phi }}_{\mathrm{0}}\right)\right]\\ & +B_{\mathrm{m}}\sin \left[m\left(\mathit{\phi }-{\mathit{\phi }}_{\mathrm{0}}\right)\right]\mathit{\}}.\end{array}\end{array}$$

3. The V_S(t) function describes the time regulation of SAPS:

   $$\begin{array}{}\text{(12)}& V_{\mathrm{S}}=\left(\mathrm{0.75}\mathrm{kV}\right){\mathrm{Kp}}^{\mathrm{2}}.\end{array}$$

The initial plasma density distribution is assumed to change as a function of the L shell ($\mathrm{2}\le L\le \mathrm{7}$) according to the model obtained from Sheeley et al. (2001) (the formula of the Sheeley model is expressed as Eq. 1) with no MLT dependence. In addition, for regions where the L shell is larger than 7, the electron density remains at 5 cm⁻³. All particles (approximately 128 000 electrons in total) emitted in the model are considered to be cold electrons and assumed to have an initial energy of 10 eV. Here, the motions of electrons are assumed to be adiabatic. We calculate the drift velocity as a combination of the velocity due to E×B drift and the bounce-averaged velocity due to gradient and curvature drifts (Roederer, 1970; Ganushkina et al., 2005; Li et al., 2021). Here, the pitch angles of the 128 000 electrons are deemed arbitrary, because the electron energy is considered to be small enough that the associated gradient-curvature drift velocity is very small (Roederer and Zhang, 2014). The motions of electrons are mainly contributed by the E×B drift. According to the simulation results, we can calculate the particle density in a certain area to reflect the evolution of the plasmasphere and plume. Notably, the actual shape of the plasmapause is too complicated to obtain its actual electron distribution function, so the above typical model electron density distribution is adopted in the research.

3.2 Evolution of plume from 19 to 20 May 2017

The geomagnetic and solar wind indices during the geomagnetic storm from 19 to 20 May 2017 are displayed in Fig. 7. As shown in Fig. 7b, at 07:02 UT, the IMF turned southward, leading to a larger negative B_Z, and the geomagnetic storm began. The beginning of the geomagnetic storm is denoted by the vertical dashed blue line. The minimum Dst index of this geomagnetic storm was −39 nT. As shown in Fig. 7c, during the main phase of this geomagnetic storm (from approximately 12:00 UT on 19 May to 09:00 UT on 20 May), the maximum E_IM index reached 0.4684 mV m⁻¹. As shown in Fig. 7d, the Kp index that drove the SAPS model reached a maximum of 4+.

The test particle simulation started at 07:02 UT on 19 May 2017 (the beginning of the geomagnetic storm). The initial distribution of particle density is shown in Fig. 8a. In the contribution of the convection electric field, the particles in the plasmasphere obviously move sunward within 5 h from 07:02 to 12:02 UT on 19 May 2017 (as shown in Fig. 8b). Meanwhile, some of the plasmaspheric particles expand to high locations with L>8 and may be lost to the magnetopause boundary (Spasojevic et al., 2005). During the interval of the next 5 h, more particles move sunward and reach the model boundary (as shown in Fig. 8c), and the L shell of the plasmapause on the nightside obviously decreases. As shown in Fig. 8d and e, from 22:02 to 03:02 UT on 20 May 2017, the width of the plasmaspheric bulge gradually shrinks, and a plume gradually forms on the afternoon side. Furthermore, under the action of a corotation electric field, the plasmaspheric plume gradually shifts toward the nightside. From 08:02 to 13:02 UT on 20  May (as shown in Fig. 8f and g), as the convection electric field and SAPS electric field increase again, a large number of particles move in the sunward direction, and the plasmaspheric plume narrows with time. For most times in the recovery phase, the convection electric field becomes weaker, and the formed plasmaspheric plume slowly rotates from the afternoon side to the nightside (as shown in Fig. 8g to h). The plasma density detected by the VAP-A from 17:36 to 22:39 UT on 20 May 2017 is shown in Fig. 8i. From 18:12 to 19:53 UT on 20 May 2017, the VAP-A operated near the apogee of its orbit on the duskside, and the plume also rotated to the duskside. The orbit of the VAP-A, while it actually observed the plume, is indicated by the black curve in Fig. 8h. We can see that the observed plume roughly coincides with the simulated plume in Fig. 8h. As the positions of the simulated plume and observed plume are roughly identical at this time, we believe that the initial distribution of particles and magnetospheric electric field models used in this paper are basically reliable.

Figure 7 The indices of geomagnetic activity and solar wind during the geomagnetic storm that occurred in the time interval of 19–20 May 2017: (a) Dst index, (b) B_Z index in GSM coordinates. The dotted red line indicates the position where B_Z is equal to 0, (c) E_IM index, and (d) Kp index. The start time of the geomagnetic storm (07:02 UT on 19 May 2017) is marked by the vertical dashed blue line.

Figure 8 (a–h) Equatorial plots of the plasmasphere and plume obtained by simulation in the time interval of 19–20 May 2017. The dotted white circles, from inside out, indicate L shell values of 4, 6, and 8. The simulation duration and the corresponding actual duration are represented in the title. The solid black line in (h) represents the orbit arc in which the VAP-A observed a plasmaspheric plume from 18:12 to 19:53 UT on 20 May. (i) The plasma density detected by the VAP-A from 17:36 to 22:39 UT on 20 May 2017. The gray shadow indicates the plasmaspheric plume detected by the VAP-A.

Figure 9 The indices of geomagnetic activity and solar wind during the geomagnetic storm that occurred in the time interval of 8–10 June 2015; the format is the same as that in Fig. 7.

Figure 10 (a–h) Equatorial plots of the plasmasphere and plume obtained by simulation in the time interval of 8–10 June 2015. The dotted white circles, from inside out, indicate L shell values of 4, 6, and 8. The simulation duration and the corresponding actual duration are represented in the title. The solid white line in (h) represents the orbit arc in which the VAP-A observed the plasmaspheric plume from 06:39 to 07:31 UT on 10 June. (i) The plasma density detected by the VAP-A from 05:00 to 11:16 UT on 10 June 2015. The gray shadow indicates the plasmaspheric plume detected by the VAP-A.

Table 2 The widths of the simulated plume at different L shells during the above two geomagnetic storms.

3.3 Evolution of plume from 8–10 June 2015

The geomagnetic and solar wind indices during the geomagnetic storm from 8 to 10 June 2015 are shown in Fig. 9. The geomagnetic storm began at 00:18 UT on 8 June 2015 (denoted by the blue line). As shown in Fig. 9a, the minimum Dst index of this geomagnetic storm was −74 nT, which was much lower than that during the storm presented in Sect. 3.2. At the beginning of the geomagnetic storm, the B_Z turned southward, and the E_IM value (calculated from E_SW) became a relatively large positive value. The maximum E_IM in this geomagnetic storm was 1.0674 mV m⁻¹, which was much greater than the 0.4684 mV m⁻¹ obtained from the geomagnetic storm presented in Sect. 3.2. Meanwhile, the maximum Kp index reached 6 in the main phase. It was also much larger than 4+ presented in the last storm. Larger E_IM and Kp indices both implied that convection during this geomagnetic storm was much more intense than that during the geomagnetic storm from 19 to 20 May 2017.

Figure 10a shows the electron density distribution during the first minute of the simulation during this geomagnetic storm. Then, under the contribution of a continuous convection electric field, the plasma in the outer part of the plasmasphere moves along the E×B-drift paths from 00:19 UT on 8 June to 06:39 UT on 10 June. As shown in Fig. 10b and c, the more intense convection electric field brings about a larger number of particles on the dayside moving sunward and extending out of the model boundary, and the particles on the nightside in the plasmasphere move faster toward the Earth. As shown in Fig. 10d, the particles in the outer part of the plasmasphere dissipate at 00:18 UT on 9 June, and a narrower plume emerges near the duskside. During the recovery phase of the geomagnetic storm, as shown from Fig. 10e to h, the formed narrow plume revolves around the Earth. Finally, the VAP-A observes the structure of the plume from 06:39 to 07:31 UT on 10 June 2015 (as shown in Fig. 10i), and the orbit of the probe during this time interval is indicated in Fig. 10h. The observations and simulations both suggest that the L shell of the plasmapause is lower than that presented in Sect. 3.2. Although there is some difference in the MLT between the simulated plume and observed plume (approximately 1 MLT, which may be due to a simulation time that is too long), the results imply that the width of the plume is much narrower than that during the geomagnetic storm from 19 to 20 May 2017.

3.4 Comparison of simulation results

To better exhibit the difference in the simulated plume width during the above two geomagnetic storms, we calculated the Cartesian distances and ΔMLTs at positions where the L shell was equal to 5, 5.5, and 6. The calculated results are shown in Table 2. As presented in Table 2, at the same L shell, regardless of Cartesian distances and ΔMLT, the simulated plume driven by stronger geomagnetic storms is always narrower.

4 Discussion and conclusion

In this study, we present statistical research on the relationship between the widths of plumes and the intensities of geomagnetic storms by analyzing the data collected by the VAP-A. Here, the widths of the detected plume are defined as the Cartesian distance and ΔMLT of the detected plume orbit. In the first step, by directly analyzing all 586 potential plume events after the minimum Dst of a geomagnetic storm, we find that there is almost no correlation between plume width and the level of storm intensity. This result is similar to the conclusion obtained from Borovsky and Denton (2008), which suggests that the linear correlation coefficient between them is almost 0. Since the plasmasphere can be eroded by the enhanced convection electric field during geomagnetic storms (Krall et al., 2017), the enhanced convection field causes low energy plasma drainage to the magnetopause (Denton et al., 2005). We consider that the level of storm intensity may affect the width of the plume in some conditions. In the second step, we define the end of the recovery phase as 5 d after the main phase finishes. Only the plumes in the recovery phase interval are analyzed. Moreover, the criterion of plumes for statistical investigation is further improved. This result suggests that there is a negative correlation between the plume width and absolute value of the minimum Dst value during the storm, although the negative correlation is not very strong.

To explain the negative correlation between them during the recovery phase of the geomagnetic storm, the group test particle simulation is adopted to reveal the dynamic evolutions of the plasmasphere and plume during two geomagnetic storms (with minimum Dst values of −39 and −74 nT, respectively). By comparing the evolutions during the two storms, we find that in the more intense geomagnetic storm, the erosion of the plasmasphere is more severe, most particles in the outer part of the plasmasphere are dissipated during the initial and main phases, and a narrower plume is exhibited during the recovery phase. Although the evolutions of plasmaspheres and plumes may be very complicated, the above simulation results provide a reasonable candidate explanation for the negative correlation between level of geomagnetic storm intensity and plume width during the recovery phase of geomagnetic storms.

As shown in Fig. 4, most plume events are mainly observed on the duskside; the relationship between the plume width and its MLT is also a meaningful work, which will be studied in our next project. Since there may be a short time delay with several minutes between the changes of E_SW and E_IM (Nishimura et al., 2009), a more precise and real-time E_IM model will be discussed and explored in the future. This paper analyzes the plume events with Dst values ranging from −90 to −15 nT during the recovery phase. In the future study, we will use other geomagnetic activity indices to analyze the relationship between plume width and magnetic storm intensity, such as Kp, and AE, thus studying the correlation between the plume width and substorm activity.