Nowcast and forecast of galactic cosmic ray (GCR) and solar energetic particle (SEP) fluxes in magnetosphere and ionosphere – Extension of WASAVIES to Earth orbit

Real-time estimation of cosmic-ray fluxes on satellite orbits is one of the greatest challenges in space weather research. Therefore, we develop a system for nowcasting and forecasting the galactic cosmic ray (GCR) and solar energetic particle (SEP) fluxes at any location in the magnetosphere and ionosphere during ground-level enhancement (GLE) events. It is an extended version of the WArning System for AVIation Exposure to SEP (WASAVIES), which can determine event profiles by using real-time data of the count rates of several neutron monitors (NMs) at the ground level and high-energy proton fluxes observed by Geostationary Operational Environmental Satellites (GOES) satellites. The extended version, called WASAVIES-EO, can calculate the GCR and SEP fluxes outside a satellite based on its two-line element (TLE) data. Moreover, organ absorbed-dose and dose-equivalent rates of astronauts in the International Space Station (ISS) can be estimated using the system, considering its shielding effect. The accuracy of WASAVIES-EO was validated based on the dose rates measured in ISS, as well as based on high-energy proton fluxes observed by POES satellites during large GLEs that have occurred in the 21st century. Agreement between the nowcast and forecast dose rates in ISS, especially in terms of their temporal structures, indicates the usefulness of the developed system for future mission operations.


