Early efforts to model the plasmaspheric electron density included effects due to solar activity and season, with the first models providing the density using simple empirical relations depending on the McIlwain parameter L in Earth radii. Carpenter and Anderson (1992) derived a “reference profile” of the plasmaspheric electron density, valid for 2.25 < L < 8, to describe the saturated plasmasphere, N L = 10 − 0.3145 L + 3.0943 , with additional dependences to include perturbations due to season and phase of the solar cycle. This model uses the International Sun-Earth Explorer (ISEE) measurements and is limited to the local time interval of 0–15 MLT omitting plasma expansion on the dusk side, which these authors had aimed to treat separately (Carpenter and Anderson, 1992). Lyons and Thorne (1973) used a form N L = 1000 4 / L 4 , that was consistent with Carpenter et al. (1964), to derive electron lifetimes that yielded equilibrium flux profiles for L = [1,5] using a Fokker-Planck radial diffusion code, even though the density model was not valid below L = 2. Albert (1999) used an exponential, N(L) = 16400e^−0.875L, instead of a power law in L. Sheeley et al. (2001) reported N(L) = 1390 (3/L)^4.8 ± 440 (3/L)^3.6 in the plasmasphere for 3 ≤ L ≤ 7, where the authors show that the standard deviation captures differences between a newly filled and saturated plasmasphere. They did not find a magnetic local time (MLT) dependence for the plasmasphere, but did model the MLT-dependence of the plasma trough. A combination of Albert (1999) within the plasmasphere and Sheeley et al. (2001) within the plasma trough is used in the wave-particle interaction simulations of Ripoll et al. (2017) when satellite observations are lacking.

Gallagher et al. (2000) developed the Global Core Plasma Model (GCPM), a single unified model of the whole plasmasphere using an ‘amalgam’ of previously developed ‘region-specific’ models. GCPM addresses the density, temperature and composition of the plasmasphere, plasmapause, trough and polar cap. It depends on solar and geomagnetic indices, but is intended to be ‘representative’ of these conditions rather than used as a dynamic model. GCPM uses a modified version of the reference profile of Carpenter and Anderson (1992), N(L) = 10^−0.79L+5.3, added to the perturbations due to the solar cycle and season. It joins the topside ionosphere model of IRI to the equatorial plasma density model by first extrapolating the slope of the IRI model above the F2 peak using an exponential function and extrapolating the slope of the equatorial model downward in altitude with another exponential function, then blending the two functions with hyperbolic tangents. At higher latitudes, the shape of the exponential function is determined from IRI above the F2 peak, but the form is shifted by a constant so that the exponential decays to the equatorial value. The plasmapause location and width depend on local time. GCPM could be considered as the best compilation of all empirical density models. However, the GCPM model has not been directly coupled to radiation belt codes or wave particle interactions codes (to the knowledge of the authors) but it has been used for the validation of other plasma density models, themselves used in radiation belt codes (e.g., Ozhogin et al., 2012).

There have been recent efforts to model the variation of electron density with magnetic latitude. Denton et al. (2006) used satellite measurements from Polar and the Combined Release and Radiation Effects Satellite (CRRES) to model the latitudinal variations as a power law of the radial distance R to any point along the field line, N L , R = N e q L R E / R α , and fit α as a function of L and equatorial density for L > 2. Denton et al. (2006) fit this model to IMAGE RPI data and found that it did not perform well at high magnetic latitudes. Reinisch et al. (2004) and Huang et al. (2004) found that the following form fits data from IMAGE RPI well: N L , λ = N e q L 1 + γ λ / λ I N V cos π / 2 α λ / λ I N V − β (1) where λ is the magnetic latitude, λ I N V is the invariant magnetic latitude, and the fitting parameters are α , β , γ . Ozhogin et al. (2012) built on the earlier work of Reinisch et al. (2004) and fixed the values of the fitting parameters to α = 1.01 ± 0.03 , β = 0.75 ± 0.08 , γ = 0   a n d   N e q L = 10 − 0.4903 ± 0.0315 L + 4.4693 ± 0.0921 (2)

