© X. Xu et al., Published by EDP Sciences 2026

1 Introduction

The inner Van Allen Belt, extending from several hundred kilometers to nearly two Earth radii, contains dense population of energetic protons. These arise from a dynamic balance of many processes, including cosmic ray albedo neutron decay (CRAND), solar energetic particle injections, collisional losses, and radial diffusion (Singer, 1958; Farley & Walt, 1971; Farley et al., 1972; Baker et al., 2018; Li & Hudson, 2019). The inner belt comes closest to Earth’s surface in the South Atlantic Anomaly (SAA), where weaker magnetic field strength results from the dipole tilt, field eccentricity, and the offset from Earth’s geodetic center (Stassinopoulos et al., 2015). Intense trapped proton fluxes pose significant radiation risks to spacecraft systems and crew in this environment, making accurate modeling essential for mission planning and risk assessment (Jun et al., 2024).

Gradual geomagnetic field changes, coupled with atmospheric interactions, can significantly alter the trapped proton environment over decades (Heynderickx, 1996). The static AP8 model, based on observations in 1960s (Vette, 1966; Sawyer & Vette, 1976), no longer provides reliable estimates for current missions (Daly et al., 1996; Heynderickx et al., 1996a; Badavi et al., 2006). In addition, solar activity influences both particle source mechanisms and atmospheric interactions. This produces flux variations over the solar cycles which are particularly pronounced at altitudes below 1000 km (Nakano & Heckman, 1968; Dragt, 1971; Huston & Pfitzer, 1998; Huston et al., 1998; Qin et al., 2014; Li et al., 2020). While the AP9-International Radiation Environment Near Earth (IRENE-AP9) model incorporates a broad observational dataset and produces statistical outputs (Ginet et al., 2013; Johnston et al., 2015), at the time of this study, it does not explicitly incorporate solar cycle modulations. This limits the applicability for certain mission radiation risk assessments. This is especially true for human missions to Low Earth Orbit, which tend to last less than a year – a small fraction of the full solar cycle.

Trapped charged particles have three basic motions: gyration around the guiding magnetic field lines, bouncing between mirror points and drifting longitudinally around the Earth. The location of a particle’s drift shell within the radiation belts is commonly characterized by the McIlwain L parameter (McIlwain, 1961). A particle’s equatorial pitch angle α follows the relation $\sqrt{B_0/B}=\sin \left(\alpha \right)$, where B is the local field strength and B₀ is the equatorial magnetic field strength of the drift shell. Particles with lower equatorial pitch angles reach lower altitudes. For particles reaching low altitudes (<1000 km), the probability of atmospheric loss increases sharply as altitude decreases. The minimum equatorial pitch angle (α_m) required for a particle to complete the longitudinal drifts around the Earth defines the drift loss cone angle. From the perspective of a given location, flux at geomagnetic coordinates (L, α) includes particles with equatorial pitch angle in the range [α_m, α].

A key parameter for understanding how the trapped environment responds to secular changes in the geomagnetic field and solar cycle modulations is H_min, which has two context-dependent but physically unified definitions. For a location, it represents the lowest altitude among its equivalents in geomagnetic coordinates (L, α); for a particle, it is the minimum mirror altitude along its bounce-drift shell. Thus, a particle with drift shell L and equatorial pitch angle α shares the same H_min (L, α) as the corresponding geomagnetic location.

The residence time of a trapped particle, i.e., the duration from being trapped to loss, is strongly influenced by the atmospheric density along its drift shell. Particles with H_min in the high-altitude region can remain trapped for extended periods (Selesnick et al., 2007), while those with H_min in the low altitude region are removed more rapidly due to increased atmospheric interactions. At the so-called atmospheric cutoff altitude (~100 km), the atmospheric density becomes so high that any particles reaching this altitude are lost through collisions.

Ongoing changes in Earth’s geomagnetic field lead to a decreasing trend in H_min(L, α) (Farley et al., 1972; Lemaire et al., 1990; Selesnick et al., 2007). As a result, particles confined to a given drift shell (L, α) reach lower altitudes at later epochs, leading to an increase in α_m and a narrowing pitch angle range [α_m, α]. Additionally, particles mirroring at low altitudes within the same drift shell in later epochs are more likely lost due to increased atmospheric interactions. This explains why coupling 1960s observation-based AP8 flux database with the geomagnetic fields of later epoch’s results in unrealistic overestimates of proton flux in the SAA.

Solar cycle modulation of the trapped proton environment is primarily through solar activity-induced atmosphere density changes. Its effect on CRAND is relatively small (Selesnick et al., 2007; Bruno et al., 2021). Heightened solar activity heats Earth’s upper atmosphere, increasing density and collisional loss. Proton flux at low altitudes typically peaks near solar minimum when the atmosphere is coolest. A notable feature at the higher altitudes within the low-altitude region is the hysteresis effect, a delayed response in flux maxima or minima (Huston et al., 1998; Qin et al., 2014). This delay is related to the residence time of particles mirroring at these locations; longer residence times enable the influence of earlier solar conditions to persist and manifest in subsequent observations.

2 Satellite observations for model development

Developing a dynamic model of Earth’s trapped proton environment requires long-term, high-coverage observational data, ideally with high spatial resolution. The Van Allen Probes, launched in September 2012, consisted of twin satellites in highly elliptical orbits, with L-shell values ranging from 1.1 to 7 and altitudes spanning from several hundred kilometers to approximately 30,000 km. They provided extensive coverage of a broad region of the inner belt, including the locations with the highest proton fluxes. In contrast, the POES satellites, orbiting at ~840 km with an inclination of 98°, provided extensive spatial and long-term observations of the SAA region (Evans et al., 2008). Together, these datasets offer valuable insight into the evolving trapped proton environment.

