Modeling Field Line Curvature Scattering Loss of 1–10 MeV Protons During Geomagnetic Storms

The proton radiation belt contains high fluxes of adiabatically trapped protons varying in energy from ∼one to hundreds of megaelectron volts (MeV). At large radial distances, magnetospheric field lines become stretched on the nightside of Earth and exhibit a small radius of curvature RC near the equator. This leads protons to undergo field line curvature (FLC) scattering, whereby changes to the first adiabatic invariant accumulate as field strength becomes nonuniform across a gyroorbit. The outer boundary of the proton belt at a given energy corresponds to the range of magnetic L shell over which this transition to nonadiabatic motion takes place, and is sensitive to the occurrence of geomagnetic storms. In this work, we first find expressions for nightside equatorial RC and field strength Be as functions of Dst and L* to fit the TS04 field model. We then apply the Tu et al. (2014, https://doi.org/10.1002/2014ja019864) condition for nonadiabatic onset to solve the outer boundary L*, and refine our expression for RC to achieve agreement with Van Allen Probes observations of 1–50 MeV proton flux over the 2014–2018 era. Finally, we implement this nonadiabatic onset condition into the British Antarctic Survey proton belt model (BAS‐PRO) to solve the temporal evolution of proton fluxes at L ≤ 4. Compared with observations, BAS‐PRO reproduces storm losses due to FLC scattering, but there is a discrepancy in mid‐2017 that suggests a ∼5 MeV proton source not accounted for. Our work sheds light on outer zone proton belt variability at 1–10 MeV and demonstrates a useful tool for real‐time forecasting.


Introduction
The proton radiation belt is characterized by high fluxes of protons varying in energy from ∼one to hundreds of megaelectron volts (MeV).Particles in this region undergo periodic gyration, bounce and drift motion through the geomagnetic field and conserve an adiabatic invariant over each timescale.On the nightside of Earth, where magnetic local time (MLT) is near zero, the field becomes stretched as radial distance increases.Eventually, field line radius of curvature becomes comparable to the proton gyroradius near the magnetic equator for a given value of the first invariant μ.This leads to field line curvature (FLC) scattering, whereby nonadiabatic changes δμ accumulate over each traversal through the region (Young et al., 2002).The nominal outer boundary of the proton belt at a given energy corresponds to the range of magnetic L shell over which this transition to nonadiabatic motion takes place in the equilibrium magnetic field, that is, L ∼ 3 at ∼20 MeV (Selesnick & Looper, 2023) and 7 ≲ L ≲ 9 at hundreds of keV (Sergeev & Tsyganenko, 1982).However, stretching of the field during geomagnetic storms can induce scattering at lower L, leading to inward re-positioning of the boundary (Hudson et al., 1997).Engel et al. (2016) studied the loss of radiation belt protons during a storm in March 2015 by simulating sample particle trajectories and succeeded in reproducing the observed post-storm flux distribution.Their experiment showed that storm-time electric fields associated with strengthening of the ring current play a crucial part in such loss events by causing the outward adiabatic expansion of drift orbits into regions of higher curvature to conserve the third adiabatic invariant.Models based on a radial diffusion framework require an empirical approach in order to include such storm-time losses, as well as the onset of nonadiabatic motion controlling the quiet-time boundary.
The ratio ϵ = ρ/R C has been used in various literature to parameterize conservation of μ, where ρ is the particle gyroradius and R C the field line radius of curvature, with ϵ ≪ 1 conducive to adiabatic motion and ϵ ∼ 1 indicating the onset of scattering (i.e., Alfven & Falthammar, 1963;Büchner & Zelenyi, 1989;West et al., 1978).The normalized change δμ/μ over each traversal through the minimum-B point has been shown to vary exponentially with 1/ϵ and sinusoidally with gyrophase ψ (Birmingham, 1984;Gray & Lee, 1982;Hastie et al., 1969), and with equatorial pitch angle α eq as shown in Figure 4 of Young et al. (2002).Young et al. (2002) showed that a model of δμ/μ should also depend on additional parameterizations of field line curvature and in order to generalize well to different regions of the magnetosphere, where S max is the maximum-ϵ point on the field line.Young et al. (2008) then developed an expression for FLC scattering-induced pitch angle diffusion coefficients D α eq α eq ϵ,α eq ,ζ 1 ,ζ 2 ) .Recently, Selesnick and Looper (2023) implemented these diffusion coefficients into a proton belt model, using the T89 external field model (N.Tsyganenko, 1989) to evaluate them through time (see also Yu et al., 2020;Zhu et al., 2021).Their work demonstrates variation in trapped proton intensity as a result of FLC scattering that compares well with observations, and demonstrates the challenge of capturing the extent of scattering observed at low L during storms.Scattering is highly sensitive to field geometry, and the authors suggest that results might be improved using a more accurate external field model.
An alternative approach to model FLC scattering has been to consider a stable trapping condition.For example, Anderson et al. (1997) show that the δ B (α eq , ϵ) parameter of Birmingham (1984) corresponds well to the onset of nonadiabatic motion at the threshold δ B ∼ 0.01 in a simulation of ring current particles.A more recent condition for stable trapping was derived by Tu et al. (2014) using particle tracing experiments at E ≥ 25 MeV in the T89c field model: ϵ max < ϵ onset (E, α eq , Kp), where ϵ max is maximum ϵ along a drift shell, and the Kp dependence takes into account variability of the T89c field.Such conditions are simple to evaluate, however, they approximate the outer boundary as a sudden transition to nonadiabatic motion.
There is currently a call from industry to extend the capabilities of radial diffusion models to include space weather forecasting (Horne et al., 2021).For protons, this requires evaluating losses due to FLC scattering using indices available in real-time.A key energy range for space weather applications is 1-10 MeV because it is the primary contributor to solar cell degradation of satellites transiting the proton belt (Lozinski et al., 2019;Messenger et al., 1997).However, this poses another challenge due to a lack of spacecraft measurements at this energy near the magnetic equator, which are required to capture the full pitch angle distribution.
In this work, we develop a method to evaluate FLC scattering loss of radiation belt protons in real-time.Our approach involves using the TS04 external magnetic field model (N. A. Tsyganenko & Sitnov, 2005) to derive expressions for nightside R C , B e as continuous functions of Dst and L*, thereby including the effect of a strengthening ring current during geomagnetic storms.We use these expressions to solve for ϵ max and determine the region subject to FLC scattering by using the Tu et al. (2014) onset condition.This region is also determined directly from Van Allen Probes observations of 1-50 MeV proton flux during several large storms over the 2014-2018 era, and we adjust our expression for R C to produce good agreement between observation and calculation.
We then use the British Antarctic Survey proton belt model (BAS-PRO; Lozinski et al., 2021) to solve the time evolution of proton flux in the region L ≤ 4 and compare with observations.The updated BAS-PRO model accurately captures the inward extent of losses following large geomagnetic storms, but there is a discrepancy near the end of the modeling period (mid-2017) that suggests a source of ∼5 MeV particles near the outer boundary not accounted for by the model.With the additional capability to model FLC scattering loss in real-time, we shed light on outer zone proton belt variability at 1-10 MeV and demonstrate a useful tool for space weather forecasting.

