The Contribution of Plasma Sheet Bubbles to Stormtime Ring Current Buildup and Evolution of Its Energy Composition

The formation of the stormtime ring current is a result of the inward transport and energization of plasma sheet ions. Previous studies have demonstrated that a significant fraction of the total inward plasma sheet transport takes place in the form of bursty bulk flows, known theoretically as flux tube entropy‐depleted “bubbles.” However, it remains an open question to what extent bubbles contribute to the buildup of the stormtime ring current. Using the Multiscale Atmosphere Geospace Environment Model, we present a case study of the 17 March 2013 storm, including a quantitative analysis of the contribution of plasma transported by bubbles to the ring current. We show that bubbles are responsible for at least 50% of the plasma energy enhancement within 6 RE during this strong geomagnetic storm. The bubbles that penetrate within 6 RE transport energy primarily in the form of enthalpy flux, followed by Poynting flux and relatively little as bulk kinetic flux. Return flows can transport outwards a significant fraction of the plasma energy being transported by inward flows, and therefore must be considered when quantifying the net contribution of bubbles to the energy buildup. Data‐model comparison with proton intensities observed by the Van Allen Probes show that the model accurately reproduces both the bulk and spectral properties of the stormtime ring current. The evolution of the ring current energy spectra throughout the modeled storm is driven by both inward transport of an evolving plasma sheet population and by charge exchange with Earth's geocorona.