The SEM2 instrument onboard the POES satellite provided proton count data for the >16 MeV, >36 MeV, >70 MeV, and >140 MeV energy channels (Evans & Greer, 2004). The RPS instruments onboard the Van Allen Probes measured protons across a broader energy range of approximately 60 MeV to 2000 MeV using 20 discrete channels (Mazur et al., 2013). As this work focuses on developing a dynamic modeling methodology for the trapped proton environment under an evolving geomagnetic field and varying solar cycle modulations, this analysis is limited to a single energy bin for the >70 MeV trapped protons. The POES dataset used in this work for model development consists of 16-second averaged Level 2 proton count data from NOAA-15 spanning from June 1, 1998, to December 31, 2013.¹ The RPS dataset comprises the 1-minute averaged level-2 proton observations from year 2013.²

Extensive data processing is performed to prepare satellite observations for model development. This includes calculating the RPS 1-minute integral flux of >70 MeV protons by summing the product of channel flux and channel energy width; computing yearly averaged >70 MeV POES measurements, binned into 1° longitude by 1° latitude grid boxes for each year from 1998 to 2013; calculating key magnetic field parameters, such as L-shell value, B, α, and H_min for each data point; removing background contributions from galactic cosmic rays (GCR); and performing cross-calibration between RPS and POES >70 MeV proton flux observations.

The twin RPS instruments (RPS-a and RPS-b), which follow nearly identical orbits, are expected to record similar fluxes. However, significant discrepancies exist between their measurements. Given their similar orbital configurations, these differences are likely due to measurement artifacts rather than the environment (O’Brien et al., 2021). It is found that equatorial fluxes, generally recognized as relatively stable over time, have good agreement between RPS-b observations and the simulations by IRENE-AP9 (V1.57, mean), whereas RPS-a observations are substantially lower. Since the IRENE-AP9 V1.57 fluxes have been extensively cross-checked and calibrated against multiple satellite datasets, RPS-a measurements are excluded from this model development. In summary, only RPS-b and POES measurements are used in this work.

In cross-calibrating the RPS-b and POES measurements for >70 MeV proton fluxes, a technique similar to that used in the IRENE-AP9 model database development (Ginet et al., 2013) is employed. Measurements from RPS-b and POES in 2013 are compared over the same geomagnetic region (L: 1.2–1.5, H_min > 500 km), where the background noise and GCR contributions are negligible relative to trapped proton fluxes. Analysis indicates that observations from the two satellites are approximately proportional, with POES measurements consistently several times larger. To align the datasets, a downscaling factor of 3.42 is applied to POES observations to match the RPS-b observation. This scale factor is consistent with the value reported in the release note and model validation report of IRENE-AP9 V1.50 (IRENE development team, 2017a, 2017b).

Analysis shows that storm-time solar energetic proton injections have very limited impact on the POES >70 MeV channel, even during highly active periods such as 2003. Therefore, no specific filtering is needed to remove the possible storm-time solar energetic proton injections. The yearly GCR background is extracted from POES observations at locations with H_min lower than 100 km and is described as a function of L. When subtracting the GCR background from RPS measurements, altitude-dependent variations due to Earth’s shadowing are taken into account.

The geomagnetic parameters for each data point are computed using the International Geomagnetic Reference Field model (IGRF-12; Thébault et al., 2015), with the epoch set to the midpoint of the observed year. The L parameter is computed with the method developed by Kluge (1972) and later modified by Bilitza (1992). The value of B₀ is computed as M/L³, where M is Earth’s dipole moment for the given epoch. H_min is calculated using a quadratic function of ln(B/B₀) for a given L-shell, following the approach suggested by Heynderickx et al. (1996b). Instructions and the coefficient tables for calculating H_min within the internal geomagnetic fields of AP8 and the IGRF-12 geomagnetic fields (covering epochs from 1995 to 2020 at 5-year intervals) are provided in Section 1 of the Supplement.

Figure 1 shows processed >70 MeV proton observations from RPS-b and POES at locations with H_min above 100 km for year 2013. The left panel uses the (L, B/B₀) coordinate system, while the right panel uses (L, H_min/H_min,_eq) to better visualize the spatial dependence of the flux distribution. H_min,eq denotes the H_min of the equatorial locations of the drift shell.

Figure 1 Flux observations (>70 MeV protons) by RPS-b and POES (cross-calibrated to RPS-b) in year 2013 at locations with Hmin above 100 km in two different coordinates: Left, (L, B/B0); Right, (L, Hmin/Hmin.eq). Dashed line: Hmin (L) = 1000 km, the reference line to separate the low-altitude region and the high-altitude region of the environment.

3 Model development

3.1 Modeling framework

The new model represents the trapped proton environment as the product of two primary components: one computes flux within a strictly adiabatic geomagnetic system, and another quantifies atmospheric collisional losses. In the strictly adiabatic system, the proton flux f at a given location (L, X) is expressed as

$$f\left(L,X\right)=N_0\left(L\right)·g\left(L,X\right),$$(1)

where the L-shell represents particle’s drift shell and X the bounce location. N₀ (L) denotes the flux of the adiabatic system on the geomagnetic equatorial plane. g(L, X) is the geomagnetic field confinement function, which gives the relative flux at the location normalized to N₀ (L). X is determined by the approach used in the g function. Two forms of the g function, g_b (L, B/B₀) and g_h(L, H_min), are independently developed in this work in Section 3.3.

Detailed in Section 3.4, the function r(H_min,s) provides a scale factor reflecting the impact of collisional loss, where s is a solar cycle activity component. The r function is built upon a scaled atmospheric density function associated with the H_min of the location. With atmospheric collisional loss accounted for, the proton flux in Earth’s trapped environment is calculated as

$$f\left(L,X,s\right)=N_0\left(L\right)·g\left(L,X\right)·\left[1-r\left(H_{\min },s\right)\right].$$(2)

3.2 N₀ (L)

Due to the lack of high-quality, cross-calibrated equatorial observations spanning multiple decades, incorporating long-term variation into N₀ (L) is currently impractical. With support from both theoretical reasoning and observational evidence, the model treats N₀ (L) as time independent.