Overview
The Van Allen Probes pair of satellites (formerly known as the Radiation Belt Storm Probes, RBSP) were launched into elliptical orbit (∼600 km pergiee to ∼5.8R E apogee) at 10°inclination on 30 August 2012 (Kessel et al., 2013).This work makes use of proton flux data recorded by the following onboard instruments (in order of increasing energy): the Radiation Belt Storm Probes Ion Composition Experiment (RBSPICE, Manweiler et al., 2022;Mitchell et al., 2013); the Magnetic Electron Ion Spectrometer (MagEIS, Blake et al., 2013); and the Relativistic Electron-Proton Telescope (REPT, Baker et al., 2012) of the RBSP-Energetic Particle, Composition, and Thermal plasma (RBSP-ECT) instrument suite (Spence et al., 2013).
We previously made use of the RBSPICE and MagEIS data in Lozinski et al. (2021), and found evidence of contamination from electrons and higher energy protons that limited its use to L ≥ 2. Consequently, we only use the data in this region.We refer to Lozinski et al. (2021) to supplement the discussion in Sections 2.2 and 2.3 with more details about the instruments and caveats.The REPT data set is also subject to interference as discussed in Section 2.4, but is used throughout the region 1.2 ≤ L ≤ 4.
The modeling period of 2014-2018 was selected based on overall data availability.The top panel of Figure 1 shows several minima in Dst index corresponding to geomagnetic storms during this period.Prior to data processing, we applied a minima-finding algorithm to the time series of Dst index to identify time windows encompassing each minimum below 65 nT.These were identified as storm periods, generally lasting ∼1-3 days, during which large changes in flux tended to occur.A 120-hr time window was also taken before and after each storm period, resulting in a set of three adjacent data average windows for each storm: the pre-storm, storm and post-storm windows respectively.The remaining time periods, spanning from the end of each post-storm window to the beginning of the next pre-storm window, were then divided into data average periods equal in length and as close to 120 hr as possible.Equatorial pitch angle distributions were fit to the flux data averaged within each window.This method highlighted pre and post-storm distributions of flux, and reduced the standard deviation of each fit by not averaging across storm periods.
Figure 1 (second panel and below) shows fitted differential flux at 90°equatorial pitch angle for selected energies across the range of all three instruments in the outer proton belt at 2 ≤ L* ≤ 4. Figure 1 shows that the occurrence of large geomagnetic storms corresponding to minima in Dst, such as the 223 nT storm in March 2015, coincides with depletions in flux levels near the proton belt outer boundary.The recovery period following each depletion is shown to be energy-dependent, with lower energies (i.e., 2.97 MeV) recovering in a matter of months between storms, whereas the March 2015 storm causes a depletion in the 46.7 MeV channel that seems to last at least 1 year.