The Ozhogin et al. (2012) model is restricted to altitudes greater than 2000 km and L > 1.5, up to L = 4, and does not address dependence on MLT, season, solar activity, or differences in density between the Northern and Southern hemispheres. This model is plotted in Figures 1G, H to illustrate the increase of density with latitude. Models of Carpenter and Anderson (1992), Gallagher et al. (2000), Denton et al. (2006), Sheeley et al. (2001) and Ozhogin et al. (2012) are compared in Figure 8 of Ozhogin et al. (2012).

Empirical and first-principles models of the cold plasma density in the plasmasphere have been in development since the 1960’s, but there is just one model (Ozhogin et al., 2012) that is valid below L = 2 and includes latitudinal dependence. To our knowledge, there was no valid empirical model below L = 1.5, nor one that includes variations due to solar activity, season, local time, and hemispheric differences, in addition to L , λ below L = 2. Recently, Hartley et al. (2023) combined Van Allen Probes data for latitudes below 20° with Arase data up to 40° for 1 < L ≤ 3 and derived a new electron density model with both a latitudinal and MLT dependence. Comparison with the L dependence of the Ozhogin model shows good agreement above L = 1.5. Below L = 1.5, a fitting form similar to the Ozhogin model is adopted with new parameters defined as α = 1.03 , β = 0.44 . An MLT dependence of the plasma density was identified, which is consistent with the diurnal variation of ionosphere. This variation is strongest at low L, but persists out to L = 3. All empirical electron density models discussed in this article are listed and succinctly synthetized in Table 1.

There exist a variety of empirical plasmapause models currently used in radiation belt simulations. They generally provide the radial distance of the plasmapause in the equatorial plane as a simple function of the geomagnetic activity level. Tu et al. (2009) used, for instance, the CRRES data driven model of O’Brien and Moldwin (2003), whereas Tu et al. (2013) implemented the Carpenter and Anderson (1992) plasmapause model (noted CA92). The CA92 plasmapause model is the most commonly used model to our knowledge. It has been largely used for radiation belt studies over the last 10 years (Subbotin and Shprits, 2009; Kim et al., 2011; Shprits et al., 2013; Ripoll et al., 2016; Ripoll et al., 2019; Wang and Shprits, 2019; Cervantes et al., 2020a; Cervantes et al., 2020b; Malaspina et al., 2020; Saikin et al., 2021).

Recently, Ripoll et al. (2022a) derived both a plasmapause and a 100 #/cc density level models based on the entire Van Allen Probes mission (2012-2019) from both the Electric and Magnetic Field Instrument Suite and Integrated Science (EMFISIS) suite’s (Kletzing et al., 2013) and the Electric Field and Waves (EFW) (Wygant et al., 2013) data. The cold plasma densities are either determined by the upper hybrid resonance (UHR) method from EMFISIS measurements (Kurth et al., 2015) or by using the spacecraft floating potential (Escoubet et al., 2007; Torkar et al., 2016; Torkar et al., 2019) measured by the Electric Field and Waves (EFW) instrument (Wygant et al., 2013).

