Simulation of Cold Ion Transport Originating From the SED Plume Into Dayside Magnetosphere

In this study, based on the Dynamic‐Kinetic Fluid (DyFK) model, we conduct simulations on the field‐aligned transport of cold ions (<1 eV) within three independent closed flux tubes from different starting points (L = 3.0 Re, L = 3.122 Re and L = 3.5 Re, respectively) during the 2015 St. Patrick's Day storm. These flux tubes magnetically connect the plasmaspheric plume with the ionospheric Storm Enhanced Density (SED) plume but drift toward the dayside magnetopause under E×B effects. Meanwhile, we can obtain the density variations of three ions species (H+, O+, and He+, respectively) around the flux tubes' apex (at the magnetic equator) as a function of the L shell. The simulations present the temporal evolution of cold plasma along these drifting flux tubes. We found that with increasing initial L shell location of the flux tube, there are more ionospheric cold light ions (H+ and He+) flowing into the equatorial region while the inflows of cold heavy ions (O+) toward the equator decrease. As the simulated flux tube drifts toward dayside magnetopause from the main plasmasphere (L = 3.0 Re), the SED sourced heavy ions pile up near the magnetic equator at low L shells (L < 4.0 Re) due to gravitational effects, while the light ions rapidly reach density equilibrium along the flux tube. With the flux tube drifting beyond L > 4.0 Re, the piled‐up cold heavy ions will later transport into the magnetic equatorial region and furnish the dayside magnetosphere, demonstrating that the SED is another supplier of magnetospheric cold heavy ions (O+).

Previous research also found that the ionospheric storm enhanced density (SED) plume if it exists, is magnetically connected to the plasmaspheric plume (Coster et al., 2003;Coster & Skone, 2009;Foster et al., 2002;Kelle et al., 2004;Yizengaw et al., 2008). Both simulations and observations implied that the underlying ionospheric SED plume should be able to supply cold plasma for the plasmaspheric plume extended into the dayside magnetosphere by field-aligned flows (Borovsky et al., 2014;Huba & Sazykin, 2014;Moldwin et al., 2016;Yuan et al., 2009). In addition to plasmaspheric cold plasma transported sunward by flux tube convection, Krall et al. (2018) suggested that the high latitudinal ionosphere could be a sufficient source for magnetospheric cold populations through supersonic field-aligned outflows. Simulation of Wang et al. (2016) using the DyFK model also showed significant mass loading in the vicinity of the magnetopause. However, our previous simulation study (Qiao et al., 2019) on field-aligned cold plasma transport from SED plume to plasmaspheric plume during a superstorm only found slight equatorial mass loading at sub-auroral latitudes (L < 4.0 Re) which is different from that at high latitudes. In that work, the temporal evolution of cold plasma within one single flux tube initialized at L = 3.122 Re was simulated. To better understand the processes of the field-aligned cold plasma transport from the ionospheric SED plume to the plasmaspheric plume, it is necessary to simulate temporal variations of the field-aligned distribution of cold plasma along these drifting flux tubes initialized at different L-shells within the SED plume. In this study, using the DyFK model we examine the contribution of cold ion populations in the SED plume to the magnetosphere in the dayside. During the known geomagnetic storm that occurred on St. Patrick's Day in 2015, we focus on temporal evolution and spatial variation of field-aligned cold plasma transport within closed sunward drifting flux tubes, which are initialized at different L shells simultaneously. The next section presents the geophysical conditions during the event studied and briefly describes the numerical model used. Section 3 presents the simulation results and discussion. Summary and conclusions are given in Section 4.

