© G. Charpentier et al., Published by EDP Sciences 2024

1 Introduction

Through the Global Exploration Roadmap (ISECG, 2018), space agencies have identified a number of challenges to expand human exploration on the surface of Mars, one of the most critical being radiation protection for astronauts on long-duration missions. The deep space radiative environment is composed of two main sources of radiation: Galactic Cosmic Rays (GCRs) and Solar Particle Events (SPEs). GCRs are mainly fully ionized atomic nuclei composed of 87% protons, 12% helium and 1% heavier ions (Simpson, 1983). They form a continuous and isotropic flux, anti-correlated with the solar cycle, over energy ranges from a few keV/A to hundreds of TeV/A. Particles generated during SPEs are ejected by the Sun and accelerated by solar flares or Coronal Mass Ejections (CMEs), and although their probability of appearance is correlated with the solar cycle, they are unpredictable. They are mainly composed of protons, with energies ranging from a few keV/A to a few GeV/A, although their flux values may be much higher than GCRs (Desai & Giacalone, 2016). The Martian atmosphere, as on Earth, is composed of several layers, a troposphere from 0 to 45 km altitude and a mesosphere up to 110 km altitude (for an average surface temperature of −63 °C). The surface atmospheric pressure depends partly on the altitude and ranges from 82 Pa at the top of Mons Olympus (22 km above the datum) to 1200 Pa within the Hellas Basin (−8 km to the datum) (Zhang et al., 2023). The Martian surface radiative environment has three components: the first being the secondary particle flux resulting from nuclear interactions between primary particles from GCRs or SPEs with the atmosphere, known as atmospheric shower (Auger et al., 1939). The second component consists of primary particles which did not collide with the atmosphere due to their low density. Such particles are recorded as secondary particles as soon as they hit the Martian surface. The third, known as albedo flux, are the particles emitted by the interaction between projectile nuclei and the Martian regolith, consisting mainly of neutrons and photons (Drake et al., 1988).

Surface flux particles, whether directly or indirectly ionizing, can then interact with habitats and the human body, ultimately depositing energy. To reduce uncertainties in dose calculations and therefore risks for long-term missions, accurate surface flux spectra are essential. In this respect, the Radiation Assessment Detector (RAD) instrument (Hassler et al., 2012) onboard the Mars Science Laboratory (MSL) Curiosity rover has been providing a number of surface flux measurements for about one solar cycle (Hassler et al., 2014; Köhler et al., 2014; Wimmer-Schweingruber et al., 2015; Ehresmann et al., 2014, 2017; Guo et al., 2021b; Martinez Sierra et al., 2023). RAD in-situ measurements constitute the only data available for developing and improving calculation codes modeling the Martian environment. Although numerous studies (Simonsen et al., 1990; Ehresmann et al., 2011) had been carried out prior to the landing of Curiosity in August 2012, Matthiä et al. (2016) first compared RAD data with particle spectra obtained with GEANT4 Monte Carlo code (Agostinelli et al., 2003; Allison et al., 2006). Other particle transport codes, such as PHITS (Niita et al., 2006) or HZETRN (Wilson et al., 1991a), have been benchmarked during a workshop attended by leading specialists in the field (Hassler et al., 2017), to validate and discuss these deterministic or probabilistic codes. A global validation of these codes (Slaba & Stoffle, 2017; Matthiä et al., 2017) emerged through comparison with RAD data, for both charged and neutral particle spectra from GCRs. Validation and methods for modeling and comparing SPE spectra (Guo et al., 2018) with RAD observations (Ehresmann et al., 2018) have also been performed.

Such general validations provide a good level of confidence in the particle transport codes, but improvements can then be made to address uncertainties and differences between models and measurements in certain specific energy domains. Indeed, since RAD does not directly provide spectra for neutral particles, convolutions between the incident spectrum and the scintillator detector response functions are necessary to obtain neutron and gamma surface particle spectra. To this end, Köhler et al. (2014) have carried out power-law inversion and full 16-bin inversion calculations to compute MSL-RAD neutron and gamma spectra, from 5 to 900 MeV. The inversion method produced fitting estimates from models such as PLANETOCOSMICS (Desorgher et al., 2005), with satisfactory accuracy. However, a deviation between the MSL-RAD measurements and the models can be observed above 100–200 MeV (MSL-RAD flux being higher than models flux). Radiation exposure can induce numerous biological effects for astronauts, such as carcinogenesis or acute effects (Durante & Cucinotta, 2011). The need to model the spectra of neutral particles as accurately as possible is especially important as such particles account for a significant fraction of exposure on the Martian surface, due to the albedo effect. To provide input for radiation protection sizing computation, a flexible tool for particle spectra assessment, accounting for the atmospheric condition under any exposure scenarios and mission parameters, turned out to be necessary for Centre National d’Etudes Spatiales (CNES). Originally adapted from the Radiation Atmospheric Model for SEE Simulation (RAMSEES) (Cintas et al., 2023),a model of Earth’s atmospheric environment, the Atmospheric RAdiation Model for Ionizing spectra on martian Surface (ARAMIS) was developed as part of this study, providing control over particle physics and enabling all particle-specific information to be accessed for database construction.

2 Materials and methods

2.1 Global modelling strategy

As ARAMIS aims to provide input spectra for radiation protection calculations on the Martian surface, it has been designed to perform calculations in different scenarios, without having to repeat all the simulation steps. The approach adopted is to establish a surface flux response function that can be applied to any GCR or SPE input spectrum, a method which has been previously investigated and proven effective in reducing computation time and complexity (Guo et al., 2018; Banjac et al., 2019; Guo et al., 2019a). To this end, a distinction has been drawn between two types of input parameters: those describing the mission scenario and those relating to the exposure scenarios. Parameters relating to the mission scenario (dates, altitude and landing site coordinates) are fixed within the Monte Carlo code. Exposure parameters, on the other hand, are included in post-Monte Carlo codes, using GCR models or the addition of SPEs with varying intensities.

Monte Carlo calculations (using the GEANT4 toolbox developed by CERN) are not conducted with a particular input flux, but with a set of monoenergetic particles launches as described in Table 1 and Figure 1. These results are then post-processed using Python 3.8.4 in order to be assembled and then convolved with any input spectrum. The main advantage of this method is to reduce GEANT4 computation times while ensuring a high level of confidence in the reliability of the results. Indeed, since GEANT4 is a Monte Carlo statistical code, obtaining a good statistical result depends on the number of primary particles launched, which can have a major impact on calculation time. Running ARAMIS represents a very substantial number of calculation hours for a single computer, regardless of how powerful it may be. Therefore, the CNES High-Performance Computing (HPC) cluster has been used as it gives access to high computing power, a large number of processors, a large amount of memory (12,288 CPUs: AMD processors divided into 128 cores per server, 1.1 TB of RAM, 3 GHz per core) as well as a large storage volume and high-performance disk space.