RBSPICE Measurements
The RBSPICE instrument measured ions from ∼20 keV to several MeV.The Ion Species High Energy Resolution Low Time Resolution (ISRHELT) data product from this instrument was used to provide data over the modeling period.ISRHELT data is contained within the level 3 Common Data Format (CDF) files obtainable online at http://rbspiceb.ftecs.com/Level_3/ISR-HELT/(Lanzerotti et al., 2024).ISRHELT CDF files provide proton differential unidirectional flux as a 3D array of values, with dimensions epoch, energy channel and telescope.There were six telescopes recording simultaneously, and each measurement of flux was taken in the instantaneous look direction of the corresponding telescope, rotating with the spacecraft.
Measurements taken within a given data average window were first resolved into bins of equatorial pitch angle using the same method as described in Section 2.1.1 of Lozinski et al. (2021), and a cadence was applied to average the flux measurements in each bin across 1 min intervals.For a given L* under consideration, data outside L* ± 0.02 were then filtered out based on the L* coordinate of the spacecraft at each epoch provided by the CDF files.As per the method in Lozinski et al. (2021), measurements from five of the six telescopes were combined but the remaining telescope data excluded.An equatorial pitch angle distribution was then fit to the resulting data points spanning the data average window for each energy channel, using the function where j is differential unidirectional flux at equatorial pitch angle α eq , and A and n are fitting parameters.In order to avoid contamination issues, usage of the ISRHELT data was limited to the energy channels spanning from 0.29 to 1.68 MeV.

MagEIS Measurements
The MagEIS "high" unit electron spectrometer on Van Allen Probe B housed an ion range telescope with three silicon detectors, and the 2500-micron detector in this arrangement was used to measure ∼2-20 MeV protons (Blake et al., 2013).This data is accessible via the "FPDU_pix2" variable in the Level 3 RBSP-B CDF files available online at the RBSP-ECT Science & Data Portal at https://rbsp-ect.newmexicoconsortium.org/data_pub/rbspb/ (Spence, Reeves, Blake, et al., 2024).The CDF files provide a 3D array of differential unidirectional flux with dimensions epoch, local pitch angle and energy channel, with local pitch angle in terms of 15 equally spaced bins spanning 0°-180°.
Within a given data average window, the local pitch angle distribution provided by the CDF files was used to derive a time series of equatorial pitch angle distributions at each CDF epoch using the same method as described in Lozinski et al. (2021), Section 2.2.1.For a given L*, data outside L* ± 0.02 were filtered out using the spacecraft L* from the CDF files, and Equation 1 was then used to fit all observations spanning the data average window.In order to avoid contamination issues, usage of the MagEIS pix2 data was limited to the six energy channels covering 2.97-12.4MeV.
An equatorial pitch angle distribution was fit for each data average window, L* bin and energy channel, using the function where A, n 1 , and n 2 are fitting parameters and j 90 is flux at 90°equatorial pitch angle.Equation 2 was chosen because in the energy range of REPT, cosmic ray albedo neutron decay is a significant source of off-equatorial protons which causes pitch angle distributions to exhibit wide flanks before falling to zero near the loss cone, which Equation 1 cannot fit well.
At 21:54 UT on 21 May 2015, the instrument equation used to determine proton counts was changed in order to reduce background noise from galactic cosmic rays (GCRs), as described in Section 4.2 of Baker et al. (2021).Following this event, a step change occurs in raw flux measurements affecting various energy channels.
In order to use REPT proton data prior to the mode change, we applied a correction via the following steps: (a) compare the fitted flux measurements immediately before and after the change; (b) find a constant factor for each energy channel, equatorial pitch angle and L* bin that describes the drop in flux; (c) multiply measurements made before 21:54 UT, 21 May 2015 by the correction factor, then re-fit the pitch angle distributions.This method eliminated the step change, however it required making the assumption that the GCR background was constant.Furthermore, the correction factors indicated that GCR background was comparable in magnitude to trapped flux.Therefore, we use the five energy channels covering 27.6-78.9MeV to initialize our simulations, but otherwise for qualitative comparisons only.

