Plasmasphere Refilling after the 1 June 2013 Geomagnetic Storm

Abstract

Plasma content and distribution are key parameters in the dynamics of the inner magnetosphere. The plasmasphere contributes, for the most part, to the plasma mass density, and its properties are very dependent on the history of the magnetosphere and geomagnetic activity. In this work, we investigated plasmasphere dynamics and plasmasphere–ionosphere coupling, focusing on the refilling process that followed the geomagnetic storm that occurred on 1 June 2013. The equatorial plasma mass density used to evaluate the refilling rates was remotely sensed by observation of the field line resonance (FLR) frequencies of the geomagnetic field, driven by ultra-low-frequency magnetic waves. The FLR frequencies were retrieved by performing an analysis of signals detected by several station pairs of the European quasi-Meridional Magnetometer Array. We estimated the rate at which the refilling process occurred, concentrating on both the diurnal and the day-to-day refilling rates. The estimated contraction rate during the main phase of the storm was higher than $\sim 3.5$ $R_E{d}^{-1}$, while the average expansion rate was $\sim 0.4$ $R_E{d}^{-1}$. We investigated the radial dependence of the refilling rates, using a novel approach based on fit plasma mass density profiles, and we related their variation to the plasmasphere boundary layer and the zero-energy Alfvén boundary. We found evidence supporting the idea that flux tubes mapping in the region between these two boundaries experience an enhanced refilling process.

1. Introduction

The plasmasphere is a toroidal region of cold plasma that encircles the Earth, and corotates with it. The source of this population is the ionosphere, the composition of which results from a balance of production and loss mechanisms, and from diffusion processes (see [1] and references therein). The ionosphere properties change drastically from the dayside to the nightside region, and only ion species having a small recombination rate are present during the nighttime. During the daytime, the net result is that the plasma streams outward, by ambipolar diffusion, along the magnetic field lines, populating the plasmasphere. Ambipolar diffusion is more efficient for lighter ions: thus, the principal constituents of the plasmasphere are H⁺ ions. After sunset, the ionosphere charge content strongly decreases, causing a downward flux of particles from the plasmasphere to the ionosphere (diurnal breathing).

The plasmasphere extends from 1000 km above the Earth’s surface to 4–6 Earth radii ($R_E$). The equatorial plasma density typically shows a decrease as a function of the geocentric distance ranging from $\sim {10}^4$ to 10 cm${}^{-3}$. The outer boundary—the plasmapause or plasmasphere boundary layer (PBL, [2])—is generally characterized by an abrupt decrease in the plasma density, and, in a steady state situation, separates the corotating plasmasphere from the plasmatrough, a region of tenuous plasma (1–100 cm${}^{-3}$) that generally streams toward the dayside magnetopause.

This differentiation is the net result of the superposition of the two main electric fields that drive the bulk motion of the cold plasma in the inner magnetosphere [3,4,5]: the first is the corotation electric field, $E_{cor}$, determined by the rotation of the Earth, surrounded by a conducting plasma, through the magnetic field [6]; the second is the convection electric field, $E_{conv}$, which is due to the penetration of the magnetosphere by part of the interplanetary electric field [7]. Representing the geomagnetic field as a rotating dipole, the corotation potential in the equatorial plane can be written as

$${\mathsf{\Phi }}_{cor}\left(L\right)=-\frac{\omega B_o{R}_E^2}{L},$$

(1)

where L is the McIlwain parameter [8], $\omega =2\pi /86,400$ s${}^{-1}$ is the angular velocity of the Earth and $B_0=30,000$ nT is the equatorial geomagnetic field at $1R_E$. The corotation potential is shown in Figure 1a: the corresponding electric potential, $E_{cor}$, is directed radially inward toward the Earth (black arrow), and the relative plasma streamlines are concentric with the Earth’s equator. The simplest equatorial convection model assumes a spatially uniform dawn-to-dusk electric field, e.g., [9]. The streamlines, therefore, are directed from the nightside to the dayside region (Figure 1b). Figure 1c shows the streamlines resulting from the superposition of the two electric fields in the equatorial plane. A separatrix (the magenta line)—also called the Alfvén layer or the zero-energy Alfvén boundary (ZEAB)—divides a region close to the Earth, in which $E_{cor}$ dominates and the streamlines are closed, and a region in which $E_{conv}$ is predominant, and the plasma streams towards the dayside magnetopause on open drift paths. In the dusk sector, there is a bulge with a stagnation point, where the total electric field vanishes (indicated by two opposite arrows): in this scenario, the ZEAB is usually interpreted as the plasmapause, and presents a pronounced dawn–dusk asymmetry. The corotating region is refilled with plasma by diffusion from the underlying ionosphere, following the diurnal breathing, while in the region external to the ZEAB, the flux tubes convect towards the magnetopause, where the plasma is lost.

The size of the ZEAB depends on solar wind conditions [9] and geomagnetic activity [10]. In particular, during geomagnetically disturbed periods, the enhanced convection causes an earthward movement of the ZEAB. Using this simple steady state model, Grebowsky [5] suggested that the plasmapause and the ZEAB do not necessarily coincide, due to the different characteristic time for changes in the dawn–dusk electric field, which is of the order of minutes [11], with respect to the time a convecting charged particle takes to reach the magnetopause, which is of the order of hours. In the few hours following a sudden enhanced convection—such as the onset of a geomagnetic storm—high-density plasma, formerly corotating with the Earth, drifts along open lines, until it is lost at the magnetopause, generating a high-density plasma region in the dayside magnetosphere (sunward surge). The nightside plasmapause position reaches the ZEAB position more rapidly, while in the afternoon sector a drainage plume forms and the plasmapause position reaches the new ZEAB position 15–20 h after the convection enhancement, assuming no further variation of the convection electric field. On the other hand, a large reduction of the convection usually happens during the storm recovery phases, possibly causing the formation of high-density corotating structures beyond the average plasmapause position.

The mechanism described above is only a simplified description of what really happens. The magnetosphere is not in a stationary configuration—it changes in shape and dimension, as the solar wind conditions change: the convection electric field changes accordingly, causing a variation of the pattern presented in Figure 1. Reinisch et al. [12] gave an overview of the main electric field empirical models widely used in plasmaspheric research, the most popular of which is probably the Volland–Stern (VS) convection potential, given by [13,14]:

$${\mathsf{\Phi }}_{VS}(L,\varphi )=-A{L}^{\lambda }sin\varphi ,$$

(2)