Figure 1 Schematic of the global ARAMIS organisation and structure.

2.2 Primary GCR or SPE models

Various GCR models are used by space agencies and scientists, having been established either empirically, such as the DLR model (German space agency) (Matthiä et al., 2013) and the CREME model (Nymmik et al., 1996; Adams et al., 2012), or numerically, such as HelMod (Bobik et al., 2013; Boschini et al., 2018), SDEMMA (Luo et al., 2019; Chen et al., 2023) or BON20 (Slaba & Whitman, 2020), the latter developed by NASA. Liu et al. (2024) have shown that among these models, BON20 appears to be considerably robust for most low-energy elements and can predict with high accuracy the fluxes of the two species contributing most to radiation dose, H and He, at high energies. Therefore, the model used for the results presented in this document is BON20, over 2 time periods: from 19 August 2012 to 17 February 2013 and from 15 November 2015 to 15 January 2016 with respective averaged solar modulation potential $\bar{\mathrm{\Phi }}$ of 599 and 590 MV. The planetary protection from GCRs provided by Mars is not initially included in free-space GCRs models. Therefore, only the flux over half the solid angle was considered to obtain the Martian boundary environment, as shown in Figure 2.

Figure 2 Free space GCRs spectrum in Mars surrounding from 19 August 2012 to 17 February 2013 (corresponding to $\bar{\mathrm{\Phi }}=599 \mathrm{MV}$). The BON20 model is shown by solid lines while DLR model is represented by dotted lines. Due to different approaches (numerical for BON20, empirical for DLR) and input parameters (Liu et al., 2024), those two models differ significantly at low energies for most particles (for energies below 100-200 MeV). As detailed in 3.3.3, most particles of these energies do not reach the surface, as the Martian atmospheric cut-off is about 160 MeV for protons, for instance. Differences also exist at higher energies for heavier ions with lower GCR abundances compared to H and He. Plotted with data extracted from On-Line Tool for the Assessment of Radiation in Space (OLTARIS).

ARAMIS is also able to consider primary SPE spectra as an input parameter. From the major events shown in Figure 3 and Table 2, the October 2003, October 1989, August 1972, November 1960 and February 1956 events have been considered for the results presented in this study. As SPE are very sudden events lasting only a few days, the most appropriate quantity to describe them is fluence. These fluences have been conservatively considered at r = 1 AU, r being the distance to the sun, in accordance with ECSS E-ST-10-04, even though Mars is located at r = 1.53 AU, which would lead to an attenuation of a factor of 2.3 if a reduction ratio of 1/r² had been applied. Moreover, SPEs being both directional and unpredictable, it is impossible to know whether it would hit Mars on the side faced by the mission location, or on the other side, with the planet acting as a shield in the latter case. Therefore, fluences have been fully taken into account without any mitigation factor. In addition to distance from the sun, another parameter makes it difficult to describe SPE flux at Mars from Earth measurements, namely the difference in longitude between the two planets. It has been shown that a longitudinal difference influences the intensity of event measurement (Richardson et al., 2014), so the longitude gap between Mars and Earth is likely to result in very different SPE flux patterns for the two planets.

Figure 3 Different SPE fluences in free space at 1 AU, corresponding to the historical major events recorded or reconstructed. Plotted with data extracted from OLTARIS and OMERE (for October 2003 data).

2.3 Physics used in Monte Carlo

The interactions of cosmic rays with the Martian atmosphere have been modelised in ARAMIS using the GEANT4 toolkit version 11.1.0. This Monte Carlo code integrates various advanced physics of particle and nuclear interactions. The so-called physics lists used for ARAMIS are detailed in Table 3, alongside their hadronic processes. The various physics lists chosen are derived from FTFP_BERT and QGSP_BERT, the two physics lists recommended for “cosmic ray applications”.¹ Each of these physics lists uses the HP option: a High Precision Neutron Model which offers greater accuracy in describing elastic and inelastic scattering, capture and fission for low-energy neutrons (<20 MeV). These physics are described by different interaction models depending on the energy range. Fritiof (FTF) (Andersson et al., 1987) and Quark Gluon String (QGS) (Allison et al., 2016) models are used for high-energy particles (≥5 GeV and ≥20 GeV respectively). For lower energies, various cascade models describe atmospheric showers: the Bertini Cascade model (BERT) (Bertini & Guthrie, 1971), the Intra-Nuclear Cascade model of Liege (INCL++) (Mancusi et al., 2014) and the Binary Cascade model (BIC). All standard electromagnetic processes are included by default in these physics lists. Elastic scattering is modelled in the same way in all four configurations: using G4ChipsElasticModel from 0 to 100 TeV for protons, and identically for neutrons, adding NeutronHPElastic from 0 eV to 20 MeV.

2.4 Atmospheric design

Geometric modelling of the Martian atmosphere was carried out using a 30 km square column base, divided into a parametric number of layers representing the density distribution defined by Mars Climate Database (MCD) version 6.1 (Forget et al., 1999; Millour et al., 2019), according to chosen altitude, coordinates and date. The geometry stack rises to 90 km above the Mars Orbiter Laser Altimeter (MOLA) reference level (MOLA₀). As density decreases exponentially with altitude, it has been decided to allow the user to select the number of layers and their thickness in the lower atmosphere, in order to maximise accuracy between the surface and 40 km above MOLA₀. In the results presented in this article, the layers are 5 km thick between MOLA₀ and 40 km, then 10 km above. As surface altitude can vary between roughly −8 km and 21 km above MOLA₀, the thickness of the first atmospheric layer depends on the chosen landing site. The composition of the atmosphere consists of 95.482% CO₂, 2.705% N₂, 1.603% Ar, 0.13% O₂ and 0.08% CO as given by OLTARIS (Singleterry Jr et al., 2011). Beneath the surface, a 5 km-thick layer of regolith has been added to take account of the albedo particles, with its specific characteristics shown in Table 4 and depending on the landing site. The Sol-averaged surface pressure used is respectively 848 Pa and 826 Pa for the Curiosity (−4.6°N, 137.4°E) and Phoenix (68.2°N, 234.3°E) landing sites over the time period from 19 August 2012 to 17 February 2013, while it is 878 Pa for the Curiosity landing site over the second period considered, between 15 November 2015 and 15 January 2016. Indeed, RAD measurements have made it possible to demonstrate the relation between GCR-induced dose and pressure variations (Rafkin et al., 2014). In particular, the surface dose was found to be anti-correlated with surface pressure for solar modulation potentials below 900–1000 MV (Guo et al., 2017).