In a slowly evolving geomagnetic system marked by a weakening dipole moment and contracting drift shells, variations in N₀ are expected to be subtle. Lemaire & Singer (2012) projected that fluxes decrease with the dipole weakening, while other studies indicated that shell contraction may increase flux (Farley et al., 1972; Selesnick et al., 2007). These opposing effects likely offset one another to some extent. Given the slow pace of geomagnetic field change relative to the model’s intended time applicability, significant variation in N₀ (L) over several decades is unlikely.

Observational data provide further support for the assumption of N₀ (L)’s temporal stability. Equatorial measurements remain relatively stable in regions where geomagnetic disturbance is minimal (He et al., 2023). In comparison, N₀ has greater stability, since it is unaffected by variations in atmospheric collisional loss driving by geomagnetic field changes or solar modulation.

From a theoretical standpoint, N₀ is not directly observable. Equatorial measurements include contributions from particles mirroring at low altitudes, where atmospheric collisional losses reduce flux intensity. As a result, the equatorial measurements inherently include the effects of atmospheric interaction, although such effects are likely minimal at many equatorial locations.

Constructing N₀ (L) from observations requires isolating it from atmospheric collisional loss effects. As the systematically available equatorial measurements (denoted as f_eq) in this work are limited to 2013, N₀ (L) in the model is given by

$$N_0\left(L\right)=\frac{f_{\mathrm{eq}, 2013}\left(L\right)}{1-r_{2013}\left(H_{\min ,\mathrm{eq}}\left(L\right), s\right)},$$(3)

where f_{eq, 2013} (L) is f_eq of year 2013 and r₂₀₁₃ (H_min,eq (L), s) is r(H_min, s) specified for the geomagnetic equatorial locations for year 2013. Extracted from the 2013 observations with B/(B₀) < 1.02, f_eq,2013 is shown in Figure 2.

Figure 2 Equatorial fluxes (>70 MeV protons) as a function of L-shell value. Discrete data points: observations in year 2013 by RPS-b for locations with B/B0 < 1.02; Solid line in red: feq,2013.

Replacing the N₀ (L) in equation (2) with the expression in equation (3) results in:

$$f\left(L,X,s\right)=\frac{f_{\mathrm{eq}, 2013}\left(L\right)}{1-r_{2013}\left(H_{\min ,\mathrm{eq}}\left(L\right), s\right)}·g\left(L,X\right)·\left[1-r\left(H_{\min },s\right)\right].$$(4)

This formulation is subsequently used to fit the observational data during the model development and serves as the basis for constructing the g and r functions.

3.3 Geomagnetic field confinement functions

This section presents the derivation and numerical formulations of the g_b (L, B/B₀) and g_h(L, H_min) functions, which compute the relative flux at a given location, normalized to N₀ (L), within a strictly adiabatic magnetic system.

A g_b (L, B/B₀) function

The distribution of trapped particles along the field line is highly dependent on the local magnetic intensity. The local flux, denoted J(L, B), can be computed by integrating contributions from particles with equatorial pitch angle in the range of [α_m, α]. A commonly used equatorial pitch angle distribution (PAD) form is [sin(α_i)]^k(L) (Selesnick et al., 2014; Ni et al., 2015; Shi et al., 2016; Roussos et al., 2018; Xu et al., 2019; Greeley et al., 2021). Accordingly, with the equatorial flux at 90° for a given L-shell denoted as j₀ (L), the local flux J(L, B) is calculated using the integral formulation provided in Roederer (1970):

$$J\left(L,B\right)=4\frac{B}{B_0}\underset{\cos \left(\alpha \right)}{\overset{\cos \left({\alpha }_m\right)}{\int }}{j}_0{\left(L\right)·\left[\sin \left({\alpha }_i\right)\right]}^{k\left(L\right)}·\frac{\cos \left({\alpha }_i\right)}{\sqrt{1-\frac{B}{B_0}\left(1-{\cos }^2\left({\alpha }_i\right)\right)}}d\left(\cos \left({\alpha }_i\right)\right).$$(5)

Function g_b (L, B/B₀), the relative flux normalized to N₀ (L), is calculated by J(L, B)/N₀ (L). Since N₀ (L) = J(L, B₀), function g_b (L, B/B₀) is given as

$$g_b\left(L,\frac{B}{B_0}\right)=\frac{J\left(L,B\right)}{J\left(L,B_0\right)}=\frac{\frac{B}{B_0}\underset{\cos \left(\alpha \right)}{\overset{\cos \left({\alpha }_m\right)}{\int }}{\left(\sin \left({\alpha }_i\right)\right)}^{k\left(L\right)}· \frac{\cos \left({\alpha }_i\right)}{\sqrt{1-\frac{B}{B_0}\left(1-{\cos }^2\left({\alpha }_i\right)\right)}}d\left(\cos \left({\alpha }_i\right)\right)}{\underset{0}{\overset{\cos \left({\alpha }_m\right)}{\int }}{\left(\sin \left({\alpha }_i\right)\right)}^{k\left(L\right)}·d\left(\cos \left({\alpha }_i\right)\right)} .$$(6)

It is worth noting that in an evolving geomagnetic field, the time dependence in g_b (L, B/B₀) is obtained through the time variation in α_m (L).

Calculations using equation (6) require the value of α_m of the shell, or equivalently, the corresponding magnetic field strength B_m. For a particle to complete the longitudinal drift around the Earth, the minimum altitude of its gyration-bounce-drift orbit is required to be above Earth’s surface. Accurate determination of B_m involves a three-dimensional search for the lowest mirroring point where the vertical component of the particle’s gyration radius equals H_min(L, α_m). Trapped protons come closest to Earth’s surface in the SAA region. To simplify the analysis, an approximation is used for the vertical component of the gyration radius near Earth’s surface within this region. Using the formulation in Armstrong et al. (1990), the gyration radius of a 100 MeV proton at the Earth’s surface in the SAA is approximately 60 km. Slightly smaller than the full gyration radius, the vertical component is approximated as 50 km. Accordingly, H_min(L, α_m) = 50 km is used to estimate the corresponding B_m.