where $\varphi$ is measured counterclockwise from noon, A is a coefficient related to the electric field intensity, and $\lambda$ can be interpreted as a shielding parameter for lower L-shells—an effect possibly related to the ring current pressure [15,16]. For $\lambda =1$, there is no shielding, and the electric field is uniform and directed from dawn to dusk, as in Figure 1b. For $\lambda >1$, the electric field near the Earth decreases. Several works [13,14,17] have suggested that $\lambda =2$ is the value that best fits the observations. Maynard and Chen [18] included a dependence on the $Kp$ index in the A coefficient, using $\lambda =2$, and obtaining (in kV/$R_E^2$)

$$A\left(Kp\right)=\frac{0.045}{{(1-0.159Kp+0.0093K{p}^2)}^3}.$$

(3)

Another widely used model is the E5D [19], which describes the $Kp$-dependent convection electric potential (in kV) at the equator as

$${\mathsf{\Phi }}_{E5D}(L,\varphi )=\{L[0.8sin(\varphi -\pi )+0.2cos(\varphi -\pi )]+3\}\left[1+0.3Kr\right]H,$$

(4)

where $R_{ar}=9.8-1.4cos(\varphi -\pi )-(0.9+0.3cos(\varphi -\pi ))Kr$, $Kr=Kp/(1+0.1Kp$) and $H=1/[1+{(0.8R_{ar}/L)}^8]$.

The dynamics of the plasmasphere have been studied for decades, since the first observations of the PBL by Carpenter [20]. Measurements by the retarding ion mass spectrometer (RIMS, [21]), onboard the Dynamics Explorer 1 satellite, have boosted our knowledge of the plasmasphere’s ion composition. The Dynamics Explorer mission also provided direct measurements of field-aligned ion flows associated with refilling, e.g., [22]. More recently, the most significant advances in our understanding of plasmasphere dynamics have been achieved thanks to the Extreme Ultraviolet Imager (EUV, [23]) onboard the IMAGE satellite, which provided global pictures of the plasmasphere. An important contribution to the development of crucial concepts—such as erosion and refilling, as well as the PBL formation and properties—has been made by ground-based observations of electron density, obtained by very-low-frequency (VLF) whistler measurements, e.g., [24,25], and plasma mass density, inferred by ultra-low-frequency (ULF) signals recorded by magnetometers, e.g., [26,27]. Nowadays, it is well-known that after a period of geomagnetic activity the outer plasmasphere is eroded, due to the enhanced convection, and that the PBL can assume very complicated shapes [28]. During the recovery phase, the depleted flux tubes slowly recover, thanks to enhanced refilling from the ionosphere, e.g., [29]. Due to the dynamic nature of the magnetosphere, the refilling process is active most of the time, and characterizing in detail its spatial and temporal variations is of great importance to understanding the physical processes involved, and to consideration of the active role of cold plasma in wave propagation and wave-particle interactions [30]: basically, magnetic shells beyond $L=3.5$–4 spend most of the time in a state of refilling [25,26,31].

The refilling process can be analyzed both on a daily basis, considering only the upward plasma flux from the ionosphere during daytime hours, and on a day-to-day basis, considering the net change in the plasmaspheric density due to the combined effects of diurnal refilling and nighttime depletion, e.g., [24,32]. While the contraction of the plasmasphere takes place very quickly, e.g., [11,33]—even in half an hour—the day-to-day refilling timescale is much longer, and continues over a period of days, e.g., [34,35,36,37], provided no other convection enhancements occur. For very prolonged quiet periods, the plasmasphere can reach a state of saturation, in which the density variations are limited to diurnal breathing [38]. Recently, using EUV images, Gallagher et al. [39] described the plasmasphere erosion and the subsequent refilling as a kind of deep breathing, estimating the variation of the total mass content as a whole.

Gallagher and Comfort [31] discussed several aspects that remain open questions in our understanding of plasmasphere refilling, including the influence of thermospheric winds, ionospheric chemistry, photoelectron heating, mass dependence and the interaction of cold plasma with more energetic populations and plasma waves. Although extensive analysis has been done, focusing on both electron and mass density variations (e.g., [24,25,32,36,39,40,41,42,43,44]), quantitative observations are still needed, to improve and validate the existing models, e.g., [45,46,47,48,49]. In this work, we analyzed the refilling process that followed the storm event of 1 June 2013, using plasma mass density estimates inferred from ULF signals detected by the European quasi-Meridional Magnetometer Array (EMMA, [44]). We based our analysis on the assumption that PBL and ZEAB do not generally coincide [35]: we tested this assumption, and studied how the refilling process changes through time and space, especially in relation to the two boundaries. In particular, after the plasmasphere erosion that follows the geomagnetic storm, a decrease in the magnetospheric convection should determine a condition in which the ZEAB lies at larger geocentric distances compared to the PBL; the outward movement of the PBL should then be driven by the refilling rate. Other works have investigated the properties of the refilling rates at different L-shells, e.g., [32,35,36], but this is the first attempt to simultaneously and quantitatively relate the radial variation of the refilling rate to the PBL and ZEAB positions, proposing a novel approach that takes advantage of the EMMA latitudinal extension. The PBL position was determined by inspecting radial plasma mass density profiles derived from EMMA observations. Ground-based magnetometer observations, to monitor the PBL, have been successfully used in the past, by other authors, e.g., [50,51,52]. We adapted the methodology proposed by Del Corpo et al. [52], to retrieve PBL information in a more systematic and quantitative way, and we used fit radial profiles of the plasma mass density, to estimate the latitudinal dependence of the refilling rate. For the evaluation of the ZEAB position, we relied on the equatorial electric field empirical models introduced in Equations (2)–(4), and we tested their consistency with the EMMA observations.

Section 2 describes the method used to derive the plasma mass densities from the ground-based magnetic measurements, the estimation of the PBL and the procedure adopted to determine the ZEAB boundary. Section 3 shows the plasmasphere dynamics, and both the diurnal and the day-to-day refilling rates evaluated for the period under investigation. Section 4 presents a discussion of our outcomes, and their comparison to results from past works. Section 5 summarizes the main results.

2. Data and Methods

2.1. Plasma Mass Density

The event is part of the data set created by Del Corpo et al. [53], used to statistically describe the average plasma mass density variations [52,53] and plasma composition [54] for different geomagnetic conditions. A detailed description of the procedure used to infer the plasma mass density can be found in the references above: here, we briefly recall the main aspects. The procedure relies on the methods developed by Baransky et al. [55] and Waters et al. [56], which allow for detecting signatures of field line resonances (FLRs), by performing cross-spectral analysis of ULF magnetic signals simultaneously recorded by two magnetometer stations slightly separated in latitude. The FLR frequency detected is associated with the field line passing through the midpoint between the stations. The plasma mass density at the magnetic equator is determined by solving numerically the magnetohydrodynamic (MHD) equation for the toroidal mode [57,58], using the FLR frequency as a known parameter. In solving the MHD equation, the TS05 model [59] is used to describe the geomagnetic field, and the following power law dependence of the plasma mass density along the magnetic field line is assumed:

$$\rho ={\rho }_{eq}{\left(\frac{r_{eq}}{r}\right)}^m,$$

(5)

where r is the geocentric distance, $\rho$ is the plasma mass density and the suffix $eq$ indicates quantities at the magnetic equator.

Station pairs, located at different latitudes, provide information on the plasma mass density at different geocentric distances. Because of the latitudinal extension of EMMA, it is possible to monitor the equatorial plasma mass density for nominal L-shells ranging from 1.6 to 6.2 (values determined using the IGRF model [60]). For this event, we inferred the plasma mass density for 18 station pairs. The TS05 model realistically reproduces the variations of the magnetospheric current system, which—especially during disturbed geomagnetic periods—causes an extension of the upper radial limit sounded by EMMA up to 8 $R_E$ and beyond [53,61,62].

The region sounded is limited in magnetic local time (MLT), because the detection rate of FLRs at nighttime is relatively low [53,63]. The plasma mass densities inferred at nighttime are also questionable, because to solve the governing MHD equation, a perfect reflection of the ULF waves is assumed at ionospheric heights: in other words, the ionospheric conductivity is assumed infinite, but this approximation is not necessarily adequate during nighttime. If the ionospheric conductivity is sufficiently low, the frequency observed could be related to a free-end mode configuration, and a finite conductivity should be used in the boundary conditions for solving the MHD equation [64]. When only one of the footprints is sunlit, the FLR frequency observed could refer to quarter-wave modes [65,66].

Nonetheless, some recent works have suggested that the nighttime values obtained by the above procedure are not necessarily wrong. Takahashi et al. [67] observed a resonance compatible with a fixed-end mode configuration across midnight. More recently, Regi et al. [68] observed a nocturnal boundary crossing of the flux tubes monitored by EMMA, which was qualitatively compatible with total electron content observations from ground-based GNSS receivers: this suggests that, although plasma mass densities can be overestimated [69], nocturnal observations provide valuable information about plasmasphere dynamics.

In this work, we used nighttime density estimates to describe the plasmasphere dynamics, but we used only diurnal values to estimate the refilling rates. By "diurnal values", we mean points at which both footprints of the field line are lighted by the Sun. The footprints were evaluated at 120 km, the height at which the Pedersen conductivity maximizes. The procedure used to estimate the refilling rates is described in Section 3 below.

The proper power law index, m, which is employed in Equation (5), has been investigated by several authors. For low latitudes, we primarily refer to the work by Vellante and Förster [70], who suggested a value of 1 in the L-range 2.3–3.4; they also suggested that a higher value would be necessary for $L<2$. In the L-range 4–6, different works [71,72,73] have suggested values between 0.5 and 2. Menk et al. [26] investigated the mass density index, using ground-based observations of FLR harmonics in the L-range 2.6–4.6: they pointed out that m is a highly variable parameter, probably reflecting the plasma redistribution processes caused by the variation in the geomagnetic activity; they found, however, that the choice of m is not critical for the L-shells considered. Vellante et al. [74] tested the range [0, 2] as possible values for m in the L-range 2–6.5: they found differences smaller than the typical plasma mass density uncertainty, which, on average, is of the order of 30%, but which can vary between 10% and 50%, depending on L and the geomagnetic activity level. For the discussion above, and considering the L-range investigated in this paper, we considered it reasonable to use a value of $m=1$.

As discussed by Menk et al. [26], the evaluation of the plasma mass density described in this section was carried out by making several assumptions. The most obvious were that the magnetic field was accurately represented by the TS05 model, and that the plasma mass density followed the power law in Equation (5). There are also some inconsistencies between the MHD equation and the magnetic field model used [75]. We used a purely toroidal-mode MHD equation [57], in which background electric currents effects are ignored; however, there may be a coupling between toroidal and poloidal modes in reality, and the TS05 model used to solve the equation actually reproduces magnetic distortions, due to the previously mentioned currents [75]. Assuming that the FLR frequency detected referred to the field line passing through the midpoint between the stations was also an approximation. The actual location depends on the way the Alfvén speed varies, and on the resonance width [76]. In this study, the uncertainty of the final estimation of $\rho$ was determined by the uncertainty of the FLR frequency, which is usually considered the main source of error, e.g., [26,51]. Uncertainties due to the other approximations were ignored.

2.2. Solar Wind Parameters and Geomagnetic Activity Indices

Determination of the geomagnetic field by TS05 requires knowledge of the solar wind parameters and the $Dst$ index. We retrieved this information from the OMNI data available from the NASA/GSFC Space Physics Data Facility’s OMNIWeb service. We also used the $Kp$ index, available from the same data set, to determine the geomagnetic activity level during the geomagnetic storm.

2.3. Fit Profiles and PBL Inner Edge Evaluation

The radial position sounded by a given station pair changes with time. During disturbed periods, the equatorial region monitored by high-latitude pairs can vary by more than 2 $R_E$ in one rotation of the Earth. Again, we took advantage of the EMMA extension to obtain interpolated points at fixed positions. Following Del Corpo et al. [53], we fitted the radial profiles with at least 2 points with a smoothing spline. Fit portions with point separation larger than 2.5 $R_E$ were excluded from the analysis.

Figure 2a shows the EMMA observations (open circles) and their uncertainty (vertical bars) on 2 June 2013 (day 153) at 07:00 UT. The blue solid curve is the fit profile, and the dashed lines represent its uncertainty.

The fit density at a given time can be used to estimate the PBL inner edge position in the MLT sector sounded by EMMA [52]. Compared to electron density, mass density profiles exhibit a less steep variation across the PBL, mainly due to the presence of heavy ions [77]. Statistical works suggest that heavy ion mass loading is more effective in the dawn and pre-midnight MLT sectors [54,78]. Spatial integration effects [79] can also contribute to smoothing the profiles observed by ground magnetometers. The presence of the PBL can remove any FLR signatures in the cross-spectra [50,51]; this circumstance might prevent its identification: for example, in the profile shown in Figure 2a, two points are missing between 2.5 and 3.3 $R_E$, where a density drop occurred. Following the procedure adopted by Del Corpo et al. [52], the inner edge of the PBL can be determined by searching for a minimum ($r_{min}$) in the second derivative of the fit radial profiles. Likewise, the outer edge ($r_{out}$) corresponds to a maximum. Figure 2b shows the second derivative for the fit profile in panel a. An estimate of the maximum width of the PBL is $r_{out}-r_{inn}$ and, to a first approximation, we can reasonably assume that the real inner edge of the PBL is somewhere between $r_{inn}$ and the midpoint ($r_{mid}$) between $r_{out}$ and $r_{inn}$. In our example, we estimated that the PBL inner edge was between 2.4 and 2.9 $R_E$.