2.5 Postprocessing methodology

As illustrated in Figure 1, ARAMIS has been built in such a way that each distributionally calculated GEANT4 simulation requires reconstruction to recover the surface flux. First, the results are collated in databases containing all the specific features of each particle detected: (Seed, RunId, ParticleId, IdParent, Type, Energy, Position, Momentum and Excitation energy). With these individual particle counts, databases of surface spectra, dN_ij/dT*, by primary particle type i of energy T and secondary particle type j of energy T*, have been constructed using Python codes. From these secondary spectra on each energy bin, relative error (White, 2010), R_kj, was calculated:

$$R_{\mathrm{kj}}={\left[\frac{\sum _{p=1}^N\mathrm{}{T_p^{\mathrm{*}}}^2}{{\left(\sum _{p=1}^N\mathrm{}{T}_p^{\mathrm{*}}\right)}^2}-\frac{1}{N}\right]}^{1/2},$$(1)with N the total number of secondary particles type j in the energy bin k, and $T_p^{\mathrm{*}}$ the energy of a p particle. Any primary GCRs or SPEs spectrum, dϕ_i/dT, has been convolved with the spectra of secondary particles per energy bin to obtain the surface flux by secondary type, dN_ij/dT* (T, z₀), following the method used by Cintas et al. (2023). Integration by energy bin was performed using trapezoidal integration method:

$${\frac{\mathrm{d}{\varphi }_j}{\mathrm{d}{T}^{\mathrm{*}}}|}_{z_0}=2\pi \sum _i\mathrm{}{\int }_T\mathrm{}\frac{\mathrm{d}{N}_{\mathrm{ij}}}{\mathrm{d}{T}^{\mathrm{*}}}\left(T,z_0\right)·\frac{\mathrm{d}{\varphi }_i}{\mathrm{d}T}\mathrm{d}T,$$(2)applying a 2π factor to account for incident particles from all azimuth φ and zenith θ angles over the entire solid angle. With this data structure, complementary spectra enabling a more accurate analysis of the ARAMIS results were constructed in databases by primary Z atomic number, primary energy, secondary particle type and upward (albedo) or downward (atmospheric) origin.

3 Results and discussion

3.1 Validation and comparison of ARAMIS with existing models and data

In order to validate and discuss the different models, statistical performance indicators have been used, namely the χ² and absolute relative difference |RD|, the latter being adapted from Liu et al. (2024):

$${\mathrm{\chi }}^2=\sum _{k=0}^N\mathrm{}\frac{{\left[F_{j,\mathrm{model}}\left(T_k^{\mathrm{*}}\right)-F_{j,\mathrm{RAD}}\left(T_k^{\mathrm{*}}\right)\right]}^2}{{\sigma }_k^2},$$(3)

$$\left|\mathrm{RD}\right|=\frac{1}{N}\sum _{k=0}^N\mathrm{}\frac{\left|F_{j,\mathrm{model}}\left(T_k^{\mathrm{*}}\right)-F_{j,\mathrm{RAD}}\left(T_k^{\mathrm{*}}\right)\right|}{F_{j,\mathrm{RAD}}\left(T_k^{\mathrm{*}}\right)}\times 100\%,$$(4)with all the χ² and |RD| values displayed in Appendix A. These parameters have been calculated using the flux $F_j={d{\varphi }_j/d{T}^{\mathrm{*}}|}_{z_0}$ for a secondary particle type j for both model and RAD measurements, with N being the number of energy bins, $T_k^{\mathrm{*}}$ the energy and σ_k the uncertainty on RAD data for the kth bin.

3.1.1 Surface neutron and photon spectra with GCR exposition

In this section, the results provided by ARAMIS are presented and compared with experimental measurements from the Curiosity rover and with other existing models (GEANT4 and PHITS simulation from Matthiä et al. (2016) OLTARIS results using the HZETRN code). Figure 4 shows the secondary neutron and photon surface fluxes for the FTFP_INCLXX_HP physics list. With regard to neutrons, ARAMIS has quite decent statistical indicators, although not the strongest, since it is the GEANT4-Matthia model that has the best fit with the MSL data over the overall and <150 MeV spectra parts (see Table A.1.a). One should also note that the OLTARIS model shows much poorer statistical indicators than the other three models for all spectral categories. Moreover, we observed in ARAMIS results for physics lists using INCL++, that fluxes are very slightly underestimated at low energies (T^* < 2 MeV) while physics lists using the BERTINI cascade seem to overestimate fluxes between 10 and 100 MeV. Furthermore, ARAMIS predicts a significantly increased flux between 200 MeV and 4 GeV. This also corresponds to an inflexion of results between the model and experimental measurements, since for energies above 200 MeV, other models have less accurate fitting parameters (both χ² and |RD|) with the data of MSL-RAD on this spectral domain. It must be noted that the χ² parameter is smaller at high energy for all models in comparison with other spectral domains due to bigger uncertainties on the MSL-RAD measurements. ARAMIS therefore slightly closes the gap between simulation and experiment. To analyze this discrepancy in results, the value of the relative error on the Monte Carlo calculation, R_kj, provides some insight. Indeed, through the architecture of ARAMIS, it is possible to differentiate fluxes and particles at each stage of the reconstruction according to their ascending or descending direction. The statistics provided in Figures 5 and 6 allow very good confidence in the estimates of ARAMIS since the relative error is less or equal to 5% for nearly all the domains over which secondary neutrons are generated. Next, it can be seen that albedo neutrons and “atmospheric” neutrons do not have the same domains of predominance on the spectrum, since albedo neutrons are present between 10⁻⁹ and 10² MeV, while atmospheric neutrons predominate between 10⁻² and 10⁵ MeV. Therefore, the difference observed between ARAMIS and other models above 200 MeV can be attributed to the atmospheric neutron. Indeed, it has been shown that the high-energy neutron flux is dependent on atmospheric conditions, since it depends directly on surface pressure (or column depth), and that when surface pressure increases, the high-energy neutron flux increases (Zhang et al., 2022). As mentioned above, the MCD model induces a surface pressure of 848 Pa, which is in the high range of the conditions tested by Zhang et al. (2022) (between 82 Pa and 1200 Pa), hence the increase in high-energy neutron flux. Moreover, as detailed in the dedicated section, ARAMIS allows parameterization of the number of atmospheric layers and their thickness to better represent the evolution of atmospheric density with altitude. Lastly, it can be seen that beyond its influence on the proportion of each species in the incident GCR fluxes, the atomic number of the incident particle influences the Monte Carlo statistics, as the heavier the ion, the more energized it has to be in order for the atmospheric shower to reach the ground. One explanation is that larger nuclei are broken up and shattered higher up in the atmosphere, leaving secondary particles to be dissipated and absorbed by the medium.