About the accuracy of these density measurements, we note the densities derived along with the corresponding spacecraft potentials are fit to a function with a non-linear least squares fit. The resulting fits typically have a Pearson (R ²) coefficient in the range of ∼0.75–0.95 and an average percent error between the selected fit derived densities and the densities used to perform the fit of ∼15%. Experiments with individual orbits show that fits of the functional form can capture the density voltage relation over a range of densities from ∼few cm^-3 up to 3,000 cm^-3 with the lower densities still agreeing with the EMFISIS UHR densities to within 10%. However, using the same fit for a longer period (larger than an orbit) the EFW and EMFISIS densities may diverge by over factor of two at densities <10 cm^-3. The reason for this is the variability of the plasma environment outside the plasmasphere. For periods during which the upper hybrid resonance line is clearly resolved in the High Frequency Receiver (HFR) spectral data, the EMFISIS density product is generally more accurate than the EFW. However, during times in which there are high levels of wave activity that make identification of the upper hybrid line difficult or impossible, resulting in increased uncertainty in the EMFISIS densities, the EFW density fits still return densities by applying the relevant fit equation to the spacecraft potential. Regarding the semi-automated process for determining the EMFISIS density from the UHR (Kurth et al., 2015), there is a 8.7% mean percentage difference between the manual process and the semiautomatic process, which is less than the ∼10% resolution available for an individual measurement. This difference is visible in Figure A2 of Goldstein et al. (2014a), where the average difference is often low (∼7%), is less than 20% in general, but can be up to 100% for a very small number of data points. Another main source of error is the spectral resolution, due to the upper hybrid resonance that can only be defined at specific values dictated by the binned frequency spectrum. This translates to a density resolution of Δne/ne of about 10%. The uncertainty increases when the spectra become difficult to interpret, as discussed in Goldstein et al. (2014a). In most instances, the spectral resolution uncertainty is estimated to be between 10% and 20% (Hartley et al., 2016). The EMFISIS and EFW densities from ∼10 cm^-3 to 3,000 cm^-3 are statistically compared in Jahn et al. (2020), who found that the EFW values predominantly fall in a range of 50%–200% of their corresponding EMFISIS measured value (e.g., 0.5 to 2 times the actual value), while most of the EFW to EMFISIS points used for comparison are ∼100% (e.g., n_EFW ∼ n_EMFISIS). Further comparisons of the 100 #/cc level are carried out in Ripoll et al. (2022a) in which Figures 1K, L confirm the good agreement between both methods, with the bulk of normalized differences below ±20%.

A comparison of the CA92 plasmapause model with Van Allen Probes measurements is performed in Ripoll et al. (2022a). These authors first recover the CA92 model using EMFISIS data and a gradient method to localize the plasmapause, showing the practical reliability of the CA92 model. However, direct comparisons of the 100 cm^-3 level deduced from Van Allen Probes EFW measurements and the CA92 model show the dense plasmasphere expands farther out than predicted by the CA92 model. Departure of the CA92 model from the 100 cm^-3 EFW data increases as the maximum value of the Kp index over the last 24 h (Kp) increases and L-shell decreases, and storm-induced erosions are less deep than predicted by the CA92 model (Ripoll et al., 2022a).

The model of O’Brien and Moldwin (2003) based on CRRES data is the first to show the MLT dependence as well as the relevance of parametrizing the plasmapause model with various indices, such as Kp, AE, and Dst (see also Moldwin et al. (2002)). Carpenter et al. (2000) states the experimental error in the CRRES density is associated with measuring the UHR or plasma frequencies on the SFR records. They estimated to be +/− 6% in spectral resolution (Δf/f), which corresponds to +/− 12% in density. Kwon et al. (2015) derived the median/mean plasmapause locations from the electron density inferred from the Time History of Events and Macroscale Interactions during Substorms (THEMIS) spacecraft potential under steady quiet conditions (Kp ≤ 1). The comparison of their plasmapause model with the estimated Lpp from models such as GCPM (Gallagher et al., 2000), Moldwin et al. (2002), and O’Brien and Moldwin (2003) with Kp = 1 shows the plasmapause is farther extended ∼1–2 L from the Earth (i.e., GCPM and CRESS based-models models underestimate the extend of the plasmapause). Ripoll et al. (2022a) show the underestimation of the plasmapause position is caused by the gradient method that fails identifying gradients particularly during quiet times and on the dusk.