Since the integration over the equatorial pitch angle range [α_m, α] cannot be performed analytically, a numerical approach is used in this work. To minimize uncertainties, the integration is carried out using a fine step size; specifically, one hundredth of the width of the equatorial pitch angle range [α_m, α] is used in the calculations presented in this work.

To determine the geomagnetic location dependence of component k, analysis of high-altitude RPS-b observations indicates that, for a given particle energy, modeling k as a function of L-shell alone is sufficient (see Fig. 3). Accordingly, constructing the g_b function for >70 MeV protons reduces to numerically specifying the component k(L).

Figure 3 y, >70 MeV proton flux observations normalized to feq(L), as a function of B/B0 in six L shells of the IGRF field for year 2013. Discrete data points: based on observations by RPS-b at high altitudes (Hmin > 1000 km) in year 2013; Solid line: gb function, computed with value of k indicated (see Eq. (6)).

B g_h (L, H_min) function

The g function has a value of 1.0 at the equatorial location. As the location moves along the field line toward Earth’s surface, it decreases gradually and reaches zero at the location where the corresponding equatorial pitch angle α equals α_m. To capture this behavior, a power-law function is deliberately chosen:

$$g_h\left(L,H_{\min }\right)={\left(\frac{H_{\min }- H_{\min }\left(L, {\alpha }_m\right)}{H_{\min ,\mathrm{eq}}\left(L\right)-H_{\min }\left(L, {\alpha }_m\right)}\right)}^{\eta (L, H_{\min })}.$$(7)

Preliminary analysis based on RPS-b observations indicates that, for a given particle energy, it is adequate to model the exponent η as a function of L-shell alone (see Fig. 4). With H_min (L, α_m) approximated as 50 km, as discussed in Section 3.3A, the g_h (L, H_min) equation simplifies to:

$$g_h\left(L,H_{\min }\right){=\left(\frac{H_{\min }- 50}{H_{\min ,\mathrm{eq}}\left(L\right)-50}\right)}^{\eta \left(L\right)}.$$(8)

Figure 4 y, >70 MeV proton flux observations normalized to feq(L), as a function of x = (Hmin − 50)/(Hmin,eq (L) − 50) in six L shells of IGRF field for year 2013. Discrete data points: based on observations by RPS-b at high altitudes (Hmin > 1000 km) in year 2013; Solid line: gh function of the shell (see Eq. (8)).

The remaining task in building the g_h function for >70 MeV protons is to numerically determine the component η(L).

3.4 Collisional loss function with solar cycle modulations

This section presents the derivation and numerical formulation of the collisional loss function r. Using collision theory (Pechukas, 1976), collisional loss rate is quantified by the following function:

$$r\left(H_{\min },s\right)={1-e}^{-\mathrm{a\rho }\left(H_{\min }, s\right)}.$$(9)

In the equation, a is a collisional interaction coefficient and ρ is a scaled atmospheric density that depends on the H_min of the location and solar cycle activity parameter s. The function r(H_min, s) increases with ρ, ranging from 0 (no collisional loss on local flux) to 1 (complete particle loss).

The scaled atmospheric density component ρ(H_min, s) is computed using a formulation similar to those used for estimating atmospheric density at satellite altitudes (Walker, 1964; Martin, 1965; Jacchia, 1970; Emmert, 2015):

$$\rho \left(H_{\min },s\right)=e^{-\left[{\left(\frac{H_{\min }}{z_1+z_2\bullet s}\right)}^b\right]}.$$(10)

In this expression, the solar modulation effect is incorporated into the atmospheric scale height term, z₁ + z₂s, where both z₁ and z₂ are positive-valued parameters.

The solar radio flux at 10.7 cm, known as F10.7 and having a long-standing history in atmospheric density modeling at satellite altitudes (Hedin, 1991; Bilitza, 2001; Picone et al., 2002; Emmert, 2015), is employed here as a proxy for tracking solar cycle activity levels.³ Particles mirroring at locations with higher H_min generally have longer residence times than those mirroring at low altitudes. As a result, the flux at locations with higher H_min in the SAA region is influenced by solar modulation effects over a longer temporal window. To incorporate this effect, the solar activity level parameter s is calculated as the average F10.7 over a period of m months preceding the target epoch, where m depends on H_min. Specifically,

$$s= \overline{F10.7\left(m\right)}, \hspace{1em}\mathrm{where} m=\frac{50}{1+e^{-0.0055\bullet \left(H_{\min }-500\right)}}$$(11)

This dependency of m on H_min is empirically derived from POES observations of >70 MeV protons. Complete evaluation of the collisional loss requires specifying the coefficients a, b, z₁, and z₂.

3.5 Completion of the model for >70 MeV protons

The final step in model development for >70 MeV protons is to specify the pitch angle distribution exponent k(L) needed to calculate g_b, the altitude dependence exponent η(L) needed to calculate g_h, and the coefficient a, b, z₁, and z₂ required for the collisional loss function r(H_min,s). This is accomplished by fitting equation (4) to the observational data from RPS-b (2013) and POES/SEM2 (1998–2013).

The optimization process involves three main steps. First, the numerical forms of k(L) and η(L), along with preliminarily coefficients estimates, are derived from RPS-b 2013 observations from locations with H_min > 1000 km, where the impact of collisional loss is low. Second, using the k(L) and η(L) from previous step, the collisional interaction coefficient a and the atmospheric scale height coefficients z₁, and z₂ are estimated by fitting to the POES/SEM2 1998–2013 observations, with the exponent b fixed at 1.0. Finally, all coefficients in the k(L), η(L), and r(H_min, s) functions are refined through fine-tuning to fit the satellite observations.

