The Badhwar‐O'Neill 2020 GCR Model

The Badhwar‐O'Neill (BON) model has been used for some time to describe the galactic cosmic ray (GCR) environment encountered in deep space by astronauts and sensitive electronics. The most recent version of the model, BON2014, was calibrated to available measurements to reduce model errors for particles and energies of significance to astronaut exposure. Although subsequent studies showed the model to be reasonably accurate for such applications, modifications to the sunspot number (SSN) classification system and a large number of new high‐precision measurements suggested the need to develop an improved and more capable model. In this work, the BON2020 model is described. The new model relies on daily integral flux from the Advanced Composition Explorer Cosmic Ray Isotope Spectrometer (ACE/CRIS) to describe solar activity. For time periods not covered by ACE/CRIS, the updated international SSN database is used. Parameters in the new model are calibrated to available data, which include the new Alpha Magnetic Spectrometer (AMS‐02) and Payload for Antimatter Matter Exploration and Light‐nuclei Astrophysics (PAMELA) high‐precision measurements. It is found that the BON2020 model is a significant improvement over BON2014. Systematic bias associated with BON2014 has been removed. The average relative error of the BON2020 model compared to all available measurements is found to be <1%, and BON2020 is found to be within ±15% of a large fraction of the available measurements (26,269 of 27,646 → 95%).