Figure 4 Martian surface flux of neutrons (left) and photons (right) at Gale Crater from 19 August 2012 to 17 February 2013, given by ARAMIS and compared to measurements and other models with FTFP_INCLXX_HP. Similar plots with other physic lists are displayed in annex.

Figure 5 Relative error per energy range in secondary spectra calculation for downward neutrons, for 105 launches using FTFP_INCLXX_HP and different primary particles: Proton (a), Helium (b), Carbon (c), Oxygen (d), Silicon (e) and Iron (f).

Figure 6 Relative error per energy range in secondary spectra calculation for upward neutrons, for 105 launches using FTFP_INCLXX_HP and different primary particles: Proton (a), Helium (b), Carbon (c), Oxygen (d), Silicon (e) and Iron (f).

Figure 7 Martian surface flux of proton, electron & positron, muon and pion at Gale Crater from 19 August 2012 to 17 February 2013, given by ARAMIS and compared to measurements and other models with FTFP_INCLXX_HP. For protons, solid lines represent cone-viewed downward spectra, while dashed lines represent direction-averaged spectra.

Figure 8 Effective dose rate Yield function Yi(T, z0) (left) and relative error on dose rate Yield function calculation δi(T, z0) (right) at Gale Crater between 8 August 2012 to 30 January 2013, for different primary particle types and energies, with FTFP_INCLXX_HP.

Figure 9 Martian surface neutron fluxes at Gale Crater (Curiosity landing site) and Vastitas Borealis (Phoenix landing site) from 19 August 2012 to 17 February 2013, given by ARAMIS with FTFP_INCLXX_HP. Neutron fluxes are multiplied by energy to enhance peak readability.

Figure 10 Evolution of the CO2 column depth (in g · cm−2) with the solar longitude, Ls, for several coordinates corresponding to the landing sites of different probes and rovers (Curiosity, Perseverance, Phoenix, Viking 2). Plotted with data from the MCD v6.1. The central point value corresponds to 15 g · cm−2, hence the curve shapes.

Figure 11 Martian surface flux of neutrons (left) and photons (right) at Gale Crater from 15 November 2015 to 15 January 2016, given by ARAMIS and compared to measurements and other models with FTFP_INCLXX_HP.

Figure 12 Martian surface fluence of neutrons (left) and photons (right) at Gale Crater from 19 August 2012 to 17 February 2013, given by ARAMIS with FTFP_INCLXX_HP and OLTARIS for different SPE historical events. Curves have been respectively scaled by factors of 10 and 100 for the November 1960 and August 1972 events for sake of clarity. The OLTARIS October 2003 event is not displayed as this event does not exist in the OLTARIS database and exporting OMERE fits to OLTARIS would imply using the fit function for only T < 300 MeV, as it is a combination of fits on the overall spectrum.

For photons, similar overall observations can be made, with a better fitting of ARAMIS flux with MSL data at high energies in terms of χ² or |RD| but also on the overall spectrum with respect to |RD|, the PHITS model performing better at energies <150 MeV. Thus, ARAMIS offers much better agreement with experimental measurements over almost the entire photon energy spectrum. These conclusions are observed for all physics lists tested, although physics lists using the BERTINI cascade seem to slightly overestimate fluxes between 1 and 10 MeV. As for neutrons, the error in the Monte Carlo calculation is very small, with a validity range of between 10⁻² and 10⁵ MeV (see Figures B.1 and B.2). Looking at the whole picture, whether for neutrons or photons, physics lists using INCL++ seem better suited for modelling the atmospheric environment of Mars. Indeed, as shown in Table 3, physics using INCL++ handles reactions between 6 GeV and 15 GeV, unlike physics using BERTINI (where this range is handled by the FTF/QGS part of the model). Furthermore, INCL++ has been shown to be particularly effective for modelling neutron production by spallation, as pointed out by Guo et al. (2019a), who reached the same conclusion regarding the choice of physics lists. The atmosphere of Mars being particularly conducive to neutron production, INCL++ physics then provides a more efficient model. Hence the results using FTFP_INCLXX_HP and QGSP_INCLXX_HP are quite similar and in the following, the analysis will be performed for the physics list FTFP_INCLXX_HP.

3.1.2 Other particle spectra