The forms of k(L) and η(L) are empirical, chosen for their effectiveness in reproducing the observation data. The resulting optimal expressions are:

$$k\left(L\right)=2.5+{2.5\left(L-1.1\right)}^{4.0}$$(12)

and

$$\eta \left(L\right)=1.39+0.2\left(L-1\right)+{\left(L-1\right)}^{5.5}.$$(13)

The optimal values for the coefficients used in the collisional loss function for >70 MeV proton flux are summarized in Table 1.

Analysis of the model simulations reveals that the optimal g_b and g_h functions for >70 MeV protons yield very similar values for shells with L < 2.0, with some noticeable differences emerging in shells with L > 2.0. A representative comparison of the g_b and g_h functions for year 2013 is presented in Figure 5. Analysis also shows that the slight difference in the optimal b values between the r functions for g_b and g_h results only minor differences between the two r functions.

Figure 5 Comparison of the optimal g functions for >70 MeV protons (red solid line: gb; dashed blue line: gh), as a function of the corresponding equatorial pitch angle α of a geomagnetic location, evaluated in six L shells of IGRF field for 2013.

While 50 km is a reasonable approximation for H_min (L, α_m), sensitivity tests using 0 km and 100 km are performed to investigate how the H_min (L, α_m) approximation impacts the model’s performance. The results indicate that the flux uncertainties introduced by the approximation are relatively small and primarily confined to low altitude locations. Since collisional loss to the atmosphere has a greater impact on fluxes at the low altitudes, the uncertainties are largely absorbed into the r function. Tuning the r function yields comparable model performance in fitting the observations (Figs. S1–S4). Overall, the model is robust to reasonable variations in the assumed value of H_min (L, α_m).

4 Results

This section compares model-calculated >70 MeV proton fluxes or the times and locations of the satellite measurements. In general, the model calculations agree well with satellite observations from RPS-b (2013) and POES (1998–2013), as shown in the scatter plots of Figure 6, where the coefficient of determination (R²) is provided as a measure of goodness of fit. Model outputs, using either the g_b or g_h function along with the respective r function, have similar performances.

Figure 6 Comparison of model calculated >70 MeV proton fluxes with observations. Left panels: model calculations with gb function. Right panels: model calculations with gh function; Upper panels: fluxes along the 1-minute RPS-b trajectory in year 2013; Lower panels: fluxes at POES satellite altitude for each year between 1998–2013; y = x: the perfect fit line (red); Coefficients of determination R2 are denoted.

Fluxes in the low-altitude region exhibit strong spatial gradients. To assess model performance, Figure 7 presents the >70 MeV trapped proton fluxes measured by POES in geographic coordinates for the years 2003 and 2010 (one year after the solar activity extremes), with H_min and L overlaid. The new dynamic model accurately locates the flux center and effectively captures the spatial gradient and solar modulation induced variation.

Figure 7 >70 MeV proton fluxes in the year following solar activity extremes at POES altitude in the geographic coordinate, with L (1.2, 1.35, 2.0, and 4.0, in black) and Hmin (100 km, 300 km, and 600 km, in red) isolines overlaid, and with the flux peaking area marked by the lines crossing at [50° W, 26° S]; Left panels: observations by POES; Middle panels: model calculations with gb function; Right Panels: model calculations with gh function.

Fluxes in the low-altitude region are particularly sensitive to changes in the geomagnetic field and solar cycle modulation due to their interaction with the neutral atmosphere. To evaluate the model’s capability in capture these flux changes, Figure 8 presents model-calculated yearly flux profiles from 1998 to 2013 at POES altitude along 26° S latitude and 50° W longitude, compared with corresponding POES observations. The new dynamic model successfully reproduces the variations over the 16-year span, and notably it accurately represents the yearly northwestward drift. This provides strong evidence of model’s ability to account for both the evolving geomagnetic field and the solar modulation.

Figure 8 Yearly >70 MeV proton fluxes crossing 26° S (upper panels) and 50° W (lower panels) at POES orbit altitude during 1998–2013 (year color coded). Left panels: POES observations; Middle panels: model calculations with gb function; Right panels: model calculations with gh function.

To examine the model’s ability to capture solar modulation-induced time variations, including the hysteresis effects, Figure 9 presents monthly F10.7 indices alongside the observed and model-simulated fluxes averaged over three H_min zones. Among the three zones, the hysteresis effect is most pronounced in the zone with the highest H_min. As the solar activity reaches a maximum in 2002 or a minimum in 2009, fluxes in the H_min >750 km zone do not reach the corresponding minimum or maximum until one or two years later. Although some uncertainties are evident in the high-altitude zone, the new model captures the time variations in the low-altitude region well, including the hysteresis effects.

Figure 9 Yearly variations in >70 MeV proton flux at POES altitude, averaged over the locations with L in the range of 1.25–1.6 for three Hmin zones; Observations by POES (black); model calculations with gb (solid red); model calculations with gh (blue dashed); Monthly F10.7 variation (green).

Although the model has demonstrated strong performance in reproducing the long-term POES observations, a thorough validation on longer timescales is constrained by the lack of high-quality observations across extended periods. To partially address this, the model’s capability to predict the northwestward drift of the low-altitude region is evaluated through its ability to reconstruct the SAA location in the AP8-MIN epoch, which corresponds to one of the earliest periods with systematic observations of the trapped proton belt. Figure 10 presents a comparison of model-simulated flux distributions at the POES altitude for the AP8-MIN epoch (1965) and for the year 2013, alongside AP8-MIN model outputs and POES 2013 observations, respectively. The results show that the new model effectively reproduces the drift of the low altitude region over the ~50-year span. For the 1965 environment simulations, the model captures the overall flux distribution well. The differences between the new model and AP8-MIN are largely due to the low spatial resolution of the AP8-MIN model and the lack of cross-calibration between RPS-b and the observing instruments used in the AP8-MIN model construction.