2.4. ZEAB Position Evaluation

To relate the refilling rate latitudinal dependence to the ZEAB position, we used the Volland–Stern–Maynard–Chen (VSMC) model (Equations (2)–(3)) and the E5D model (Equation (4)), which are both dependent on radial distance L and local time (LT). Figure 3 shows the equatorial electric potential maps obtained by the two models on 2 June 2013 at 01:24 UT, when $Kp$ was equal to 4.

The magnetopause position was estimated by the Shue et al. [80] model. In the VSMC model, the stagnation point $L_S$ is given by [81]

$$L_S^{VSMC}={\left(\frac{\omega B_o{R}_E^2}{2A}\right)}^{\frac{1}{3}}.$$

(6)

The ZEAB position for any LT can then be inferred by retrieving the position of the equipotential line passing through $L_S^{VSMC}$, which takes the value

$${\mathsf{\Phi }}_S^{VSMC}=-\frac{3}{2}\frac{\omega B_o{R}_E^2}{L_S^{VSMC}}.$$

(7)

We obtained a second estimate of the ZEAB from the E5D model. The stagnation point $L_S^{E5D}$ can be determined by calculating the coordinates for which the electric field $E=-\nabla ({\mathsf{\Phi }}_{cor}+{\mathsf{\Phi }}_{E5D})$ is null, where ${\mathsf{\Phi }}_{cor}$ and ${\mathsf{\Phi }}_{E5D}$ are given by Equations (1) and (4), respectively. The ZEAB for any LT is the equipotential line, having magnitude

$${\mathsf{\Phi }}_S^{E5D}={\mathsf{\Phi }}_{cor}\left(L_S^{E5D}\right)+{\mathsf{\Phi }}_{E5D}(L_S^{E5D},\varphi ).$$

(8)

In Figure 3, the estimated ZEABs are drawn as magenta lines.

3. Results

3.1. Observing the Plasmasphere Dynamics with EMMA

In this section, we discuss the plasmasphere dynamics in response to changes in the geomagnetic conditions, as observed by EMMA. The fit radial profiles obtained at different local times during a single Earth rotation can be merged in a single image, in polar coordinates, to obtain a two-dimensional representation of the density in the inner magnetosphere. Figure 4 shows a sequence of magnetospheric plasma configurations during the storm event under investigation.

Color is used to indicate plasma mass density on a logarithmic scale, as indicated by the color bar at the bottom. Magnetic noon is to the left, and shadowed areas indicate the plasma mass density values obtained for field lines having at least one footprint not lighted by the Sun. The $Dst$ and $Kp$ geomagnetic indices during this period are shown in the top panels. Before the storm, and after two days of very low geomagnetic activity ($Kp<1$), the plasmasphere was in a state of saturation (Figure 4c). Figure 4d shows the main phase of the storm and the beginning of the recovery phase. One can clearly see the formation of a drainage plume on the afternoon side. At the subsolar point, the magnetopause (magenta line) retreats to 8 $R_E$, in response to the interplanetary disturbance that impacted the magnetosphere. The magnetopause position was estimated using the Shue et al. [80] model. Early in the recovery phase (Figure 4e), the plasmasphere appears eroded up to 2–3 $R_E$. Panel f shows the plasmasphere configuration after three days: the replenishment seems to be effective only up to 4–5 $R_E$, and there is evidence of a bulge in the dusk sector.

Figure 5 shows in more detail the temporal variation of the equatorial plasma mass density, ${\rho }_{eq}$, at 2.5, 3.5, 4.5 and 5.5 $R_E$, over the period 30 May–7 June 2013 (Figure 5c–f).

The density at fixed geocentric distances was derived using the fitting procedure described in Section 2. Flux tubes having one or both footprints not lighted by the Sun are indicated by light and dark shadowed areas, respectively.

The storm started at the very beginning of day 152, marked by a sudden decrease of $Dst$ up to −119 nT. Through days 150–151, before the storm onset, a slight daily increase in the density profiles suggests that saturation condition was reached, at least for $r_{eq}\lesssim 4.5R_E$, while at 5.5 $R_E$, this trend is more pronounced. Such a small daytime increase is probably the signature of the gradual refilling by the ionosphere of the magnetospheric flux tubes that were partially depleted during the nighttime hours (diurnal breathing). Note also a day-to-day increase at 5.5 $R_E$, indicating that the flux tubes were still in a phase of recovery, following a previous event of high geomagnetic activity.

On day 152, during the main phase and the first stage of the recovery phase of the storm, the pattern is more confused, because of the rapid change in the magnetospheric configuration, and the appearance of complex structures, such as the drainage plume (see Figure 4d). Before the plume crossing—visible as a sudden increase in the density for $r_{eq}\gtrsim 3.5R_E$ near the terminator—a general density decrease is observed, which is more remarkable at 3.5 $R_E$.

On day 153, the effects of the erosion are visible at $r_{eq}\gtrsim 3.5R_E$. At these geocentric distances, the density had generally decreased by a factor of ∼4–5, with respect to day 151, while at 2.5 $R_E$, it had returned to the typical pre-storm level. The significant plasma depletion also gave rise to a more pronounced daytime refilling process at all radial distances.

In the following days, the geomagnetic activity progressively decreased, and the plasmasphere slowly recovered towards a saturated condition. At $r_{eq}=3.5R_E$, the day-to-day refilling started on day 153, while at $r_{eq}\gtrsim 4.5R_E$ the erosion proceeded, at least until day 154. The pre-storm conditions were approximately reached on day 156, at 3.5 $R_E$, while at higher geocentric distances, both diurnal and day-to-day refilling processes appear to have been still in progress on day 157.

At the end of day 157, a new storm event interrupted the plasmasphere recovery, and gave rise to a new erosion of the plasmasphere, which is clearly visible on day 158: although the geomagnetic indices would suggest that it was a minor disturbance, compared to that of day 152, a more intense erosion apparently took place. Inspection of the OMNI data highlights that the interplanetary magnetic field (IMF) was directed southward for almost 24 h across days 157 and 158, which was much longer than the 7 h interval characterizing the main storm on day 152: this could explain the larger plasma depletion observed on day 158. In addition, the MLT sector sounded by EMMA when the IMF turned southward may have influenced the observed dynamics [82,83].

3.2. Day-to-Day Refilling

Further information on the temporal–spatial variation of the plasma density can be derived from the analysis of the radial profiles for a specific local time sector. Figure 6 shows the radial profiles at 07:00 UT, approximately corresponding to 09:00 MLT (on average, the MLT sector sounded by EMMA corresponds to UT + 2) on days 151, 153, 155 and 157.

