This study analyzes the intense geomagnetic storm that was triggered by the arrival of a coronal mass ejection (CME) and began with a sudden commencement on 1 June 2013, around 0300 UT, with a main phase that lasted for about 6 h until ∼ 0900 UT, followed by a long recovery phase that extended over at least 4 days. This storm has been selected due to the coincident availability of data from TWINS/LAD, which is needed to estimate exospheric H distributions, and Van Allen Probes instruments used to perform data-model comparisons. Figure 1 shows an overview of the geomagnetic activity, global thermospheric/exospheric conditions, and interplanetary environment during the storm’s development. Panels from top to bottom show the time series of (a) the horizontal component asymmetry index (SYM-H), (b) the globally-averaged exobase temperature ( T exo ) , (c) the globally-averaged exobase density ( n Hexo ) , and (d) the planetary K-index (Kp). Values from panels (b) and (c) were obtained from the National Center for Atmospheric Research (NCAR) Whole Atmosphere Community Climate Model with thermosphere and ionosphere extension (WACCM-X). The gray-shaded interval corresponds to the period of study in which the ring current is simulated and the terrestrial exosphere is estimated from NASA’s TWINS radiance data.

The variability of exospheric densities at ring current altitudes during a geomagnetic storm is typically associated with the global distribution and temporal variability of the temperature and densities at the exobase ( ∼ 500 km) Chamberlain (1963). During the storm’s development, ion precipitation slowly increases the thermospheric temperature via the Joule heating process, propagating heat from the poles to the equator (Zesta and Oliveira, 2019). Light H atoms are highly sensitive to these temperature variations, resulting in a density re-distribution where atomic H density decreases near the exobase ( ∼ 500 km altitude) and increases at high altitudes (ring current region) (Qin et al., 2017). This anti-correlation between density and temperature at the exobase is shown in panels (b) and (c). Furthermore, the energy partition of exospheric H atoms is continuously altered by plasmaspheric ions via charge exchange interactions, especially by protons ( H + ) with an energy range of ∼ 1 − 2 eV. During geomagnetic storms, the plasmaspheric density decreases, and its external boundary, the plasmapause, moves inward, thus reducing the rate of charge exchange and varying the 3-D structure of the exosphere (Kuwabara et al., 2017). Panel (d) shows the Kp index, which is linearly correlated to the plasmapause location (Moldwin et al., 2002) and serves as another indicator of exospheric density variability.

We estimate global exospheric densities using radiance data from the TWINS mission. NASA’s TWINS is comprised of two spacecraft, TWINS1 and TWINS2, each of them featuring a highly elliptical orbit, an apogee in the Northern ecliptic hemisphere of 7.2 R E geocentric distance, and an orbital precession of 1 year. Each satellite has two Lyman-Alpha Detectors, LAD1 and LAD2, able to acquire FUV photon flux with a central wavelength of 121.6 nm (Ly- α ) (McComas et al., 2009). The LADs are single-pixel optical sensors with a 4 ° field-of-view (FOV) assembled on a rotating platform in the spacecraft. Due to its geometrical configuration, the LADs acquire column-integrated photon flux measurements above ∼ 2 R E geocentric distance every ∼ 0.66 s.

To analyze the storm that occurred on 1 June 2013, we use radiance data acquired by LAD1/2 onboard the TWINS1 spacecraft for 3 days from May 31 to June 2. The LADs onboard TWINS2 were out of commission after 2012.

In this study, we estimate time-dependent, global H density distributions from Ly- α radiance data for the “optically thin region” of the exosphere (beyond 3 R E geocentric distance), where the density is sufficiently low that the scattering of Ly- α photons by H atoms occurs only once. This condition establishes a linear relationship between the detected photon flux and the volumetric atomic H density according to the following equation: I r , n ̂ , t = g * 1 0 6 ∫ l = 0 l = L max n H l Ψ β d l + I IPH , (1) where I is the measured Ly- α photon flux along a line-of-sight (LOS) in units of Rayleigh (1 R = 1 0 6 / 4 π photons/ cm 2 /sec/str), r indicates the 3-D spatial location of the optical detector, n ̂ is the direction of the detector LOS, and t is the time of the acquisition. Also, the factor g * ( sec − 1 ) is the local scattering rate that includes the solar Ly- α flux at time t and the cross-section of the atom-photon interaction, n H is the volumetric H density along the LOS in units of [atoms/cc], Ψ is the scattering phase function that handles the anisotropic nature of the photon re-emission (Brandt and Chamberlain, 1959), and β is the angle between n ̂ and the direction of the solar Ly- α flux (- x ̂ GSE ). The line integral is evaluated from the satellite position ( l = 0 ) to an upper boundary where scattered Ly- α photon contribution to I is considered negligible ( l = L max ) . Further, the term I IPH represents the interplanetary Ly- α background radiance emitted by H atoms in the boundary of the heliosphere.