Functional Forms for R C and B e
The exponential dependence of δμ/μ on 1/ϵ during a particle's traversal of the minimum-B point shows how FLC scattering timescales are sensitive to R C and B e encountered along a drift orbit.These quantities, along with other parameterizations of curvature such as ζ 1 and ζ 2 , can be determined using a magnetic field model then used to evaluate empirical scattering models (i.e., Selesnick & Looper, 2023;Yu et al., 2020).The Dst dependence of the TS04 model suggests it would be a good choice for this, as Dst has been shown to correlate with the inward L extent of storm-driven loss in the proton belt (Selesnick et al., 2010).However, TS04 also depends on other parameterizations of geomagnetic and solar wind conditions that make it more complex to evaluate, especially in real time.We have therefore taken the empirical approach to FLC scattering further by constructing a model for nightside radius of curvature R C (Dst, L*) and equatorial field strength B e (Dst, L*) based on average variation of the TS04 field with Dst.These expressions can be used to evaluate ϵ max , the maximum ϵ encountered by a particle during its drift orbit, and in turn evaluate empirical models of FLC scattering onset or timescales.The dependence on Dst alone cuts out the need to evaluate an external field model, and allows the study of any time period overlapping with historical records of Dst (i.e., Nose et al., 2015).
To begin, a set of Dst bins were defined, with central values as shown in the left column of Table 1.A minimafinding algorithm was then applied to a time series of Dst spanning from 1983 to mid-2023, provided by NASA omni data and accessed via the Python spacepy library (Morley et al., 2011).Events corresponding to a minimum in Dst separated by at least 0.5 days from other minima were identified and collected in the corresponding bin.The Dst = 50 nT bin was skewed as shown in row 2, column 2 of Table 1 so that the average Dst of events in this bin was close to 50 nT.For the Dst ∼ 0 bin, events were collected when Dst crossed from a negative to positive value.
The remaining TS04 driving parameters P dyn , B yIMF , B zIMF , W 1 , W 2 , W 3 , W 4 , W 5 , and W 6 (see N. A. Tsyganenko & Sitnov, 2005) were then sampled from NASA omni data at all event epochs and the average of each parameter was  1.
taken across the events within each Dst bin. Figure 2 plots the average value of each TS04 driving parameter (blue dots) against the Dst bin centers, along with the standard deviation of each average and the original values (light gray dots).Figure 2 shows an approximately linear correlation between Dst and most TS04 driving parameters.This suggests that the dependence of TS04 can be reduced to Dst alone with the caveat that doing so will not capture the whole variety of magnetic field geometry possible.
For each Dst bin, the set of average driving parameters shown in Figure 2 was then used to drive TS04, with the goal of producing an average field configuration for that Dst.TS04 was evaluated using the IRBEM package as implemented in spacepy.For each resultant magnetic field configuration, R C and B e were determined along the nightside magnetic equator as a function of L* for a particle with 90°equatorial pitch angle.The choice of L* as a coordinate, as opposed to L, was necessary to take into account the Dst effect, since L is not conserved during the outward adiabatic expansion of drift orbits accompanying storm-time electric fields.The results of this calculation are shown in Figure 3 (dots).To perform the calculation, the IRBEM package was first used to locate the magnetic equator at MLT = 0 for each L* queried, giving a location where the equatorial plane was intersected by a field line of the given drift shell.Field strength B e was taken at this location, and a radius of curvature was optimized to fit a circle to the locus of points given by tracing a short distance up and down the field line in small steps away from the equator.This process was repeated for different L* to build up each profile as shown.
The radius of curvature data shown in Figure 3 (top panel, dots) was fit using the function R C0 = 0.47723L * 0.128 exp (0.4309L * ) in units of Earth Radii where 1R E = 6.3712 × 10 6 m, with constants c 1 to c 5 listed in Table 2 and calibration parameter X equal to zero by default, but utilized later to adjust this fit.The expression for R C0 represents a fit to the radius of curvature when Dst = 0. Figure 3 shows Equation 3 evaluated for each Dst using the coefficients in Table 2 (solid line), as well as Equation 3 evaluated with a different set of coefficients for each Dst (not listed) to achieve a fit as close as possible to each set of original data points (dashed lines).The latter set of curves demonstrates the flexibility of the functional form in Equation 3 and highlights the difference between a general fit versus a fit for a specific Dst. Figure 3 shows R C → 0 at high L*, which corresponds to an unrealistically thin current sheet.Usage of Equation 3 should therefore be limited to the domain where R C is above some minimum value at high L*, but this limit could not be constrained using TS04.
The equatorial field strength data shown in Figure 3 (bottom panel, dots) was fit using the function where B 0 is mean field strength along the equator, taken as 3 × 10 5 T at present day, and constants d 1 -d 6 as listed in Table 2.As with R C , a general fit using the coefficients in Table 2 is shown (solid line), as well as a fit to each Dst using several different sets of coefficients not listed (dashed line).