Geophysical Conditions and DyFK Model Description
A major purpose of this work is to investigate field-aligned flows of cold plasma originating from the SED plume into the plasmaspheric plume which might intersect with dayside magnetopause. We chose the known geomagnetic superstorm that occurred in 2015 on St. Patrick's Day. The observed geomagnetic indices such as By/ Bz components in Geocentric Solar Magnetospheric coordinate systems, the solar wind dynamic pressure, Kp, Symmetric disturbance index of the horizontal component of the geomagnetic field (SYM-H), and F 10.7 indexes are shown in Figure 1, together with the input Kp index to the simulation model (the dashed line in panel c).
The Dynamic Fluid-Kinetic (DyFK model simulates the field-aligned plasma transport within a single closed flux tube subjected to ExB drift (Tu et al., 2004(Tu et al., , 2005Wang et al., 2016). The DyFK model couples Field Line Interhemispheric Plasma (FLIP) model and Generalized Semi-Kinetic (GSK) model with an overlapped region in 800-1,100 km altitude range based on a semi-kinetic approach (Estep et al., 1999;Wang et al., 2015). The FLIP model, which has been extensively used in ionospheric and plasmaspheric research for decades, solves the time-dependent mass, momentum, and energy transport equations including comprehensive chemistry of multiple ion species (O + , H + , He + and N + ) along a closed flux tube above 120 km altitude (Richards & Torr, 1990;Richards et al., 2000). In the DyFK model, the FLIP model is truncated at 1,100 km altitude and treated as the ionospheric part in both hemispheres. Above 800 km, ions are treated as simulation particles and electrons serve as fluids maintaining quasi-neutrality (Wilson et al., 1992). The field-aligned motion of simulation particles is subjected to ambipolar electric field (Mitchell et al., 1992), gravity, magnetic mirror force, centrifugal force, and Coulomb collisions (Nanbu, 1997;Nanbu & Yonemura, 1998). Meanwhile, the density and velocity of electrons are derived from the quasi-neutrality condition and parallel current-free assumptions, respectively. The electron 3 of 15 temperature is calculated from an empirical temperature model (Titheridge, 1998). In both hemispheres, the ionospheric fluids are passed to the kinetic domain at 800 km on the basis of gyro-averaged drift Maxwellian velocity distributions of different ion species (H + , O + , and He + ) through particle injection (Estep et al., 1999;Wang et al., 2015). While the bulk parameters (densities, bulk velocities, temperatures, and heat fluxes) in the kinetic domain are fed back to the fluid domain as upper boundary conditions at 1,100 km. The parameters in the two simulation domains reach a smooth connection for the coupled FLIP and GSK models in the overlapped region. Based on the convection trajectory of a single flux tube together with Kp and F 10.7 indexes as input parameters, the DyFK model will provide field-aligned distributions of ion density, velocity, anisotropic temperatures, as well as velocity distributions of multiple ion species (O + , H + , He + ), including the effects of geometrical change associated with flux tube E×B drift.
Since we focus on the cold (<10 eV) ion populations in this study, we only consider the E×B drift of the flux tubes, neglecting curvature and gradient drifts. We combine the E5D magnetospheric convection electric potential model (McIlwain, 1986), corotation electric field, and G05 magnetospheric SAPS electric potential model  to calculate the drift trajectories of flux tubes as done in Qiao et al. (2019). In addition, in Qiao et al. (2019), we have shown the projection of the open passage of the flux tube on the global ionospheric total electron content (TEC) map. To check if the chosen flux tubes follow the plasmaspheric plume, instead, we compared the calculated E×B drift trajectories of these flux tubes' apex (at the magnetic equator) with the magnetic equatorial projection of the global positioning system (GPS) TEC map using the magnetospheric magnetic field model T04s (Tsyganenko, 2005) and the plasmapause location derived from plasmapause test particle (PTP) simulation (Goldstein et al., 2014). Determinations of vertical TEC have been binned in

Simulation Results and Discussion
In this study, our purpose is to determine the contribution of SED plumes to magnetospheric cold populations.
Considering that the Earth's magnetosphere would be compressed under extreme solar wind conditions, we stop the flux tube drift trajectory at geosynchronous orbit (L ∼ 6.6 Re). In our previous study (Qiao et al., 2019), the simulation results showed SED plume filling plasmaspheric plume at sub-auroral latitudes instead of at high latitudes, in other words, there is no significant mass loading due to the ionospheric plume to magnetosphere at high L shells (L > 4.0 Re). In that simulation, the single flux tube spent quite a while (nearly 3 hr) drifting into the high latitudinal region (L > 4.0 Re) from its initial location (L = 3.122 Re), and a large amount of SED sourced cold Patrick's Day mapped to magnetic equator plane in SM coordinate system based on T04_s magnetic field model (Tsyganenko, 2005). In each panel, the sun is to the right, the black dotted line shows the location of plasmapause derived from plasmapause test particle simulation results (Goldstein et al., 2014), the black solid line is the whole convection path of the flux tube drifting from L = 3.0 Re, and the magenta dot refers to the footprint of single flux tube in magnetospheric equatorial plane at current time epoch.
plasma passed through 900 km altitude and flowed upward along the flux tube. The initial location of flux tube likely affects the drift path as well as the upward plasma flows, which are considered in this study.
In the following, we present simulation results of three independent flux tubes starting from L = 3.0 Re, 3.122 Re, and L = 3.5 Re at 1200 MLT drifting toward magnetopause, respectively. Note that the results for the flux tube initialized from 3.122 Re is adapted from Qiao et al. (2019).   of L-shell of each flux tube shown in Figure 5a. As shown in Figure 5b, when the start point of the flux tube is located deeper inside the plasmasphere, the flux tube will experience a longer duration of convection in the plasmasphere before reaching the stagnation point.

The Field-Aligned Transport of O + , H + , and He + Ions
Based on the calculated convection trajectories, we conduct DyFK simulation on these three flux tubes. Verified by comparing with TEC observations as mentioned above, the flux tubes all magnetically connect SED plume and plasmaspheric plume but start sunward drifting from different initial L shells around MLT noon. Figures 6a and 6c present the density of three ion species (O + , H + , and He + ) varying with L shell at about 900 km altitude and at the flux tube's apex as a function of L shell, respectively. Note that the densities of O + and He + in Figure 6c are multiplied by 2000 and 10, respectively, so that the densities of different ion species could be shown in the same range. In the DyFK model, the division of grid cells depends on their latitudinal angles along the flux tube. so that the volume and cross section of grid cells expand with outward drifting flux tubes, which would dilute the particle fluxes passing through. Here in Figure 6b, we show the number flux of ions (units in s −1 ) passing through the cross section of ionospheric and apex grid points along the flux tube derived from bulk parameters, corresponding to the data points in 6a and 6c, respectively. The black dashed lines in each panel of Figure   As mentioned above, in order to avoid the expansion effects of cross section in grid cells when the flux tube is drifting sunward, Figure 6b displays field-aligned number flux of different ion species (units in s −1 ). It can be clearly seen that the ionospheric field-aligned flux of light ions (H + and He + ) is in sync with variations of their ionospheric density (Figure 6a). In contrast, the changes in ionospheric O + flux show a significant difference from ionospheric O + density variations. Particularly in simulations starting from L = 3.0 Re (panel 1 of Figure 6b) and L = 3.122 Re (panel 2 of Figure 6b), the ionospheric O + flux goes up and down while the ionospheric O + density is constantly increasing, indicating prominent gravitational impacts on heavy particles. The gravitational force to a large extent alters the field-aligned transport of O + particles while having little effect on the H + and He + ions, which would cause a significant discrepancy in equatorial densities between light and heavy ions (Figure 6c). For the same ion specie, the trends of number flux both in the ionosphere and around the flux tube apex are similar, but with different magnitudes. For example, in panel 1, 2 and 3 of Figure 6b, the number flux of ionospheric O + (red lines) shows similar variations, and so does the number flux near the apex (gray lines). For light ions (H + and He + ), the magnitudes of number flux around flux tube apex are likely to increase with increasing initial L shell location of flux tubes, while for heavy ions (O + ) the trend tends to be descending.
On the other hand, the red lines in panel 1, 2, and 3 of Figure 6b show a declining trend of total ionospheric supply of heavy ions (O + ) toward increasing initial L shell of the flux tubes, while for light ions (H + and He + ) the trend is not obvious. Around the flux tube apex, the total O + inflow (gray lines in panel 1, 2, and 3 of Figure 6a) tends to decrease with increasing initial L shell location of the flux tubes as well, while the total inflow of light ions (gray lines in panel 7, 8 and 9 for H + , gray lines in panel 4, 5, and 6 for He + ) presents a growing trend. The time integral of upward ionospheric number flux also reveals the different behaviors between heavy and light ions. As listed in Table 1, for the simulation starting at L = 3.0 Re, roughly 3.85 × 10 14 O + particles are flowing upward from the ionosphere in total by the time the flux tube reached geosynchronous orbit, while the numbers for simulations starting at L = 3.122 Re and at L = 3.5 Re are 2.68 × 10 14 and 1.48 × 10 13 , respectively, implying a declining trend as noted above. The total ionospheric upward fluxes of He + and H + basically follow the same trend as well but with lower magnitudes at the order of 10 11~1 0 12 . Meanwhile, at the apex of the flux tube, there are about 8.13 × 10 9 , 2.58 × 10 9 , and 7.50 × 10 8 total O + flowing inward for simulations starting at L = 3.0 Re, L = 3.122 Re, and L = 3.5 Re, respectively. For light ions (H + and He + ), their total inflow around the flux tube apex shows a similar trend but with much larger magnitudes.
Maintained by field-aligned fluxes, the density of heavy ions (O + ) around the flux tube apex (equatorial density) exhibits distinct characteristics different from that of light ions (H + and He + ). As shown in Figure 6c, the equatorial densities of light ions are roughly in line with the expansion of flux tube volume (black dashed lines), meaning that the inflows and outflows are almost even and that the density distribution of light ions rapidly reaches equilibrium along the flux tube. However, the equatorial O + density shows a growing tendency of particle loss around the flux tube apex toward simulations initialized at larger L shells, judging by the separations between simulated O + density and the other density decaying with volume expansion. What is more interesting is that in the simulation starting from L = 3.0 Re (panel 1 of Figure 6c), the equatorial O + density shows prominent enhancement twice, outstripping the expected values of adiabatic expansion.

The SED Sourced Mass Loading in the Dayside Magnetosphere
In our previous work (Qiao et al., 2019), as the flux tube drifts sunward, there is slight equatorial mass loading via field-aligned plasma transport sourcing from the SED region at sub-auroral latitudes instead of at high latitudes. In this study, the results display the same sub-auroral enhancements of equatorial O + density both in simulations starting from L = 3.0 Re (panel 1 of Figure 6c) and from L = 3.5 Re (panel 3 of Figure 6c). However, in the simulation starting from L = 3.0 Re, we see two instances of equatorial mass loading: the equatorial O + density at first increases then hits its minimum around L = 5.5 Re after gradual decreases, and later the equatorial mass loading occurs for the second time at large L shells.  As we discussed in Qiao et al., 2019, after the sudden increase of geomagnetic activity, the ionospheric plasma density first increases due to enhanced ionization. Then the plasma density continues to increase under positive storm effects caused by traveling amospheric disturbances (e.g., Proelss et al., 1991;Yuan et al., 2003) such as uplifting of the F layer. In addition, various other mechanisms, such as soft electron precipitation, E×B convection flows, and thermosphere wind can also contribute to density increase within SED (e.g., Heelis et al., 2009;Liu et al., 2016;Lu et al., 2012;Zou et al., 2014Zou et al., , 2017. In this paper, we focus on the single drifting flux tube condition. The intense plasma concentration triggers the growing ionospheric upward fluxes and leads to the rise in equatorial densities at sub-auroral latitudes (Figure 7). With the foot of the flux tube entering into a high-latitudinal region (L > 4.0 Re), the heating effects together with positive storm effects will again strengthen ionospheric upward supplies, even further prominently affect the tendency of equatorial densities, for example, the equatorial O + density in the simulation starting from L = 3.5 Re (panel 3 of Figure 6c). But in the simulation starting from L = 3.0 Re, it seems that neither the disturbed geomagnetic activity nor the high latitude heating effects alone are able to account for such an increase of equatorial O + density at L shells larger than 5.5 Re (panel 2 of Figure 6c).
The simulation starting from L = 3.0 Re shares a number of similarities with the simulation starting from L = 3.122 Re, such as the very close initial location of flux tubes, the simulated ionospheric O + densities (panel 1 and 2 of Figure 6a) and the characteristics of O + number fluxes both in the ionosphere and around the apex of flux tube (panel 1 and 2 of Figure 6b). However, both the total ionospheric O + upflow and equatorial O + inflow in the simulation starting from L = 3.0 Re are much higher than those in the simulation starting from L = 3.122 Re, not to mention the distinct difference of equatorial O + densities (panel 1 and 2 of Figure 6c). As shown in Figure 5b, in simulation starting from L = 3.0 Re, the flux tube takes 2 hr 45 min to reach L = 4.05 Re, while in simulation starting from L = 3.122 Re the time duration is 1 hr 50 min. And in simulation starting from L = 3.5 Re, the flux tube drifts to L = 4.05 Re after 1 hr 35 min. From that point (L = 4.05 Re), both flux tubes reach geosynchronous orbit in 50 min at roughly the same speed, while in the simulation starting from L = 3.5 Re it takes 70 min. The results demonstrate that a longer duration of convection leads to a more prominent enhancement of equatorial O + density in flux tubes at sub-auroral latitudes (L < 4.0 Re), together with disturbed geomagnetic activity. However, under enhanced magnetospheric convection effects at high latitudes (L > 4.0 Re), the flux tube drifts faster toward dayside magnetopause with rapid volume expansion, which further causes decreasing equatorial density of ions. Then, the light ions (H + and He + ) quickly reach density equilibrium along the magnetic field line, maintaining their equatorial density, while the heavy ions (O + ) take a longer time to react against gravity or fall back to the underlying ionosphere. Although heating effects will drive ionospheric particles upward when the foot of the flux tube moves into high latitudes, it is not enough to compensate for early particle loss. In fact, besides flowing into the apex grid point of the flux tube, there are a large number of ionospheric O + particles trapped near the flux tube apex in the northern hemisphere, which probably contributes to later enhancement in equatorial O + density as shown in panel 1 of Figure 6c.
Figures 7a-7c present parallel velocity distributions for O + particles along the flux tube in simulations starting from L = 3.0 Re, L = 3.122 Re, and L = 3.5 Re, respectively. The parallel velocity distributions are calculated by integrating the velocity distribution over all the gyrospeeds by binning particles according to their positions along the flux tube. In this section, we focus on the dynamics in the northern hemisphere, aiming to discuss the SED supply of cold ions for plasmaspheric plume.
In the simulation starting from L = 3.0 Re (Figure 7a), the underlying ionosphere supplies O + particles for the plasmaspheric region and fills the flux tube for a few hours at the beginning of the simulation until the flux tube arrives at L = 4.0 Re (panel 1-6 of Figure 7a). Then though the equatorial concentration of O + particles tends to decrease due to rapid volume expansion of the flux tube, the snapshots show clear transport of O + particles along the flux tube toward the flux tube apex denoted by black lines particularly when L > 5.0 Re (panel 10-15 of Figure 7a), manifested as the change of binned density distributions at latitudes lower than 20°. On the other hand, in the simulation starting from L = 3.122 Re (Figure 7b), before the flux tube arrives at L = 4.0 Re, the initial density distribution around the flux tube apex is lower (compared with the snapshots in Figure 7a)  initial L shell of the flux tube, and the filling process seems weaker as well because of shorter convection time.
The snapshots do show a clear pile-up of particles near 30° magnetic latitude (e.g., panel 6-10 of Figure 7b), but there is no sign of field-aligned transport of particles toward the flux tube apex at L > 5.0 Re. It seems that the O + particles tend to move toward the flux tube apex but mostly get trapped around 20° magnetic latitude (black vertical solid line in each panel). In the simulation starting from L = 3.5 Re (Figure 7c), the flux tube drifts to L = 4.0 Re in a short time with a much lower density distribution at plasmaspheric altitudes along the flux tube, resulting in a continuous decrease of equatorial O + density from the beginning. The underlying ionosphere supplies numerous O + particles for the plasmasphere so that the density distributions in the ionosphere are much lower than in the other two flux tubes. Although the snapshots show the transport of O + particles toward flux tube apex at high L shells (third and fourth rows), it has minor impacts on equatorial O + density as shown in the right panel in Figure 6c. The comparison of parallel velocity distribution for O + particles reveals that the filling of dayside magnetospheric cold populations might also depend on middle-latitudinal ionospheric SED, in addition to plasmaspheric plume evolution. During a geomagnetic storm, there will be a fraction of upward ionospheric O + ions getting trapped near the apex of the flux tube which drifts sunward/poleward from the inner plasmasphere. Under later strengthened extrinsic force effects, these trapped particles have the chance to flow into the magnetic equatorial region and blend into the dayside magnetosphere. In this case, a single closed flux tube spent about 3 hrs drifting from L = 3.0 Re to L = 6.6 Re and the simulation presents a distinct equatorial filling of cold O + particles in the magnetosphere where L > 5.0 Re, which originated from the middle latitudinal SED plume.
Among the three simulated flux tubes, the magnitude of equatorial plasma density remains at the order of 10 2 cm −3 at high L shells, higher than the typical plasmaspheric densities calculated by the empirical model of Sheeley et al., 2001 (shown in the Supporting Information S1), which is consistent with being inside the plasmasphere. Through the plasmaspheric drainage plume as a convection channel, the flux tubes that originally drift on closed orbits inside the main plasmasphere jump to open passages, drifting sunward while bringing plasmaspheric plume sourced cold ions and additional ionospheric SED supplied cold heavy ions to dayside magnetosphere. The density distribution of light ions (H + and He + ) rapidly reaches equilibrium along the flux tube and their equatorial density approximately declines proportionally to the volume expansion of the flux tube (Figure 6c), while distinct density gradients of heavy ions (O + ) along the flux tube exist the whole time (Figures 7a-7c). It is no doubt that the ionospheric cold light ions (H + and He + ) can move along the flux tube and reach the flux tube apex in a short time. But the ionospheric cold heavy ions tend to be trapped at low latitudes along the flux tube at low L shells (<4.0 Re), with the potential to reach flux tube apex at high L shells (∼5.0 Re). The results demonstrate that the SED plume is able to supply cold O + particles to the dayside magnetosphere along a sunward drifting flux tube. Furthermore, if the SED plume was long-lived, it could be an additional plasma source for persistent plasmaspheric plume. In the simulated flux tube starting from L = 3.0 Re, the SED plume supplies cold O + particles for the dayside magnetospheric equatorial region and raises the equatorial O + density by a factor of 20% compared to the expected value. However, the cold light ions quickly arrive at the flux tube apex and reach density equilibrium along the flux tube so that their equatorial densities decrease with the volume expansion of the flux tube as a function of L −1 . The additional contribution from the mid-latitudinal ionosphere implies that there may be a denser concentration of cold O + particles in the dayside magnetosphere than previously thought, as these particles are difficult to be detected. Though light particles (H + and He + ) dominate magnetospheric ion populations, the cold heavy ions control the dynamics influenced by inertia such as dayside reconnection in the magnetosphere. Our simulations present the middle latitudinal SED plume to be an additional source region of magnetospheric cold heavy ions, which help to complete the known mechanisms of cold plasma filling in the dayside magnetosphere. The suggestion of a higher density of cold heavy ions would also give a hint about the explanation for the source of warm particles in the dayside magnetosphere.

Summary
Based on the DyFK model, we simulate the field-aligned dynamics within three independent closed flux tubes which drift toward dayside geosynchronous orbit from different initial L shells around MLT noon during the 2015 St. Patrick's Day storm. All the three simulated flux tubes have been verified to magnetically connect the SED plume with the plasmaspheric plume. The main results are listed as below: 1. The simulations show decreasing supply of cold ions (H + , He + , and O + ) from the ionosphere with increasing initial L shell location of the flux tube. However, the inflows of cold ions at the apex of the flux tube show 10.1029/2021JA030140 13 of 15 different trends depending on ion weights: the inflows of light ions (H + and He + ) increase with a larger initial L shell location of the flux tube while the inflows of heavy ions (O + ) present a decreasing trend. We attribute this discrepancy between light ions and heavy ions to synergistic effects of external forces including gravity acting on the ions 2. Despite the equatorial density of light ions (H + and He + ) as well as electrons decreasing as a function of L −1 , which agrees with the volume expansion of flux tube implying field-aligned density equilibrium, the equatorial density of heavy ions (O + ) shows different variations. The simulation displays that the ionospheric sourced cold light ions rapidly reach density equilibrium along the flux tube while cold heavy ions (O + ) do not present a similar behavior. In the simulation starting from L = 3.0 Re, the equatorial O + density even increases by a factor of 20% compared to the expected value based on the adiabatic expansion of the flux tube.
In combination with the simulated parallel velocity distribution of O + particles along the flux tube, we found that the cold O + particles tend to "pile up" near the flux tube apex in a single hemisphere while the flux tube is drifting at sub-auroral latitudes (L < 4.0 Re). Moreover, compared to the other two flux tubes, the flux tube starting from L = 3.0 Re spent a much longer time drifting at L < 4.0 Re, which causes many more ionospheric cold O + particles to pile up near the flux tube apex under gravitational effects. As the flux tube drifts into high latitudes, these "trapped" ions start flowing toward the flux tube apex because the gravitational force effects have been overcome by other enhancing forces. The simulations show a particular process of cold O + particles sourcing from mid-latitudinal ionosphere but flowing into dayside magnetosphere, indicating SED plume to be an additional source of magnetospheric heavy populations at dayside. This simulation study provides a new view of O + filling mechanisms in the dayside magnetosphere which suggests a denser concentration of cold O + distributed in the dayside magnetosphere As the "hidden particle source" in the magnetosphere, the cold ions play an important role in solar wind-magnetosphere interaction dynamics. Due to their low energy range which is below the detection threshold of instruments onboard spacecraft, magnetospheric cold populations present many questions that remain unsolved. In this study, we explored the filling mechanisms of magnetospheric cold ions from the Earth's plumes based on DyFK simulations. Beyond the impacts from the initial location of the flux tube before drifting toward dayside magnetopause from a closed orbit, there are more factors affecting the field-aligned transport of cold ions within plumes such as the geomagnetic activity, solar activity, etc. which will be discussed in the future.