We have omitted the profile of day 152, which corresponded to the highly dynamic main phase of the storm. The profiles of days 154 and 156 have also been omitted, as they were very similar to the profiles of days 155 and 157, respectively. The solid curves are smoothing spline curves drawn through the measurements (open circles).

The density profiles show a dramatic change from day 151 to 153, with the formation of a plasmapause around ∼2.4–$2.9$$R_E$, as determined by the procedure described in Section 2.3. On the following days, the plasmapause moved outward as the recovery proceeded; however, as already stated above, the saturation condition tentatively indicated by the quiet-day density profile of day 151 (black line) was reached on day 157 only below 3.5 $R_E$. The assumption that day 151 was representative of a saturated plasmasphere is justified by the prolonged (>2 d) quiet conditions of the magnetosphere in the preceding days, as indicated by the magnetic indices. For at least two days, $Kp$ was less than 1, and no remarkable variation of the $Dst$ index was observed. A further confirmation comes from the EMMA plasmasphere model [52], which is a good representation of a saturated plasmasphere, and is reported as a magenta line in Figure 6. The model is LT-dependent, and derived between 2.3 and 4.5 $R_E$: extrapolated up to 7 $R_E$, it seems to match, reasonably well, the experimental points.

Using the procedure described in Section 2.3, we estimated the PBL inner edge position from each radial profile in the interval 03:00–13:00 UT (roughly 05:00–15:00 MLT) for the entire recovery phase (days 153–157). The average values assumed by $r_{inn}$ and $r_{mid}$ in each day are reported in

Table 1. Average PBL inner edge position in the interval 03:00–13:00 UT (05:00–15:00 MLT).

	Day 153	Day 154	Day 155	Day 156	Day 157
${\mathit{r}}_{\mathit{i}\mathit{n}\mathit{n}}$ [${\mathit{R}}_{\mathit{E}}$]	2.5	2.9	3.0	3.4	4.0
${\mathit{r}}_{\mathit{m}\mathit{i}\mathit{d}}$ [${\mathit{R}}_{\mathit{E}}$]	3.0	3.7	3.7	4.4	4.9

.

These values can be used to obtain information about the timescale of the contraction and expansion of the plasmasphere. We have no information about the plasmapause position before the storm, but we can argue that it was beyond $7R_E$, if any was present—at least until 17:30 UT, the time of the last saturated radial profile observed by EMMA. After about 30 h, a plasmapause appears in the radial profiles near 2.5 $R_E$, and remains clearly visible for the entirety of day 153 (e.g., the blue curve in Figure 6): on that day, the PBL was, on average, in the range ∼2.5–$3.0$ $R_E$, which implies that the contraction took place at an average rate of at least ∼3.5 $R_E{\mathrm{d}}^{-1}$. The recovery took place over several days. On day 157, the day of maximum expansion before the occurrence of the new storm, the plasmapause was at ∼4.0–$4.9$ $R_E$: this suggests an average rate of expansion of approximately ∼0.4$\pm 0.2$$R_E{\mathrm{d}}^{-1}$, which is an order of magnitude slower than the contraction rate. Inspection of Table 1 suggests that the expansion did not take place in a linear way: for example, the PBL remained stationary from day 154 to 155, while on the other days it moved outward at a rate of ∼0.6 $R_E{\mathrm{d}}^{-1}$.

Figure 7 shows the day-to-day density variation at 09:00 MLT (the red open circles) and 12:00 MLT (the blue open circles) at different geocentric distances (Figure 7d–g).

The error bars represent the uncertainty due to the fitting process described in Section 2. Figure 7a,b re-propose the geomagnetic indices, and Figure 7c shows the ZEAB position, determined from the E5D (green) and VSMC (brown) models in the MLT sector sounded by EMMA (∼UT + 2). Both models predict an earthward motion of the ZEAB during the main phase, and an outward motion during the recovery phase. The VSMC model predicts a bigger reduction of the ZEAB radial distance in the morning–noon sector, while in the dusk sector—visible as diurnal peaks—the models are quite similar.

At 2.5 $R_E$, the difference between the 09:00 MLT and 12:00 MLT density values was significant only on days 152 (during the main phase) and 153, when the plasmapause was at ∼2.5 $R_E$. To allow for an erosion at 2.5 $R_E$, the ZEAB needed to be, for some time, below this geocentric distance. This condition is satisfied only by the VSMC model, which approaches 2 $R_E$ during the main phase, while the E5D model never drops below 4 $R_E$ in the entire period: this is in contrast to the EMMA observations that clearly suggest an erosion below that radial distance. On the contrary, the radial profiles in the noon sector during the main phase (not shown) revealed a density drop at around 1.7–1.9 $R_E$, compatible with a PBL, although the plasmasphere erosion was less evident, due to other effects, such as the presence of the plume.

Accordingly, we can assume that the VSMC better describes the ZEAB motion, at least for the event under investigation. The flux tube mapping at 2.5 $R_E$ was between the plasmapause and the ZEAB, and the recovery took place in, at most, one day. At 12:00 MLT, the day-to-day refilling can be evaluated as the difference between the density on day 153 and the density on day 152, giving $\Gamma =430\pm 150$ amu cm${}^{-3}$ d${}^{-1}$.

In the following days, the plasmapause moved progressively outward, reaching ∼3.1–3.5 $R_E$ on day 156. At 3.5 $R_E$, the density at both local times reached the minimum value on day 153, and then approximately recovered to the pre-storm conditions on day 156, when the flux tube was very close to the plasmapause. In the same period, the ZEAB predicted by the VSMC model during daytime (the non-shaded area in Figure 5) moves outward, from 3.5–4.5 $R_E$ to 4.5–5.5 $R_E$.

At 4.5 $R_E$ (Figure 7e), recovery took place between days 154 and 157, for both local times. The recovery was only partial, because it was interrupted by a new disturbance on day 158.

The flux tube mapping at 5.5 $R_E$ (Figure 7d) shows only small signs of recovery, mainly on days 156 and 157—behavior compatible with the ZEAB derived from the VSMC model, which, as observed above, is below that distance most of the time.

We evaluated the day-to-day average refilling rates, $\Gamma$, by linear-fitting the points from day 153 to day 156, for both local time sectors, at 3.5 $R_E$, and from day 154 to day 157, for both local time sectors, at 4.5 and 5.5 $R_E$. The linear fit is a simple but effective way to obtain information on the average refilling rate, and has also been used by other authors in the past, e.g., [41]. However, as pointed out in this section, the recovery might not happen in a linear way. The results are reported in

Table 2. Day-to-day refilling rates.