Introduction
Estimation of the fluxes of high-energy particles in the Earth's magnetosphere is very important for designing space missions because they can adversely affect the health of astronauts and cause single-event upsets of semi-conductor devices used in satellites. Three sources must be considered in the estimation, namely, galactic cosmic rays (GCRs), trapped particles (TPs), and solar energetic particles (SEPs). The fluxes of GCR and TP are relatively stable and predictable compared to that of SEP, and the procedures for calculating GCR (Nymmik et al., 1996;Matthiä et al., 2013;O'Neill et al., 2014;Slaba & Blattnig, 2014) and TP (Ginet et al., 2013) fluxes are rather well established. In addition, software for calculating their mean fluxes in certain Earth orbits, such as CREME96 (Tylka et al., 1997) and SPENVIS (Heynderickx et al., 2004), were developed and released to public.
By contrast, SEP fluxes during a space mission are unpredictable because they increase suddenly when a large solar particle event (SPE) occurs. Thus, the worst-case scenarios are generally considered in mission design (Xapsos et al., 1999;Aghara et al., 2015;Jiggens et al., 2018). On the other hand, estimation of real-time SEP fluxes during SPE is beneficial for mission operation to take adequate actions, such as sheltering astronauts in well-shielded locations in their spacecraft (Townsend et al., 2018). A few studies (Hu et al., 2016;Matthiä et al., 2018) have been conducted to reproduce the SEP fluxes in spacecraft during certain SPEs, but these studies were mostly dedicated to post-exposure evaluation. Therefore, it is desirable to develop a system that can nowcast or forecast the SEP fluxes during a SPE. An active dosimeter-based method for estimating astronaut acute radiation risk during SPE in real time was recently proposed (Mertens et al., 2018).
With these situations in mind, we set out to develop a new computational method that can nowcast and forecast SEP fluxes at any location in the magnetosphere and ionosphere during a large SPE associated with a ground-level enhancement (GLE). It is an extended version of the WArning System for AVIation Exposure to SEP (WASAVIES) Kataoka et al., 2018;Sato et al., 2018b), which can nowcast and forecast radiation doses in the atmosphere by using real-time data of the count rates of several neutron monitors (NMs) at the ground level and high-energy proton fluxes observed by Geostationary Operational Environmental Satellites (GOES). The contribution of GCR can be also calculated by the system by using the PARMA model (Sato, 2015(Sato, , 2016. The most important feature of WASAVIES is that it is fully based on physics models of SEP transport from the Sun to the ground level of the Earth, and this feature enables us to smoothly extend the system applicable to the magnetosphere and ionosphere. The extended system is called WASAVIES-EO, where EO represents Earth Orbit. The two-line element (TLE) data of a satellite must be supplied to the system for calculating the GCR and SEP fluxes on its orbit. The organ doses and dose equivalents of astronauts in International Space Station (ISS) can be also estimated by the system considering its shielding effect, which has been evaluated using the Particle and Heavy Ion Transport code System (PHITS, Sato et al., 2018a) coupled with a virtual ISS model. Detailed calculation procedures of the original WASAVIES can be found in our previous paper (Sato et al., 2018b). Thus, the present study focuses on describing the extended part of the calculation procedures, together with its validation results, based on the experimental data measured by ISS and Polar Orbiting Environmental Satellites (POES) during large GLEs that have occurred in the 21st century.
2 Calculation procedures Figure 1 shows a flowchart of the calculation procedures of the system developed in this study. The algorithm employed in the original WASAVIES is summarized briefly in Section 2.1, while that of the extended part is described in detail in Sections 2.2 and 2.3.

Basic algorithm of WASAVIES
In WASAVIES, the profiles of each SPE are characterized by four parameters, namely, injection profile (IP), power index c of the primary SEP around the Sun, north-south tilt angle of Fig. 1. Flowchart of calculation procedures of WASAVIES-EO. The pastel green, cyan, and yellow boxes indicate the real-time analysis procedure, model or code, and data or database, respectively. The boxes in the blue frame are used in the original WASAVIES, while those with red frame are newly implemented in WASAVIES-EO. All databases included in the gray box are used for calculating the best-fit IP, h t , c, and N 0 , while only the database of SEP fluxes at 1 AU is used for calculating SEP flux outside satellite. The Kp index is automatically updated at intervals of 3 h.
T. Sato et al.: J. Space Weather Space Clim. 2019, 9, A9 ambient magnetic field incident to the magnetosphere h t , and normalization coefficient of the SEP fluence N 0 . These parameters are determined automatically from the count rates of several NMs at the ground level and the proton fluxes measured by GOES, which are downloaded continuously at intervals of 5 min. Note that we adopted the power law spectrum in rigidity of the primary SEP around the Sun.
Four databases must be prepared before executing the realtime and automatic analysis program of WASAVIES, and they are as follows: 1. Time profile of primary SEP fluxes at 1 astronomical unit (AU) for a certain IP and c, U 1AU,IP,c (t, E 0 , u), where t is the time after flare onset, E 0 is energy of the primary SEP, and u is pitch angle incident to the magnetosphere from the ambient magnetic field, the so-called Parker spiral's interplanetary magnetic field defined here as 45°declined from the Sun-Earth line. 2. The probability densities of the pitch angle u of protons with energy E 0 that can penetrate through the magnetosphere, P G;O;Kp;ht (E 0 , u), where G indicates the geographic coordinate of the arrival location at the top of the atmosphere, O denotes Earth's orbital condition characterized by month and hour, Kp is Kp index. 3. The response function of the fluences of particle i with energy E i in the atmosphere generated through extensive airshowers (EAS) induced by mono-energetic protons with energy E 0 , R EAS,i (E 0 , E i , h), where h denotes altitude. 4. The response functions of the standard NM64 NM for the incidence of particle type i with energy E i , R NM,i (E i ).
The U 1AU,IP,c (t, E 0 , u) database was developed by solving one-dimensional focused transport equations by using the method developed by Kubo et al. (2015), while P G;O;Kp;h t (E 0 , u) was prepared by tracing the trajectories of anti-protons emitted from the top of the atmosphere by using the proton trace model (Miyake et al., 2017) based on the empirical geomagnetic field model T89 (Tsyganenko, 1989). The other databases were constructed by performing Monte Carlo simulations by using the PHITS version 2.88 (Sato et al., 2018a). Details of the PHITS simulation procedures can be found in Sato et al. (2014).
In the real-time analysis program, GOES proton fluxes above 100 MeV and the count rates of 13 selected NMs are downloaded automatically from the Space Weather Prediction Center of NOAA (ftp://ftp.swpc.noaa.gov/pub/lists/particle/) and the Neutron Monitor Database NMDB (http://www.nmdb. eu/), respectively, at 5-min intervals. Information about the selected NM stations can be found in our previous paper (Sato et al., 2018b). At the end of each day, the daily count rates of each NM are calculated and used for determining the force field potential of the day. The GCR fluxes and the corresponding dose rates in the atmosphere are then calculated using the PARMA model (Sato, 2015(Sato, , 2016 coupled with the evaluated force field potential. The occurrence of GLE is checked by comparing the background and the current data of the NM count rates, as well as the GOES proton fluxes (Sato et al., 2018b). When the program detects GLE, it activates the SEP mode and automatically determines the best-fit values of the four parameters, namely, IP, c, h t , and N 0 , to reproduce the NM count rates and GOES proton fluxes. The SEP fluxes and the corresponding dose rates in the atmosphere are then calculated in real time by using the best-fit parameters. In parameter determination and dose evaluation, the aforementioned four databases, as well as the database of the effective dose conversion coefficients for isotropic irradiation (ICRP, 2010), are used. In addition, WASAVIES can forecast the SEP dose rates up to 24 h after flare onset, assuming that the evaluated parameters are time-independent, i.e., the time variation of the future SEP fluxes at 1 AU exactly follows their pre-calculated data contained in the database of U 1AU,IP,c (t, E 0 , u). However, this assumption is occasionally not adequate because the power index of the primary SEP, c, varies with time for some events, as discussed in our previous paper (Sato et al., 2018b).

Preparation of databases used in WASAVIES-EO
Two databases were additionally prepared for extending WASAVIES for application to Earth orbits, and they are as follows: 1. Response function of the fluences of particle i with energy E i at certain locations in a satellite for isotropic irradiation of mono-energetic particles j with energy E j , R S,ji (E j , E i ). 2. Response functions for converting the fluence of particle i with energy E i to various types of doses, R D,i (E i ).
For preparing R S,ji (E j , E i ), we performed a cosmic-ray transport simulation by using PHITS coupled with a virtual ISS model developed by Japan Aerospace Exploration Agency (JAXA) (Sato et al. 2018c). Figure 2 shows a three-dimensional (3D) view of the virtual ISS model drawn using PHITS. Only the Kibo, Columbus, Harmony, and Destiny modules have been modeled thus far. The models are composed mainly of the walls and the payload racks of these modules, which are made of Al alloys and Al, respectively, and their masses are approximated to those of the real ones. Protons and ions with energies up to 1 TeV/u and charges up to 28 (Ni) were considered as source particles j, and the fluences of those particles, as well as neutrons, pions, electrons, positrons, and photons, at several locations in the Kibo module were scored in the simulation. In this study, R S,ji (E j , E i ) inside the pressurized module of Kibo was used to validate the accuracy of WASAVIES-EO. The mean shielding thickness of the module was estimated to be 35.6 g cm À2 .
For R D,i (E i ), we prepared the response functions of two silicon-based active detectors-Liulin-5 in the MATROSHKA-R phantom and DOSTEL (DOSimetry TELescope)-to reproduce the experimental data presented by Reitz et al. (2005), Semkova et al. (2014), and Berger et al. (2018). The response function of Liulin-5 was evaluated by performing the PHITS simulation for mono-energetic particle irradiation incident on the MATROSHKA-R phantom-a spherical tissue-equivalent material with a radius of 17.5 cm-and scoring the absorbed dose around the location of Liulin-5 used for our validation, that is, at a depth of 4 cm depth from the surface. By contrast, the response function of DOSTEL was assumed to be the same as the stopping power of silicon because of the thin thickness of the detector-315 lm. In addition, fluence to dose and Q(L)-based dose-equivalent conversion coefficients for red bone marrow (RBM) of the reference adult male phantom were included in the database, where their numerical values were taken from the International Committee on Radiological Protection (ICRP) Publications 116 and 123 (ICRP, 2010(ICRP, , 2013.

Extension to real-time analysis program
The extended part of the real-time analysis program in WASAVIES-EO is executed after all analyses processes in T. Sato et al.: J. Space Weather Space Clim. 2019, 9, A9 the original WASAVIES are completed. During the solar quiet period, at first, the GCR fluxes of particle j with energy E j outside the magnetosphere, U GCR,1AU,j (E j ), are calculated using the GCR model developed by Matthiä et al. (2013) coupled with the daily force field potential evaluated from the count rates of the selected NM monitors. Note that U GCR,1AU,j (E j ) is assumed to be constant for a day, even during GLE, although the actual GCR fluxes fluctuate slightly if Forbush decrease occurs. This assumption indicates that the SEP contributions to the NM count rates and the GOES proton fluxes during GLE can be simply determined from the difference between their real-time and background data, where the background data were evaluated before the occurrence of GLE. Then, the geomagnetic transmission function of particle j with energy E j at the satellite location at the current time t c , T S,Kp,j (t c , E j ), is calculated using the particle trace model (Miyake et al., 2017) based on the T89 model (Tsyganenko, 1989). The TLE data of a satellite should be supplied by a user for this calculation, where the data for most satellites can be downloaded from the website of Space Track (https://www.space-track.org). Note that TLE is a de facto standard data format encoding a list of orbital elements. The GCR fluxes of particle j outside the satellite at t c , U GCR;Sout;j (t c , E j ), can be determined simply by multiplying U GCR,1AU,j (E j ) with T S,Kp,j (t c , E j ).
During GLE, the probability densities of the pitch angle u of protons with energy E 0 that can penetrate through the magnetosphere to the satellite location at t c for the best-fit h t , P S;Kp;ht (t c , E 0 , u), are calculated using the same particle trace model and TLE data. In this calculation, the program traces the trajectory of anti-protons with energy E 0 emitted from the satellite location in 288 directions, covering a solid angle of 4p, and estimates their pitch angle u escaping from the magnetosphere, where u is determined from the direction of the ambient magnetic field and the asymptotic direction. The current SEP fluxes outside the satellite, U SEP;Sout (t c , E 0 ), are calculated using the following equation: where N 0 (t c ) is the normalization coefficient of the SEP fluxes calculated using WASAVIES, U 1AU,IP,c (t, E 0 , u) is the time profile of the primary SEP fluxes at 1 AU for the best-fit IP and c, and t ct 0 is the time since flare onset. Note that the lowest value of energy E 0 included in the U 1AU,IP,c (t, E 0 , u) database is 80 MeV, which is sufficiently low to estimate aircrew doses but excessive for estimating astronaut doses. Thus, the program extrapolates the SEP fluxes down to 10 MeV by using the power index of U 1AU,IP,c (t c , E 0 , u) at 80 MeV. Both GCR and SEP fluxes outside the satellite are converted to the corresponding data inside the satellite by using the following equation: where E j should be replaced with E 0 for the SEP case. Then, the dose rates whose response function is included in the database of R D,i (E i ) can be estimated in real time by using the following equation: For forecasting the dose rates at a certain time t, T S,Kp,j (t, E j ) and P S;Kp;ht (t, E 0 , u) are calculated based on the expected satellite location at that time deduced from the TLE data. Then, Equation (1) is replaced by where U SEP;Sout (t, t c , E 0 ) denotes the forecasted SEP fluxes outside the satellite at time t, assuming the evaluated four parameters, IP, h t , c and N 0 , remain stable after current time t c . Finally, the forecasted GCR and SEP fluxes inside the satellite, as well as corresponding dose rates, are determined in the same manner as in Equations (2) and (3).

Results and discussion
The accuracy of WASAVIES-EO was examined using four sets of experimental data, which are listed in Table 1. The absorbed doses in silicon measured by DOSTEL and Liulin-5 in MATROSHKA-R were calculated using their response functions, R D,i (E i ), as described before. However, information about the shielding configurations around those detectors is not available. Thus, in this study, R S,ji (E j , E i ) of the pressurized module of Kibo was substituted with that of the modules on which the detectors were actually loaded during the GLEs. For reproducing the high-energy proton fluxes measured by Medium Energy Proton and Electron Detector (MEPED) loaded on POES, we simply integrated the calculated proton fluxes on their orbits, U GCR;Sout (t c , E 0 ) and U SEP;Sout (t c , E 0 ), for energies higher than 100 MeV. Figures 3-5 show a comparison of the absorbed dose rates measured in ISS with the corresponding data obtained from WASAVIES-EO. The irregular peaks observed in the measured dose rates can be ascribed to the passage of ISS through the South Atlantic Anomaly (SAA), which cannot be reproduced by our calculation because we did not consider the contribution of TP. Except for those peaks, the measured dose rates oscillated with a period of approximately 50 min because ISS approaches the polar regions closely with that interval.
The measured dose rates were mostly between 1 and 8 lGy h À1 , but they were occasionally considerably higher during GLEs, especially GLE 60 and 72.
Agreements between the measured and calculated dose rates are generally satisfactory, except for the SAA peaks, indicating the validity of WASAVIES-EO for calculating the GCR and SEP dose rates. However, WASAVIES-EO underestimated and overestimated the measured dose rates at the peaks of GLE 60 and 72, respectively. These discrepancies can mostly be ascribed to the use of R S,ji (E j , E i ) of the pressurized module of Kibo instead of that of the modules on which the detectors were loaded. The use of an empirical geomagnetic field model T89 in calculating P S;Kp;ht (t, E 0 , u) might be another cause of the discrepancies because they are more apparently observed for a certain hemisphereduring the northern and southern passages for GLE 60 and 72, respectively. In addition, extrapolation of the SEP fluxes down to 10 MeV by using the same power index at 80 MeV may have led to the overestimation observed at the peaks of GLE 72 because the SEP spectrum generally hardens with decreasing energy (Mertens et al., 2010;Mewaldt et al., 2005). Note that the assumption of the softer spectrum results in overestimation of the lower-energy SEP fluxes in our system because the model parameters were determined to reproduce the observables related to the high-energy SEP fluxes, that is, the count  Fig. 3. Comparison between dose rates measured using DOSTEL in US Laboratory of ISS during GLE 60 and corresponding data calculated using WASAVIES-EO. The peaks are alternatively attributable to the northern (N) and southern (S) passages of ISS except for the SAA peaks.
T. Sato et al.: J. Space Weather Space Clim. 2019, 9, A9 rates of the NMs and the GOES proton fluxes higher than 100 MeV. Figure 6 shows a comparison between the proton fluxes higher than 100 MeV measured by POES during GLE 72 and the corresponding data obtained using WASAVIES-EO. The data measured by the three satellites,  were downloaded from the Solar-Terrestrial Physics site of NOAA (https://www.ngdc.noaa.gov/stp/satellite/poes/). The corresponding data measured by GOES are also plotted in each graph. Similar to the absorbed dose rates shown in Figures 3-5, the GCR and SEP fluxes oscillate with a period of approximately 50 min, while irregular TP peaks are observed owing to the complex structure of the radiation belt. Except for the TP peaks, WASAVIES-EO can reproduce the measured data fairly well, indicating the validity of WASAVIES-EO not only to ISS but also other satellites. However, the WASAVIES-EO overestimated the peak SEP fluxes in the early stages of the GLE, where the calculated fluxes agreed with the data observed by GOES instead of POES. The reason for causing this discrepancy is currently under investigation. Figure 7 shows a comparison between the organ dose and the dose-equivalent rates for RBM of the reference male phantom inside ISS during GLE 69 and 70, as calculated using WASAVIES-EO, and the effective dose rates at the conventional flight altitude (12 km) above the South Pole, as calculated using WASAVIES. It is evident from the graphs that the dose  T. Sato et al.: J. Space Weather Space Clim. 2019, 9, A9 rates inside ISS increased in the late and early stages of GLE 69 and 70, respectively. This is because the orbital inclination of ISS is approximately 52°, and the SEP fluxes do not always increase when ISS approaches the polar regions closely because of lower L values on the orbit. Thus, the doses to ISS astronauts depend significantly on the timing of flare onset; GLE 70 occurred nearly at the worst timing in terms of the SEP exposure to astronauts because ISS passed through higher L-value locations up to around six during the GLE peak. Table 2 summarizes the total doses calculated by integrating the data shown in Figure 7 over 24 h from the flare onset. The mean quality factors are approximately 2.3 and 1.6 for the GCR and SEP doses, respectively. It should be mentioned that the calculated RBM doses as well as their mean quality factors are lower than the corresponding data measured by small detectors such as DOSTEL because of the self-shielding of human body. In addition, both GLEs occurred during Forbush decreases of earlier events, and thus, the GCR dose rates during the GLEs were slightly suppressed in comparison to the solar quiet condition.
Owing to the worst timing of occurrence of GLE 70 as discussed before, the total doses inside ISS during the event are higher than those during GLE 69, despite GLE 69 has been the largest event that has occurred in the 21st century. The dose equivalents in ISS are comparable to the corresponding effective doses at the flight altitude instead of the lower shielding thickness of ISS than that of the atmosphere, owing to the protection provided by the magnetosphere. If ISS were outside the magnetosphere during the entire span of GLE 69, the total dose equivalent would have been approximately 8 mSv. Note that the . Comparison between proton fluxes higher than 100 MeV measured using POES during GLE 72 and corresponding data calculated using WASAVIES-EO. The data measured using GOES are plotted in each panel.
T. Sato et al.: J. Space Weather Space Clim. 2019, 9, A9 total effective doses shown in the table are expected to be considerably higher than the values actually received by an aircrew because it is unrealistic to stay at the flight altitude for 24 h. Figure 8 shows the forecasted dose rates by WASAVIES-EO and WASAVIES for the same condition as in Figure 7, where the forecast started 20 min after GLE detection, that is,  Comparison between organ dose and dose equivalent rates for red-bone marrow (RBM) of the reference male phantom inside ISS during GLE 69 and 70 calculated using WASAVIES-EO and effective dose rates at conventional flight altitude (12 km) above South Pole calculated using WASAVIES. L-value at the ISS location is also plotted. adequate actions against SEP exposure, such as sheltering astronauts in well-shielded locations in the spacecraft. However, the absolute values of the forecasted dose rates at the peaks of SEP exposure are smaller and larger than the corresponding data shown in Figure 7 for GLE 69 and 70, respectively. This discrepancy can be ascribed to the gradual changes in the best-fit parameters evaluated using WASAVIES, as discussed in our previous paper (Sato et al., 2018b).

Conclusions
A physics-based warning system for aviation exposure to SEP, WASAVIES, was extended to be capable of nowcasting and forecasting the GCR and SEP fluxes not only in the atmosphere but also in the magnetosphere and ionosphere during a large SPE associated with a GLE. The extended version, called WASAVIES-EO, can calculate the GCR and SEP fluxes outside a satellite based on the TLE data of the satellite. In addition, the system can estimate the organ dose and dose-equivalent rates of astronauts in ISS, considering the shielding effect of ISS, as evaluated by PHITS coupled with a virtual ISS model. The accuracy of WASAVIES-EO was well validated by the dose rates measured in ISS as well as by the high-energy proton fluxes observed by POES during large GLEs that have occurred in the 21st century. Agreement between the nowcast and the forecast dose rates in ISS, especially in terms of their temporal structures, suggested the usefulness of the system for future mission operations. A function to calculate the TP fluxes should be implemented before developing an operational WASAVIES-EO system for practical use.
Two issues should be addressed for the further improvements of WASAVIES-EO. One is the limitation of 1-D simulation to solve the focused transport equation because the  T. Sato et al.: J. Space Weather Space Clim. 2019, 9, A9 assumption of well-connected magnetic field geometry is not always the case for GLEs. Energetic particles occasionally pass through complex 3-D dynamic structures as caused by the multiple occurrences of CMEs (Meyer et al., 1956;Richardson et al., 1991). We therefore need to properly incorporate the results from 3-D focused transport simulations, allowing the complexity of the background interplanetary magnetic field. The other issue is the inaccuracy in dose estimation due to the simple power-law extrapolation of the primary SEP fluxes below 80 MeV. It is known that the time evolution of CMEs associated interplanetary shocks and the magnetic field connectivity to the shocks make the energy spectra of SEP below 100 MeV very complicated (Luhmann et al., 2017). For nowcast, we can utilize the real-time proton data measured by satellites, but we need to study and model the complex behavior of such low-energy SEPs for forecast. For tackling both issues, realistic 3-D magnetic field structures of the inner heliosphere, as obtained from the 3-D magnetohydrodynamic simulations of CMEs such as SUSANOO-CME model (Shiota & Kataoka, 2016), will play essential role. The replacement of the geomagnetic field model T89 by its updated model (Tsyganenko & Andreeva, 2015) is also desirable.