Introduction
The galactic cosmic ray (GCR) environment in space is a complex spectrum of ions and energies that pose a significant risk to astronauts and sensitive electronics. Accurately and efficiently describing this radiation environment is therefore vital for mission planning, shield design, and astronaut safety. The Badhwar-O'Neill (BON) model was created to meet this need and is periodically updated as new data become available (Adriani et al., 2011(Adriani et al., , 2013Aguilar et al., 2017Aguilar et al., , 2018aAguilar et al., , 2018bAguilar et al., , 2018c and improved methods for calibration and uncertainty quantification are developed (Slaba et al., 2014;Slaba & Blattnig, 2014a, 2014b. Although differences between past versions of the BON model exist, all of them rely on the same basic formalism (O'Neill et al., 2015). First, the local interstellar spectrum (LIS) for each ion is described by a simple parametric function. Next, the state of the heliosphere at any time and radial distance from the Sun is described by a diffusion coefficient, which depends, in part, on the time-dependent solar modulation potential, ϕ(t) (MV). The modulation potential is obtained through empirical relationships to solar indices such as neutron monitor counts or sunspot number (SSN). Finally, the LIS is transported to 1 AU using a numerical solution to the steady-state, spherically symmetric Fokker-Planck equation with the aforementioned diffusion coefficient. The main differences between the various versions of the BON model are found in the solar indices used to derive modulation potential and the data and methods used to calibrate free parameters in the LIS function for each ion. The most recent version, BON2014 (O'Neill et al., 2015), used SSN to obtain the modulation potential, and parameter calibration was guided by sensitivity studies highlighting the energies and ions of interest to space radiation protection (Slaba & Blattnig, 2014a).
In~2015, the Royal Observatory of Belgium revised their SSN classification system (Clette et al., 2016). The resulting SSN database, referred to as version 2 (V2), contained markedly different values than the original database (V1). The change in SSN values had a significant impact on ϕ(t) and the subsequent behavior of the GCR energy spectrum predicted by BON2014. The observatory abandoned the V1 database after May 2015, thereby leaving no means of determining ϕ(t) from the original SSN classification system beyond the termination date. The V2 values could also not be used directly in BON2014 without completely rederiving the ϕ(t) relationship. Further analysis revealed that between January 2014 and May 2015, the ratio of V1 to V2 values approached an average value of~0.7. To enable the use of BON2014 beyond the May 2015 termination date of the V1 database, the V2 values were simply scaled by 0.7. This provided a temporary, albeit suboptimal, extension of the V1 database to current dates. Although this approach continues to be used for the BON2014 model, it is clear that more robust methods of determining ϕ(t) are needed.
Along with the SSN classification system changing, there has been a significant number of new and important measurements from the Alpha Magnetic Spectrometer (AMS-02; Aguilar, et al., 2015aAguilar, et al., , 2015bAguilar, et al., , 2017Aguilar, et al., , 2018aAguilar, et al., , 2018bAguilar, et al., , 2018c and the Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics (PAMELA; Adriani et al., 2011Adriani et al., , 2013Adriani et al., , 2017Martucci et al., 2018). Some of these data have been used in comparative studies with available GCR models (Norbury et al., 2018;Whitman et al., 2019), revealing areas where spectral characteristics and solar modulation effects can be improved. Both experiments have also published GCR measurements on monthly timescales covering the descending phase of solar cycle 23 and the ascending phase of cycle 24. The AMS-02 collaboration released monthly proton and helium data from May 2011 to May 2017 (Aguilar et al., 2018c), and the PAMELA group released monthly proton data from October 2006 to January 2014 (Adriani et al., 2013;Martucci et al., 2018). Since GCR protons and helium account for~70% of the astronaut exposure behind shielding (Slaba & Blattnig, 2014a), it is clear that these new datasets fill an important measurement gap for space radiation protection and create a new opportunity to improve models.
In this work, the BON2020 model is presented. This updated model contains new methods to compute ϕ(t), and the new data from AMS-02 and PAMELA are utilized in LIS parameter calibration. Using previously described methods for uncertainty quantification (Slaba et al., 2014), it is shown that BON2020 is significantly better than its predecessor, BON2014. Uncertainty propagation methods are used to show that these updates reduce propagated errors in effective dose equivalent as compared to BON2014. We refrain from making direct comparisons to other GCR models, as these models have not yet been updated as a result of the new measurements. It is expected that all of the models will eventually account for the new data, at which time a fair comparison may be made.

Model Description
The BON model is based on a numerical solution to the Fokker-Planck equation accounting for diffusion, convection, and adiabatic deceleration under the assumptions of a quasi-steady state and spherically symmetric interplanetary medium (O'Neill et al., 2015). The boundary condition for the Fokker-Planck equation is the LIS, which is described within the model for each ion as where T is the ion kinetic energy (MeV/n), m is the proton rest mass (MeV), β is the ion velocity relative to the speed of light, and j 0 , γ, and δ are free parameters for each ion. The units of f LIS are particles/(m 2 -MeV/ n-sr-sec). The values of j 0 and γ fix the magnitude and high-energy slope of the GCR spectrum for each ion. The δ parameter controls the lower energy behavior in conjunction with the solar modulation potential.
The state of the heliosphere at any radial distance from the Sun, r, and time, t, is described by the diffusion coefficient, κ, evaluated as where R is the ion rigidity (MV), ϕ(t) the solar modulation potential, V sw = 400 km/s is the solar wind speed (specified as a constant), r 0 = 4 AU, and k 0 = 8.8 × 10 20 cm 2 /s. The values selected for V sw , r 0 , and k 0 have been discussed elsewhere (O'Neill, 2006(O'Neill, , 2010 and were left unchanged for BON2020. Mertens et al. (2013) describe the theoretical basis for the functional form of equation 2.
The overall mechanism of particle transport is also unchanged in BON2020. Mertens et al. (2013) provide a complete derivation and discussion of the transport equation, and the numerical methods used to solve the Fokker-Planck equation are described in detail elsewhere . A summary of the transport equation and relevant quantities is provided here.
The following details apply to BON2020, BON2014, and all previous versions of the Badhwar-O'Neill model. In order to produce modulated spectra at Earth (1 AU), an LIS is generated for all ions at the edge of the heliosphere, assumed to be located at 100 AU. This spectrum is transported to 1 AU according to the particle transport equation, first proposed by Parker (1965) and written by Fisk (1970) in the form of a Fokker-Planck equation, given by where r is the distance from the Sun, U(r,T) is the omnidirectional particle density distribution, and α = (T + 2 m)/(T + m). As mentioned above, this equation assumes a quasi-steady, spherically symmetric heliosphere and includes the effects of convection of particles with the solar wind, adiabatic cooling due to the expansion of the solar wind, and diffusion of particles due to scattering off irregularities in the magnetic field.
The level of modulation is controlled by the single parameter, ϕ(t), in the diffusion coefficient. As described in the introduction, various solar indices have been used to determine ϕ(t), and BON2014 relied on SSN with an empirical delay function (Nymmik, 2000) that accounts for the time-lag between observed sunspots and magnetic field disturbances capable of modulating GCR spectra near Earth. Because the BON model relies on a one-dimensional transport equation, the time delay also serves as a proxy for the three-dimensional drift effects that cause positive particles to drift into the heliosphere over the poles when the sun's polarity is positive and along the heliospheric current sheet when the sun's polarity is negative. This drift effect, typically called the charge-sign drift effect, is responsible for the alternating peaked and plateau shapes observed in neutron monitor and other GCR data from one solar cycle to the next (Lockwood et al., 2001;Potgieter, 2014). The two different paths that GCR particles take into the heliosphere (over the solar poles or along the heliospheric current sheet) affect the timescales on which particles will start to modulate at Earth following changes in the solar magnetic field. Because this is a 3D effect, it cannot be modeled directly by BON; however, the time-lag function within the model is constructed to treat positive and negative solar cycles differently, serving to approximate the drift effect and produce alternating peaked and plateau solar cycles.
Comparison against data from the Advanced Composition Explorer Cosmic Ray Isotope Spectrometer (ACE/CRIS) as a function of time revealed this approach to be reasonably accurate, although moderate uncertainties in excess of 20% were found in some time periods (O'Neill et al., 2015). It is important to note that in most circumstances, the BON model is evaluated by specifying mission start and end dates. The function relating SSN to ϕ is integrated over this time period to obtain an average value of ϕ to be used in equation 2. Alternatively, a value of ϕ may be input directly into the model, thereby bypassing the dependence on SSN or any other solar index. One must be aware in such cases that externally derived models for ϕ should not be used directly in the BON model.

Parameter Calibration
For BON2020, the two primary goals were to develop a more descriptive and robust solar modulation function, ϕ(t), and to utilize the new AMS-02 (Aguilar et al., 2015a(Aguilar et al., , 2015b(Aguilar et al., , 2017(Aguilar et al., , 2018a(Aguilar et al., , 2018b(Aguilar et al., , 2018c and PAMELA (Adriani et al., 2011(Adriani et al., , 2013(Adriani et al., , 2017Martucci et al., 2018) data sources to recalibrate the LIS parameters. More specifically, for the solar modulation potential, we wanted to extend the capability of the model to utilize either ACE/CRIS data or SSN. Although it is expected ACE/CRIS data should provide a better description of solar activity than SSN for the applications in which the BON model is used, there are still advantages in using SSN. For example, space radiation protection studies often consider historical solar epochs (e.g., 1977 solar minimum) over which ACE/CRIS data are unavailable.
This section is organized as follows. First, the computational tools used to determine ϕ from ACE/CRIS integral flux and SSN are described. These tools are utilized in the LIS parameter calibration and are therefore described first. Next, the multi-step parameter calibration procedure developed for ions between Z = 5 and Z = 28 is discussed. Finally, the modified procedure for Z = 1-4 ions is defined. Specific results are shown and discussed for Z = 8 throughout this section for illustrative purposes.
3.1. ϕ Z (t) From ACE/CRIS ACE/CRIS measures low energy (<500 MeV/n) ions between Z = 5 and Z = 28 that may be used as a proxy for solar modulation. For example, Matthia et al. (2013) used the ACE/CRIS carbon flux in a specified energy range to define solar modulation in their model. Corti et al. (2016Corti et al. ( , 2019 analyzed BESS and AMS-02 data and found that protons and helium have distinct modulation characteristics in some epochs due to differing mass-to-charge ratios (which therefore affects the β term in the diffusion coefficient). This suggests that a single data source (e.g., ACE/CRIS data for a single ion) may not precisely specify solar modulation for all ions.
To account for this, Z-dependent data from ACE/CRIS may be used to define a time-dependent solar modulation potential function for each ion, referred to as ϕ Z (t). Oxygen (Z = 8) is used in this section as an example to help describe the computational procedures.
First, we utilize the BON2014 LIS parameters for oxygen (j 0 = 1.50 × 10 −6 , γ = 2.725, δ = −1.9) and evaluate the model over a range of input ϕ values in the interval [100, 2,000]. This interval extends over the range of ϕ values estimated for the past 250 years from SSN. The resulting model fluxes are integrated over the ACE/CRIS energy bands to provide an integral flux as a function of ϕ. Likewise, the 27 day average fluxes measured by ACE/CRIS are integrated across the same energy bins. This integration is performed by multiplying the measured differential fluxes in each energy bin of the data set by their bin widths and summing across all bins.
Note that ACE/CRIS energy bands only cover the low-energy portion of the spectrum for heavy ions (50-500 MeV/n), and the integral fluxes being discussed in this section are therefore only a representative component of the total integral flux for each ion. Nevertheless, ACE/CRIS covers an energy region that is highly sensitive to solar activity and therefore very useful in quantifying solar modulation effects.
A comparison of the measured and model-generated integral flux values is provided in Figure 1. From this comparison, a value ϕ can be determined that provides the best agreement between the model and each of the ACE/CRIS measurements. Figure 2 shows the optimal ϕ values obtained from this analysis plotted as a function of ACE/CRIS integral flux.
The data in Figure 2 were found to be well-described by the simple function where F ACE is the integral ACE/CRIS flux, and a and b are fit parameters.
For the oxygen data shown in Figure 2, it was found that a = 52.51, b = 6.50, and the fitting produced a coefficient of determination value of r 2 = 0.99999126 (perfect agreement corresponds to r 2 = 1). The computational procedure for determining the optimal ϕ may be repeated for each ion between Z = 5 and Z = 28 resulting in a unique solar modulation function for each Z.
Since F ACE is easily obtained from the ACE/CRIS data set as a function of time for a given ion, the time-dependent solar modulation potential may be more generally written as where the Z-dependence of the a and b parameters has been made explicit for clarity. This approach produces a distinct value of ϕ for each ion and is markedly different than past versions of the BON model where a common ϕ was defined for all ions.

ϕ Z (t) From ACE/CRIS Oxygen Flux
Driving the modulation of each ion directly by its associated ACE/CRIS integral flux has the benefit of linking the model results directly with real measurements of the GCR environment for a specified mission duration. This approach becomes difficult to interpret, however, if the model needs to be evaluated with a prescribed level of solar modulation (i.e., input ϕ) instead of mission dates. In this case, a user would need to know the value of ϕ for all GCR ions, which adds significant complexity and possible confusion when using the model.
A simpler approach is to link the behavior of all ions directly to ACE/CRIS oxygen. We first consider the correlation between ACE/CRIS integral oxygen flux and all other ions. An example is shown in Figure 3 for oxygen and silicon where a high degree of correlation between the ions is observed. The solid line in the plot was obtained from linear regression (r 2 = 0.951), and the dashed lines show that a majority of the silicon flux values are within ±15% of the simple linear fit. The regression coefficients between ACE/CRIS integral oxygen flux and all other ions between Z = 5 and Z = 28 were computed and are given in Appendix A.
The computational procedure for determining ϕ Z (t) from only ACE/CRIS integral oxygen flux is described as follows. For a given set of input mission dates, the average ACE/CRIS integral oxygen flux is determined. The regression coefficients from Appendix A are used to estimate all other ACE/CRIS integral ion fluxes between Z = 5 and Z = 28. These approximate integral flux values are inserted into equation 5 to obtain ϕ Z (t). In cases where the solar modulation potential, instead of mission dates, are specified for model evaluation, only a single input ϕ value is required. In this case, the value is interpreted by the model as the oxygen solar modulation potential. The inverse of equation 4, given as is evaluated to determine the associated ACE/CRIS oxygen integral flux value. The regression coefficients from Appendix A are employed to estimate the ACE/CRIS integral fluxes for all other ions. Finally, these approximate integral flux values are inserted back into equation 5 to obtain ϕ Z (t). This computational procedure guarantees self-consistent results whether mission dates or a specific level of solar modulation are considered. Figure 4 shows the relationship between V2 SSN and ACE/CRIS integral oxygen flux as a function of time. Using the methods from BON2014, the monthly SSN values have been smoothed over a 12 month forward window and are plotted as a function of the delayed-time determined from the empirical model of Nymmik (2000). Also shown in the plot is a parametric fit that relates time-delayed SSN values to ACE/CRIS integral oxygen flux. The function is given as
To obtain the oxygen solar modulation function from time-delayed SSN, equation 7 is evaluated to estimate the ACE/CRIS oxygen flux, and the solar modulation potential is determined from equation 5. Note that the parametric function in equation 7 depends only on the time-delay function of Nymmik (2000), but not on the LIS parameters or theoretical transport formalism of the BON model.
The computational procedure for determining ϕ Z (t) from SSN for each ion is described as follows. For a given set of input mission dates, the time-delayed average SSN is computed (exactly the same as in BON2014). The average SSN value and equation 7 are used to estimate the associated ACE/CRIS integral oxygen flux value, F 8 . The procedure described in the previous section for obtaining ϕ Z (t) from ACE/CRIS integral oxygen flux is then carried out.

ϕ Z (t) Summary
For ions between Z = 5 and Z = 28, the solar modulation potential is obtained from ACE/CRIS data, if available. If the input mission dates fall outside the time-period covered by ACE/CRIS, the solar modulation potential is determined from the V2 SSN values which extend from the year~1750 to present day. The oxygen solar modulation potential is used for ions between Z = 1 and Z = 4.

Parameter Calibration (Z = 5-28)
The computational procedure outlined in this section was used to calibrate the three LIS parameters (j 0 , γ, δ) for each ion. A flowchart of the process is provided in Figure 5. We begin by selecting the BON2014 LIS parameters as the initial values and label them as j (1) , γ (1) , and δ (1) . Using these parameters, an initial solar modulation function, ϕ Z (t), is determined with methods outlined in the previous subsections. This step provides a solar modulation function that ensures good agreement between the model and measurement data that are sensitive to solar activity.
The parameters j (1) and γ (1) set the magnitude and high-energy slope of the LIS and ideally would be calibrated against measurements taken at sufficiently high energy to be insensitive to solar activity; however, gaps in the measurement database make this difficult in practice. These parameters are therefore calibrated using data extending down to 1.5 GeV/n, thereby increasing the amount of available data but introducing moderate dependence on solar activity. Fortunately, the first step in this process yields a solar modulation potential ensuring a match between the low energy part of the model spectrum and ACE/CRIS measurements. This means that the initial LIS parameter δ (1) coupled with ϕ Z (t) well-represents the level of modulation, allowing j (1) and γ (1) to be calibrated against available data >1.5 GeV/n while maintaining a reliable description of lower energy regions that are heavily influenced by solar modulation effects.
To perform the calibration for j (1) and γ (1) , the model is evaluated over a broad range of j (1) and γ (1) values and compared to selected measurements. The following error metrics are considered (9) where f mod and f meas are the model and measured fluxes, respectively, and N denotes the total number of measurements in energy bins found in a prescribed energy range. The error metrics were computed separately for measurements covering the energy range 1.5-4 GeV/n and measurements covering the energy range above 4 GeV/n. This yielded six error metrics for each j (1) and γ (1) pair. The average errors associated with each pair were sorted by rank, and an average rank was computed by giving equal weight to the six metrics. The combination of parameters giving rise to the best average rank was identified as the optimal solution. It was found that separating the energy domain consistent with prior sensitivity studies (Slaba & Blattnig, 2014a) and considering three error metrics instead of one provided better results that were less likely to introduce systematic model bias. The optimal parameter pair is labeled as j (2) and γ (2) . A comparison of the oxygen flux computed by BON2014 and the corrected model to AMS-02 measurements (Aguilar et al., 2017) is shown in Figure 6. It can be seen that the corrected model is in excellent agreement with AMS-02 and is clearly an improvement from the BON2014 result.
We now have an intermediate set of corrected parameters, denoted as j (2) , γ (2) , and δ (1) . Before optimizing δ (1) , the correlation between LIS parameters and the solar modulation function needs to be addressed. To accomplish this, the solar modulation function, ϕ Z (t), is recalculated many times with the corrected parameters j (2) and γ (2) and a range of δ (1) values in the interval [−5, 5].
The final step in the calibration process is to correct δ (1) . In this part, the model is evaluated using j (2) and γ (2) and a range of δ (1) values in the interval [−5, 5] with the corresponding optimal solar modulation function for each δ (1) . Model results are then compared to all available measurements. The following error metrics are considered for this step (Slaba & Blattnig, 2014b) where w i are the weights defined over specific energy bins identified in Table 2 of Slaba and Blattnig (2014a), and the error metrics R D and R |D| have already been defined. Equations 11 and 12 are simply the average relative errors computed over the predefined energy bins. Equations 13 and 14 quantify the expected relative error propagated into effective dose equivalent behind 20 g/cm 2 shielding as a result of errors in the GCR model. It was found that consideration of all four error metrics provided an improved global fit compared to the outcomes achieved with just using one of the metrics. The average errors associated with each δ (1) value were sorted by rank, and an average rank was computed by giving equal weight to the four metrics. The δ (1) value giving rise to the best average rank was identified as the optimal parameter, δ (2) .
The flow chart in Figure 5 suggests a repeatable loop wherein LIS parameters can be further refined. It was found that repeating the process beyond the first iteration provided negligible corrections to the LIS parameters for many of the ions. This was especially true for ions characterized by high precision measurements from AMS-02 or PAMELA. For low abundance ions that are characterized by limited measurement data with moderate uncertainty, repeated loops of the optimization scheme converged to distinct results that depended on the combination of error metrics chosen (see equations 8-14), as might be expected. For consistency, LIS parameters for all ions were determined by using only one iteration of the optimization scheme in Figure 5. Comparisons of the final model to measurement data in the next section provide ample justification for this overall approach.

Modified Process for Z = 1-4
For ions with Z = 1-4, there are no ACE/CRIS measurements to find an optimal solar modulation potential. The available measurements for these ions also do not contain sufficient time resolution to directly quantify the solar modulation function. As a result, the oxygen solar modulation function is used for these ions. The remainder of the calibration process outlined in the previous section is then repeated to obtain corrected LIS parameters for these ions. The parameters obtained for all ions are given in the Appendix B.
• BON2020 ACE-O -Same as BON2020 ACE , except that only daily ACE/CRIS integral oxygen flux is used to described solar activity. The linear regression coefficients are used to estimate ACE/CRIS integral fluxes for all other ions. • BON2020 SSN -Same as BON2020 ACE , except that only V2 SSN values are used to obtain the solar modulation potential.
These extra comparisons are being included to demonstrate that the options for describing solar activity within the model are self-consistent and yield similar, if not identical, results in most cases. This consistency is particularly important for historical analyses wherein the model may be driven directly by ϕ, SSN, or ACE/CRIS data.
A summary of the measurements used for LIS parameter calibration and uncertainty quantification is given in Table 1. Visualizing the measurement data demonstrates the significant coverage for some ions in energy and time and extremely sparse coverage for others. Figure 7 shows the large number of observations for H (left) and He (right). The most recent high-precision AMS-02 fluxes better define the H and He spectral shapes at high energies. Both AMS-02 (Aguilar et al., 2018c) and PAMELA (Martucci et al., 2018) recently published H fluxes on monthly timescales. AMS-02 also published 78 He fluxes on monthly timescales (Aguilar et al., 2018c).   Figure 9, have been observed by multiple experiments across a wide energy range. However, measurements at high energies have large errors bars. The multitude of ACE/CRIS data below 500 MeV/n for Z > 5 are not shown in the plots but represent a significant source of data for model calibration and validation, as will be seen later in this section. The χ 2 metric also provides a measure of dispersion between model and measurement, but because the differences are squared, extreme errors are amplified in the summary value.
It can be seen in Table 2 that all variants of BON2020 are markedly better than BON2014. The systematic under prediction associated with BON2014 has been corrected, and the overall disparity has been reduced by more than a factor of two. One can also see that the average errors associated with BON2020 ACE and BON2020 ACE-O are nearly identical, suggesting that driving the model with all ACE/CRIS ion data or just ACE/CRIS oxygen data yields similar results. It can be seen that errors associated with BON2020 SSN are larger than BON2020 ACE or BON2020 ACE-O , as expected since the ACE/CRIS data provide a better description of solar modulation than SSN.
Further information is gained by examining Figure 10, which shows the distribution of model errors giving rise to the average values shown in Table 2. The distributions in the figure for the BON2020 variants are centered near 0% and appear to be symmetric, suggesting no systematic bias in the models. The BON2014 distribution is shifted towards negative values and is also noticeably broader than any of the BON2020 variants.
The width of the distributions is related to the level of disparity between model and measurements. For the BON2014 model, the fraction of R D values that fell between ±15% was 0.73, while for BON2020 ACE , BON2020 ACE-O , and BON2020 SSN , these fractions were increased to 0.95, 0.95, and 0.83, respectively. Phrased differently, the 95% confidence level (CL) of the BON2020 ACE and BON2020 ACE-O relative errors is ±15%.  Specific ions and energy groups account for large fractions of the exposure received behind shielding in space. It is therefore instructive to consider model errors in this way. Figure 11 shows the average model relative differences (R D ) separated into energy groups shown along the horizontal axes and charge groups indicated by different colored symbols. The error bars in the plot reflect the interval associated with measurement uncertainty. The figure shows negative relative errors for BON2014 for Z = 1 and Z = 2 in the three highest energy bins which have been shown to be significant for space radiation protection applications. The BON2020 variants all show minimal error for these same ions and energy groups.
Figures 12-14 provide more detailed information regarding the time-dependent behavior of the model for protons and helium. Figure 14 is of particular interest since the Cosmic Ray Telescope for the Effects of Radiation (CRaTER) data (Zeitlin et al., 2019) were not used in the calibration process and therefore provide an opportunity for independent validation.
The BON2020 ACE and BON2020 ACE-O results in these figures are identical since the proton and helium solar modulation functions in both versions of the model are obtained from the oxygen solar modulation function (i.e., ACE/CRIS oxygen data). It can be seen that these models track very well with the measurements, suggesting that the short term heliospheric perturbations reflected in ACE/CRIS oxygen data are a suitable proxy for protons and helium. The BON2020 SSN model lacks the same level of detail since it is driven by SSN, but it still shows an overall improvement compared to the BON2014 results as the systematic underestimation of the proton flux in solar cycle 23 (Figures 12 (left) and 14) and helium flux (Figure 12 [right]) by BON2014 have been corrected.
A quantitative assessment of the results shown in these figures is given in Table 3. The r 2 values for BON2020 ACE and BON2020 ACE-O are identical, as expected, and are all greater than 0.92, indicating a high degree of correlation between model results and measurement data. The R |D| values are also given in the table to highlight that not only is the updated model more closely correlated with time-dependent  measurement data but also provides flux values that are closer in overall magnitude to measurement data than BON2014. Figure 15 shows the average model relative differences compared to ACE/CRIS data for Z = 5-28 as a function of time. It can be seen that again BON2020 ACE and BON2020 ACE-O yield almost identical results. This is a significant result since it establishes that the linear regression coefficients used to relate integral ACE/CRIS oxygen flux to all other ions are a suitable proxy for describing solar activity for all ions between Z = 5 and Z = 28. Also evident in the plot is that the method outlined in the previous section used to generate an optimal ϕ allows the model to accurately reproduce the low energy ACE/CRIS data throughout the solar cycle. BON2020 SSN lacks the same level of fidelity but is still an improvement over BON2014, especially near the 2000 solar maximum.
Finally, in Figure 16, the LIS from BON2014 and BON2020 (all variants of BON2020 use the same LIS) are compared to the models of Corti et al. (2016) and HELMOD (Boschini et al., 2018). These additional models were developed specifically to represent the shape and magnitude of the LIS and utilize over 10 free parameters each to represent available data. The BON LIS, on the other hand, is a consequence of only three parameters that are calibrated only by near Earth measurements. The comparison in Figure 16 therefore provides  an indirect test of the physics included in BON describing solar modulation and heliospheric transport and also provides independent verification for the specified LIS parameters.
When compared to HELMOD or Corti et al., the overall BON2020 LIS spectral shapes show a significant improvement compared to BON2014. The BON2020 proton LIS shows a more realistic rise, rollover, and high energy spectral slope. The BON2020 helium LIS corrects the systematic underestimate at high energies that was present in BON2014 and is in excellent agreement with the other models across all energies up to 2 × 10 5 MeV/n. Discrepancies at larger energies are of little consequence to space radiation protection applications since the flux of particles at these energies is exceedingly small.

Propagated Errors
Model uncertainties described in the previous section may be propagated into effective dose behind shielding (Slaba & Blattnig, 2014b), as shown in Figure 17. Results in the figure are for the 1 June 2009-1 June 2010 solar minimum. The methods used to compute effective dose are described elsewhere (Slaba et al., 2010a(Slaba et al., , 2010b(Slaba et al., , 2013. The green dashed line with symbols in each figure corresponds to the nominal model result. These are the values one would obtain if the mission dates were input into the GCR model to provide a boundary condition for the remainder of the analysis tools. The shaded region corresponds to the distribution of effective dose values attributable to the uncertainties outlined in the previous section. The standard deviation interval (1-σ) and 95% CL of the distributions are explicitly shown. For each model, the line appearing closest to the nominal result is the median of the distribution. Percent differences between the 1-σ and 95% CL value relative to the median are explicitly shown at depths of 20, 50, and 100 g/cm 2 .
In the case of the BON2020 variants, the nominal and median model results are nearly identical. This is consistent with the results of Figure 10 where it was shown that these models have negligible systematic bias. For BON2014, it is seen that the nominal model result is below the median, suggesting that the nominal model result is underestimating what would be expected from measurements (also consistent with Figure 10). The BON2020 ACE and BON2020 ACE-O error intervals are all smaller than those found for BON2014. Moderate reductions in uncertainty compared to BON2014 are also found for BON2020 SSN . As in all other results shown thus far, the BON2020 ACE and BON2020 ACE-O are almost indistinguishable and were found to be within 0.5% of each other across all depths.
A final comparison of the models is given in Figure 18, which shows effective dose as a function of time for 20 g/cm 2 aluminum shielding. The BON2020 ACE and BON2020 ACE-O results were within 3% of each other across all mission times. Prior to~1998, all of the BON2020 variants revert back to SSN to obtain ϕ and therefore yield identical results. After 1998, the BON2020 SSN values appear consistent with BON2020 ACE and BON2020 ACE-O . The average R |D| between BON2020 ACE and BON2020 SSN over all mission times was found Figure 14. Model integral flux compared to PAMELA measurements. The shaded region represents total measurement uncertainty. BON2020 ACE and BON2020 ACE-O are nearly overlapping.

10.1029/2020SW002456
Space Weather to be 2%, although R D values as high as ±24% were found. The results from BON2014 appear consistently lower than the others.

Summary
The BON2020 model has been fully described and its improved performance demonstrated. The updated model has been calibrated to available measurement data, including the high precision measurements from AMS-02 and PAMELA. Solar activity is quantified through a simple parametric function relating daily ACE/CRIS integral flux to the solar modulation potential for each ion. For time periods outside of the region covered by ACE/CRIS, the model utilizes empirical relationships to relate time-delayed SSN values to ACE/CRIS integral oxygen flux. The estimated oxygen flux is combined with linear regression coefficients to estimate integral fluxes for all other ions between Z = 5 and Z = 28. These estimated fluxes are then used to compute the solar modulation potential for each ion. In all cases, the oxygen solar modulation potential is used for Z = 1-4 ions.
Three variants of the BON2020 model were identified: • BON2020 ACE -The solar modulation potential is obtained for each ion between Z = 5 and Z = 28 from daily ACE/CRIS integral flux if

10.1029/2020SW002456
Space Weather available and from V2 SSN elsewhere. The Z = 8 solar modulation function is used for ions between Z = 1 and Z = 4. • BON2020 ACE-O -Same as BON2020 ACE , except that only daily ACE/CRIS integral oxygen flux is used to describe solar activity. The linear regression coefficients are used to estimate ACE/CRIS integral fluxes for all other ions. • BON2020 SSN -Same as BON2020 ACE , except that only V2 SSN values are used to obtain the solar modulation potential.
Along with BON2014, the BON2020 variants were compared to available measurements, and uncertainty was quantified. It was found that the BON2020 is a significant improvement over BON2014. The systematic tendency of BON2014 to underestimate measurements was corrected. The dispersion between model and measurement, expressed through the R |D| and χ 2 metrics, was reduced by more than a factor of two. Comparison of the models against time-resolved AMS-02, PAMELA, CRaTER, and ACE/CRIS data showed that the BON2020 is able to accurately describe short and long-term solar activity effects.
It is suggested that the BON2020 ACE-O variant be used for space radiation protection applications. It was clearly shown that the BON2020 ACE and BON2020 ACE-O models yield similar, if not identical, results in terms of uncertainty quantification and radiation protection quantities such as effective dose. The