${\mathit{r}}_{\mathit{e}\mathit{q}}$ [${\mathit{R}}_{\mathit{E}}$]	$\Gamma$ [amu cm${}^{-3}$ d${}^{-1}$]	MLT Sector	Source
3.5	$92\pm 28$	09:00	our observations
4.5	$18.6\pm 4.8$	09:00	our observations
5.5	$1.27\pm 0.43$	09:00	our observations
2.5	$430\pm 150$	12:00	our observations
3.5	$97\pm 28$	12:00	our observations
4.5	$25.7\pm 7.1$	12:00	our observations
5.5	$3.99\pm 0.86$	12:00	our observations
2	650	15:00	Chi et al. [43]
2.3	$426\pm 162$	11:00	Obana et al. [32]
2.6	$423\pm 104$	11:00	Obana et al. [32]
3.7	$83\pm 20$	12:00	Lichtenberger et al. [44]
4.1	$79.9\pm 7.4$	11:00–12:00	Dent et al. [36]

, together with estimates from previous studies.

3.3. Diurnal Refilling

The day-to-day refilling process discussed in Section 3.2 is directly linked to the refilling that occurs during daytime hours. In this section, we discuss the variation of the diurnal refilling, according to latitude and geomagnetic activity.

The average diurnal refilling rate, $\gamma$, was obtained by linear-fitting the diurnal profiles of the plasma mass density at a given geocentric distance. We used only data where flux tubes had both footprints getting the sunlight at ionospheric heights: those data lie outside the shadowed areas in Figure 5. Only daily profiles with at least five points were fit.

Figure 8 shows the values of $\gamma$ obtained for the entire period, and for different geocentric distances.

The error bars represent the uncertainty due to the fitting process. The values assumed by $\gamma$ on day 151 can be interpreted as the typical refilling rates for saturated conditions, at least for $r_{eq}\lesssim 4.5R_E$.

On day 152, $\gamma$ reflects the high variability of the main phase. The negative values observed for $r_{eq}\le 3.5$$R_E$ were not due to a flux directed towards the ionosphere. They were, rather, the result of the action of the convection during the main phase of the storm and the early stage of the recovery phase, in which the plasmasphere was progressively eroded during the day, leading to a decrease in the plasma density over time. Conversely, the crossing of the plume in the afternoon sector affected the estimates for $r_{eq}\ge 4.5$$R_E$ (the full red dots in Figure 8).

The general behavior of $\gamma$ during the recovery phase strictly reflects the behavior of $\Gamma$. As noted in Section 3.2, the picture of a diurnal refilling that is enhanced in flux tubes mapping to the region between the PBL and the ZEAB seems to be confirmed: $\gamma$ enhancements are clearly visible on day 153 at 2.5 $R_E$, on days 154–155 at 3.5 $R_E$ and on days 155–157 at 4.5 $R_E$. After the enhancements at 2.5 and 3.5 $R_E$, the refilling rate assumes approximately the pre-storm values, confirming that conditions close to saturation were reached again. At 4.5 $R_E$, the plasmasphere recovery was interrupted by the disturbance on day 158.

The positions of the enhancements are strictly related to the PBL and the ZEAB: on day 153, the peak is present only at 2.5 $R_E$, when the ZEAB was close to 3.5 $R_E$; on day 154, the ZEAB was between 3.5 and 4.5 $R_E$; in the following days, it moved between 4.5 and 5.5 $R_E$, until day 157.

At 5.5 $R_E$, an increase in $\gamma$ is observed mainly on days 156 and 157. At this geocentric distance, the corresponding flux tubes were beyond the ZEAB for almost all the recovery phase, and approached the ZEAB only on days 155–156. This was not true for days 150 and 151, when the PBL and the ZEAB probably coincided, and were well beyond 5.5 $R_E$.

4. Discussion

The strong L dependence of the plasmasphere recovery after a period of enhanced convection is a well-known aspect, e.g., [25,26,31,41], and is confirmed also by our results. In the storm presented, starting from a saturated condition, the plasmasphere was eroded up to $2.5R_E$, and after four days it recovered to pre-storm conditions only below $3.5R_E$. An accepted explanation is that flux tubes closer to the Earth, i.e., those mapping to low latitudes, need less time to reach diffusive equilibrium than those mapping at high latitudes: this is mostly due to the volume of the flux tubes, which approximately increases as $L^4$ [32]. Recovery at high L-shells is also slowed by minor enhancements in the magnetospheric convection, driven by variations in the solar wind conditions during the recovery phase [26,31]: this latter aspect is probably what caused the non-linear expansion of the plasmasphere during the recovery phase (see Table 1). We found that the plasmapause remained stationary from day 154 to 155. During the night, between days 154 and 155, the IMF turned intermittently southward for about 10 h. Even if the GSM-z component of the IMF never fell below −3 nT, an enhancement of the magnetospheric convection probably occurred, visible as a small increase in the $Kp$ index during the same time interval: this could explain the pause observed in the plasmasphere recovery.

We focused on the idea that the ZEAB and the PBL do not coincide during the recovery phase of a storm [5,35]. Indeed, the ZEAB responds very rapidly to changes in the convection electric field, while the plasmapause moves outward on timescales determined by the refilling process.

Looking at the day-to-day average refilling rates summarized in Table 2, no local time difference appears at 3.5 $R_E$, while at 4.5 $R_E$ the refilling rate at 12:00 LT is 40% higher than the one obtained at 09:00 MLT. At 5.5 $R_E$, the increase is 310%. This confirms that the value of the refilling can be significantly dependent on the choice of the local time used for the calculation [36].

A few authors have evaluated day-to-day refilling rates, using ground-based FLR measurements: their results are also summarized in Table 2. The value of Dent et al. [36] was evaluated by us, by linear-fitting the data reported in Table 2 of their paper. Except for the work of Chi et al. [43], all the other rates were estimated around the noon local time sector, and can be directly compared with our estimates at 12:00 MLT. The comparison is shown in Figure 9, where also the measure by Chi et al. [43] is reported.

We can see good agreement between the data. The general trend is very well fit by the expression (black line in Figure 9):

$${\Gamma }_{\rho }\left(r_{eq}\right)={10}^{(4.09\pm 0.16)-(0.592\pm 0.045)r_{eq}}.$$

(9)

Figure 9 also shows similar trends obtained by other authors from derived measurements of the equatorial electron number density, n; therefore, the refilling rates have the dimension of cm${}^{-3}$ d${}^{-1}$. Denton et al. [41], using observations from the IMAGE RPI instrument [85], derived the model

$${\Gamma }_n\left(L\right)={10}^{2.22-0.006L-0.0347L^2},$$

(10)

based on many events in the year interval 2001–2005 (the magenta line). Krall et al. [84], by averaging refilling rates at dawn and dusk for a single event, obtained the following exponential dependence (the green line):