In the previous section, focus has been given to neutrons and photons, but secondary fluxes are determined for any type of elementary particle (mesons, baryons, leptons). ARAMIS also provides surface spectra of secondary heavy ions and associated isotopes, which are not presented in this article but may be the subject of a future publication. Figure 7 presents the surface fluxes for proton, e⁺, e⁻, π⁺, π⁻, μ⁺ and μ⁻, using identical mission and exposition parameters as before. Whereas RAD-reconstructed neutral particle measurements are averaged in all directions, charged particle measurements follow a conical view of the sky, with an opening angle of 60° (Ehresmann et al., 2014). Using the ARAMIS parameterized data structure, separate databases can be reconstructed for downward particles entering a 60° cone of view, so that RAD measurements can be compared with their simulation equivalent. In this way, models can be compared with each other for downward proton spectra seen from a cone, which provides a valuable indication of model validity for their spectra computation in all directions. A much better agreement of ARAMIS with the data can then be observed both visually and statistically (see Table A.1.c) for cone-viewed downward protons over the overall energy range of the RAD instrument. This spectrum is less accurate at lower energies (<20 MeV), where the convergence of results is slightly poorer than for the rest of the spectrum. Improving this convergence is difficult to achieve in terms of computation time. One solution would be to consider only H, He, Li, C, O, Fe and Si ions (representing about 98% of GCR abundance) as primary particles and then increase the number of particles launched. Indeed, such a solution would reduce ARAMIS surface fluxes by around 1% (in the case of the test parameters used in this study) and reduce computation times by 80%. This is a possible approach for the future, if its validity were to be verified for all atmospheric and GCR modulation conditions. The e⁻/e⁺ fluxes then show very good agreement with those of OLTARIS, with very slightly higher values, that can once more be attributed to the parametric choice of atmospheric modelling, which for a higher number of layers gives a better approximation of density evolution with altitude. It must be noted that a method enabling to obtain electron and positron counts measured by RAD with an inversion approach has been established (Köhler et al., 2016). These count rates have then been used to validate integrated fluxes obtained by PLANETOCOSMICS simulations of the RAD instrument, see Figure 7of Köhler et al. (2016). This enables a first indirect comparison to be performed between the e⁻/e⁺ fluxes obtained with the RAD measurement and the ARAMIS/OLTARIS models, as shown in Table 5. It emerges that over the RAD range 4–20 MeV, the ARAMIS integrated fluxes are the most accurate relative to the RAD measurements, while the fluxes obtained with OLTARIS slightly underestimate the measurements by a factor of 2, although remaining within the same order of magnitude. However, it is important to be cautious, as these values are integrated fluxes and not differentials, and no more precise comparative analysis is possible at this stage. Concerning pions and muons, no measurements are currently available to validate the codes, although a comparative analysis between OLTARIS and ARAMIS is worthwhile. Unfortunately, pion and muon fluxes are not in good agreement between the ARAMIS and OLTARIS models, despite being closer above 10 GeV for pions. The underestimation of pion fluxes by HZETRN can be explained by the non-inclusion of kaons (Slaba et al., 2020), which can decay into pions from high-energy nuclear reactions. Since muons result from charged pion decay, a similar conclusion can be made regarding the underestimate of HZETRN muons spectra as compared with ARAMIS.

3.2 Yield effective dose rate

Beyond presenting surface fluxes by secondary particle type, ARAMIS can also be used to perform an advance dosimetry analysis. For this, the effective dose rate yield function, Y_i, used by Larsen & Mishev (2023) has been implemented. This function (different from the effective dose rate) provides the dose distribution on the primary spectrum, at surface altitude z₀. It is calculated using the flux, F_ij, of secondary particles j (neutron, proton, γ, e⁺, e⁻, π⁺, π⁻, μ⁺, μ⁻) with energy T^* and primary particles i with energy T, and with the flux-to-effective-dose conversion coefficients, C_j(T^*), given by ICRP Publication 116 (Petoussi-Henss et al., 2010):