Radiance data from LAD1/2 onboard TWINS1 are used to generate a dataset that should be pre-processed and filtered according to the following constraints:

• Re-emission from the lower exosphere known as Earth’s Albedo is a secondary source of Ly- α photons that violate the linearity between photon flux and H density expressed in Equation 1, especially near the 3 R E boundary. Therefore, in our study, LAD measurements whose LOSs have an impact distance (perpendicular distance from the center of Earth to the LOS) lower than 3.75 R E are discarded from the analysis,

We use a tomographic approach to convert 1-D radiance measurements from LADs into 3-D hydrogen density distributions beyond 3 R E , following the procedure described in Cucho-Padin and Waldrop (2018), Cucho-Padin et al. (2022). First, we select the solution domain as a spherical region spanning radii from 3 to 20 R E . Then, we divide this region into regular spherical voxels with dimensions (radial, azimuthal, and latitudinal) Δ r = 0.25 R E , Δ ϕ = Δ θ = 12 ° which provide sufficient spatial resolution to capture gradients in the exosphere. The total number of voxels is N = 39,600 . Also, we define the [ N × 1 ] vector x that contains the H density values for each voxel. Second, a [ M × 1 ] vector y of background-free measurements is created such that each element follows y m ( r , n ̂ , t ) = I ( r , n ̂ , t ) − I IPH ( r , n ̂ , t ) where I is the radiance acquired by a LAD and I IPH is the corresponding background Ly- α radiance estimated using the model developed by Zoennchen et al. (2013). Third, a [ M × N ] observation matrix L is generated by intersecting the m t h LAD’s LOS with the solution domain, producing line sectors in some spherical voxels. The length of a line sector within the n t h voxel multiplied by its corresponding factor g * ( t ) Ψ ( β ) / 1 0 6 is used to populate the element ( m , n ) of L . Hence, this procedure discretizes Equation 1 to yield a simple algebraic system y = L x where y is formed with LAD/TWINS data, the L matrix can be created a priori with the satellite position, LOS direction, and voxels sizes, and the vector of H densities x is the only unknown [see more technical details in Cucho-Padin et al. (2022)].

The time-dependent reconstruction of the exosphere requires solving the dynamic system y k = L k x k , where radiance measurements and the acquisition geometry vary each time step k . To solve this dynamic inverse problem, we closely follow the approach thoroughly described in Cucho-Padin et al. (2024) in the context of soft X-ray tomography. This methodology is based on Kalman filtering (KF) theory, which uses a set of observations of a variable ( y k ) as well as the previous state of that variable ( x k − 1 ) to predict its next state ( x k ) . The numerical implementation of this process is expressed as: x ̂ k = L k T R − 1 L k + Q k − 1 − 1 − 1 L k T R − 1 y k + Q k − 1 − 1 x ̂ k − 1 (2) where x ̂ k is the estimate of H densities at time k and R is the [ M × M ] covariance matrix of measurements defined as a diagonal matrix whose elements are identical to y k (i.e., R k = diag ( y k ) ). Also, the term Q k − 1 − 1 is known as the [ N × N ] precision matrix, which is formed following the Gaussian Markov Random Field (GMRF) theory described in Cucho-Padin et al. (2022) and serves (a) to provide statistical information (covariance among voxels) from the previous estimate x k − 1 to the next one and (b) to impose smoothness to the solution through minimization of the first and second derivatives of x k − 1 . The reader is referred to Cucho-Padin et al. (2022), Cucho-Padin et al. (2024), Zoennchen et al. (2024) for technical details in the implementation of Equation 2.