$${\Gamma }_n\left(L\right)=3.81{\left(\frac{6.8}{L}\right)}^{4.94}.$$

(11)

A non-local perspective was adopted by Gallagher et al. [39], who, by investigating seven events from IMAGE EUV observations, modeled with two different trends the radial dependence below and above $L=4.5$ (cyan lines):

$$\begin{array}{ccc}\hfill {\Gamma }_n\left(L\right)& =& {10}^{-(0.62\pm 0.29)+(0.66\pm 0.10)L}(L<4.5),\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill {\Gamma }_n\left(L\right)& =& {10}^{(4.50\pm 0.40)-(0.51\pm 0.08)L}(L>4.5).\hfill \end{array}$$

(13)

Our model best agrees with the Krall et al. [84] profile, which was obtained by a procedure similar to the one adopted by Denton et al [41]: the main difference between the two is that Denton et al. [41] used 34 time periods (with two MLT sectors considered for each of them), whereas Krall et al. [84] considered only one event. Denton et al. [41] reported high variability in the observed refilling rate, highlighting also a solar cycle dependence. In particular, at the solar maximum, when the event analyzed by Krall et al. [84] occurred, they found a larger value for $L<4$. The highest refilling rates observed by Denton et al. [41] between 2 and 3 $R_E$ were comparable with the ones observed by Krall et al. [84] in the same range: thus, it is not surprising that the radial variation observed for a single event can differ substantially from the average trend.

Considering all the curves, they show comparable values only around 3.5–4 $R_E$, which is also the average PBL position observed during the recovery phase of a storm (see, e.g., Carpenter and Anderson [38] and Del Corpo et al. [52]). Below 3.5 $R_E$, the model by Denton et al. [41] and, more evidently, the one by Gallagher et al. [39] show lower values compared to our results. In the derivation of both models, the refilling rates at each L were evaluated over the same time interval, without considering whether the shell spent most of the time inside the plasmasphere or not. On the other hand, our estimates were made over an L–dependent time interval, determined considering the minimum and maximum density observed at a fixed radial distance during the recovery. Moreover, from the plasmasphere dynamics observed, all the distances analyzed were supposed to be outside or in the proximity of the PBL during such time intervals. A similar approach was also adopted to estimate the other rates reported in Table 2, that were used to derive Equation (9).

Even though a density decrease inside the eroded PBL has been reported in several studies (see Gallagher et al. [39] and references therein), the magnitude of the depletion is much lower with respect to the regions external to the PBL, and certainly negligible for shells located deep inside the plasmasphere. Shells inside the PBL presumably recover faster than the ones located well outside the boundary. An evaluation of the refilling rate over an interval longer than the one needed to reach saturation inevitably results in an underestimation, and this effect becomes more pronounced when considering inner shells. Similarly, the refilling rates evaluated for shells located outside but very close to the PBL, for which the recovery proceeds faster than more external shells, might be underestimated: for example, if we had used the day-interval 152–157 to evaluate the refilling rate at 2.5 $R_E$, we would have found the value $0\pm 60$ amu cm${}^{-3}$ d${}^{-1}$. By comparison, at $L=2.5$$R_E$, the model by Gallagher et al. [39] gives ∼11 cm${}^{-3}$ d${}^{-1}$, but they reported a case in which the rate was as low as ∼1 cm${}^{-3}$ d${}^{-1}$ in the range 2.0–2.5.

The effect described above may explain the discrepancy between our observations and the models by Denton et al. [41] and Gallagher et al. [39] for $L<3.5$. Another aspect to be considered is that the ratio of mass density to electron number density is generally close to 1 in the plasmasphere, but outside the plasmapause higher values are expected, especially during the recovery phase, e.g., [54]: this would suggest higher values of ${\Gamma }_{\rho }$ with respect to ${\Gamma }_n$, as is the case in this L-range. The higher values found by Denton et al. [41] with respect to Gallagher et al. [39] might be related to the higher number of events considered, but also to the different approaches used to derive the equatorial densities. Equations (12) and (13) were derived by evaluating the integrated plasmaspheric content across all MLT sectors for different L-ranges. The mass refilling rates evaluated were eventually translated to average equatorial electron density refilling rates. This approach provided a non-local view of the refilling process that should have prevented any azimuthal transport variation affecting the refilling rates estimation: for example, the not-exact corotation of the plasmasphere with the Earth, which has often been reported (e.g., [28,86]) could have increased the indeterminacy of the method presented in this work, because of the limited ability of EMMA to distinguish between density variations due to ionospheric refilling or transport. The local approach is also a limitation for the other models presented, which rely on in situ satellite observations that are generally unable to separate temporal and spatial variation of the plasma.

From the discussion above, it is reasonable to expect that the refilling rate decreases with L for shells that spend most of the time outside the PBL: basically, all the observations used to derive our model are in this category; accordingly, a decreasing trend is found. Furthermore, the other models show similar trends beyond $L=4$ (i.e., beyond the expected average PBL position), but the rates found by Gallagher et al. [39] are one order of magnitude greater than the others. Future works should focus on understanding whether the differences in the methods used to derive the models can account for this and for the other differences discussed above.

Further information on the refilling process and its relation with the PBL and ZEAB is provided by the analysis of the radial variation of the diurnal refilling rate, as summarized in the left panels of Figure 10, where all the recovery phase is presented.

The average PBL position, as derived in Table 1, is represented as light blue stripes. The average diurnal ZEAB positions, determined in the same time interval and on the same days, from the VSMC and E5D electric potential, are represented as red and green solid lines, respectively. From the comparison of the ZEAB position with the observed plasmasphere dynamics, we argue that, for the present study, the ZEAB position was better described by the VSMC model, but we reported results from both models, for comparison. The refilling rates obtained at different geocentric distances refer to flux tubes with very different volumes: hence, identical upward plasma fluxes from the ionosphere lead to very different equatorial refilling rates at different L-values. To take into account this effect, we also evaluated the normalized diurnal refilling rate, ${\gamma }_n$. For each day, the data were first normalized to the plasma density value at 12:00 MLT, and then were fit linearly, to obtain ${\gamma }_n$. This procedure was repeated for all the geocentric distances in the range 2–7 $R_E$ at 0.1 $R_E$ step.

The radial variation of ${\gamma }_n$ is shown in the right panels of Figure 10. Focusing on the left panels, one can see a general decrease in the observed refilling rate with increasing geocentric distance, although a clear peak appears near 2.6 $R_E$ on day 153, and moves toward higher geocentric distances as the recovery proceeds. The peak value was $\gamma \simeq 100$ amu cm${}^{-3}$ h${}^{-1}$ at $r_{eq}\simeq 2.6R_E$ on day 153, and it gradually decreased in the following days, up to reach $\simeq 20$ amu cm${}^{-3}$ h${}^{-1}$ on day 157 at $\simeq 3.9R_E$.