Fitting to Observations of Loss Onset
It can be determined whether or not a proton drift orbit is adiabatic by evaluating the condition for nonadiabatic onset ϵ max ≥ ϵ onset .Here, ϵ max is given by evaluating ϵ = p sin α eq /qB e R C using Equations 3 and 4, and ϵ onset is given by  However, we found that this approach tended to overestimate the apparent maximum trapped L* when compared with the observations shown in Figure 1.One potential reason could be that R C is overestimated by TS04 during times of minimum Dst, so that using the fit of Equation 3 leads to underestimates of ϵ max .Another source of error may be attributed to our usage of Equation 5 from Tu et al. (2014), which was developed from particle tracing experiments in the T89c field at E ≥ 25 MeV.The good fit demonstrated by the original authors may therefore be contingent on applying this condition with T89c field geometry, whereas our expressions for R C and B e are based on TS04, and the energy dependence may be incorrect when extended to 1-10 MeV.However, we found that a good fit to observations of the apparent trapping boundary could be achieved by applying a modification to Equation 3 to reduce R C .
To derive this modification, we first selected six storm events corresponding to the six largest minima in Dst, and analyzed their pre and post-storm flux profiles to infer the change in maximum trapped L* that occurred during each storm's peak.Figure 4 (left panels) shows the pre and post-storm flux profiles (solid and dashed lines respectively) for two example storms at selected energies (rows).We denote the L* at which each set of profiles diverges as L max , and detect it algorithmically as the lowest L* where log 10 j 1 /log 10 j 0 ≥ 0.955, where j 0 and j 1 are pre and post-storm flux respectively.Using this criteria, we were able to infer L max at multiple energies and times, with examples shown in Figure 4 by the vertical red lines (left panel) and red crosses (right panel).
This method for inferring L max was limited to observations of large storms which provided the clearest evidence of changes to the trapping boundary.During small (∼ 50 nT) storms, L max seems to correspond to regions where measurements of flux contain high levels of uncertainty, especially at ∼10 MeV, making it difficult to infer precisely.Our criteria was also unable to identify L max during certain storms in the low energy (∼1 MeV) channels, as indicated by the missing red crosses in Figure 4 (right, top panel), because the change in post-storm flux was too subtle.
Finally, we optimized a calibration parameter X in Equation 3 so that the L* of nonadiabatic onset, where ϵ max = ϵ onset , was calculated to match the set of L max inferred from observations (red crosses in Figure 4, right panels).The optimized calculation for L* of nonadiabatic onset is plotted in Figure 4, right panels, as a white line.Figure 4 demonstrates that this calculation seems to match the inward extent of losses observed following minima in Dst.The effect of the calibration parameter X is to reduce R C as shown in Figure 5, which compares R C calculated using X = 0 and the optimized value of X = 4.7469 × 10 3 .This modification to R C is in line with that demonstrated in Figure 11 of Selesnick et al. (2010), which was required to better model the inward extent of losses observed by SAMPEX following a 288 nT storm on 5th April 2000.As those authors state, reducing R C but leaving B e unchanged requires assuming a thinner current sheet but same total ring current.This calibration decreases the L* at which R C from Equation 3 becomes unrealistically small.However, the critical value of R C to determine adiabaticity occurs where ϵ max = ϵ onset .Considering a 7.5 MeV proton with ϵ onset (Kp = 6, α eq = 90°, E = 7.5 MeV) ∼ 0.34 in a Dst = 200 nT storm, ϵ max = ϵ onset occurs at L* ∼ 2.8, where Figure 5 (purple dotted line) shows R C ∼ 0.1.For the same particle in a 350 nT storm, ϵ max = ϵ onset occurs at L* ∼ 2.5, where Figure 5 (gray dotted line) shows R C ∼ 0.1 also.Therefore in practice, the values of R C near zero are not utilized by this method.In the next section, we apply our expressions for B e and modified R C to model FLC scattering and simulate the observation period shown in Figure 1.