Tomographic reconstructions of exospheric density distribution are conducted “hourly” (i.e., the cadence of k is 1 h) using available LAD radiance data in this period. The solar Ly- α flux needed to calculate the scattering rate g * ( t ) is obtained daily from the Solar Extreme Experiment (SEE) instrument onboard NASA’s Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED) mission.

Since ring current modeling needs H densities below 3.75 R E geocentric distance, we extrapolate our data-based model to low altitudes using H density estimates from the WACCM-X model. Specifically, we run the WACCM-X during the period from May 31 to 2 June 2013. Then, for each hour, we select H densities spanning altitudes from ∼ 100 to ∼ 500 km at latitudes and longitudes corresponding to our high-altitude grid (i.e., 12° resolution in the GSE coordinate system). Hydrogen density values from our reconstruction from 3.75 to 10 R E and those from the WACMM-X model corresponding to 100–500 km altitude are fit together using an order four spherical harmonic function, which facilitates their inclusion into the ring current model.

The first column of Figure 2 shows the time-dependent, global tomographic reconstruction of the exospheric n H during the 1 June 2013 storm (see Section 2.3). It displays unfolded radial shells of H density with a geocentric radius of 5.5 R E (well into the ring current location) in geocentric solar ecliptic (GSE) coordinates. We selected six periods of time that cover the pre-storm conditions (31 May 1200 UT), the main phase of the storm (1 June 0500 UT/1 June 0900 UT), and the recovery phase (1 June 2300 UT/2 June 0400 UT). All radial shells show two strong H-density structures, a dayside nose near the ecliptic plane and a nightside tail, which are mainly produced by solar radiation pressure, a force exerted over H atoms due to their constant interactions with solar Ly- α photons (Beth et al., 2014). The second column of Figure 2 shows the H density maps obtained from the WACMM-X simulation during the same periods at the exobase (500 km altitude). It is noteworthy that estimated H-density values from the WACMM-X model are ∼ 2 to 3 times greater than those obtained by the NRLMSIS model during the quiet time period.

During the geomagnetic storm, there was an evident global enhancement of H densities beyond 3 R E starting on May 31, ∼ 1800 UT, which lasted for almost 1 day until June 1, ∼ 1600 UT, followed by a slow recovery to quiet-time conditions. The radial shell in the third row of Figure 2 (left side) shows a total increase of exospheric n H of ∼ 22 % with respect to quiet-time densities. The increase is significantly higher in the dayside nose (including the dawn sector) and the nightside tail compared with other regions of the exosphere. It is noteworthy that exospheric H densities derived from TWINS/LAD instruments during these days are greater than those described by the (Rairden et al., 1986) model during both quiet- and storm-time, having a constant value of 81 atoms/ cm − 3 . On the other hand, H density values in the exobase (right column) globally decrease during the storm with respect to quiet time conditions as a response to the increase in temperature.In this work, we use two exospheric models: the spherically symmetric and temporally static exosphere implemented by Rairden et al. (1986), hereafter referred to as the “static H” (StH) model, and the hybrid model derived from TWINS/LAD radiance data and WACCM-X predictions, referred to as the “dynamic H” (DynH) model. Figure 3 shows a numerical comparison between both models during the 1 June 2013 geomagnetic storm. Panel A shows the difference in H density distributions between both models expressed in a logarithmical scale for better visualization. Panel B shows the DynH/StH ratio as an alternative method to observe the variability of density distributions among the models.

Through the analysis of these plots, we found that the DynH model exhibits larger H densities than the StH model for all periods and at all altitudes (500 km–10 R E ). Although Figure 2 only shows the differences and ratios within the XY S M plane, we verified that H density is also greater at all latitudes. Panel A shows the density difference between models on a global scale. Since StH is fixed for all periods, an expansion of the DynH exosphere is evident (mostly in the dayside dawn sector) starting on 1 June 0500 UT, and lasting until 1 June 2300 UT. By the morning of June 2, the exospheric density distributions are similar to pre-storm conditions (e.g., see the contour line with level 1.75 at different time periods). Furthermore, panel B shows that within the ring current region in the equatorial plane (3–6 R E ), the DynH model provides densities ∼ 2 times greater than those from the StH model.