Other plasmapause models include Bandić et al. (2016) based on CRRES data, Cho et al. (2015), Liu and Liu (2014) and Liu et al. (2015) based on THEMIS data, and Larsen et al. (2007) based on IMAGE data. Verbanac et al. (2015) plasmapause model is based on CLUSTER data and analytical relationships obtained from geomagnetic and solar wind observations. Bandic et al. (2017) derived a plasmapause model from a large dataset including multiple sources. A comparison of these models is provided in Pierrard et al. (2021c) showing a great variability of mean plasmapause empirical models (see also Guo et al., 2021). All empirical plasmapause models discussed in this article are listed and succinctly synthetized in Table 2 (see also Table 1 in He et al. (2017) listing the model dependences).

The large variability of the measurements underlying the mean empirical plasmapause models (more generally mean plasma density empirical models) is one major limitation of this type of models that calls for the use of either physic-based models or machine learning technics.

The 3D global ionosphere/plasmasphere fluid model SAMI3 (Huba and Krall, 2013; Krall and Huba, 2013) of the Naval Research Laboratory (NRL) solves the continuity and momentum fluid equations for seven ion species (H⁺, He⁺, N⁺, O⁺, N⁺ ₂, NO⁺ and O⁺ ₂) and includes the thermospheric wind-driven dynamo electric field. It is based on SAMI2 (Huba et al., 2000). SAMI3 uses the partial donor cell method (Hain, 1987; Huba, 2003) and a newly implemented 4-order flux-corrected transport scheme for ExB transport perpendicular to the magnetic field (Huba and Liu, 2020). The temperature equation is solved for three atomic ion species and electrons. The model has a co-rotation potential, a neutral wind dynamo potential (with winds from HWM93 (Hedin, 1987)), and a time-dependent Volland-Stern-Maynard-Chen potential. In Huba and Krall (2013), SAMI3 density results are compared at the equator for 4 MLT sectors with the quiet time empirical electron density of Berube et al. (2005) defined as n_eq = 10^0.51L+4.56 for L in 2–5 from IMAGE RPI data between May 2000 and May 2001. They find the SAMI3 electron density is lower by a factor 2 attributed to a lower F10.7 index used in the simulation.

The SAMI3 model has been recently modified to support the NASA ICON mission and provide ionosphere and thermosphere properties during this mission (Huba et al., 2017). SAMI3 recently integrated an improved model of counterstreaming H⁺ outflows from the two hemispheres during storm using a two fluid species for H+ (Krall and Huba, 2019) in order to avoid non-physical high-altitude ‘top-down refilling’ density peaks (Krall and Huba, 2021). SAMI3 is currently used to try to reproduce the formation of density ducts in the plasmasphere (e.g., Jacobson and Erickson, 1993; Loi et al., 2015) caused by the thermosphere composition and winds on plasmaspheric refilling outflows (Krall et al., 2018) as observed from the Murchison Widefield Array (MWA) interferometric radio telescope in Australia (Helmboldt and Hurley-Walker, 2020). SAMI3 recently coupled to the atmosphere/thermosphere code WACCM-X (Whole Atmosphere Community Climate Model with thermosphere and ionosphere extension) provided the first high-resolution global simulation using realistic thermospheric conditions of the formation and penetration of plasma bubbles into the topside F layer (Huba and Liu, 2020). These structures will further propagate to higher altitudes and introduce longitudinal and seasonal dependence structures into the plasmasphere. Further coupling of SAMI3 and applications are discussed in Huba (2023).