Modeling
The BAS-PRO model, detailed in Lozinski et al. (2021), solves a distribution function proportional to proton phase space density f(μ, K, L) ∝ j/p 2 as a function of the three adiabatic invariants μ, K, and L, where p is momentum and j is differential, unidirectional flux.The full model equation is where J = 2pK/ ̅̅̅̅̅̅ B m √ in terms of mirror field strength B m , but approximated here using J ≈ 2pLaY(y) in terms of the sine of equatorial pitch angle y, radius of Earth a, and Y given by Equation 161in Schulz (1991).Equation 6includes the effects of: Coulomb collisions with various constituents of the atmosphere, ionosphere and plasmasphere (second and third terms on the left side), modeled as frictional change dμ/dt fric ; scattering of the third invariant L due to electromagnetic perturbations modeled as radial diffusion (first term on the right side); the cosmic ray albedo neutron decay source S n ; loss due to inelastic nuclear scattering with inverse timescale Λ; and FLC scattering loss (last term on the right side).The geomagnetic field is modeled as a dipole, with B 0 = 2.9867 × 10 5 T (calculated for the year 2015 using coefficients of the IGRF magnetic field model, Alken et al., 2021).We used the following expression for D LL because it produced optimal results, but we note that there is much uncertainty surrounding this parameter (see for example Figure 5 of Holzworth & Mozer, 1979): There are two differences in the new version of model compared with the version presented in Lozinski et al. (2021): one is the addition of FLC scattering loss; the other is in modeling of loss cone equatorial pitch angle α ecc.
h (L, h), corresponding to K = K lc , which is now based on a fit to particle tracing calculations of protons in an offset dipole field at various atmospheric scale heights h as opposed to a centered dipole.A loss cone based on an offset dipole field leads to significantly more anisotropic pitch angle distributions at L < 2 compared with a centered dipole, and the new loss cone equatorial pitch angle model is detailed in Appendix A.
At each model timestep, the Tu et al. ( 2014) condition for nonadiabatic onset ϵ max ≥ ϵ onset (E, α eq , Kp) was evaluated across the model grid using Equation 3for R C (Dst, L*), modified with X = 4.7469 × 10 3 , and Equation 4for B e (Dst, L*).Since BAS-PRO uses a dipole field, the model L coordinate was used in place of L*; the effect of FLC scattering on the model distribution at a given L is thereby approximated by assuming the particles forming this distribution have L* ∼ L when on the nightside, where Equations 3 and 4 apply.Maximum adiabatic L was then interpolated and any higher L was considered subject to FLC scattering.However, observations of post-storm 1-10 MeV flux in Figure 1 and the example profiles in Figure 4 (left panels) show that measured flux was not necessarily zero in scattering regions.This implies either that loss following nonadiabatic onset is not instantaneous, or that trapped flux is restored quickly following a storm.We tested the latter hypothesis by setting flux equal to zero throughout all nonadiabatic regions, and found that post-storm flux did not recover quickly enough to match observations even with enhanced rates of radial diffusion.Therefore, we did not assume instantaneous loss due to FLC scattering, but rather apply a timescale for exponential decay τ flc like so: where τ d is the drift period.Our simple approximation of τ flc = 100τ d in the nonadiabatic region was made because 1 MeV (10 MeV) protons have a drift period of ∼800s (80s) at L ∼ 3, therefore 100τ d corresponds to ∼1 day (2 hr).This gives a qualitative agreement with the decrease in post-storm ∼1 MeV flux at L = 3 following an average storm period spanning 1-3 days, yet at higher energies leads to a faster decay in flux as observed.In reality, we expect that the loss timescale at a given adiabatic coordinate changes with time as the field distorts over multiple drift periods during buildup of the ring current.Figure 4 (right panels) together with Equation 8shows how deeply our model predicts FLC scattering loss to occur through time.
Similar to our previous study, four boundary conditions were applied to the model: (a) f(μ, K, L i ) = 0 where where L o = 4 and f b (μ, K, t) is the time-varying outer boundary spectrum specified from the Van Allen Probe measurements and extrapolated across the range in μ; (c)f(μ, K lc , L) = 0 where K lc (L) is the lowest K on the model grid that is inside the loss cone at L; and (d) f(μ max , K, L) = 0 where μ max is chosen such that it corresponds to an energy of 80 MeV at L max , which is observed to be higher than the adiabatic trapping limit during quiet times at all L. The initial distribution of phase space density at the simulation start time is solved by fixing the model to phase space density interpolated from Van Allen Probes observations everywhere available, then allowing the model to reach steady state to fill gaps in the spectrum between these fixed observations.