To determine the ring current development during the 1 June 2013 storm, we use the Comprehensive Inner Magnetosphere-Ionosphere (CIMI) model designed and implemented by Fok et al. (2014), Fok et al. (2021), which simulates the ring current ion fluxes and the magnetospheric electric field. CIMI is a kinetic model that calculates ion (0.1–500 keV) and electron (1 keV–6 MeV) distributions, the subauroral Region-2 field-aligned currents, the sub-auroral ionospheric potentials, and the plasmasphere distributions by solving (a) the bounce-averaged Boltzmann equation for the distribution functions of energetic ions and electrons, (b) the conservation equation of plasmasphere particles, and (c) the ionospheric current conservation equation for the ionospheric potential. The spatial domain of CIMI is circumscribed to the closed magnetic field region bounded by the dayside magnetopause located at a typical radial geocentric distance of 10 R E . The model grid resolution is defined by 53 and 48 grid bins in magnetic latitude and local time of the ionospheric foot-points of the magnetic field lines, respectively. Current species supported by CIMI include H + , O + , He + , and high-energy electrons. Also, CIMI accounts for several particle loss mechanisms such as (a) particle diffusion in energy and pitch angle due to interactions with plasma waves, (b) particle precipitation into the loss cone, and, of particular interest in this study, (c) charge exchange between ring current ions and exospheric neutral atoms. The charge exchange interaction is calculated in CIMI using the formula − v σ s H ⟨ n H ⟩ f s where v is the ion velocity, σ s H is the energy-dependent charge exchange cross-section between ion species “s” and atomic hydrogen, and the terms ⟨ n H ⟩ and f s are the averaged H density and averaged ion distribution function along a given field line, respectively. Specific settings in the CIMI model used in this study are provided in the next section. The reader is referred to Fok et al. (2014), Fok et al. (2021) for further details of the CIMI model.

To assess the role of a realistic terrestrial exosphere in the ring current decay, we use the CIMI model to simulate the ring current ion fluxes and their dynamic response to two distinct exospheric models (StH and DynH) during 3 days from May 31 to 2 June 2013. Then, the CIMI model is configured to calculate the fluxes of H + , O + , and high-energetic electrons ( e − ). We use the plasma sheet density and temperature from Tsyganenko et al. (2003) as the outer boundary condition, which is controlled by solar wind parameters. The ion composition of the plasma sheet is determined by the (Young et al., 1982; Pandya et al., 2018) formulations based on Kp, F10.7, and solar wind dynamic pressure (Psw). The ionospheric electric potential at the poleward boundary (71.4 deg) is specified by the Weimer 2K model (Weimer, 2001), and the electric field is calculated self-consistently within the CIMI domain (Fok et al., 2001). The magnetic field used is based on the TS04 model (Tsyganenko and Sitnov, 2005). Also, we only consider electron interactions with chorus waves, and besides the plasma loss to the magnetopause, charge exchange interactions with atomic H are the only ion loss mechanism included. Hence, we conducted two ring current simulations with these plasma parameters and magnetic field configuration, wherein only the exospheric H model was changed.

In this section, we calculated the temporal evolution of the total ring current energy (RCE) in units of [ 1 0 31 eV] during the 1 June 2013 storm. The left panel in Figure 10 shows the total RCE from both runs using DynH (solid red line) and StH (dashed red line), along with the Dst index (dotted blue line) and its pressure-corrected version Dst* (dashed blue line). In general, the time-dependent structure of the total RCE from both runs exhibits good agreement with the trends in Dst*. Also, while the simulated ring current in both runs has nearly equal peak strength, the DynH run exhibits a faster recovery (after 1 June 0900UT), which is associated with the larger H densities in the DynH model at ring current altitudes and consequent stronger charge exchange loss during this period. To further quantify the relationship between RCEs, the ratio RCE-DynH/RCE-StH is calculated and displayed in the right panel of Figure 10. Before the storm onset, the total RCE in the DynH run is ∼ < 5 % smaller than in the StH run. During the recovery phase, the total RCE in the DynH run constantly decays and becomes ∼ 32 % smaller than in the StH run on 2 June 2300UT. A simple analysis of the slope of the RCE-DynH/RCE-StH ratio during the recovery phase (after 1 June 0900UT) indicates that the “energy loss per day” in the ring current using the DynH model is ∼ 12.5 % faster than using the StH model.