Figure 10 Northwestward drift of the low-altitude region of the SAA from 1965 to 2013. The plots show the contour lines of >70 MeV proton flux intensities at 840 km from: a) AP8-MIN; b) POES 2013 observations; c) model calculations (solid red line: with gb; dashed blue: with gh) for year 1965; d) model calculations for year 2013 (solid red line: with gb; dashed blue line: with gh).

To evaluate the model’s predictive capability beyond its development period, Figure 11 compares model-projected >70 MeV proton fluxes at POES altitude with yearly averaged POES/SEM2 observations for 2020 (solar minimum) and 2024 (solar maximum). The POES observations, available as raw Level 1a data,⁴ are processed using the same method outlined in Section 2, including a consistent scaling factor. The model reproduces the geographic distribution of >70 MeV proton fluxes with high fidelity. Direct comparisons between model simulations and POES observations yield correlation coefficients of 0.998 for 2020 and 0.985 for 2024, indicating excellent spatial agreement. Despite this strong correlation, the model underestimates by approximately 10%. The origin of this bias remains uncertain and may stem from limitations in the model or from the observational dataset. For example, differences in data processing levels or degradation of instrument hardware, as NOAA-15 approaches the end of its operational life, could distort the measurement data. For reference, AP8-MIN and AP8-MAX simulations (shifted 5° northward and 15° westward to match the observed flux peak), as well as IRENE-AP9 simulations (version 1.57, mean mode), are also compared with POES observations for both years.

Figure 11 >70 MeV proton fluxes in year 2020 (solar minimum) and year 2024 (approximately solar maximum) at POES altitude. upper row: POES observations; second row: simulations by the new method; third row: simulations by IRENE-AP9 (V1.57, mean mode); forth row: simulation by AP8-MIN and AP8-MAX, with geomagnetic field shifting (5° N, 15° W) to match the flux peak location; bottom row: scatter plots to compare model simulations with POES observations.

5 Discussion

High-performance trapped proton environment models are essential for mission radiation risk analysis. Modeling the low-altitude region is particularly challenging due to its strong spatial and temporal dependencies, driven by geomagnetic field changes and solar modulation, and further complicated by the particle-atmosphere interactions.

This work introduces a framework for a dynamic model composed of two primary components. The first calculates flux within a strictly adiabatic geomagnetic system, constructed from two key elements: N₀ (L), representing the equatorial fluxes of a strictly adiabatic geomagnetic system, and the g function, which provides the relative flux normalized to N₀ (L). The second primary component quantifies atmospheric collisional loss effects via the r function. The resulting flux is expressed as the product of N₀ (L), the g function, and the attenuation factor 1 − r. The IGRF-12 model is used for the geomagnetic field description. Observations from POES (1998–2013) and RPS-b (2013) are used in model construction. N₀ (L), derived from RPS-b data for the year 2013, is treated as time-invariant.

The model successfully reproduces the key observational features, such as the northwestward drift, solar cycle variations, and hysteresis effects. This indicates that the model framework and the formulation of the g and r functions effectively capture the underlying mechanisms. The model’s ability to predict the flux peak location and overall flux distributions in the low-altitude region for the AP8 epoch further supports its validity and suggests its applicability for predictions for future decades.

The trapped proton environment consists of protons over a wide range of energies and is highly directional, particularly at low altitudes. The applicability of the current modeling approach to a wider energy range, as well as its possible extension to a directional model, is discussed here.

Particle angle distribution in a strictly adiabatic system is governed by the radial and pitch angle diffusion processes (Roederer, 1970). Based on a widely used formulation for PAD, the term j₀ (L)[sin(α_i)]^k(L) is expected to remain valid across a broad energy range, while k(L) is expected to vary with energy. Such expectations are supported by the analysis presented in Figure S5. Moreover, with the PAD function embedded, the g_b function is well-suited for future extension from an omnidirectional model to a directional model.

Compared to the g_b function, the g_h function has a very simple structure. Preliminary analysis of the RPS-b observations confirms that the structure of the g_h function performs well across a wider range of particle energies, with energy dependence encapsulated in the component η (see Fig. S6 for an example). The g_h function does not incorporate flux directionality and that constrains its applicability in a directional model. However, its simpler structure offers a computational advantage. The g_b function requires numerical integration over the pitch angle range [α_m, α] with small step sizes and needs the value of α_m (or equivalently B_m) calculated, which is time-dependent and varies with L-shell. In contrast, the g_h is analytical and does not require α_m or B_m, as the relevant information is encoded in H_min (L, α_m), which can be approximated as 50 km and yields comparable modeling performance. Whether this performance persists across a wider energy range remains to be explored in the development of a full-spectrum model. If confirmed, g_h could serve as a fast and effective replacement for g_b when an omni-directional model is sufficient.

Collisional loss significantly contributes to flux directionality, with stronger attenuation in directions where particles mirror at lower altitudes. In a directional environment model, this directionality can be introduced via direction-associated H_min. Energy dependence in the r function is expected, as the gyroradius of trapped protons increases with proton energy, and further analysis is needed to quantify and incorporate it. Nonetheless, the general form of the r function, grounded in collisional theory, is expected to perform well across a broad energy range and for differential fluxes.

This modeling approach has certain limitations due to its simplifying assumptions, and several temporal variations are not accounted for. First, the Earth’s magnetic field is a sum of the internal and external fields. Substantial field changes during geomagnetic storms can lead to flux variations, most pronounced in shells near the outer edge of the radiation belt. Second, some solar activity-related processes are not included in the model. These include solar proton injections, which primarily affect fluxes below a few tens of MeV in the outer shells, and solar modulation of CRAND, whose effects are largely confined to low altitudes in very low-L shells (Selesnick et al., 2007; Bruno et al., 2021). Additionally, the slow change in N₀ induced by the geomagnetic field is not considered.