The Ionosphere-Plasmasphere-Electrodynamics (IPE) model is derived in Maruyama et al. (2016) to investigate the connection between terrestrial and space weather (e.g., Fuller-Rowell et al., 2008). IPE provides 3D thermal plasma densities for nine ion species, electron and ion temperatures, and parallel and perpendicular velocities of the ionosphere and plasmasphere. The parallel plasma transport is based on the Field Line Interhemispheric Plasma (FLIP) Model (Richards et al., 2010). There is a detailed model of the Earth’s magnetic field using Apex coordinates (Richmond, 1995) and the International Geomagnetic Reference Field IGRF (as in SAMI3). The transport is computed with the same solver all the way from the equator to the pole on a global static grid with a semi-Lagrangian scheme that allows for the global plasma transport perpendicular to magnetic field lines. There is a self-consistent photoelectron calculation enabling more accurate studies of the longitudinal/UT dependence of the ionospheric mass loading process. IPE is generally defined from 90 km to approximately 10,000 km. The spatial resolution of the radial direction in the plasmasphere varies from 0.05 R_E (L = 1.5) to 0.46 R_E (L = 5). IPE has used to reproduce the Weddle Sea Anomaly (Sun et al., 2015) and for studying extreme plasmaspheric erosion as low as L∼1.7 (Obana et al., 2019). Current applications of IPE include plasmaspheric drainage plumes, ionospheric storm enhanced density (SED) plumes, plasmaspheric refilling, and plasmaspheric composition. The Whole Atmosphere Model (WAM)(e.g., Akmaev and Juang, 2008) has been coupled with IPE (WAM-IPE) and provides today space weather forecast 24/7 at NOAA SWPC (https://www.swpc.noaa.gov/products/wam-ipe). WAM-IPE has recently be used to simulate ESF irregularities (Hysell et al., 2022).

The IRAP Plasmasphere Ionosphere Model (IPIM) uses a 16-moment approach for strong temperature anisotropy at high altitude and for accurately modeling the transition between collision dominated at low altitude and collisionless media at high altitude (Marchaudon and Blelly, 2015). IPIM solves the interhemispheric hydrodynamics convection and corotation of six ions and thermal electrons along flux tubes at different distances from Earth. IPIM has a kinetic model for suprathermal electrons and solves for the chemical reactions in the ionosphere. IPIM has been used to simulate the depletion of the ionospheric F₂ layer by a high-speed stream for short-term behavior on the scale of a few hours Simulations were found to be consistent with EISCAT radar and the ionosonde measurements (Marchaudon et al., 2018). For the long-term evolution of the plasmasphere-ionosphere system and during quiet conditions, IPIM simulations indicate that the plasmasphere in not stable in MLT and that no real dynamic equilibrium can be reached (Marchaudon and Blelly, 2020).

Different plasmasphere models combining semi-empirical relations and physics-based backgrounds have been developed to reproduce the inner magnetosphere, the plasmapause, and even the plasma trough above the plasmasphere limit (see Pierrard et al., 2009 for a review of the plasmasphere models before 2009).

Pierrard and Stegen (2008) have developed the Belgian SWIFF Plasmasphere Model (BSPM), a 3D dynamic kinetic model of the plasmasphere. The BSPM model is based on physical mechanisms, including the interchange instability for the formation of the plasmapause (Pierrard and Lemaire, 2004), and provides the density and the temperature of the electrons, protons and other ions, both inside and outside the plasmasphere in the plasma trough. It has been coupled to the ionosphere (Pierrard and Voiculescu, 2011) using the IRI model as a boundary condition and is continuously improved by including other physical processes like plasmapause thickness and plasmaspheric wind (Pierrard et al., 2021b). The input of the model is the date that determines the geomagnetic indices Kp and Dst. The plasmapause position does strongly correlate with the Bartels geomagnetic index, Kp index, which is retained as the main parameter used in the model to determine the plasmapause position. These indices may be predicted values when forecasting is required, or observed values when past events are simulated. They determine also the convection electric field. As BSPM uses the IRI model, it also depends on IRI parameters listed in Section 2. The BSPM model includes plasmapause erosion during geomagnetic storms as well as refilling, and is able to reproduce the plumes generated during storms and other structures like shoulders. It uses the kinetic approach that allows for the inclusion of non-Maxwellian distributions (Pierrard and Lemaire, 2001). The last version of the BSPM model is shown in Figures 1D–F for quiet, substorm, and storm activity. On 16 March 2015 1H (UT) with a quiet period with Kp∼2 (and almost constant during several hours), the plasmasphere is quite extended and almost circular (to compare with Figure 1A). A few hours later after a substorm injection on 16 March 2015 19 h (UT) with Kp∼4, there is formation of a plume in the dusk sector rotating with the Earth (to compare with Figure 1B). On 17 March 2015 21 h (UT) during an intense storm with Kp = 8^-, the model shows a strong erosion of the plasmasphere and formation of a long plume rotating with Earth (to relate with the statistics of Figure 1C).

The kinetic approach based on particle-in-cell simulations has also been combined with the fluid approach in Wang et al. (2015) to develop a dynamic fluid-kinetic model for plasma transport within the plasmasphere. A semi-kinetic model of plasmasphere refilling following geomagnetic storms has also been recently developed by Chatterjee and Schunk (2020b) and compared with hydrodynamic models to explore their differences. In hydrodynamic plasmasphere models, the non-linear inertial terms in the plasma transport equations are retained (Chatterjee, 2018; Chatterjee and Schunk, 2019; Chatterjee and Schunk, 2020a; Chatterjee and Schunk, 2020b). Limitations of such models are generally related to the difficulty to reproduce the mechanisms implicated in the formation of the plasmapause and the refilling process that is a key physics-based problem to solve to obtain a fully coupled plasmasphere-ionosphere model.

A two-dimensional physics-based plasmasphere model called Cold Plasma (CPL) (Jordanova et al., 2006; Jordanova et al., 2014) is used in a ring current-atmosphere interactions model of the source and loss processes of refilling and erosion driven by empirical inputs to simulate equatorial plasmaspheric electron densities. The performance of CPL has been evaluated against in situ measurements by the Van Allen Probes (Radiation Belt Storm Probes) for two events (De Pascuale et al., 2018). This study finds that severe erosion is best captured by an effective Kp-index for scaling the inner-magnetospheric potential governing E x B flows while refilling subsequent to moderate activity requires a solar wind parameterization of the quiet time background after the onset of a geomagnetic storm. Empirical models driving plasmasphere dynamics can be improved by capturing localized enhancements in electric field measurements and asymmetric profiles in electron density observations. More specific simulations were dedicated to comparisons with Van Allen Probes plasmapause observations (Goldstein et al., 2014a; Goldstein et al., 2016).

Physics-based models also provide the plasmapause location (e.g., Pierrard et al. (2021b)), with some models integrating Van Allen Probes measurements and plasma trough densities (e.g., Botek et al., 2021). Goldstein et al. (2003), Goldstein et al. (2005) developed a plasmapause test particle (PTP) dynamic model that represents the plasmaspheric boundary as an ensemble of E × B-drifting particles. The PTP model uses an electric field which is driven by the solar wind E field and Kp. The evolution of the plasmapause is modeled by the changing shape of the curve defined by the aggregate of the test particles evolving in a time-varying convection E-field. PTP simulation for the moderately disturbed interval 18–20 January 2000 shows a narrow drainage plume followed by significant plasmaspheric erosion, forming a second plume that coexists with the residue of the first plume (Goldstein et al., 2014b). Observations from three of the Los Alamos National Laboratory geostationary satellites are globally consistent with this PTP simulation in terms of the durations of plume sector transits while the MLT widths and timings of the simulated plumes do not precisely agree (Goldstein et al., 2014b). Goldstein et al. (2019) further generated a plasmapause statistical model from the simulations of 60 storms with Dst,_PEAK ≤ −60nT based on Van Allen probes data yielding over 7 million model plasmapause locations. The epoch-binned PTP simulation results are combined in order to create an analytical plasmapause model for moderate storms (−120nT ≤ Dst,_PEAK ≤ −60nT) and strong storms (Dst,_PEAK ≤ −120nT) that explicitly includes plumes. This model depends on the duskside plasmapause radius and two fitted coefficients, all three depend on epoch time (from −24 h to 36 h).

A new promising approach is to couple a global geospace model of the magnetosphere with a physics-based density model. Figure 2 provides an example of the global geospace model, GAMERA (Zhang et al., 2017; Sorathia et al., 2020; Sorathia et al., 2021) coupled to RCM (Toffoletto et al., 2003). With a two-way coupling, these models are subparts of the Multiscale Atmosphere-Geospace Environment (MAGE) (e.g., Chen et al., 2021; Pham et al., 2021; Lin et al., 2022). The details of GAMERA’s core MHD numerics and its verification are presented in Zhang et al. (2019). GAMERA uses high-order spatial reconstruction for the preservation of sharp structures. For typical MHD problems, Zhang et al. (2019) showed lower-order reconstruction (e.g., Second-order) requires four to eight times finer grid resolution (corresponding to a 250-4,000 factor increase of the cost resolution in 3D) as the higher-order (seventh- or eighth-order) reconstruction to reach the same accuracy. In addition to coupling the global MHD model to the inner magnetosphere model via ring current pressure ingestion (e.g., Pembroke et al., 2012), here the RCM is additionally evolving a cold fluid to model the evolution of the plasmaspheric density. In this coupling, the plasmasphere density is initialized using an empirical model (Gallagher et al., 2000) and refilling rate (Denton et al., 2012), and evolved using the same dynamically-calculated electrostatic potential as in the MHD simulation (e.g., Merkin and Lyon, 2010). Note that RCM can further be coupled with SAMI3 as done by Huba et al. (2017b) to study the ionosphere-plasmasphere system response to the 17 March 2015 geomagnetic storm. The coupling occurs through the electrostatic potential equation (Huba et al., 2005; Huba and Sazykin, 2014) in which the conductance is defined by the sum of the conductance associated with solar activity computed by SAMI3 and the auroral enhanced conductance provided by RCM.

Figure 2A depicts localized dipolarizations on the nightside and the formation of a dayside plasmaspheric plume during the 17 March 2013 geomagnetic storm (Sorathia et al., 2018). There is a complex interacting mesoscale process with nightside flows, boundary Kelvin-Helmholtz on the dayside and flanks, and rolling dense plasmaspheric plume and structures. The plume is shown at 12 UT, i.e., 6 h after the CME impacted the Earth, with a typical expansion in the dusk-day sector that reaches L∼6 and has started to roll around Earth.

The Kelvin-Helmoltz instability we see forming on the magnetopause and rolling side way of the magnetosphere (Merkin et al., 2013) may contribute to transfer shear and turbulence to the plume as it expands and removes pockets of dense plume plasma. In this way, the plume may potentially inherit a complex shape that is here captured by the global MHD model. The dense plasmasphere has a circular aspect for levels above 1,000 #/cc and there are structured plasma pockets of low density from 1-10 #/cc on the nightside beyond the main plasmapause gradient at L∼3. On the night side, the magnetic field (and similarly the electric field) has a fine scale structure (with finger-like regions of higher field value) that reach the L∼6 region and imprint a fluctuating profile to the dense pockets down to the plasmapause layer. As simulation resolution increases, some aspects of these structures become finer and more torturous. However, understanding the full cascade of energies down to the smallest scales requires global kinetic modeling.

These kinds of mesoscale structures play a critical role in shaping both the global-scale and micro-scale processes of the magnetosphere. Localized injections are believed to be an important part of the transport of magnetic flux and energetic particles into the inner magnetosphere (e.g., Gkioulidou et al., 2014; Merkin et al., 2019), thus building the large-scale ring current and affecting global dipolarization of the inner magnetosphere, as well as resulting in density enhancements at the dayside magnetopause that will alter local reconnection rates (e.g., Zhang et al., 2017) with potentially global consequences. Additionally, these mesoscale processes shape the different wave populations of the inner magnetosphere: anisotropic ion injections provide free energy for the ElectroMagnetic Ion Cyclotron (EMIC) wave population and the evolving plasmapause boundary correlates with the relative distribution of hiss and chorus waves, with important consequence on flux enhancements of energetic trapped particles in the radiation belts. Physics-based models discussed in this article are listed and succinctly synthetized in Table 3.