Results
We present simulated 90°differential proton flux from 2014 to February 2018 in Figure 6 for comparison with the observations shown in Figure 1. Figure 7 shows the corresponding ratio of modeled to observed flux where data was available, with areas of red representing an overestimate and blue representing an underestimate.
Model accuracy can be examined in more detail using Figure 8.It shows modeled differential flux at 90°equatorial pitch angle before (blue) and after (red) selected storms (solid lines) compared with Van Allen Probe observations (dots) at selected energy channels (one per row).Figure 8 shows a different storm in each column, with minimum Dst during the storm shown at the top of each column alongside the dates of the pre and post-storm flux profiles.The standard deviation is shown for each observation (vertical lines), and the observation is taken from the fitted equatorial pitch angle distributions as described in Section 2. Each observation plotted is therefore also the result of time averaging spacecraft measurements inside each pre and post-storm window before this fit is applied.Figure 9 shows the same results for 60°e quatorial pitch angle, which has high levels of uncertainty.

Discussion
Figure 6 shows losses occurring due to FLC scattering during several large geomagnetic storms throughout the modeling period.The recovery of flux levels following storm-driven loss is due to a combination of outward radial diffusion from lower L unaffected by the storm, and inward radial diffusion from the data-driven outer boundary.To investigate this, Figure 10 (top panel) plots the distribution function f at a fixed value of the first invariant μ = 500 MeV/G, for two values of second invariant K before and after the 223 nT geomagnetic storm in March 2015.This panel shows that prior to the storm (blue line), the gradient in phase space density is relatively flat at L ≳ 2.5, which is typical during quiet times throughout the modeling period.Shortly after the storm (yellow line), there is a minimum in phase space density at L ∼ 3.5, which is gradually filled in over several weeks as shown by the subsequent profiles (orange lines) spaced 8 days apart.Figure 10 (bottom panel) plots the corresponding timescales for FLC scattering loss and radial diffusion at the time of minimum Dst, with τ flc given by Equation 8. Figure 10 thus shows the source of particles from inward radial diffusion following the storm is greater than from outward radial diffusion due to the sharp dependence of D LL on L. The bottom panel also shows that scattering occurs at lower L for particles with larger K/smaller α eq .Figure 6 shows an increase in flux at (second and fourth panels respectively).One possible reason could be that the rate of radial diffusion was overestimated, however the better agreement with nearby surrounding energy channels suggests some energy channels could also be mislabeled.
An area for improvement is the model's underestimation of flux observed by energy channels spanning from 4.92 to 12.4 MeV, highlighted by the blue area near the end of the modeling period in Figure 7 (fifth panel).Figure 1 does show observations of loss corresponding to the geomagnetic storms of minimum Dst ∼ 125 nT and Dst ∼ 142 nT in May 2017 and September 2017 respectively, which are shown in more detail by columns four and five of Figures 8 and 9.However, flux in both the pre and post-storm distributions for these two storms were underestimated by the simulation at ∼5-10 MeV.Some possible explanations are that: the rate of radial diffusion increases near the end of the modeling period rather than remaining constant; there is a source of protons below the outer boundary unaccounted for by the model; or that a Dst dependence alone is not enough to incorporate nondipolar effects on scattering.
One limitation of the method demonstrated is the approximate timescale τ flc = 100τ d used to model scattering loss in the nonadiabatic region.Examining Figures 8 and 9, it can be seen that this timescale leads to an underestimation of post-storm flux at L ∼ 3 and above at 1.52 and 4.92 MeV (first and third rows).At higher energies however, this approximation is less important because particles scatter faster in nonadiabatic regions and the outer boundary tends to follow that predicted by the onset condition more closely.To improve this method, one could develop functional forms for the ζ 1 and ζ 2 parameters as done here for R C and B e .This would enable the evaluation of the Young et al. (2008) pitch angle diffusion coefficients to use as inverse scattering timescales, and we leave this as future work.