The connection of such enhancements with the dynamic of the plasmasphere can be better explained by comparing the radial profiles of the normalized diurnal refilling rate (the right panels of Figure 10): ${\gamma }_n$ generally assumed small values at small geocentric distances (saturated condition), and increased with increasing radial distance, showing a peak at a geocentric distance that increased with increasing time; as opposed to the $\gamma$ profiles, the magnitude of the ${\gamma }_n$ peak is comparable from one day to another; the enhancements clearly took place beyond the plasmapause, and generally before the predicted ZEAB position. This last condition is less clear for day 156, and was definitely not satisfied on day 157, when the ZEAB moved to lower L-shells. A possible explanation for this discrepancy is that the enhancement of the convection took place when EMMA was already beyond the midnight sector, as confirmed by the $Kp$, which did not exceed $1^{+}$ until 03:00 UT of day 157, when EMMA was in the 05:00 MLT sector. As described by Chappell [82], the plasmapause observed during daytime hours is dependent on the magnetic activity changes that occurred even hours before, when the same flux tube transited the midnight sector. According to this, EMMA is able to promptly see variations due to a convection enhancement only if it occurs when EMMA is near the midnight sector, as happened for the main disturbance of this event study, or for the small enhancement during the night between days 154 and 155. The 8 h average values of the ZEAB across midnight, between days 156 and 157, are reported in Figure 10 as dashed lines. The two models predict basically the same position, correctly placing the ZEAB beyond the PBL.

Past estimates of the diurnal refilling rate from FLR detection were conducted by Chi et al. [43], who found a value of 200 amu cm${}^{-3}$ h${}^{-1}$ at $L=2$; they found no difference in the value of $\gamma$ evaluated before and after a strong geomagnetic storm, for which the PBL was supposed to move earthward up to 1.8 $R_E$. The latitudinal dependence of $\gamma$ for a selected day during the recovery phase of a moderate storm was studied by Obana et al. [32], who found $\gamma =248$, 110, 39 and 13 amu cm${}^{-3}$ h${}^{-1}$ for the L-shells 2.3, 2.6, 3.3 and 3.8, respectively. The decreasing trend beyond 2.5 is compatible with the radial dependence showed in Figure 10a. We did not observe, however, the refilling rates found at $L=2$–$2.3$ in both studies, which were 2–2.5 times higher than the highest value observed in our analysis, which was $\sim 100$ amu cm${}^{-3}$ h${}^{-1}$ at 2.4 $R_E$ (Figure 10a): a possible explanation for this discrepancy could be the different state of the magnetospheric plasma after the erosion, which is, in general, different from storm to storm; the time elapsed from the start of the refilling process and the time of the refilling rate evaluation may also have affected the final estimate. Indeed, as is clear from Figure 8, $\gamma$ strongly depends on the phase of the storm, and can change by a factor of $\sim 10$ at 2.5 $R_E$ (see, for example, the change in the day-interval 153–155 in panel g). The difference in magnitude form one day to another becomes less important for external L-shells, which have higher volumes to refill and generally low plasma densities at the equator.

It is worth noting that the upward flux from the ionosphere is not the only source that can increase the refilling rates. For example, Gallagher and Comfort [31] suggested that cold plasma observations, like the one presented in this study, may actually be altered by the presence of the warm plasma cloak [87], a population reach of O⁺ ions with energies of 10 eV to 3 keV, typically observed in the morning sector, outside the plasmasphere. The scenario of a plasmapause determining an enhancement of the refilling process, as observed from plasma mass density, could confirm the idea of the significant contribution of heavy ions during the recovery phase of a storm [54]. From this point of view, the diurnal variation of the plasma mass density, and the related refilling rate, could be the result of two concurrent processes: the upward plasma flux provided by the particle diffusion from the ionosphere directly through the flux tube (the proper refilling process), and the injection of O⁺ ions that originate in the polar regions and accumulate near the plasmapause [78].

5. Conclusions

We have presented a study of the plasmasphere refilling in the aftermath of the geomagnetic storm that occurred on 1 June 2013. The storm was preceded by a prolonged period of extremely low geomagnetic activity, in which the plasmasphere was in equilibrium with the ionosphere.

We evaluated the diurnal plasma refilling rate and the day-to-day refilling rate for different geocentric distances and different phases of the storm. We have presented a novel approach, based on fit plasma mass density radial profiles, to evaluating the refilling rates, and to relating their radial variations with the PBL and ZEAB positions.

The enhanced diurnal refilling directly caused the day-to-day refilling, and ultimately the plasmasphere recovery. The main aspects can be summarized as follows:

• Inside the PBL, the flux tubes were not significantly depleted during the main phase, resulting in small refilling rates in all of the recovery phase.
• The region between the PBL and the ZEAB consisted of depleted flux tubes that corotated with the Earth. As suggested by Obana et al. [32], the resulting upward plasma flux during daytime was enhanced. As a result, the recovery of the plasmasphere took place principally in the region between the PBL and the ZEAB, and proceeded at an ever-increasing distance, until a new disturbance occurred.
• Outside the ZEAB, the plasma supplied by the ionosphere could not be sufficiently trapped, and was lost through convection towards the magnetopause. As pointed out by Denton et al. [42], the observed density variation may not be the result of refilling of a particular flux tube. This is especially true well outside the ZEAB, where the plasma typically does not corotate with the Earth, and what we see is the plasma density sampled at different moments on different drift paths. Strictly speaking, the refilling concept is not applicable to the region outside the ZEAB, with the exception of the post-dawn sector, when the flux tubes convecting from the nightside might corotate with the Earth [35], and be filled with ionospheric plasma for a few hours before they are lost to the magnetopause, e.g., [88]. However, it is worth noting that because of the plasma loss, an analysis of this region will result in an apparently smaller diurnal refilling rate.

The contraction rate due to convection, and the expansion rate that followed in the recovery phase, were also evaluated: the contraction rate was higher than ∼3.5 $R_E{d}^{-1}$ while the expansion rate was estimated to be ∼0.4 $R_E{d}^{-1}$.

We compared our observation of the day-to-day refilling rates with past estimates made of plasma mass density by several authors, finding good agreement. The radial variation of the day-to-day refilling rate is well-fit by the following expression:

$$\Gamma \left(r_{eq}\right)={10}^{(4.09\pm 0.16)-(0.592\pm 0.045)r_{eq}}.$$

Comparison to similar models obtained from equatorial electron number density observations shows only partial agreement, underlining the necessity to cover more extensively this topic, both experimentally and theoretically.