$$\{\begin{array}{l}{Y}_i(T,z_0)=\sum _j\mathrm{}{\int }_{T^{\mathrm{*}}}\mathrm{}{F}_{\mathrm{ij}}(T,T^{\mathrm{*}},z_0)C_j\left(T^{\mathrm{*}}\right)\mathrm{d}{T}^{\mathrm{*}},\\ F_{\mathrm{ij}}(T,T^{\mathrm{*}},z_0)=2\pi \left[\frac{\mathrm{d}{N}_{\mathrm{ij}}}{\mathrm{d}{T}^{\mathrm{*}}}(T,T^{\mathrm{*}},z_0)·\frac{\mathrm{d}{\varphi }_i}{\mathrm{d}T}\right].\end{array}$$(5)

Yield function matrices can only be constructed as a result of the architecture of ARAMIS, which, with its intermediary databases, enables such data to be accessed. In order to validate those dose rate values, the relative error of the Yield function, δ_i, is also determined by propagating the relative errors of the Monte Carlo calculation R_kj:

$${\delta }_i(T,z_0)=\frac{{\left[\sum _j\mathrm{}\sum _k\mathrm{}{\left(F_{\mathrm{ij}}{R}_{\mathrm{kj}}\right)}^2{C}_j(T^{\mathrm{*}}{)}^2\right]}^{1/2}}{Y_i(T,z_0)}.$$(6)

Figure 8 shows Y_i and δ_i matrices, from which several conclusions can be drawn. The main contributors to the dose are protons and helium, in line with the abundance of the GCR spectrum, although other ions are not negligible. Moreover, for all species of the primary spectrum, energies of about 10 GeV contribute most to the dose. The uncertainty of the Monte Carlo calculation then propagates in such a way that the relative error on the Yield function is less than 2% for almost the entire spectrum, but for the lowest energies, the delta value exceeds 15%, a rather significant level. This domain also corresponds to very low Yield dose rate values, which may be attributed to the low energies at which atmospheric particle showers do not produce enough secondary particles to make a significant contribution to the dose. This could be addressed by significantly increasing the number of incident particles, although this would result in unrealistically long computation times – for comparison, simulations with 2.52 × 10⁷ incident particles required 3 weeks per configuration, with CNES HPC performances. Furthermore, as the Yield dose rate is very low in this particular area, it is acceptable to have a slightly higher relative error there, without altering the overall Yield dose rate, which therefore benefits from a very high level of confidence for these results.

3.3 Influence of exposition and mission scenario

3.3.1 Landing site choice influence

Landing site selection has a significant influence on the surface flux characterization. Indeed, areas close to the poles such as Vastitas Borealis (landing site of the Phoenix probe) are estimated to contain up to 50% water in the soil composition (Table 4). Röstel et al. (2020) have shown that the presence of hydrogen in the form of water in soils (andesite in this case) can modify the secondary neutron spectrum by decreasing it by up to 25% in the <10 MeV spectral range for a composition of 50% water (or a 5% reduction for a 10% water content). ARAMIS underlines this point by showing a reduction in secondary neutron fluxes (up to a factor of 2.5) in the <100 MeV energy range, as illustrated in Figure 9. This reduction is entirely due to a lower albedo effect: as mentioned in section 3.1.1, albedo neutrons dominate the spectrum up to 200 MeV, after which downward neutrons predominate. Indeed, as detailed by Goldhagen et al. (2004) and highlighted by Matthiä et al. (2017), neutrons generated by GCRs show two peaks in the spectrum, the first being located around 1 MeV and corresponding to neutrons evaporating from excited nuclei, while the second sitting around 100 MeV corresponds to neutrons created by spallation processes with primary and atmospheric nuclei. The results obtained with ARAMIS slightly shift those peaks to 2.3 MeV and 240 MeV, as shown in Figure 9. Indeed, these peaks may differ slightly depending on the modeling used, see Figure 2 of Matthiä et al. (2017). It should be noted that the first peak is not visible on RAD measurements (energy below the range of the instrument) but that a very slight peak forms around 250–300 MeV, in a zone of higher measurement uncertainty. As the density of Martian regolith is notably higher than the atmospheric density, the number of excited nuclei is also much higher in the soil, which explains the predominance of low-energy albedo neutrons between 1 and 10 MeV in ARAMIS results. Therefore, ARAMIS spectrum suggests that for similar altitudes (−4488 m for Curiosity and −4125 m for Phoenix), surface pressure (848 Pa for Curiosity and 826 Pa for Phoenix) and dates, the effect of atmospheric variations is negligible compared with variations in soil composition and density. Moreover, this finding supports the Martian poles (for non-extreme altitudes) as a suitable destination for future human missions, from a radiation protection perspective, although a compromise must be made to find the best landing coordinates. Indeed, a higher water content reduces the albedo neutron spectrum, while, as pointed out by Zhang et al. (2022), certain extreme altitudes would induce an increase in the downward neutron spectrum due to higher surface pressure.

3.3.2 Seasonal and solar cycle variation

Fluxes observed by MSL-RAD or determined with ARAMIS tend to vary over time. The two main factors influencing these variations are the solar cycle, since the primary GCRs are anti-correlated with it, and the solar longitude L_s, which represents the seasonal variations of a Martian year. As shown in Figure 10, the evolution of CO₂ column depth differs from one point to another, and as a function of L_s varies by as much as a factor of two (in the case of the Viking 2 lander coordinates, for instance). Furthermore, the variation in intensity of GCR spectra depends on the energy range: as an illustration, the flux of 100 MeV protons is multiplied by 10 at solar minimum compared with solar maximum, while the flux of 10 GeV protons hardly varies at all (Slaba & Whitman, 2020).

It is therefore appropriate to validate the spectra calculated by ARAMIS over a period of solar minimum, between November 15, 2015, and January 16, 2016, shown in Figure 11. As for previous conditions, ARAMIS follows the pattern of RAD measurements at higher energies for photons (and therefore differs from the other models) but not for neutrons. It must be pointed out that for this configuration, the χ² and |RD| values indicate a significant overestimation of neutron fluxes by OLTARIS and photon fluxes by the GEANT4-Matthia model. Since the latter uses an old and obsolete version of the MCD (v4.3), such differences can be generated by modifying atmospheric conditions at the desired dates. Between MCD v4.3 and v5.2, a major update was carried out with the implementation of a new Global Climate Model (GCM). It included updates in mean atmospheric composition to match measurement made by Curiosity SAM (Sample Analysis at Mars), as well as improved CO₂ cycling, convective boundary layer and ionosphere. Later on, from v5.2 to v6.1, some corrections in density extrapolation were made with the implementation of a high-resolution mode. Therefore, use of version v6.1 of the MCD seems to be the most appropriate modelling option, although the statistical indicators of PHITS (using v5.2 of the MCD) are equally reliable, especially at low energies. Even though the spectrum shapes are similar, the flux appears slightly lower than in Figure 4, as the CO₂ column depth of Gale Crater (Curiosity landing site) barely varies around its mean value with L_s (see Fig. 10), which explains the marginal decrease in flux, with the primary GCR flux becoming the only time variable (with a solar modulation of 590 MV in this configuration, compared with 599 MV in the previous scenario). So, in the case of a point where the variation in mass density is significant over the Martian year, such as Phoenix landing site, significant flux variations in the ARAMIS spectrum can be expected.

3.3.3 Solar event impact

Regardless of the mission parameters, events influencing exposure, such as SPEs, may occur. ARAMIS therefore supplies particle spectra from the Martian surface for such events. Figure 12 represents the secondary neutron and photon fluences at the landing site of Curiosity, taking into account the atmospheric state for the same conditions as in Section 3.1.1, for different historical events.

At low energies (below 70 MeV for neutrons and below 15 MeV for photons) the August 1972 event presents extremely strong photon and neutron surface spectra, even though it is only stronger than the others events in the primary [30 MeV, 150 MeV] range. Indeed, low-energy primary protons (below 100 MeV) experience difficulties to be transported through the atmosphere and generate no secondary spectrum, which has been verified by ARAMIS simulations on the low energy range. This finding further validates an observation made by Guo et al. (2019a), which find that at Gale Crater the atmosphere stopped primary protons below the energy of about 160 MeV for an average pressure of 830 Pa. According to the ARAMIS results, the atmospheric cut-off is therefore slightly lower, around 100 MeV. Due to this effect, the October 1989 event does not provide the highest surface fluences, even though it has a stronger primary spectrum than any other event, for energies <20 MeV (see Fig. 3). At very high energies (above 150 MeV), the November 1960 event has the highest flux, dominating the primary flux above 200 MeV. However, this same event also shows a relatively elevated surface flux in the lower energy ranges for neutrons and photons, although it does not predominate over those ranges of the primary spectrum. This is explained by the tendency of high-energy primary particles to generate more low-energy secondary particles through spallation reactions in the atmosphere. Indeed, it has been previously shown that Martian atmosphere shielding has a major influence on primary particles of low energies (Zhang et al., 2023). An approach to approximate the surface absorbed dose on Mars generated by SPE was previously found by using a pivot energy of the solar event spectrum before entering the Martian atmosphere (Guo et al., 2019a). The method, using a pivot energy of 300 MeV, is valid for large primary spectra where the proton flux extends to energies above 500 MeV. It defines the surface absorbed dose in mGy · h⁻¹ as D = 4.45 · I₃₀₀ (±9.8%), where I₃₀₀ is the SEP intensity at 300 MeV in s⁻¹ · sr⁻¹ · cm⁻² · MeV⁻¹ and 4.45 the conversion coefficient in mGy · h⁻¹/(s⁻¹ · sr⁻¹ · cm⁻² · MeV⁻¹). An extended method, including surface pressure in the calculation of the pivot energy based on more realistic SEP spectra which are double power-law functions fitted to the actual measurement was also developed by Zhang et al. (2023). The latter method defines the absorbed surface dose, D_surf(p), based on the following formula:

$$D_{\mathrm{surf}}\left(p\right) \left[\mathrm{mGy}\right]=c_D\left(p\right)·F_{E_0\left(p\right)},$$(7)where the fitting coefficient, c_D(p), and the pivot energy, E₀ – used to compute original GLE fluence, $F_{E_0\left(p\right)}$ – are defined based on the pressure:

$$\{\begin{array}{l}{E}_0\left(p\right)=-0.000114·p^2+0.350·p+86.9,\\ \\ \frac{c_D\left(p\right)}{10^{-7}}=-0.000614·p^2+1.08·p+388.\end{array}$$(8)

To enable comparison, the surface fluences, F_j, provided by ARAMIS and OLTARIS were then multiplied by the flux-to-skin absorbed dose conversion coefficients, c_skin(T*) given by ICRP Publication 116 (Petoussi-Henss et al., 2010), for male and female, to obtain the mean surface absorbed dose on the skin, D:

$$D \left[\mathrm{mGy}\right]=\sum _j\mathrm{}{\int }_{T^{\mathrm{*}}}\mathrm{}{F}_j\left(T^{\mathrm{*}}\right)c_{\mathrm{skin}}\left(T^{\mathrm{*}}\right)\mathrm{d}{T}^{\mathrm{*}}.$$(9)

For all events investigated in this paper, both pivot methods give relatively close results. Indeed, the method of Zhang et al. (2023) redefines the pivot energy as a function of surface pressure, and for this application (P = 848 Pa) the pivot energy (301.72 MeV) is relatively close to the one of Guo et al. (2019b). It should be pointed out that the method used by Zhang et al. (2023) brings the calculated doses closer to those observed with OLTARIS and ARAMIS, which both provide equivalent doses for all events with the exception of August 1972, for which OLTARIS overestimates the dose by a factor of 3 (compared with a factor of 2.1 for ARAMIS). Variations in surface dose between the different models can be explained by both differences in Martian transport processes and the phantoms used to compute the dose. Firstly, the different models do not apply the same high-energy cut-off energies, which can cause the dose to vary slightly. Hence, the pivot method high-energy cut-off value stands at 1 GeV, compared with 7.5 GeV for ARAMIS, while OLTARIS uses no cut-off but rather the full primary spectrum over its continuous range for analytical resolution. The cut-off energy of ARAMIS being relatively high, there is very little difference between ARAMIS and OLTARIS surface dose for all events except August 1972. For August 1972, such deviation can be explained with the event spectral shape, which is very high at low energy and drops drastically above 100 MeV. Therefore, for such an event using the pivot method can lead to a dose slightly underestimated. Similarly, a significant discrepancy between the models and the pivot dose is observed for November 1960. As shown in Figure 3 and Table 2, this event consists of a sum of 2 exponentials, with the pivot energy applied at ∼300 MeV along the second exponential, much more energetic than the first but smaller. This tends to underestimate the dose value relative to the actual spectral pattern at low energies. For combinations of events, the pivot method is therefore not the most appropriate, although it provides extremely accurate results for other spectral shapes. The second source for variation relates to the type of phantom used, where the doses calculated with OLTARIS and ARAMIS use ICRP absorbed skin dose coefficients, calculated using the Adult Reference Computational Phantoms of ICRP (2009), the pivot method used a 0.5 mm water slab (Guo et al., 2019b) and a 15 cm radius water sphere (Zhang et al., 2023), to approximate respectively the fine skin structures and the torso. Such disparities in phantom compositions explain for the remaining discrepancies, in addition to the accuracy of the transport process. In particular, the August 1972 event is very gentle, with numerous low-energy particles on the surface of Mars (as shown in Fig. 12) which are sensitive to phantom self-shielding, and therefore cause dose variations from one phantom to another. Finally, as many studies have already demonstrated (Guo et al., 2018), there is a need for modelisation accuracy (as well as forecasts) when it comes to the dose delivered by high-energy SPEs, as their effects can generate major risks for astronauts.

4 Summary and perspectives

The Monte Carlo simulation-based ARAMIS has been presented and validated against measurements made by the RAD instrument on board the Curiosity rover and against previous models. ARAMIS uses a response-function-based Monte Carlo simulation reconstruction to obtain the surface fluxes, a method which, as demonstrated, not only maintains the structure and consistency of the results but also allows the implementation of a multitude of exposure scenarios without the need to recompile the GEANT4 code. We observed that among the various physic lists tested, FTFP_INCLXX_HP has the strongest physics components to represent the Martian atmosphere, a conclusion also shared by Guo et al. (2019a), and that the latest version of the MCD (v6.1) is the most suitable atmospheric model to compute surface radiative spectra. In addition, the results show a dependence of the low-energy flux intensity on the hydrogen composition of the Martian soil (in the form of water), particularly for energies below 100 MeV, where the flux intensity can be reduced by a factor of up to 2.5 when the water mass ratio increases from 7.4% to 50%. Alongside a dependence on soil composition, GCR-induced flux at the surface is also highly dependent on atmospheric conditions, materialized by column depth or surface pressure Guo et al. (2017), although variations in GCR-induced dose rate vary by around 10% with pressure conditions, as RAD measurements and simulations have shown (Hassler et al., 2014; Zhang et al., 2022). In addition, the surface radiative environment depends on diurnal variations (Rafkin et al., 2014), local topography (Guo et al., 2021a), solar modulation and zenith angle (Khaksarighiri et al., 2023). Taking these elements into account in the modelling, in addition to a parametrizable geometry, enables ARAMIS to present better statistical indicators of fit with MSL-RAD measurements, for photons and neutrons, while remaining very close to the data on the rest of the spectrum, where the PHITS code is very well suited to the lowest energies Matthiä et al. (2016). Such high-energy agreement is mostly attributable to an improved estimate of the atmospheric neutrons that drive the high-energy secondary flux, which also increases with increasing surface pressure (Zhang et al., 2022). Furthermore, ARAMIS presents a strong fit with RAD data, as demonstrated by the better overall accuracy over different time periods of the |RD| performance indicator, less problematic to analyse than χ², which may vary as a function of RAD measurement uncertainties. Accurate modelling of secondary protons may prove crucial, especially for dose estimation due to SPEs, which are also dependent on surface pressure variations. The surface doses absorbed by the skin induced by OLTARIS and ARAMIS spectra have then been compared with the dose provided by the pivot methods (Guo et al., 2019b; Zhang et al., 2023). We found very good consistency between the pivot method and the doses induced by ARAMIS and OLTARIS. However, it is necessary to use results provided by ARAMIS or OLTARIS to estimate the dose of events with spectral shape combination, as the pivot method is not suitable for these events. Once validated against measurements, ARAMIS enabled us to observe the contribution of the various GCR spectrum components to the calculation of effective dose through the associated efficiency function, as well as the propagation of relative Monte Carlo errors in this function. The remaining limitations of this model to date are its computation time consumption and its poor approximation of surface heavy ion fluxes, which pave the way for further improvements, as does the integration of near-surface topography (craters, hills, etc.) into the model. Additionally, to reduce ARAMIS user time, other approaches could be taken, such as reducing the number of Z primaries to only take into account H, He, Li, C, O, Si and Fe (the most abundant particles in the GCR cycle). As a preliminary result of such a study, we found that for the atmospheric and exposure conditions presented in this study, taking only these particles into account would reduce neutron, photon, electron and proton integrated fluxes by around 1%, for a reduction in user time of around 60%. In addition, it was shown that primary heavy ions (Z > 2) contribute to surface dose by around 10%, with this contribution depending on surface pressure. A global validation for all pressure and exposure conditions is now required prior to an ARAMIS use with GCR sources reduced to these 7 ions. However, ARAMIS already offers a wide range of potential applications and studies for optimizing mission parameters or mapping surface doses and fluxes at different points for different exposures.

Acknowledgments

Monte Carlo and postprocessing calculations have been performed using HPC resources from the CNES Computing Center. This work, bearing the reference EUR CARe N°ANR-18-EURE-0003, has benefited from support managed by the Agence Nationale de la Recherche under the Programme Investissements d’Avenir. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

Data availability statement

The data associated with this article are available in Zenodo at https://doi.org/10.5281/zenodo.14035143.

Appendix A.

χ² and |RD| fitting values

Appendix B.

R_kj of secondary photons using FTFP_INCLXX_HP

Figure B1 Relative error per energy range in secondary spectra calculation for downward photons, for 105 launches using FTFP_INCLXX_HP and different primary particle: Proton (a), Helium (b), Carbon (c), Oxygen (d), Silicon (e) and Iron (f).

Figure B2 Relative error per energy range in secondary spectra calculation for upward photons, for 105 launches using FTFP_INCLXX_HP and different primary particle: Proton (a), Helium (b), Carbon (c), Oxygen (d), Silicon (e) and Iron (f).

Appendix C.

Values of effective dose yield function

References

All Tables

All Figures

	Figure 1 Schematic of the global ARAMIS organisation and structure.
In the text

Figure 2 Free space GCRs spectrum in Mars surrounding from 19 August 2012 to 17 February 2013 (corresponding to $\bar{\mathrm{\Phi }}=599 \mathrm{MV}$). The BON20 model is shown by solid lines while DLR model is represented by dotted lines. Due to different approaches (numerical for BON20, empirical for DLR) and input parameters (Liu et al., 2024), those two models differ significantly at low energies for most particles (for energies below 100-200 MeV). As detailed in 3.3.3, most particles of these energies do not reach the surface, as the Martian atmospheric cut-off is about 160 MeV for protons, for instance. Differences also exist at higher energies for heavier ions with lower GCR abundances compared to H and He. Plotted with data extracted from On-Line Tool for the Assessment of Radiation in Space (OLTARIS).

In the text

	Figure 3 Different SPE fluences in free space at 1 AU, corresponding to the historical major events recorded or reconstructed. Plotted with data extracted from OLTARIS and OMERE (for October 2003 data).
In the text

	Figure 4 Martian surface flux of neutrons (left) and photons (right) at Gale Crater from 19 August 2012 to 17 February 2013, given by ARAMIS and compared to measurements and other models with FTFP_INCLXX_HP. Similar plots with other physic lists are displayed in annex.
In the text

	Figure 5 Relative error per energy range in secondary spectra calculation for downward neutrons, for 105 launches using FTFP_INCLXX_HP and different primary particles: Proton (a), Helium (b), Carbon (c), Oxygen (d), Silicon (e) and Iron (f).
In the text

	Figure 6 Relative error per energy range in secondary spectra calculation for upward neutrons, for 105 launches using FTFP_INCLXX_HP and different primary particles: Proton (a), Helium (b), Carbon (c), Oxygen (d), Silicon (e) and Iron (f).
In the text

	Figure 7 Martian surface flux of proton, electron & positron, muon and pion at Gale Crater from 19 August 2012 to 17 February 2013, given by ARAMIS and compared to measurements and other models with FTFP_INCLXX_HP. For protons, solid lines represent cone-viewed downward spectra, while dashed lines represent direction-averaged spectra.
In the text

	Figure 8 Effective dose rate Yield function Yi(T, z0) (left) and relative error on dose rate Yield function calculation δi(T, z0) (right) at Gale Crater between 8 August 2012 to 30 January 2013, for different primary particle types and energies, with FTFP_INCLXX_HP.
In the text

	Figure 9 Martian surface neutron fluxes at Gale Crater (Curiosity landing site) and Vastitas Borealis (Phoenix landing site) from 19 August 2012 to 17 February 2013, given by ARAMIS with FTFP_INCLXX_HP. Neutron fluxes are multiplied by energy to enhance peak readability.
In the text

	Figure 10 Evolution of the CO2 column depth (in g · cm−2) with the solar longitude, Ls, for several coordinates corresponding to the landing sites of different probes and rovers (Curiosity, Perseverance, Phoenix, Viking 2). Plotted with data from the MCD v6.1. The central point value corresponds to 15 g · cm−2, hence the curve shapes.
In the text

	Figure 11 Martian surface flux of neutrons (left) and photons (right) at Gale Crater from 15 November 2015 to 15 January 2016, given by ARAMIS and compared to measurements and other models with FTFP_INCLXX_HP.
In the text

Figure 12 Martian surface fluence of neutrons (left) and photons (right) at Gale Crater from 19 August 2012 to 17 February 2013, given by ARAMIS with FTFP_INCLXX_HP and OLTARIS for different SPE historical events. Curves have been respectively scaled by factors of 10 and 100 for the November 1960 and August 1972 events for sake of clarity. The OLTARIS October 2003 event is not displayed as this event does not exist in the OLTARIS database and exporting OMERE fits to OLTARIS would imply using the fit function for only T < 300 MeV, as it is a combination of fits on the overall spectrum.

In the text

	Figure B1 Relative error per energy range in secondary spectra calculation for downward photons, for 105 launches using FTFP_INCLXX_HP and different primary particle: Proton (a), Helium (b), Carbon (c), Oxygen (d), Silicon (e) and Iron (f).
In the text

	Figure B2 Relative error per energy range in secondary spectra calculation for upward photons, for 105 launches using FTFP_INCLXX_HP and different primary particle: Proton (a), Helium (b), Carbon (c), Oxygen (d), Silicon (e) and Iron (f).
In the text