It is well-accepted that the ring current is created via inward transport of the plasma sheet particles (Kistler, 2020;Liemohn & Khazanov, 2005;Wing et al., 2014).However, to what extent nightside transport at different spatial and temporal scales contributes to ring current buildup remains an outstanding question of magnetospheric physics.Transport by quasi-steady global-scale convection enforced by a dawn-to-dusk electric field that becomes enhanced during storms was proposed originally (Axford, 1969).Erickson and Wolf (1980), however, noted that steady, global convection through the tail would result in unrealistically high pressures in the inner magnetosphere, and hypothesized that time-dependent ejection of excess plasma in the form of substorms could provide a solution to this "pressure crisis."The model of time-variable but global-scale convection has since been the most widely accepted explanation for stormtime ring current formation (e.g., Daglis, 2006).
However, the existence of a stormtime enhancement in the global-scale cross-tail electric field was challenged by statistical measurements of the electric field in the plasma sheet.It was found from CRRES measurements between 7.5 and 8.5 R E (Rowland & Wygant, 1998) and Geotail measurements between 5 and 15 R E (Hori et al., 2005) that the large-scale field remained relatively weak during heightened activity.There were instead many instances of large-amplitude (several to tens of mV/m) electric fields highly fluctuating on minutes-long timescales.Hori et al. (2005) noted that the azimuthal width of the convection channels corresponding to the fluctuating field should be considered, as these enhancements may be caused by narrow-width bursty flows (Angelopoulos et al., 1997).Such "bursty bulk flows" (BBFs; Angelopoulos et al., 1992;Baumjohann, 1990) are indeed estimated to be responsible for a significant amount (>60%) of the total inward transport of magnetic flux and plasma mass and energy transport in the magnetotail (15-22 R E ; Angelopoulos et al., 1994).BBFs are observed as short-lived (∼10 s), high-speed (>400 km/s) Earthward flows that typically contain a significant enhancement in magnetic flux, referred to as a "dipolarizing flux bundle" (DFB; Liu et al., 2013), and an enhancement of hot plasma, observed as an injection (Gabrielse et al., 2012(Gabrielse et al., , 2014;;Liu et al., 2016;Runov et al., 2011;Sergeev et al., 2005).With a cross-tail size of ≲3 R E (Angelopoulos et al., 1997;Nakamura et al., 2004), they are often referred to as "mesoscale" flows since their scale size is much smaller than the global magnetotail size yet much larger than the kinetic scale.C. X. Chen and Wolf (1993) suggested that BBFs provide an alternative solution to the pressure crisis as described by Erickson and Wolf (1980) if they account for at least 40%-55% of the total inward plasma transport in the magnetotail, introducing the possibility that BBFs play a significant role in ring current buildup.
However, BBF transport within the plasma sheet does not directly correspond to contribution to the ring current, since the majority of these flows do not penetrate to within geosynchronous orbit (GEO; Ohtani et al., 2006;Takada et al., 2006).The best predictor of BBF penetration depth is not the flow velocity or magnetic field strength within the structure, but rather the flux tube entropy (FTE; Dubyagin et al., 2011), defined as S = pV 5/3 , where p is the average plasma pressure along the flux tube and V = ∫1/Bds is the flux tube volume (Birn et al., 2009;Pontius & Wolf, 1990).Plasma stability results in a distribution of magnetospheric flux tubes such that the FTE increases with radial distance from the Earth.High-magnetic flux, low-entropy "bubbles" introduced in the magnetotail are interchange unstable and are transported inwards (where they may be observed as BBFs) until they reach a location with comparable background FTE (Pontius & Wolf, 1990;Wolf et al., 2009).
The question of the relative contribution of BBFs/bubbles to stormtime ring current buildup has been notoriously difficult to answer conclusively because of the lack of multipoint measurements with sufficient coverage and resolution.Gkioulidou et al. (2014), using sub-GEO injections observed by the Van Allen Probes (RBSP), estimated that for the 17 March 2013 storm BBFs transported around 30% of the total plasma energy to the ring current.However, this result required assumptions about the azimuthal distribution of GEO-penetrating injections and the occurrence rate of injections throughout the storm extrapolated from the number observed by RBSP alone.While observations from multiple spacecraft somewhat reduce this uncertainty, studies that have utilized large conjunctions of spacecraft still have difficulty determining whether spacecraft at differing radial and azimuthal locations are observing the same BBF (e.g., Gabrielse et al., 2014;D. L. Turner et al., 2017).They also may not be capturing all BBFs that might be traversing the plasma sheet across different MLTs at any given time.Hence, more complete in situ coverage is required to conclusively address the problem.Energetic neutral atom (ENA) imaging has recently been used to detect temperature enhancements associated with BBFs (Keesee et al., 2021) and has the prospect of helping to quantify plasma energy transported by BBFs with the advantage of observing the entire tail at once.However, due to low sensitivity of the ENA camera onboard the TWINS spacecraft (due to a very small geometric factor and duty cycle; McComas et al., 2009), only limited energies can be used (2-32 keV), assuming a Maxwellian distribution, to retrieve those temperatures.

10.1029/2023JA031693
3 of 23 Because of the current observational limitations, numerical modeling has been used extensively to study inner magnetosphere transport and ring current formation, and has increasingly found that bubbles play an important role in both cases.Inner magnetosphere models that have an outer boundary around 6-9 R E , and apply a global-scale electric field there, have found that such convection can build up a ring current and reproduce observed trends in Dst (e.g., M. W. Chen et al., 1993;Ebihara & Ejiri, 2000;Fok et al., 2014;Jordanova et al., 1996Jordanova et al., , 2014;;Kozyra et al., 1998).Using the Rice Convection Model-Equilibrium (RCM-E) with a boundary at GEO, Yang et al. (2016) showed that driving their simulation with an electric field both with and without the presence of bubbles produced similar ring currents.However, both Yang et al. (2016) and Lemon et al. (2004) found that when the nightside boundary was instead placed at 15 R E , driving their simulation with a strong global-scale electric field produced a pressure enhancement outside of 6 R E , which prevented further convection inwards and subsequent buildup of the ring current.In both studies, when entropy-depleted bubbles were introduced along the 15 R E boundary, a substantial ring current formed, indicating that the presence of bubbles in the plasma sheet plays a role in ring current buildup.Cramer et al. (2017), utilizing the Open Geospace General Circulation Model (OpenGGCM) magnetohydrodynamic (MHD) code coupled with RCM, found that bubbles were responsible for 65%-85% of the total inward plasma transport to below GEO, depending on storm intensity.Sorathia et al. (2021), using the stand-alone version of the Grid Agnostic MHD for Extended Research Application (GAMERA) global MHD model (Sorathia et al., 2020;Zhang et al., 2019) and test particle tracing (Sorathia et al., 2018;Ukhorskiy et al., 2018), concluded that Earthward-propagating bubbles in the plasma sheet accounted for roughly half of the total plasma energy transport while making up only 15% of the plasma population.The studies by Cramer et al. (2017) and Sorathia et al. (2021) provided quantitative evidence of the importance of bubbles in plasma mass and energy transport in the magnetosphere.However, the contribution of bubbles to the ring current energy content was not directly examined.
The study presented in this paper is the first to quantify the contribution of bubbles to the stormtime ring current using a coupled global and inner magnetosphere model.This work builds off of that of Yang et al. (2015), which investigated this question using RCM-E and found that up to 60% of the total inner magnetosphere pressure came from plasma introduced by bubbles along the 15 R E nightside boundary.We expand on this topic through use of GAMERA two-way coupled with RCM, which has the advantage of capturing the dynamics of the global system and generating bubbles more self-consistently.GAMERA, as well as its predecessor, the LFM global MHD model (Lyon et al., 2004), has a history of accurately reproducing the statistical properties of BBFs (Sorathia et al., 2021;Wiltberger et al., 2015), and with the addition of RCM coupling is overall well-suited to address this fundamental question of magnetospheric physics.
The goal of understanding ring current buildup must also consider the evolution of its energy composition over time.Observations have found that during quiet times, protons above ∼100 keV are the dominant contributors to the ring current pressure, while during storms the main contributors to the enhanced pressure are ions with energies between several to tens of keV (e.g., Gkioulidou et al., 2016;Keika et al., 2018;Korth et al., 2000;Krimigis et al., 1985;Zhao et al., 2015).There have been a number of proposed mechanisms for this energy evolution, including the penetration of different plasma sheet populations throughout the storm (Keika et al., 2005;Lui, 1993), adiabatic effects due to localized changes in the magnetic field strength (Lyons, 1977;Lyons & Williams, 1976), radial diffusion (Jordanova & Miyoshi, 2005;Lyons & Schulz, 1989;Sheldon & Hamilton, 1993) and charge exchange (e.g., Fok et al., 1993;Jordanova et al., 1998;Kistler et al., 1989).However, similarly to investigations of bubble contribution to the ring current, the nature of in situ observations makes it challenging to make conclusive statements about which processes are ultimately responsible for the observed evolution.Reproducing this behavior in first-principle geospace models is necessary in order to self-consistently model the ring current's coupling with the ionosphere and exosphere, and as a result more accurately capture the recovery phase.In turn, such model capabilities help to inform to what extent certain processes influence the evolution of the energy spectra during storms.
In this study, we present a comprehensive picture of the contribution of bubbles to the buildup of the stormtime ring current, investigate the evolution of the ring current energy spectra throughout the storm, and examine the contributors to this evolution.This is done via a case study of the 17 March 2013 geomagnetic storm simulated with the Multiscale Atmosphere-Geospace Environment (Multiscale Atmosphere Geospace Environment (MAGE)) model, which includes two-way coupling between GAMERA and RCM.Section 2 provides an overview of the storm, the numerical model, and the data used for data-model comparisons.Section 3 presents our results, including a quantitative analysis of BBF/bubble contribution to the inner magnetosphere energy content, 10.1029/2023JA031693 4 of 23 data-model comparison of proton intensities, and an analysis of the contributors to the evolving ring current energy spectra.In Section 4, we discuss the findings and caveats of the study, and in Section 5 we present our conclusions.

Event Description
For this study we simulate the 17 March 2013 storm, one of the major storms of the Van Allen Probes era.This storm was caused by a coronal mass ejection, and had the sudden storm commencement (SSC) at about 06:00 UT on 17 March 2013 and a minimum Dst of approximately −140 nT.In this study we begin the simulation at 00:00 UT that day and, since we are focusing on ring current buildup, we only examine the main phase out to 19:00 UT, which saw a minimum Dst of about −100 nT (see the Sym-H panel in Figure 1).Additional solar wind parameters are provided in Figure S1 in Supporting Information S1.This storm featured many injections observed by the Van Allen Probes within GEO (Gkioulidou et al., 2014), providing an ideal opportunity to study the role of injections in the buildup of Earth's stormtime ring current.This storm has been studied extensively, including modeling (e.g., M. W. Chen et al., 2019;Lin et al., 2021;Raeder et al., 2016;Yu et al., 2014Yu et al., , 2015;;Wiltberger et al., 2017) and observational studies (e.g., Dang et al., 2019;Gkioulidou et al., 2014;Lyons et al., 2016;Menz et al., 2017;Tang et al., 2016;Yue et al., 2016).

Data
The solar wind data source is primarily Modified High Resolution OMNI (HRO) data.For this time period, the HRO data is based on ACE observations, which contains a 45 min data gap around the SSC.We fill any 5-min or longer data gaps with the propagated and KP-despiked Wind data set available on OMNIWeb.This method produces a significantly sharper SSC than interpolating across the 45-min data gap.
Proton data from the RBSPICE instrument aboard Van Allen Probe B (RBSP-B) is used for comparison with our model.During this storm, the RBSPICE instrument aboard Probe A was not fully operational and so its data is not used here.Because spacecraft charging can lead to bad data below 10 keV (Menz et al., 2017), we choose this as our lower energy limit.We choose an upper limit of 200 keV as this is about the maximum energy in this simulation that our inner magnetosphere model (see Section 2.3) covers throughout the entirety of the RBSP-B orbit.This energy range is fully covered using the combination of the Time-of-Flight by Pulse Height (TOFxPH) data product for proton energies between 10 and 50 keV and the Time-of-flight by Energy (TOFxE) data product for 50 keV protons and upwards (Mitchell et al., 2013).The Level-3 omnidirectional flux from both of these data products was used.This data was retrieved from cdaweb (https://cdaweb.gsfc.nasa.gov)via the cdasws python package (https://pypi.org/project/cdasws/).

Model Description
The MAGE model, developed at the NASA DRIVE Science Center for Geospace Storms (CGS), is a coupled geospace model designed to enable the investigation of multiscale phenomena, from the global magnetosphere to the ionosphere and thermosphere, during geomagnetic storms.A visual overview of the simulation at 09:28 UT on 17 March 2013 is given in Figure 1 which highlights the components of MAGE used in this study.The primary physics module is the GAMERA 3D global MHD model (Sorathia et al., 2020;Zhang et al., 2019).GAMERA is a reinvention of the well-known LFM model (Lyon et al., 2004), preserving the original high-order, non-diffusive numerics.Minimizing numerical diffusion between the low-entropy bubble with the surrounding background plasma is critical, as this enables better representation of the bubble penetration depth.GAMERA is coupled with the REMIX ionospheric potential solver, which is a rewrite of the Magnetosphere-Ionosphere Coupler/Solver (MIX) model (Merkin & Lyon, 2010).GAMERA and REMIX are also two-way coupled with the RCM, a flux tube-averaged inner magnetosphere model that resolves energy-dependent gradient-curvature drifts (Toffoletto et al., 2003;Wolf et al., 2016).The GAMERA inner boundary for this simulation is at a radius of 2 R E , where field-aligned currents and the ionospheric electric fields are mapped to/from the ionosphere using a dipole field.Self-consistently calculated mono-energetic and diffuse electron precipitation is included in the calculation of ionospheric conductance (Lin et al., 2021).
There is a "spin-up" phase, lasting several hours, where only GAMERA and REMIX are advanced using the solar wind conditions at time T = 0 (00:00 UT on 17 March 2013, in this case), with the interplanetary magnetic field (IMF) set to zero, as a constant driver.At T = 0, the full solar wind time-series is used, and the RCM starter ring current is initialized using an analytic pressure distribution described in Liemohn (2003), specifically the L-dependent profile in the nominal case where the pressure peak is centered at L = 4 and ΔL = 0.625.The pressure distribution is scaled such that its total contribution to the Dst, as calculated via the DPS-Dst relation, makes up the difference between the observed Dst and that calculated by GAMERA at T = 0.An ad-hoc temperature of 30 keV is used in combination with the pressure profile to calculate the density profile, assuming a Maxwellian distribution.Beyond 10 R E , the plasma sheet is initialized using the empirical model described in Tsyganenko and Mukai (2003).For the solar wind conditions preceding this storm, this resulted in an average plasma sheet temperature of around 10 keV.
The coupling between GAMERA and RCM enables the model to resolve energy-dependent drifts, precipitation, and charge exchange within the inner magnetosphere.The high-level coupling procedure is similar to that implemented in LFM-RCM (Pembroke et al., 2012), with several important improvements.In both LFM-RCM and GAMERA-RCM, an RCM "active domain" is defined, only within which RCM applies its advection including gradient-curvature drifts and preserves its energy invariant distribution between coupling time steps, and is where the combined density and pressure of all RCM species are passed back to the MHD model for ingestion.In LFM-RCM, the active RCM domain was not only confined to the region of closed field lines, but also to an ellipse that encompassed as much of the closed domain as possible, and was further confined to regions where the plasma β ≤ 1.This ensured the active domain was limited to regions that were magnetically dominated, consistent with the assumptions within RCM.The plasma β limiter was introduced to exclude fast flows, which may break the validity of the flux-tube averaged approximations.In GAMERA-RCM, the ellipse is more conservatively confined to be no closer than 4 radial grid cells to the open/closed boundary, and is further confined to regions 10.1029/2023JA031693 6 of 23 where the Alfven speed is sufficiently fast to communicate information along the field line relative to its evolution time-scale, such that the field-line averaged approximation is reasonably valid.We found this limiting of the active RCM domain to be more effective in excluding fast flows.Its implications in the context of the trapping/ leaking of plasma from bubbles are discussed in Section 4. For this study, the MHD density and pressure are mapped to RCM energy invariant channels along its outer boundary via a Kappa distribution, defined as in Equation 3.12 of Livadiotis and McComas (2013), with a kappa value of 6.This resulted in better agreement between the simulated and observed proton intensities compared to mapping using a Maxwellian distribution.The MHD pressure is partitioned between protons and electrons with a temperature ratio of 4. A dynamic plasmasphere is also modeled via a zero-energy channel in RCM.This channel is initialized using the Kp-dependent plasmasphere model given in Gallagher et al. (2000) and evolves self-consistently using the combined electrostatic potential from REMIX and corotation.Energy-dependent charge exchange losses are applied to the RCM protons using the empirical geocorona density profile from Østgaard et al. (2003) and cross-sections provided by Lindsay and Stebbings (2005).At each coupling time step (every 15 s) the plasma moments from the RCM hot plasma and plasmasphere channels are ingested by GAMERA as a single fluid.

Simulation Overview
Figure 2a shows, in addition to the observed and simulated Dst values, the Dessler-Parker-Sckopke Dst (DPS-Dst; Dessler and Parker, 1959;Sckopke, 1966) evaluated for all closed field lines with an equatorial extent within 6 R E .The DPS relation, assuming a dipole field, relates the total plasma energy within a defined volume (E) with the magnetic perturbation it causes at Earth's surface (ΔB z ), given by where B 0 is Earth's surface magnetic field and E m is the total energy of the magnetic field external to the Earth.Thus, a decrease in the plotted DPS-Dst directly corresponds to the increase in plasma energy within 6 R E .It can be seen in Figure 2a that the DPS-Dst stays relatively constant through the SSC, then decreases at a fairly steady rate starting around 8:30 UT, and then levels off again around 9:30 UT and reaches its absolute minimum around 10:30 UT.Therefore most of the ring current buildup in this simulation occurred within just 2 hours (08:30-10:30 UT).The DPS-Dst is roughly half the value of the total simulated Dst, in agreement with the estimates of Hamilton et al. (1988) and N. E. Turner et al. (2001).We note here that throughout the duration of the storm, the contribution of electrons to the total RCM plasma energy varied between 3% and 14%, and thus contributed minimally to the total ring current plasma energy in this simulation.Their contribution is still included through the coupling between GAMERA and RCM.However, a deeper analysis of the electron population is better suited for a future study.
The panels in Figure 2b show the plasma pressure within the inner magnetosphere in the equatorial plane, as calculated by RCM, for 4 times between 8:00 and 10:08 UT.The (c) and (d) rows show the radial velocity (V r ) and vertical component of the magnetic field with the dipole subtracted out (ΔB z ) from GAMERA.Evident in the V r panels are highly localized inward flow bursts (FBs), which are accompanied by dipolarizations (DFBs, visible in the ΔB z panels) and particle flux enhancements (not shown).Per the coupling algorithm described in Section 2.3, the RCM domain does not include the region containing these structures until the fast flow has subsided, causing the finger-like incursions of the outer boundary on the RCM active domain (row b).
The ring current pressure reaches a peak of 140 nPa around 9:15 UT and is located at a radial distance of just under 4 R E .Over time, the pressure peak is smeared out azimuthally as plasmas of different energies gradient-curvature drift at varying velocities.However, the total plasma energy within the inner magnetosphere remains fairly constant, as indicated by the DPS-Dst curve in Figure 2a.
Within MHD, the change in total plasma energy within a volume enclosed by a surface is given by the enthalpy flux through that surface (Birn & Hesse, 2005).Figure 3 presents a keogram-like view of several quantities evaluated along a 180° arc across the nightside with a radius of 6 R E from the origin at Earth. Figure 3a shows the radial component of the enthalpy flux (H) perpendicular to the magnetic field direction, defined as where p is the plasma pressure and v ⊥,r is the radial component of the bulk velocity perpendicular to the magnetic field vector.We denote inward/Earthward flux as positive, so that it corresponds to an increase in the plasma energy within 6 R E .The structure of the inward enthalpy flux is highly localized in MLT extent.Keeping in mind that at 6 R E , an arc spanning 1 MLT is roughly 1.6 R E in length, the azimuthal extent of these flows are similar to the scale-size of BBFs observed further in the plasma sheet (15-19 R E ; Nakamura et al., 2004).Nearly each instance of enhanced inward enthalpy flux is accompanied by magnetic dipolarizations (Figure 3b) and a depletion in FTE (Figure 3c), indicating that these flows are indeed bubbles penetrating below 6 R E .There are several notable exceptions in the 03-06 MLT range between 10:00 and 10:30 UT where the entropy instead appears to increase with the arrival of the enhanced enthalpy flux.These features are indicators of bubbles that did not reach 6 R E before reaching their equilibrium position, but compressed the magnetic flux Earthward of them, increasing the concentration of magnetic flux and pushing higher-entropy flux tubes to a lower L shell.
The bubbles in Figure 3 originate from the middle to distant plasma sheet and propagate to below GEO.Indeed, Figure 4 visualizes the flows associated with several bubble features during the 09:27:00-09:31:30 UT period.
The figure shows several spherical slices through the magnetotail, centered in Z on the equatorial plane, with the radial component of the bulk flow perpendicular to the magnetic field (v ⊥,r ) displayed as color contours on each.
The inner slice has a radius of 6 R E , equal to that of the nightside arc shown in Figure 3.While not shown in the figure, field lines have been traced from each grid point at each point in time, so the magnetic topology of the domain is known.Regions of each slice that are connected to non-closed field lines (i.e., open or IMF lines) have been darkened to convey the 3D topology and help to visualize the stretching of the tail.An animation of this figure is provided in Movie S1. Figure 4a shows a snapshot at 09:27 UT, where it can be seen on the furthest 2 slices, at 12 and 14 R E , that the tail is stretched and that fast flows are emerging within the current sheet.A large FB has an origin beyond 14 R E and penetrates through the 6 R E slice; this corresponds to the bubble feature shown in Figure 3 at 20 MLT at the same time.We note here that in this snapshot it may appear that there is a single large-scale inward flow spanning roughly 21-24 MLT at 12 R E , however it is clearer in Movie S1 that these are two distinct inward flows.Return flows can also be seen on either side of the most duskward bubble, spanning at least 2 R E in radial extent.In Figure 4b, a FB that does not penetrate 6 R E is highlighted.Figure 4c shows the strongest flow/bubble captured by the simulation, which reached 6 R E around 09:31 UT and also generated strong return flows.The shaded regions behind the flow's origin suggest this bubble was initiated by reconnection inside of 14 R E .This example demonstrates that not all bubbles in the simulation propagate in from the more distant central plasma sheet, and that they may not all be strictly BBFs by definition (Angelopoulos et al., 1992, which observed them between 9 and 19 R E ), but this is a byproduct of the strong storm conditions limiting the closed tail region within relatively close distances to Earth.In order to avoid any confusion due to terminology, going forward we will refer to these structures simply as bubbles, with the understanding that they all possess an inward flow, transport of plasma mass and energy as well as dipolarizing flux, and have a lower FTE relative to the local background, as shown in Figure 3.
Tailward vortical flows are observed on either side of almost every Earthward flow.Since the sign of the enthalpy flux is indicative of the flow direction, the return flows are evident from the regions of negative flux shown in Figure 3a.Such tailward flows have been identified previously in both numerical models (Birn et al., 2011;Merkin et al., 2019;Yang et al., 2019) and observationally (Panov et al., 2010).These tailward flows clearly transport a significant amount of plasma energy outwards and thus must be considered when quantifying the net contribution of bubbles to the buildup of the ring current energy, as is done in Section 3.2.flows that reach 6 R E transport energy largely in the form of plasma energy rather than magnetic or bulk kinetic energy.
A prominent feature in the enthalpy flux panel is the standing structure of multiple inward/outward flows between 09:30 and 10:00 UT.The structure is initiated by the previously-examined strong flow that arrives at 22 MLT, which kicks off a standing vortical pattern across several MLT that slowly drifts westward.The vortices are localized to just a couple of R E in the radial direction and do not contribute to continuous pumping of new energy from the tail.While this is an interesting feature, because it does not contribute to the net transport of plasma (seen more quantitatively later in Figure 5), we do not investigate it further as a part of this study.

Quantification of Plasma Energy Transport to the Ring Current
In order to quantify the contribution of bubbles to the buildup of the ring current energy content, we first construct a simple method of systematically detecting them.We define the quantity ΔB τ , which is calculated by subtracting the average B z value for the preceding τ minutes from the instantaneous B z value: where τ is a chosen time interval and l is the azumithal position along the nightside arc.With a τ on the order of minutes, this metric effectively picks out rapid enhancements of the magnetic field strength while filtering out longer timescale variations due to global dipolarization events.At the same time, this method also performs better at identifying dipolarization fronts than selection via a dB z /dt threshold, as the latter is more susceptible to fluctuations in B z that are not associated with bubbles.We find that the combined criteria of ΔB τ > 10 nT with τ = 2 min and V r < 0 (i.e., inward radial velocity), similar to the criteria used by Runov et al. (2021) to identify many dipolarization events within THEMIS observations, performs well in systematically selecting bubbles along the 6 R E arc within our simulation.For every point that meets this criteria, we also flag all points within l ± 1 R E arc-length of the point as being influenced by bubbles in order to account for the return flows brought on by the Earthward flow.Because we are triggering off of enhancements in B z , we also extend this flagged area out to 1 min before the time of the triggered feature in order to include any enthalpy flux enhancements due to compression ahead of the bubble.
The black boxed regions in Figure 5a mark the regions selected as being influenced by bubbles given the above criteria, overplotted on the enthalpy flux keogram identical to that in Figure 3a.This algorithm successfully picks out the most obvious instances of bubbles and their accompanying return flows.It is also made clearer here that a bubble at 09:30 UT and 22 MLT initiated the standing vortical structure, and that the standing structure is not a result of a continuous flow from the tail.In Figure 5b is the line integral of the enthalpy flux within regions flagged as being influenced by bubbles versus those that were not.Interestingly, while there are several bubbles throughout the initial period between 08:15 and 08:45 UT, inward plasma energy transport occurred mostly outside of these flows.There is even a surge of large-scale transport starting at 8:44 UT that leads to the greatest enthalpy flux enhancement at 08:55 UT.At this point, however, several bubbles arrive at 6 R E simultaneously and are the major contributors to this enhancement.Beyond this period, there are many instances where net inward flux is dominated by mesoscale flows rather than by global-scale convection.We also note here that transport outside of bubbles is not uniformly inward across the tail.
Figure 5c shows the value (blue) and time derivative (orange) of the total plasma energy on all closed field lines within 6 R E across all MLTs.The greatest peaks in the time derivative correlate well with the peaks in net enthalpy flux Figure 5b), indicating that the flows through the nightside arc are the main contributors to the increase of plasma energy within 6 R E , and thus the buildup of the ring current.
Separately summing up the net enthalpy flux within and outside of the boxed regions throughout the primary period of ring current buildup in the model (08:15-10:45 UT), we find that 50% of the net inward enthalpy flux through the 6 R E arc is caused by bubbles.We emphasize here that this is not the same as saying "bubbles contribute 50% of the total plasma energy to the ring current," as this analysis is only quantifying bubbles that penetrated as deeply as 6 R E .Excluded from this calculation is the consideration of plasma that ultimately made it to below 6 R E , but was only transported part of the way by bubbles.In other words, the 50% estimate serves as a lower bound on the total contribution of bubbles to the ring current energy for this simulated storm.This point is discussed in further detail in Section 4.
This result is not directly comparable to, but complements, those in Cramer et al. (2017).First, Cramer et al. (2017) quantified the number density flux whereas we are examining the transport of plasma energy, which is more directly linked to the ring current pressure.Second, Cramer et al. (2017) examined only the inward flux, and therefore did not account for the net transport due to return flows.In the presented simulation, for all of the plasma energy transported inwards by BBFs through the nightside arc, the return flows transported, on average, 40% of that amount outwards.Therefore we also conclude that the effects of return flows must be accounted for when quantifying the net contribution of bubbles to ring current buildup.However, we emphasize that this fraction was derived from this specific storm simulation, and should not be interpreted as a statement about the effect of return flows in general.

Evolution of the Ring Current Energy Spectra
We have considered, above, the processes that contribute to the buildup of the ring current through nightside transport.However, a deeper understanding of the stormtime ring current requires also an understanding of the buildup and evolution of the constituent particles across a broad range of energies.In this section we present an analysis of the proton energy spectra within the ring current region, beginning with a data-model comparison between the RCM proton intensities and that observed by RBSP-B, with the purpose of providing an in-depth understanding of which physical processes the model is accurately reproducing and those it is not.This is followed by an analysis of the evolution of the simulated proton energy spectra and the drivers of this evolution, including contribution from the source population, localized adiabatic heating, and losses due to charge exchange.

Data-Model Comparison: Proton Intensities
In this section, the proton intensities as observed by RBSP-B and calculated by the model are compared along the satellite trajectory.This is presented in Figures 6 and 7. Figure 6 shows the comparison for 10 energies between 10 and 100 keV, the energy range that constitutes the majority of the ring current pressure (Krimigis et al., 1985;Smith & Hoffman, 1973;Williams, 1985).Note that the range of the vertical axes of each panel are different in order to better resolve the comparison for each plot.A 5-min smoothing window has been applied to the RCM intensities to remove noise.The comparison begins when RBSP-B is at apogee of the orbit during which the SSC occurs.We refer to this orbit as Orbit 1, and the following orbits as Orbits 2 and 3. Included in Figure 6 are annotations denoting the phase of the orbit and azimuthal location of the spacecraft at that time.Figures 7a and 7b shows a comparison of the observed and simulated proton intensities for the continuous energy spectra between 10 and 200 keV.Figure 7c is a Dst plot equivalent to that in Figure 1a. Figure 7d shows the RCM pressure in the equatorial plane, as well as the spacecraft trajectory over the duration of the simulation, and its current location at 10:09 UT.Note that the trajectory is not a perfectly elliptic orbit; this is because the trajectory is represented by the spacecraft's 3D position mapped to the equator along the field lines traced through the model magnetic field.Vertical lines in the other panels correspond to the time shown in the top right corner of the Figure 7d.This figure is a snapshot of the animation provided as Movie S2.
Starting with Orbit 1, we observe that while the initial RCM population is relatively similar to that measured by RBSP-B, there are significant differences.SSC occurs during the inbound pass of this orbit (highlighted in Figure 7a).The intensity enhancement due to SSC compression is underestimated in the model, which can be attributed to a number of factors, including the initial energy spectrum in the model (constructed as described in Section 2.3) being at a variance with the observed spectrum, limitations in the solar wind reconstruction (Section 2.1) at producing the SSC shock (Figure 7c), and the uncertainty in the propagation of the upstream solar wind.As will be clear in the following analysis of the evolution of the energy spectra (Section 3.3.2),the pre-storm plasma in this region contributed relatively little energy to the final ring current.However, this initial population is still important because it influences the degree of stretching in the tail which in turn affects where the initial reconnection will occur, and also sets the background entropy profile that the initial bubbles must penetrate through.This is discussed further in Section 4. As discussed above, the initial ring current buildup occurred between 08:15 and 10:45 UT. Figure 7b indicates that, for almost the entirety of this period, RBSP-B was in the perigee phase of its orbit, where the spacecraft was below 3 R E and RBSPICE was not taking time-of-flight measurements due to the high voltage mode being turned off.Highlighted in Figure 7d are two Earthward penetrating bubbles that manifest in this plot as two deformations of the RCM active domain due to their fast flow speeds (see the discussion of the coupling algorithm in Section 2.3).The radial velocity and magnetic field dipolarization within these bubbles are visible in Figures 2c  and 2d, and the enthalpy flux and depleted entropy components are visible in Figures 3a and 3c.These are the only bubbles during the initial buildup phase that the virtual spacecraft has a chance of observing as an injection, however it occurs just as the spacecraft encounters the inner edge of the ring current, masking their signatures.
However, throughout the outbound pass of Orbit 2 (between 10 and 14 UT as the spacecraft passes from ∼21 MLT to just past 00 MLT), the observed and simulated proton intensities are in strong agreement.This is manifested in the timing of RCM and RBSP-B's observation of the ring current along the outbound pass and indicates that the model is capturing the location of the peak ring current quite well for this storm.The intensities at this leading edge are in good agreement in general (Figure 6), though there is an overestimate at 20 and 30 keV, and an underestimate above 80 keV.Throughout the remainder of the Orbit 2 outbound pass, the intensities agree especially well between 30 and 70 keV, the best agreement being at 50 keV with a maximum ratio between the modeled and observed intensity of 5.2, and an average ratio of 1.3.Qualitatively, injection-like features in the observations are also present in the simulated spectra throughout the 10-14 UT period, although inspection of Movie S2 shows that not all of these features are unambiguously caused by bubbles.Visible in the 60-90 keV range in Figure 6 is a drop-out in intensity near Orbit 2 apogee.It is clearer in Figure 7 that there is a general lack of a ∼5 × 10 4 cm −2 sr −1 s −1 keV −1 intensity background above 50 keV as seen in the data.A number of factors may contribute to the model missing this population, such as the initial plasma sheet being too cold, the model not including non-adiabatic processes, or a "boundary shadowing" effect, whereby particles of the relevant energies drift out of the RCM domain, which is located a few R E inward of the magnetopause.These particles are lost from the simulated ring current, whereas in reality they would have stayed trapped within the inner magnetosphere.
While the intensities at individual energies are still in good agreement for the inbound half of Orbit 2, there is a major change in the structure of the modeled energy spectra starting just before 14:00 UT (highlighted in Figure 7b).At this time, reconnection occurs close to −10 R E across the tail in the simulation, restricting the active RCM domain to below 6 R E (visible in Movie S2).As the RCM active domain moves back out, it re-populates its plasma population with a fresh Kappa distribution using the MHD moments at these locations (Section 2.3), resulting in the smoothed energy spectra seen for most of the inbound pass.
At 16:00, there are two strong injection signatures in the data, whereas the model shows two depletions in intensity at this time.It can bee seen in Movie S2 that these depletion signatures are from the virtual spacecraft passing through the wake of two bubbles that have already penetrated to a lower radial distance.There is also a stripelike feature in the modeled intensities around 60-150 keV near the end of Orbit 2. This is due to the injection related to the bubble that passed 6 R E at 10:30 UT around 20 MLT (see Figure 3) that has become significantly dispersed, so much so that it has wrapped around the Earth several times by 16:00 UT.In nature, this would not remain as such a coherent structure, as drift shell splitting would separate the different pitch angles over time (Roederer, 1967), whereas the isotropic approximation in RCM does not capture this effect.
At the end of the inbound pass of Orbit 2, the spacecraft passes through the ring current again around 3 MLT.
(Figure 6).The model computes a greater local maximum in intensity for 40-80 keV protons in the 15 to 17 UT range, but also a dropoff in intensity earlier in the orbit than is observed.This implies that the modeled ring current within the post-midnight sector is stronger but also further from Earth than in reality.The proton population in this region is dependent on a number of interacting factors, including: the location and energy of the source population of ring current protons; the competition between Eastward E × B (i.e., convection and co-rotation) drifts and energy-dependent Westward gradient-curvature drifts; and the properties of the ring current itself that contribute to the ionospheric potential through shielding via Region 2 currents and the conductance profile via electron precipitation.Because of these interdependent complexities, the divergence of the model from observations in this region is not unexpected, but the comparison yields key insights about the model's ability to capture these different factors.

Energy Spectra Evolution
Having discussed the data-model comparisons, we now turn to a more detailed investigation of the ring current energy spectra in the model.Figure 8 shows the cumulative pressure as a function of energy and time at three stationary points indicated in panels (d-g).The first point (at 3.8 R E and 21 MLT) marks the peak ring current pressure which occurred at 09:15 UT.The second and third points are located at 4.3 and 4.8 R E at the same MLT in order to sample the outer extent of the ring current.The white lines in Figures 8a-8c denote the median and quartile energies of the cumulative pressure fraction, that is, the proton energies lower than the 50% line collectively contribute half of the total pressure at that point and time.For ease of reference, we refer to this energy as the "median energy."The cumulative pressure value at the top of Figures 8a-8c is the total pressure for each point in time.Figures 8d-8g show the total RCM pressure in the equatorial plane for four select times throughout the plotted period indicated by cyan vertical lines in Figures 8a-8c.
Preceding the main period of ring current buildup (before 9:00 UT, see Figure 5), the median energy at 4.3 and 4.8 R E drops from 60 to 30 keV, respectively, to 20 keV.This is best explained by convection of the pre-existing plasma sheet (recall from Figure 5 that prior to about 08:50 UT, inward transport was dominated by convection with only several bubbles present).However, only at 09:00 UT, just after significant transport via bubbles commences, does new plasma reach the 3.8 R E probe point, where the drop to a median energy of 20 keV is accompanied by the greatest ring current pressure captured by the model.Over the course of many hours, the median energy increases to 70-85 keV, depending on the location, while the total pressure at each location decreases on average.However, the relatively constant DPS-Dst value throughout this time indicates that the total energy content within 6 R E is not appreciably changing.
We now focus our attention on the evolution of the median energy from <20 to 85 keV between 09:00 and 19:00 UT at the 3.8 R E probe location.This evolution may be caused by a localized heating of the already-present plasma population via adiabatic energization (i.e., betatron acceleration due to magnetic field enhancements), the introduction of hotter plasma from the tail later into the storm, and/or by energy-dependent loss processes.Figure 9 examines each of these processes.Figure 9a shows a time-series of the value V −2/3 , where V is the flux tube volume, evaluated at the 3.8 R E probe location.Given the equation W = λV −2/3 , where W is the kinetic energy and λ is the isotropic energy invariant (Wolf et al., 2009), changes in V −2/3 are directly proportional to energization via localized adiabatic heating.Figure 9b shows the energy flux through the nightside 6 R E arc, identical to that in Figure 3d but extended in time.Figure 9d shows the proton intensity as a function of time and energy, and Figure 9e shows a time-series of the proton intensity for the select energies of 10 and 50 keV, both evaluated at the 3.8 R E probe location.Note the log scale on the y axis.
After the initial ring current buildup, the value of V −2/3 (Figure 9a) is fairly noisy but gradually increases by less than 10%.Therefore it contributes minimally to energization over time.The initial buildup of the ring current (cumulative pressure in Figure 9c at 09:00 UT) expectedly correlates well with the period of enhanced inward plasma transport examined in Section 3.2.Just after 10:00 UT there is an enhancement in the total pressure as well as a jump in the median energy from 20 to 35 keV, and at 10:30 UT there is another jump to 40 keV (denoted with cyan arrows in Figure 9c).There are two other notable periods of enhanced transport.One is around 13:40 UT and is related to the close-in reconnection event discussed in Section 3.3.1,that results in an increase in pressure but does not notably alter the energy spectra.The other, starting around 15:00 UT, results in a lesser increase in pressure but more appreciable increase to the median energy.These instances of enhancements to the median energy and/or pressure at the 3.8 R E probe point all correlate well with periods of enhanced inward transport, implying that these changes in the energy spectra are related to transport of hotter plasma from the plasma sheet.
The intensity of the 50 keV protons (Figures 9d and 9e) remained relatively constant after the initial buildup at 09:00 UT, whereas the 10 keV population steadily decreases over time.Using the Østgaard et al. (2003) geocorona model and Lindsay and Stebbings (2005) proton cross sections (as implemented in RCM, see Section 2.3), 10 keV protons with a 45° pitch angle at 3.8 R E have a decay timescale of roughly 1.6 hr.Thus, charge exchange is capable of depleting the 10 keV protons by a factor of more than 250 over the nine hour period after the initial ring current buildup, whereas the 50 keV protons are expected to decay by a factor of 13 over the same period.
According to Figure 9e, the 10 keV protons decreased by a factor of 200 and the 50 keV protons by a factor of 13 between their max value near 09 UT and their minimum value after 18 UT, consistent with the charge exchange decay rates.In total, the evolution of the median energy from 20 to 85 keV over the course of 9 hr is best explained by a combination of impulsive adjustments due to the transport of hotter plasma from the plasma sheet and a steady increase in the median energy due to energy-dependent charge exchange with the exosphere neutrals.

Discussion
While previous modeling studies have provided important insights in transport and ring current buildup via bubbles (e.g., Cramer et al., 2017;Sorathia et al., 2021;Yang et al., 2015), this paper is the first to target directly the quantification of the buildup of the ring current plasma energy content via bubbles generated self-consistently within a global geospace model.A major difference between the simulation in this study and those of Yang et al. (2015) is that in RCM-E, bubbles are generated via a reduction in pressure along a fixed outer boundary, with properties constrained by observations, whereas in this study the bubbles' generation location, pressure, flux tube volume, and flow speed are produced self-consistently by reconnection in the magnetotail (e.g., Wiltberger et al., 2015).Nevertheless, the 60% contribution from bubbles found by Yang et al. (2015) is quite similar to the lower bound of 50% calculated in this study.We emphasize again that the inclusion of the return flows in the calculation of the net transport via bubbles is crucial, as they transport a considerable amount of plasma energy outwards.A caveat to the results presented here is that the occurrence rate of highly-depleted bubbles is known to increase as grid resolution increases, up to a resolution higher than that used in this study (2× finer resolution in each dimension; Sorathia et al., 2021).It is possible that at a higher resolution, our estimate of a lower bound of 50% contribution from bubbles would be even higher.Also, this estimate is derived for this specific storm, and how this estimate changes based on specific storm characteristics is a matter for future multi-event study (cf., Cramer et al., 2017).
Since fast flows are excluded from the RCM domain, and because GAMERA does not include heat flux on its own, the FTE of bubbles remain constant while they are being transported rapidly inwards.Particle tracing of ions (Ukhorskiy et al., 2018) and electrons (Gabrielse et al., 2017;Sorathia et al., 2018) within bubbles has demonstrated that their enhanced magnetic flux is capable of trapping relatively high-energy particles via gradient-curvature drifts whereas lower-energy particles may drift in and out of the bubble.This means that in reality the FTE is not a conserved quantity as the bubble is transported inwards.While the statistical properties of bubbles produced in GAMERA are in good agreement with observations (Sorathia et al., 2021), the interaction between the bubble and plasma at different energies is currently not captured by the model.In addition to this, the fairly conservative constraint on the RCM active domain based on the distance to the open/closed field boundary also led to the boundary incursion on the RBSP orbit, resetting the energy spectrum at that location to a smooth Kappa distribution (Figure 7b).While this did not greatly affect the intensities across the spectrum (inbound pass of Orbit 2 in Figure 6), it is unclear what effects this removal of structure may play on energy-dependent dynamics afterward.These limitations arise from the fact that inner magnetosphere models do not appropriately handle fast flows, while the MHD model does not include gradient-curvature drifts in the bulk velocity.A model that is capable of handling both of these simultaneously is needed in order to understand how the inclusion of these effects might alter bubble penetration depth and its implications on particle injections.
The peak ring current in the simulation (at 3.8 R E , Figure 8a) was composed of protons with a median energy of about 20 keV, consistent with the median energy seen earlier at 4.8 and 4.3 R E , whose origin is primarily the pre-storm plasma sheet population that has been transported inwards.After the peak ring current pressure around 09:15 UT followed a series of impulsive increases to the median energy caused by injections of hotter plasma.A multi-phase evolution of the energy spectra has been observed for a number of storms (e.g., Goldstein et al., 2017;Keika et al., 2018;Zhao et al., 2015).Keika et al. (2018), in their analysis of the 17 March 2015 storm, similarly attributed the initial stormtime ring current population to the inward transport of the pre-storm plasma sheet, with the later increases in median energy being caused by the transport of a hot, dense plasma sheet population.Given the importance of the pre-storm plasma sheet population to the initial ring current energy spectra, it is clear that the proper "pre-conditioning" of the simulation, that is, setting the pre-storm magnetosphere state including the configuration of the magnetotail and the population within the magnetosphere, is essential to capturing observed stormtime dynamics.As briefly mentioned above, properly modeling bubble penetration relies on properly capturing both the entropy profile of the background plasma and the entropy within the bubble.The pre-conditioning method employed in this simulation has performed well in reproducing the observed intensities in general, but the overestimate in the 20-30 keV range and underestimate in the 80-100 keV range (Figure 6) suggests that the method could be improved.This can be achieved for example, by constraining with spacecraft that traverse the outer plasma sheet such as THEMIS or by running the simulation for an even longer period before the storm onset to allow for a more self-consistent plasma sheet population to form.The question of the origin of the hotter plasma that was transported inwards at the later stage of the storm is beyond the scope of this paper but can be investigated with our model in the future (cf., Keika et al., 2018, for the March 2015 storm).The median energy was also strongly influenced by charge exchange with the exosphere.This study used the Østgaard et al. (2003) geocorona model to calculate charge exchange rates.However, it is known that the decay rate is highly dependent on not only the choice of exosphere model (Ilie et al., 2013), but also on the variation of the neutral density profile over the duration of geomagnetic storms (Cucho-Padin & Waldrop, 2019).The relative role of the transport of an evolving plasma sheet versus charge exchange on the evolution of the ring current energy distribution is highly dependent on these factors, and is likely to change on a storm-by-storm basis.
Using the HOPE and RBSPICE instruments aboard the Van Allen Probes, Menz et al. (2017)  and bubbles, and the resulting impact it has on ring current buildup and evolution, would most appropriately be explored with a global magnetosphere model that is not only coupled with an ionospheric outflow model (e.g., Glocer et al., 2009;Varney et al., 2016), but one that also includes feedback from non-adiabatic O + .
With the above caveats in mind, our answer to the question, "how much of the plasma energy transport through 6 R E is in the form of bubbles?" is ∼50% for this simulation.However, we note that this question is different from asking, "what is the relative contribution of bubbles to ring current buildup compared to large-scale convection?"In a quantitative sense, it is up for interpretation how much a bubble "contributed" to ring current buildup in the case where a particle ultimately makes it to the ring current but was only transported part of the way by a bubble.
A hypothetical example is the case where a bubble transports plasma up to its corresponding penetration depth outside of 6 R E , at which point global-scale convection, radial diffusion, or even compression ahead of a subsequent bubble is responsible of the remainder of the transport.Yang et al. (2015) seeded test particles in the final state of each of their simulations, back-traced them to the boundary, and determined if it entered the simulation within a bubble or through the global-scale flow.In this analysis, a bubble was considered to have contributed to the ring current even if it transported plasma only a short distance into the simulation domain, so long as the plasma remained in the domain by the end of the simulation.In this study, we examined the net enthalpy flux due to bubbles through 6 R E , including the effects of compression ahead of the bubbles, which may have resulted in plasma at certain energies being transported below their respective Alfven layer, where they otherwise would have remained outside of it in the absence of the bubble.Because the bubbles examined in this study had a guaranteed and measurable contribution to the ring current energy content, and because it did not include the complexities of bubble contribution outside of 6 R E , this analysis definitively calculates a lower bound to bubble contribution, in the context of this specific simulation.The overall role of bubbles in ring current buildup remains only partially understood, and further investigation requires not only more complex analysis methods, such as fluid-element and particle tracing, but also a more deep understanding of the physics of particle transport from the plasma sheet to the ring current and the development of the appropriate physics-based models.

Conclusion
Via a case study of the 17 March 2013 storm, we examined the contribution of BBFs/entropy-depleted bubbles to the buildup of the stormtime ring current, and the evolution of the ring current energy spectra throughout the main phase of the storm.The primary conclusions of this study are as follows: 1. Energy transport through 6 R E in the simulation is primarily in the form of enthalpy flux, followed by magnetic flux and finally by the bulk kinetic flux.In other words, by 6 R E the energy transported by bubbles is primarily in the form of plasma energy.
2. At 6 R E , 50% of the net inward enthalpy flux is contained within bubbles.This serves as a lower bound to the total contribution of plasma transported by bubbles to the buildup of the stormtime ring current, as it only considers bubbles that have penetrated to below 6 R E , which constitutes a fraction of the total number of bubbles produced by the model.3. The return flows that accompany bubbles as a result of interchange transport outwards an average of 40% of the plasma energy transported through 6 R E by the bubbles, thus their effect must be considered when quantifying the net contribution of bubbles to inward energy transport.4. The data-model comparison shows that the model is accurately capturing both the radial location and magnitude, as well as the spectral properties of the initial ring current well, implying that the model is accurately reproducing the process of ring current buildup.Possible explanations for discrepancies between the data and model include inaccuracy in the initial plasma sheet conditions, boundary shadowing of higher-energy protons as a result of a conservatively constrained RCM active domain, and the absence of O + outflow and its effects on plasma sheet composition and bubble penetration depth. 5.The evolution of the ring current energy spectra is caused by a combination of impulsive adjustments due to transport of a gradually hotter plasma sheet population and a steady increase due to energy-dependent charge exchange.
The use of state-of-the-art geospace models in combination with in-situ data continues to be a powerful tool in gaining a deeper understanding of the role of bubbles in plasma and magnetic flux transport during periods of high geomagnetic activity.While this study contributes to better understanding the drivers of stormtime ring current formation, there are many outstanding questions, such as how these results are altered depending on the solar wind driver, the more general relationship between mesoscale and global scale transport in the magnetotail, and the evolution of the plasma sheet source throughout the storm and where and how heating occurs.Deep analysis of model results enables more thorough validation via comparison to observations, and in turn increases the physical insight that may be gained through modeling.

Data Availability Statement
The  Sciola et al., 2023), which includes all information necessary to generate all Figures, except for the 3D field line data in Figure 1.The repository also includes a conda (https:// docs.conda.io)requirements file which may be used to recreate the python environment used to make the Figures generated using Matplotlib, as well as an example plotting script.

Figure 1 .
Figure 1.3D view of the simulation at 09:28 UT.The radial velocity is shown as red and blue color contours in the equatorial plane.On the same surface, ΔB z is shown as colored contour lines in 5 nT intervals, and the open/closed field boundary is shown as a black line.Overplotted in the equatorial plane is the Rice Convection Model (RCM) pressure within the RCM active domain.Field lines, seeded along constant latitudes near Earth, are shown as translucent tubes, with a cross-section scaled by |B| −1 .The field lines terminate at the Grid Agnostic MHD for Extended Research Application inner boundary at 2 R E .On the Earth's surface, downward and upward ionospheric field-aligned currents (FACs) are shown as purple and green color countours, respectively.Overlaid to the right of the Earth is another view of the northern ionosphere FACs from a perspective along SM-Z.Overlaid in the bottom left is a timeseries of the observed Sym-H, where the period of the storm that is simulated is shaded in orange, and the black vertical line is 09:28 UT, the time shown in the 3D view.

Figure 2 .
Figure 2. (a) Time-series of the observed Sym-H, the simulation Dst calculated using Biot-Savart evaluated at the center of Earth, and the DPS-Dst for all closed field lines within 6 R E .Equatorial view of the inner magnetosphere pressure (row b), v r (row c), and ΔB z (row d) at four select times corresponding to the purple dashed lines in (a).The gray dashed circle with a radius of 6 R E is included for reference.

Figure 3 .
Figure 3. Simulated quantities evaluated along a 180°arc across the nightside with a radius of 6 R E , including: (a) The radial component of the enthalpy flux, (b) ΔB z , and (c) flux tube entropy.In the enthalpy flux panel, positive values indicate radially inward flux, and negative values indicate radially outward flux.(d) Line-integral of the radial component of the enthalpy, Poynting, bulk kinetic, and total energy flux across the entire arc.

Figure 4 .
Figure 4. Spherical slices through the magnetohydrodynamic solution from 18 to 3 MLT at three select times spanning 4.5 min.The slices span 6-14 R E in radial distance, spaced 2 R E apart, and are centered in Z on the XY-plane.The color contours on the slices show the radial component of the bulk velocity perpendicular to the magnetic field (v ⊥,r ).The darker regions on the slices denote regions that are on non-closed field lines (either open or fully interplanetary magnetic field field).Each panel is marked to show flow paths through the slices.Flow vectors have been included in the animated version of this figure (see Movie S1).The Earth is also shown, with the inward/outward field-aligned currents shown in purple and green, respectively.The red line in the ionosphere marks the open-closed boundary determined via field line tracing.

Figure 5 .
Figure 5. Time-series analysis of energy flux through a 180°arc across the nightside hemisphere at 6 R E : (a) Radial component of the enthalpy flux through the nightside arc, the same as in Figure 3a.Black boxes outline the regions determined to contain bubbles; (b) Line-integrated enthalpy flux from bubbles, outside of bubbles, and in total; (c) The instantaneous value (blue line) and time-derivative (orange line) of the total plasma energy contained on all closed flux tubes within 6 R E .

Figure 6 .
Figure 6.Comparison of the proton intensities during the 17 March 2013 storm between the RBSP-B's RBSPICE instrument measurements (orange) and the Multiscale Atmosphere Geospace Environment simulation, as calculated by Rice Convection Model (RCM) (purple).The 10 panels show a comparison, each at a fixed energy denoted in the top right corner of the corresponding panel, between 10 and 100 keV.The comparison begins at 4:00 UT and ends at 19:00 UT.The RCM intensities have been smoothed over 5 min.

Figure 7 .
Figure 7.Comparison between (a) RBSPICE and (b) Rice Convection Model (RCM) proton intensities for the full energy spectra between 10 and 200 keV.(c) Comparison between the Observed Sym-H (blue), that calculated from the simulation (orange), and the DPS-Dst evaluated in the model within 6 R E (green), the same as in Figure 1.(d) RCM pressure and RBSP-B's full trajectory (orange line) and current position (white circle with orange border) mapped along field lines to the equatorial plane.The vertical lines in panels (a-c), and pressure and spacecraft location in panel (d), all correspond to the time indicated in panel (d).The thin colored bar below panel (b) shows the RCM pressure at the spacecraft location for a given time, with the same colorbar as panel (d).

Figure 8 .
Figure 8. (a-c) Rice Convection Model (RCM) cumulative pressure as a function of energy and time at three stationary radial distances, denoted in the top right in each panel, along 21 MLT.The white lines denote the median and quartile values of the cumulative pressure fraction.The cyan lines correspond to the four times shown in panels (d-g).(d-g) Equatorial views of the RCM total pressure at four select times.The three dots denote the locations at which the pressure in panels (a-c) is evaluated.

Figure 9 .
Figure 9. (a) Time-series of the quantity V −2/3 .(b) Time-series of the net energy flux through the nightside 6 R E arc, identical to Figure 3d but over a longer period of time.(c) Cumulative pressure as a function of energy and time, identical to that in Figure 8a.(d) Proton intensity as a function of energy and time.(e) Time-series of proton intensities at 10 and 50 keV.Cyan arrows denote instances of sharp increases in the median energy and/or the total pressure.
Sorathia et al. (2021)2005;Smith & Bewtra, 1978) 65% to the total ring current pressure in the L = 3-4 range, and decreased to a 50% contribution by L = 5.Inclusion of O+ in our simulation could affect both the buildup and decay of the ring current.For the latter, we would not expect the model without O+ to reproduce the observed recovery rate due to the much faster charge exchange rate of O + relative to H +(Lindsay & Stebbings, 2005;Smith & Bewtra, 1978).It is less clear how O + would influence the relative contribution of bubbles to ring current buildup, which was the focus of this study.The effect of O + outflow on ring current buildup has been studied previously using global MHD models coupled to ionospheric O + outflow and inner magnetosphere models (e.g.,Welling et al., 2011).However,Sorathia et al. (2021)used test-particle simulations to demonstrate that plasma sheet flows do not transport O + as efficiently as they do lighter ions because of significant non-adiabatic effects for the heavy species.This resulted in shallower penetration depths of O + than if these particles were constrained to adiabatic, guiding-center motion.Because the study did not include feedback from the test particles on the MHD solution, it remains unclear to what extent non-adiabatic O + may in turn affect transport by bubbles in general.Therefore the interaction between O + solar wind data source is primarily Modified High Resolution OMNI (HRO) data (https://omniweb.gsfc.nasa.gov/form/omni_min_def.html).Data gaps were filled with KP-despiked Wind data from OMNIWeb (https:// omniweb.gsfc.nasa.gov).Proton data from RBSPICE instrument aboard the Van Allen Probe B spacecraft was retrieved from cdaweb (https://cdaweb.gsfc.nasa.gov)via the cdasws python package (https://pypi.org/project/cdasws/).The model data was produced with the Multiscale Atmosphere Geospace Environment (MAGE) model, under development by the NASA DRIVE Science Center for Geospace Storms (https://cgs.jhuapl.edu).Figures1 and 3were made using the ParaView visualization software (https://www.paraview.org).All other figures were made using Matplotlib (https://matplotlib.org/).A subset of the model output data has been made available (https://doi.org/10.5281/zenodo.7921979;