6 Conclusion

Accurate simulation of Earth’s trapped proton environment, with variations driven by geomagnetic field evolution and solar modulation incorporated, is essential for effective radiation risk management. Existing static models, such as AP8 and IRENE-AP9, do not fully capture these time-dependent dynamics. To address this limitation, this work introduces a dynamic modeling framework that incorporates effects of the geomagnetic field evolution and solar cycle modulation.

The model integrates two key components: one describing proton fluxes confined within a strictly adiabatic geomagnetic system, and another quantifying atmospheric collisional losses. Solar cycle modulation is embedded directly into the loss formulation. The model is currently parameterized for >70 MeV protons, using cleaned and cross-calibrated satellite observations from RPS-b (2013) and POES (1998–2013). While some uncertainties remain at higher altitudes within the SAA region, the model effectively reconstructs the observations, including the northwestward drift of the SAA, solar modulation effects, and hysteresis. It also successfully reconstructs the 1965 location of the SAA.

These results demonstrate that the model’s physics-based structure effectively captures the influence of underlying physical mechanisms on the trapped proton environment. Although the current fitting is limited to one single energy channel, the framework is designed to be extensible to a broader energy range and directional flux configurations, providing a robust foundation for future development.

Acknowledgments

We gratefully acknowledge assistance from the IRENE-AP9 development team, especially William Johnston of the Air Force Research Laboratory, who helped with the IRENE-AP9 model installation by providing valuable and technical input. The editor thanks Daniel Heynderickx and an anonymous reviewer for their assistance in evaluating this paper.

Declaration of AI-assisted technologies in the writing process

During the preparation of this work the authors used Microsoft Copilot to review the grammar and language. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.

Funding

A significant part of this work was performed under NASA’s Advanced Exploration System RadWorks project with support from the Durability, Damage Tolerance, and Reliability Branch of NASA Langley Research Center.

Data availability statement

F10.7 solar indices used in this study are available from https://lasp.colorado.edu/lisird/data/noaa_radio_flux, hosted at Interactive Solar Irradiance Datacenter, Laboratory for Atmospheric and Space Physics, University of Colorado Boulder.

POES/SEM2 data used in this study are available from https://www.ngdc.noaa.gov/stp/satellite/poes/dataaccess.html, hosted at National Geophysical Data Center of National Oceanic and Atmospheric Administration (NOAA).

Van Allen Probes RPS 1-minute averaged Level‐2 data of year 2013 are available from https://cdaweb.gsfc.nasa.gov/pub/data/rbsp/rbspb/l2/rps/psbr-rps-1min/2013/, hosted at the Coordinated Data Analysis Web, NASA Goddard Space Flight Center.

Supplementary material

Section 1. H_min calculations. Access Supplementary Material

Section 2. Figures.

Figure S1. Model calculated >70 MeV proton fluxes (gb function, $H_{\min }\left(L,{\alpha }_m\right)$ = 0 km; optimal collisional loss coefficients: a = 18.0, z1 = −35 km and z2 = 0.362, b = 1.170), and their comparisons with the observations by RPS-b (2013) and SEM2/POES (1998–2013). The figures are similar to Figures 6–9.

Figure S2. Model calculated >70 MeV proton fluxes (with gb function, $H_{\min }\left(L,{\alpha }_m\right)$ = 100 km, optimal collisional loss coefficients: a = 17.5; z1 = −50 km, z2 = 0.351, b =1.170), and their comparisons with the observations by RPS-b (2013) and SEM2/POES (1998–2013). The figures are similar to Figures 6–9.

Figure S3: Model calculated >70 MeV proton fluxes (with gh function, $H_{\min }\left(L,{\alpha }_m\right)$ = 0 km; optimal collisional loss coefficients:a = 18.0; z1= −35 km, z2= 0.362, b = 1.142), and their comparisons with the observations by RPS-b (2013) and SEM2/POES (1998–2013). The figures are similar to Figures 6–9.

Figure S4. Model calculated >70 MeV proton fluxes (with gh function and Hmin(L, αm)=100 km, optimal collisional loss coefficients: a = 17.5, z1 = −50 km, z2 = 0.351, b= 1.142), and their comparisons with the observations by RPS-b (2013) and SEM2/POES (1998–2013). The figures are similar to Figures 6–9.