Conclusion
In this work, we developed a method to model field line curvature scattering loss of the proton radiation belt that can be utilized in real time without needing to evaluate an external magnetic field model.Functional forms for radius of curvature R C and equatorial field strength B e at magnetic midnight were developed based on geometry of the TS04 field, then used to empirically determine ϵ max as a function of L* and Dst.The onset condition of Tu et al. (2014), derived from particle tracing experiments, was then applied to calculate variations in L max of the proton belt outer boundary during geomagnetic storms.Calculations achieved good agreement with observational evidence of changes to the range of nonadiabatic L*, after calibrating the fit to R C via a constant parameter X.
Finally, we implemented the onset condition along with a simple loss timescale into the BAS-PRO proton belt This was achieved by computing the trajectories of ∼1 MeV test protons in an offset dipole field.For a given L shell, protons were initialized at the magnetic equator with equatorial pitch angles 1-89°in increments Δα eq = 0.5°.Each particle was initialized at the same magnetic longitude where the bounce loss cone was most restricted and therefore equal to the drift loss cone, which is in the direction of the vector pointing from the origin to the offset center of Earth.Earth's surface was modeled using the WGS84 ellipsoid (Decker, 1986).For each L, the largest pitch angle that collided with Earth's surface was recorded and taken as the numerical loss cone boundary, with the results shown in Figure A1 (red crosses).Furthermore, the largest pitch angles that came within altitudes of 150, 300, 450, 600, and 750 km were recorded, with the results also shown in Figure A1 (crosses).
The numerically computed loss cone as a function of L representing collision with Earth's surface was then fit using the following expression, plotted in Figure A1 (red line): α ecc.0 (L) = α dip.0 (L) + 9.7 × 10 2 x 8.12 × 10 4 x 2 + 3.83 × 10 6 x 3 (A2) where x = e 5.54/L .Using this expression as a base fit, the loss cone profiles for different altitudes could then be fit using the following expression, plotted in Figure A1 for each altitude h: α ecc.h (L,h) = α ecc.0 (L) + 10 0.0653h 5.315 z + 10 0.402h 11.7 z 2 + 10 0.725h 18 z 3 (A3) where z = h/(L 1).The centered dipole loss cone given by Equation A1 is also plotted in Figure A1 (black line), showing the significant difference particularly at low L. The numerically computed loss cone model given by Equation A3 advanced our modeling by resulting in a more realistic loss cone, but also by allowing small changes to the atmospheric scale height h controlling the loss cone boundary.Increasing h resulted in better numerical stability, allowing fine control to be exerted when choosing an optimal timestep.

Figure 1 .
Figure 1.Van Allen Probe B observations of 90°differential flux at selected energy channels from the RBSPICE instrument's ISRHELT data product (1.52 MeV), MagEIS (2.97-12.4MeV) and REPT (27.6 and 46.7 MeV), over the 2014-2018 period, alongside Dst index (top panel).The color scale has been set to a minimum flux level of 0.1 cm 2 s 1 MeV 1 sr 1 .Periods of data unavailability are colored black.

Figure 2 .
Figure 2. TS04 driving parameters sampled during minimum Dst events (gray) then averaged (blue) within bins of event Dst, plotted against the Dst bin center.Standard deviations of each average are shown in black and the number of events per Dst bin is shown in Table1.

Figure 3 .
Figure 3. Local radius of curvature R C and equatorial field strength B e given by the TS04 magnetic field model as a function of L* for a particle with 90°e quatorial pitch angle crossing magnetic midnight (dots), along with fits of Equations 3 and 4 to individual Dst (dashed lines) and a single fit for all Dst (solid lines).

Figure 4 .
Figure 4.The highest adiabatic L max during minimum Dst inferred from observations across six geomagnetic storms.Left panels show examples from two storms; the observed pre and post-storm 90°differential flux is plotted against L* for a selection of energy channels, with L max represented by vertical red bars.Right panels show L max (red crosses) for all six storms.Right panels also show L max calculated empirically using the optimized expressions for R C and B e through time (white line).

Figure 5 .
Figure 5. Local radius of curvature R C as a function of L* for a particle with 90°equatorial pitch angle crossing magnetic midnight, given by Equation 3 fitted to the TS04 magnetic field model (solid lines) and modified to fit observations of loss onset (dotted lines).

Figure 7 .
Figure 7. Ratio between simulation results shown in Figure 6 and corresponding observations shown in Figure 1, with red indicating an overestimate by the model and blue indicating an underestimate, alongside Dst index (top panel).Black regions indicate where data was unavailable.

Figure 8 .
Figure 8. Pre and post-storm profiles of 90°differential proton flux (blue and red respectively) observed by Van Allen Probe B (dots) and simulated using BAS-PRO (lines).

Figure 9 .
Figure 9. Pre and post-storm profiles of 60°differential proton flux (blue and red respectively) observed by Van Allen Probe B (dots) and simulated using BAS-PRO (lines).

Figure A1 .
Figure A1.Loss cone equatorial pitch angle α lc plotted against L. Shown is the centered dipole solution (black line) compared with a numerically derived offset dipole model given by Equation A3 for various loss altitudes h (colored lines).

Table 1
Minimum Dst Events Studied Between 1983 and 2023, and the Average Occurrence Rate per Year