Figure 11 compares simulated ring current proton fluxes (using DynH and StH) with in situ proton flux measurements acquired by the Helium, Oxygen, Proton, and Electron (HOPE) and the Radiation Belt Storm Probes Ion Composition Experiment (RBSPICE) mass spectrometers onboard NASA’s Van Allen Probes mission. Panels from top to bottom show (a) the proton (energy vs. time) flux spectrogram acquired by the HOPE + RBSPICE instruments during the period under study (May 31 to June 2) and covering ion energies from 1 × 10 2 to 5 × 10 5 eV, (b) the simulated proton flux spectrogram from the CIMI + StH run for the actual location of the Van Allen Probes A spacecraft, (c) the simulated proton flux for CIMI + DynH model for identical locations, (d) the ratio Log 10 ( Q ( E , t ) ) , where Q = j StH ( E , t ) / j vap ( E , t ) , and E represents the energy channel, (e) the ratio Log 10 ( Q ( E , r ) ), where Q = j DynH ( E , r ) / j vap ( E , r ) , (f) the L shell location of the Van Allen Probes A, and (g) the Sym-H index for the 1 June 2013 storm. With an identical format, Figure 12 shows oxygen ion spectrograms from CIMI simulations and in situ measurements from Van Allen Probes.

The analysis of CIMI + StH/RBSP (Figure 11d) and CIMI + DynH/RBSP (Figure 11e) ratios indicate that the simulated proton fluxes with both exospheric models have a reasonably good agreement with in situ observations (i.e., | Log 10 ( Q ) | < 1) when RBSP-A is located at L ∈ [ 4,6.5 ] R E . On the other hand, when the spacecraft is located at L ∼ < 2 R E , flux simulations with StH and DynH models exhibit a large ( | Log 10 ( Q ) | > 3 ) to moderate ( | Log 10 ( Q ) | ∈ [ 1,3 ] ) over- and under-estimation at all energies. These features are similar in orbits 3, 4, 5, and 6 during the recovery phase of the storm and are highlighted with black line rectangles. The overestimation of proton fluxes at low L near spacecraft perigee seems to be related to an undershielding effect in the CIMI model, which yields strong electric fields at low L pushing ions, especially protons, deeper than observed. A simple comparison of the fluxes within the black line rectangles (panels e and d) indicates that using DynH slightly reduces the overestimation owing to its high H-density at ring current altitudes, e.g., on 02 June 0300 UT (between orbits 4 and 5), the ratio for proton fluxes, whose energy channel is 2 × 1 0 4 eV, changes from Log 10 (CIMI + StH/RBSP) = + 2.0 to Log 10 (CIMI + DynH/RBSP) = + 1.5 . A similar trend is found for O + ion fluxes, as shown in Figure 12. In this case, there is a lower overestimation of fluxes near the spacecraft perigee than those observed for protons as the undershielding effect in the CIMI model (mentioned above) is less efficient for heavy ions. Similar to the proton case, the use of DynH over the O + fluxes slightly reduces the overestimation near perigee, e.g., on 02 June 0300 UT (between orbits 4 and 5), the ratio for O + fluxes, whose energy channel is 2 × 1 0 4 eV, changes from Log 10 (CIMI + StH/RBSP) = + 2.5 to Log 10 (CIMI + DynH/RBSP) = + 2 .

Figure 13 shows the integrated ion energy flux (for both species) over the energy range of 0.1–500 keV along the path of RBSP-A during the 1 June 2013 storm. For comparison, we also calculated the integrated ion flux for CIMI + StH and CIMI + DynH simulations.

The integrated proton energy flux obtained from CIMI simulations is significantly larger than those measured by RBSP-A at L < 5 R E . This trend can be attributed to (1) the simulation settings, which only include limited ion loss processes since our goal in this study is to identify the impact of the exosphere in the ring current, and (2) an undershielding effect in the CIMI model, which allows ions deeper access than observed. With this in mind, we observe that using the DynH model in the simulations certainly reduces the integrated flux up to 25% compared to the results yielded by using StH (e.g., in the top panel on 01 June 1900 UT). The results for the O + integrated fluxes exhibit a better agreement with those obtained from CIMI simulations except near the storm’s peak. Similarly, the use of DynH in the ring current yields a reduction of integrated energy flux of up to ∼ 25 % with respect to the values obtained from simulations with StH model (e.g., in the middle panel on 02 June 1400 UT).