Figure S5. Proton flux observations (normalized to the equatorial flux intensity of the shell),y, as a function of B/B0 in four L shells. Discrete data points: based on observations within the L-shell interval by RPS-b in year 2013 (red: >192 MeV; green: >300 MeV); Solid line: gb function, computed with value of k indicated (see Equation (6); Vertical dashed line: $H_{\min }\left( B/B_0\right)$= 1000 km, separating the high-altitude region and the low-altitude region of the environment.

Figure S6. y, proton flux observations (normalized to the equatorial flux intensity of the shell), as a function of $x=(H_{\min }-50)/(H_{\min ,\mathrm{eq}}\left(L\right)-50)$in four L shells. Discrete data points: based on observations within the L-shell interval by RPS-b in year 2013 (red: >192 MeV; green: >300 MeV); Solid line: gh function for the shell (see Equation (8), with the value of η indicated; Vertical dashed line: $H_{\min }\left( B/B_0\right)$ = 1000 km, separating the high-altitude region and the low-altitude region of the environment.

References

All Tables

All Figures

	Figure 1 Flux observations (>70 MeV protons) by RPS-b and POES (cross-calibrated to RPS-b) in year 2013 at locations with Hmin above 100 km in two different coordinates: Left, (L, B/B0); Right, (L, Hmin/Hmin.eq). Dashed line: Hmin (L) = 1000 km, the reference line to separate the low-altitude region and the high-altitude region of the environment.
In the text

	Figure 2 Equatorial fluxes (>70 MeV protons) as a function of L-shell value. Discrete data points: observations in year 2013 by RPS-b for locations with B/B0 < 1.02; Solid line in red: feq,2013.
In the text

	Figure 3 y, >70 MeV proton flux observations normalized to feq(L), as a function of B/B0 in six L shells of the IGRF field for year 2013. Discrete data points: based on observations by RPS-b at high altitudes (Hmin > 1000 km) in year 2013; Solid line: gb function, computed with value of k indicated (see Eq. (6)).
In the text

	Figure 4 y, >70 MeV proton flux observations normalized to feq(L), as a function of x = (Hmin − 50)/(Hmin,eq (L) − 50) in six L shells of IGRF field for year 2013. Discrete data points: based on observations by RPS-b at high altitudes (Hmin > 1000 km) in year 2013; Solid line: gh function of the shell (see Eq. (8)).
In the text

	Figure 5 Comparison of the optimal g functions for >70 MeV protons (red solid line: gb; dashed blue line: gh), as a function of the corresponding equatorial pitch angle α of a geomagnetic location, evaluated in six L shells of IGRF field for 2013.
In the text

Figure 6 Comparison of model calculated >70 MeV proton fluxes with observations. Left panels: model calculations with gb function. Right panels: model calculations with gh function; Upper panels: fluxes along the 1-minute RPS-b trajectory in year 2013; Lower panels: fluxes at POES satellite altitude for each year between 1998–2013; y = x: the perfect fit line (red); Coefficients of determination R2 are denoted.

In the text

Figure 7 >70 MeV proton fluxes in the year following solar activity extremes at POES altitude in the geographic coordinate, with L (1.2, 1.35, 2.0, and 4.0, in black) and Hmin (100 km, 300 km, and 600 km, in red) isolines overlaid, and with the flux peaking area marked by the lines crossing at [50° W, 26° S]; Left panels: observations by POES; Middle panels: model calculations with gb function; Right Panels: model calculations with gh function.

In the text

	Figure 8 Yearly >70 MeV proton fluxes crossing 26° S (upper panels) and 50° W (lower panels) at POES orbit altitude during 1998–2013 (year color coded). Left panels: POES observations; Middle panels: model calculations with gb function; Right panels: model calculations with gh function.
In the text

	Figure 9 Yearly variations in >70 MeV proton flux at POES altitude, averaged over the locations with L in the range of 1.25–1.6 for three Hmin zones; Observations by POES (black); model calculations with gb (solid red); model calculations with gh (blue dashed); Monthly F10.7 variation (green).
In the text

	Figure 10 Northwestward drift of the low-altitude region of the SAA from 1965 to 2013. The plots show the contour lines of >70 MeV proton flux intensities at 840 km from: a) AP8-MIN; b) POES 2013 observations; c) model calculations (solid red line: with gb; dashed blue: with gh) for year 1965; d) model calculations for year 2013 (solid red line: with gb; dashed blue line: with gh).
In the text

Figure 11 >70 MeV proton fluxes in year 2020 (solar minimum) and year 2024 (approximately solar maximum) at POES altitude. upper row: POES observations; second row: simulations by the new method; third row: simulations by IRENE-AP9 (V1.57, mean mode); forth row: simulation by AP8-MIN and AP8-MAX, with geomagnetic field shifting (5° N, 15° W) to match the flux peak location; bottom row: scatter plots to compare model simulations with POES observations.

In the text

	Figure S1. Model calculated >70 MeV proton fluxes (gb function, $H_{\min }\left(L,{\alpha }_m\right)$ = 0 km; optimal collisional loss coefficients: a = 18.0, z1 = −35 km and z2 = 0.362, b = 1.170), and their comparisons with the observations by RPS-b (2013) and SEM2/POES (1998–2013). The figures are similar to Figures 6–9.
In the text

	Figure S2. Model calculated >70 MeV proton fluxes (with gb function, $H_{\min }\left(L,{\alpha }_m\right)$ = 100 km, optimal collisional loss coefficients: a = 17.5; z1 = −50 km, z2 = 0.351, b =1.170), and their comparisons with the observations by RPS-b (2013) and SEM2/POES (1998–2013). The figures are similar to Figures 6–9.
In the text

	Figure S3: Model calculated >70 MeV proton fluxes (with gh function, $H_{\min }\left(L,{\alpha }_m\right)$ = 0 km; optimal collisional loss coefficients:a = 18.0; z1= −35 km, z2= 0.362, b = 1.142), and their comparisons with the observations by RPS-b (2013) and SEM2/POES (1998–2013). The figures are similar to Figures 6–9.
In the text

	Figure S4. Model calculated >70 MeV proton fluxes (with gh function and Hmin(L, αm)=100 km, optimal collisional loss coefficients: a = 17.5, z1 = −50 km, z2 = 0.351, b= 1.142), and their comparisons with the observations by RPS-b (2013) and SEM2/POES (1998–2013). The figures are similar to Figures 6–9.
In the text

Figure S5. Proton flux observations (normalized to the equatorial flux intensity of the shell),y, as a function of B/B0 in four L shells. Discrete data points: based on observations within the L-shell interval by RPS-b in year 2013 (red: >192 MeV; green: >300 MeV); Solid line: gb function, computed with value of k indicated (see Equation (6); Vertical dashed line: $H_{\min }\left( B/B_0\right)$= 1000 km, separating the high-altitude region and the low-altitude region of the environment.

In the text

Figure S6. y, proton flux observations (normalized to the equatorial flux intensity of the shell), as a function of $x=(H_{\min }-50)/(H_{\min ,\mathrm{eq}}\left(L\right)-50)$in four L shells. Discrete data points: based on observations within the L-shell interval by RPS-b in year 2013 (red: >192 MeV; green: >300 MeV); Solid line: gh function for the shell (see Equation (8), with the value of η indicated; Vertical dashed line: $H_{\min }\left( B/B_0\right)$ = 1000 km, separating the high-altitude region and the low-altitude region of the environment.

In the text