1$6$73

'6SDFH5DGLDWLRQ7UDQVSRUWLQD6KLHOGHG

,&587LVVXH6SKHUH

-RKQ::LOVRQ

2OG'RPLQLRQ8QLYHUVLW\1RUIRON9LUJLQLD

7RQ\&6ODED

/DQJOH\5HVHDUFK&HQWHU+DPSWRQ9LUJLQLD

)UDQFLV)%DGDYL

2OG'RPLQLRQ8QLYHUVLW\1RUIRON9LUJLQLD

%UDQGRQ'5HGGHOODQG$PLU$%DKDGRUL

-RKQVRQ6SDFH&HQWHU+RXVWRQ7H[DV

  

6HSWHPEHU
7KH1$6$67,3URJUDP2IILFHLQ3URILOH

6LQFH LWV IRXQGLQJ 1$6$ KDV EHHQ GHGLFDWHG WR WKH

DGYDQFHPHQW RI DHURQDXWLFV DQG VSDFH VFLHQFH 7KH

1$6$ 6FLHQWLILF DQG 7HFKQLFDO ,QIRUPDWLRQ 67,

3URJUDP 2IILFH SOD\V D NH\ SDUW LQ KHOSLQJ 1$6$

PDLQWDLQWKLVLPSRUWDQWUROH

7KH1$6$67,3URJUDP2IILFHLVRSHUDWHGE\/DQJOH\

5HVHDUFK&HQWHUWKHOHDGFHQWHUIRU1$6$¶VVFLHQWLILF

DQGWHFKQLFDOLQIRUPDWLRQ7KH1$6$67,3URJUDP

2IILFHSURYLGHVDFFHVVWRWKH1$6$67,'DWDEDVHWKH

ODUJHVWFROOHFWLRQRIDHURQDXWLFDODQGVSDFHVFLHQFH67,

LQWKHZRUOG7KH3URJUDP2IILFHLVDOVR1$6$¶V

LQVWLWXWLRQDOPHFKDQLVPIRUGLVVHPLQDWLQJWKHUHVXOWVRI

LWVUHVHDUFKDQGGHYHORSPHQWDFWLYLWLHV7KHVHUHVXOWVDUH

SXEOLVKHGE\1$6$LQWKH1$6$67,5HSRUW6HULHV

ZKLFKLQFOXGHVWKHIROORZLQJUHSRUWW\SHV

x 7(&+1,&$/38%/,&$7,215HSRUWVRI

FRPSOHWHGUHVHDUFKRUDPDMRUVLJQLILFDQWSKDVHRI

UHVHDUFKWKDWSUHVHQWWKHUHVXOWVRI1$6$

SURJUDPVDQGLQFOXGHH[WHQVLYHGDWDRUWKHRUHWLFDO

DQDO\VLV,QFOXGHVFRPSLODWLRQVRIVLJQLILFDQW

VFLHQWLILFDQGWHFKQLFDOGDWDDQGLQIRUPDWLRQ

GHHPHGWREHRIFRQWLQXLQJUHIHUHQFHYDOXH

1$6$FRXQWHUSDUWRISHHUUHYLHZHGIRUPDO

SURIHVVLRQDOSDSHUVEXWKDYLQJOHVVVWULQJHQW

OLPLWDWLRQVRQPDQXVFULSWOHQJWKDQGH[WHQWRI

JUDSKLFSUHVHQWDWLRQV

x 7(&+1,&$/0(025$1'806FLHQWLILFDQG

WHFKQLFDOILQGLQJVWKDWDUHSUHOLPLQDU\RURI

VSHFLDOL]HGLQWHUHVWHJTXLFNUHOHDVHUHSRUWV

ZRUNLQJSDSHUVDQGELEOLRJUDSKLHVWKDWFRQWDLQ

PLQLPDODQQRWDWLRQ'RHVQRWFRQWDLQH[WHQVLYH

DQDO\VLV

x &2175$&7255(32576FLHQWLILFDQG

WHFKQLFDOILQGLQJVE\1$6$VSRQVRUHG

FRQWUDFWRUVDQGJUDQWHHV

x &21)(5(1&(38%/,&$7,21&ROOHFWHG

SDSHUVIURPVFLHQWLILFDQGWHFKQLFDOFRQIHUHQFHV

V\PSRVLDVHPLQDUVRURWKHUPHHWLQJVVSRQVRUHG

RUFRVSRQVRUHGE\1$6$

x 63(&,$/38%/,&$7,216FLHQWLILFWHFKQLFDO

RUKLVWRULFDOLQIRUPDWLRQIURP1$6$SURJUDPV

SURMHFWVDQGPLVVLRQVRIWHQFRQFHUQHGZLWK

VXEMHFWVKDYLQJVXEVWDQWLDOSXEOLFLQWHUHVW

x 7(&+1,&$/75$16/$7,21(QJOLVK

ODQJXDJHWUDQVODWLRQVRIIRUHLJQVFLHQWLILFDQG

WHFKQLFDOPDWHULDOSHUWLQHQWWR1$6$¶VPLVVLRQ

6SHFLDOL]HGVHUYLFHVWKDWFRPSOHPHQWWKH67,3URJUDP

2IILFH¶VGLYHUVHRIIHULQJVLQFOXGHFUHDWLQJFXVWRP

WKHVDXULEXLOGLQJFXVWRPL]HGGDWDEDVHVRUJDQL]LQJ

DQGSXEOLVKLQJUHVHDUFKUHVXOWVHYHQSURYLGLQJ

YLGHRV

)RUPRUHLQIRUPDWLRQDERXWWKH1$6$67,3URJUDP

2IILFHVHHWKHIROORZLQJ

x $FFHVVWKH1$6$67,3URJUDP+RPH3DJHDW

http://www.sti.nasa.gov

x (PDLO\RXUTXHVWLRQYLDWKH,QWHUQHWWR

KHOS#VWLQDVDJRY

x )D[\RXUTXHVWLRQWRWKH1$6$67,+HOS'HVNDW



x 3KRQHWKH1$6$67,+HOS'HVNDW



x :ULWHWR

1$6$67,+HOS'HVN

1$6$&HQWHUIRU$HUR6SDFH,QIRUPDWLRQ

6WDQGDUG'ULYH

+DQRYHU0'
1$6$73

'6SDFH5DGLDWLRQ7UDQVSRUWLQD6KLHOGHG

,&587LVVXH6SKHUH

-RKQ::LOVRQ

2OG'RPLQLRQ8QLYHUVLW\1RUIRON9LUJLQLD

7RQ\&6ODED

/DQJOH\5HVHDUFK&HQWHU+DPSWRQ9LUJLQLD

)UDQFLV)%DGDYL

2OG'RPLQLRQ8QLYHUVLW\1RUIRON9LUJLQLD

%UDQGRQ'5HGGHOODQG$PLU$%DKDGRUL

-RKQVRQ6SDFH&HQWHU+RXVWRQ7H[DV

1DWLRQDO$HURQDXWLFVDQG

6SDFH$GPLQLVWUDWLRQ

/DQJOH\5HVHDUFK&HQWHU 

+DPSWRQ9LUJLQLD 



6HSWHPEHU
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 

443-757-5802 
Contents

Abstract .......................................................................................................................................................... 1

Introduction .................................................................................................................................................... 1

Deterministic Code Development .................................................................................................................. 1

3D Marching Procedures...............................................................................................................................11

A Simple Inhomogeneous Geometry.............................................................................................................13

Monte Carlo Benchmarks..............................................................................................................................20

Conclusions ...................................................................................................................................................24

Acknowledgements .......................................................................................................................................24

Appendix A: Monte Carlo simulations..........................................................................................................25

Appendix B: Evaluation of :i for N = 6, 10, 14, 18, 22................................................................................26

References .....................................................................................................................................................28

iii
Figures

1. Geometric relations of quantities useful in solving equation (1). The symbol G

n is a unit normal vector.

................................................................................................................................................................. 2

2. The JSC phantom torso as flown on STS-91 [Cucinotta 2004]............................................................... 5

3. Apparatus for shield attenuation studies during shuttle flight [Cucinotta 2004]. .................................... 5

4. Experimentally measured and calculated attenuation curves in aluminum and polyethylene in shuttle 

flight [Cucinotta 2004]. ........................................................................................................................... 6

5. Isotropic (iso) and forward (for) neutron spectra produced by 500 MeV protons in aluminum.............. 7

6. Particle fluence induced by the Webber SPE in 30 g/cm2 ICRU tissue slab shielded and backed by 20 

g/cm2 aluminum slabs in the vicinity of the forward (left pane) and backward (right pane) interface. 

The depths of 0 and 30 correspond to the front and back interfaces, respectively. Fluences are shown at 

the biologically significant depths (solid) in tissue defined by ICRU 51 [1993] and similar depths 

(dashed) in the aluminum shield near the interfaces. .............................................................................10

7. Alpha fluences from Fig. 6 with the vertical axis on a linear scale and the horizontal axis restricted 

below 10 MeV/n.....................................................................................................................................10

8. Comparison of fluence spectra at 0 g/cm2 (left pane) and 30 g/cm2 (right pane) in a 30 g/cm2 ICRU 

tissue slab shielded and backed by 20 g/cm2 aluminum slabs exposed to the Webber SPE spectrum. ..11

9. Geometry of 3D marching procedure.....................................................................................................12

10. Geometric relations for a sphere of tissue shielded by a spherical shell of aluminum. The left pane 

shows the geometry with variables related to ray-tracing description. The right pane shows dimensions 

of geometry (in g/cm2

) and locations of detectors where results will be shown and compared. Both 

geometries drawn in cm scale.................................................................................................................14

11. Integral (E > 1 MeV) isotropic neutron source (particles/(g-event)) induced by the Webber SPE 

incident on an ICRU sphere and aluminum spherical shell....................................................................15

12. Integral (E > 1 MeV) isotropic neutron source (particles/(g-day)) induced by 1977 solar minimum ions 

of hydrogen, helium, carbon, and iron incident on an ICRU sphere and aluminum spherical shell.......15

13. Directional components for N = 1, 2, 6, 10, 14, 18, and 22....................................................................17

14. Webber SPE convergence test in 30 g/cm2 diameter tissue sphere within a 20 g/cm2 spherical shell of 

aluminum at the two interfaces...............................................................................................................18

15. Webber SPE induced fluence within the ICRU sphere/aluminum spherical shell neighboring the near 

(left panel) and distal (right panel) tissue/shield interface with N = 22..................................................19

16. Webber induced alpha fluence within the ICRU sphere/aluminum spherical shell neighboring the near 

(left panel) and distal (right panel) tissue/shield interface with N = 22..................................................19

17. Particle spectra at detector locations 0 g/cm2 (left pane) and 30 g/cm2 from top of tissue sphere (radius 

15 g/cm2

) surrounded by aluminum shell (thickness 20 g/cm2) exposed to the Webber SPE spectrum.20

18. Same as Fig. 17 but for incident solar minimum GCR hydrogen spectrum. ..........................................21

19. Same as Fig. 17 but for incident solar minimum GCR helium spectrum. ..............................................22

20. Same as Fig. 17 but for incident solar minimum GCR carbon spectrum. ..............................................22

21. Same as Fig. 17 but for incident solar minimum GCR iron spectrum....................................................23

22. Particle spectra for 'Z = 0, 1, 2 at detector locations 0 g/cm2 (left pane) and 30 g/cm2 from top of tissue 

sphere (radius 15 g/cm2

) surrounded by aluminum shell (thickness 20 g/cm2) exposed to the 1977 solar 

minimum GCR carbon spectrum............................................................................................................23

23. Same as Fig. 22 but for incident solar minimum GCR iron spectrum....................................................24

iv
Tables

1. Comparison of straight-ahead HZETRN transport model with shuttle flight data. [Wilson et al. 

2003b]. .................................................................................................................................................... 4

2. Calculated and “measured” (STS-91) quality factors. All dose equivalents were within 12% except for 

stomach which was -19.7%. .................................................................................................................... 5

3. Webber SPE dose as a function of number of rays, N, at various depths including those defined by 

ICRU. .....................................................................................................................................................18

4. Webber SPE dose equivalent as a function of number of rays, N, at various depths including those 

defined by ICRU. ...................................................................................................................................18

5. Total CPU seconds required for benchmark calculations. Total number of histories ran in Monte Carlo 

simulations are shown in parentheses.....................................................................................................21

v
Abstract

A computationally efficient 3DHZETRN code capable of simulating High Charge (Z) and Energy 

(HZE) and light ions (including neutrons) under space-like boundary conditions with enhanced 

neutron and light ion propagation was recently developed for a simple homogeneous shield object. 

Monte Carlo benchmarks were used to verify the methodology in slab and spherical geometry, and the 

3D corrections were shown to provide significant improvement over the straight-ahead approximation 

in some cases. In the present report, the new algorithms with well-defined convergence criteria are 

extended to inhomogeneous media within a shielded tissue slab and a shielded tissue sphere and tested 

against Monte Carlo simulation to verify the solution methods. The 3D corrections are again found to 

more accurately describe the neutron and light ion fluence spectra as compared to the straight-ahead 

approximation. These computationally efficient methods provide a basis for software capable of space 

shield analysis and optimization.

Introduction

Early space radiation shield code development relied on Monte Carlo (MC) methods [Alsmiller 1967, 

Lambiotte et al. 1971] and made important contributions to the space program. Due to intensive computational 

requirements, MC methods utilized restricted one dimensional problems leading to imperfect representation of 

appropriate boundary conditions [Alsmiller et al. 1972, Pinsky et al. 2001, Armstrong and Colburn 2001, Foelsche et 

al. 1974]. Even so, intensive computational requirements for MC codes remained, and shield evaluation was made

near the end of the design process in greatly simplified geometry to enhance computer efficiency [Armstrong and 

Colburn 2001, Wilson et al. 2002]. Resolving shielding issues at the end of the design cycle had a negative impact 

on the design, since resolving issues early could have minimized shield augmentation requirements and associated 

launch costs. This is especially true in post-launch augmentation, as was done for the International Space Station 

(ISS) [Shavers et al. 2004]. Furthermore, added shielding at the end of the design process could require de-scoping 

mission objectives, as instruments and equipment may be removed to meet launch requirements, as was done on the 

Viking Project.

Improved spacecraft shield design requires early entry of radiation constraints into the design process to 

maximize performance and minimize costs. As a result, NASA has been investigating high-speed computational 

procedures to allow shield analysis to be part of the preliminary design concepts following through to the final 

design allowing shield optimization procedures [Wilson et al. 2003a, 2004a,b]. For the last several decades, NASA

has pursued deterministic solutions of the Boltzmann equation allowing field mapping within the ISS in tens of 

minutes [Wilson et al. 2007] using standard finite element method geometry common to modern engineering design 

practice [Qualls et al. 2001, Wilson et al. 2003a, 2004a]. Ray tracing procedures in complicated geometry hinders 

the application of MC methods to such engineering models. Even so, deterministic methods have relied on the 

straight-ahead approximation, resulting in the HZETRN code with loosely defined impact on model uncertainty

[Wilson and Khandelwal 1974] but it handily facilitated the ISS augmented design for which no other code was 

capable. Recently, a 3D version of HZETRN has been developed in simple geometry (homogeneous sphere) with 

improved convergence criteria and verification of 3D effects using MC methods [Wilson et al. 2014].

Herein, the homogeneous media limitation will be removed from the solution methodology, and MC codes 

are used again to verify the solution methods in simple geometry. The restriction to slab and spherical objects will 

be maintained to allow efficient MC simulations for verification purposes. The next step in this development will be 

to use more complex geometric models so that simple mapping of the present methodology into more realistic 

applications can be studied. In the present report, the current status of transport code development will be briefly 

reviewed with emphasis on extending these developments into a generalized 3D version of the HZETRN code. 

These advances will use available MC codes, Geant4 [Agostinelli et al. 2003], FLUKA [Fasso et al. 2005, Battistoni 

et al. 2007], and PHITS [Sato et al. 2006, 2013] to judge the veracity of these developments, especially with regard 

to their 3D aspects. 

Deterministic Code Development

The relevant transport equations are the linear Boltzmann equations derived on the basis of conservation 

principles [Wilson et al. 1991] for the flux (or fluence) density, (, , ) Gj x 8 E , of a j type particle in the continuous 

1
slowing down approximation (CSDA) in which atomic processes are described by the stopping power Sj(E) for each 

ion type j (vanishes for neutrons, j = n) as

4 ( , , )] ( , ', , ') ( , ', ') ' ' j jk k E k

E E E E d dE Q G TG d#< 8 8 8 8 8   ̈  ̈ x x , (1)

where the differential operator on the left hand side is defined as 

1 ( , , )] ( , , ) ( ) ( , , ) ( ) ( , , ) j j jj jj

j

E E SE E E E

A E G G G TG s w      ̄ ¢ ± s #< 8 8 8 8 8 xx xx < , (2)

and solved subject to a boundary condition over the enclosure of the solution domain as shown in Fig. 1. 

Fig. 1. Geometric relations of quantities useful in solving equation (1). The symbol G

n is a unit normal vector.

Application of the CSDA in both laboratory and space shielding has been wide and the resulting errors 

discussed elsewhere [Wilson et al. 1984, 1987a, 1994, Shinn et al. 1998, Tweed et al. 2006a, 2006b]. Equation (1)

can be rewritten as a Fredholm equation [Wilson 1977] given by

' 4

( )( ) ( , , ) , ), , () ()

( ') ( ', '', , ') ( , ', '') ' '' ', () ()

j j

j j

j j

E j

j jk k E E

j j k

S E PE

E E

S EP E

A

P E E E E d dE dE

S EP E

H

H H

H

Q

G G

T G d

    ̄ ¡ ° ¢ ±

  ̈  ̈ ̈

x x

x8

8 (	8 8

88 8 8

(3)

with 

1

[ ( )] E R d RE H j j S    , (4)

[ ( ) ( ')]   RE RE j j x x 8 8 . (5)

In equations (1) and (3), ī(:,x) is the vector connected to x along -:; U is the projection of x onto : (Fig. 1), and d

is the projection of ī(:,x) onto :. The quantity Rj(E) is the distance a type j ion with energy E will travel before 

losing all of its energy to excitation/ionization of atomic electrons, and Pj(E) is the probability a type j ion with 

energy E will not have a nuclear interaction in coming to rest in the media. The zero charge limit of equations (1)

and (3) is discussed elsewhere [Wilson et al. 1991]. The range-energy relation is given by

2
0

' ( ) ( ')

E

j j

j

dE RE A

S E   ̈ , (6)

and the nuclear attenuation function is given by

0

( ') ( ) exp ' ( ')

E j

j j

j

E

P E A dE

S E

   ̄ T  ¡ °

¢ ±  ̈ , (7)

where ıj(E ') includes nuclear elastic and reactive processes. One obstacle to solving either equation (1) or (3) is the 

need to evaluate the integral d:' at arbitrary locations within the media and development of computational methods 

to efficiently handle this limitation.

Atomic interactions limit the contributions of charged particles in the transport process. For example, the 

protons and alphas produced in aluminum below 100 MeV/n contribute to the fluence only within a few centimeters 

of their collisional source, and the heavier ions are even more restricted. This is an important factor in that the 

transported secondary charged particle flux tends to be small at low energies, and the role of additional nuclear 

reactions by these low energy ions are likewise limited, an important fact that is implemented in the current 

development.

As noted above, a prime limitation in evaluating equation (1) or (3) is the evaluation of the integral over 

d:' at arbitrary locations within the media. The approach to a practical solution of equation (1) or (3) is to develop a 

progression of solutions from the simple to increasingly complex allowing early implementation of highperformance computational procedures and establishing a converging sequence of approximations with established 

accuracy criteria and means of verification and validation. The first step leading to the lowest order solution reduces 

the evaluation by introducing the straight-ahead approximation as guided by the nucleon transport studies of 

Alsmiller et al. [1965] using MC methods in which the differential cross sections were approximated as 

( , ', , ') ( , ') ( ') T TE jk jk EE EE 88 88  . (8)

This approximation produced dose and dose equivalent results to be within the statistical uncertainty of the 

MC result obtained using the fully angle dependent cross sections in slab geometry [Alsmiller et al. 1965, Wilson et 

al. 1991]. The relation of angular dependent cross sections to spacecraft geometry was further examined by Wilson 

and Kandelwal [1974] using asymptotic expansions about angular divergence parameters, demonstrating errors in 

the straight-ahead approximation to be on the order of the square of the ratio of distance of divergence (few 

centimeters) to radius of curvature of the shield (few to several m), resulting in a small relative error a10-4 in most 

human rated space systems [Wilson and Khandelwal 1974]. Introducing the straight-ahead approximation, equation

(8), into equation (3) reduces the Fredholm equation to a Volterra equation that can be solved using marching 

procedures [Lamkin 1974, Wilson and Lamkin 1974, 1975, Wilson 1977, Wilson and Badavi 1986]. Considering the 

success of the straight-ahead approximation for nucleons encourages its use in HZE transport where the greater mass 

of heavier ions reduces even further the angle of scatter [Wilson 1977]. The success of equation (8) results from the 

fact that lateral diffusion from a given ray is compensated in a flat plate by lateral diffusion along adjacent rays. 

Numerical marching procedures were developed to solve the transport equation under the straight-ahead 

approximation, resulting in the HZETRN code. A corresponding nuclear fragmentation model, NUCFRG2, was also 

developed for HZETRN, and the verification and validation processes utilizing the NUCFRG2 database are 

described elsewhere [Wilson et al. 1987b, 2005, 2006]. 

Space flight validation of HZETRN has been limited almost exclusively to low Earth orbit (LEO), 

containing both trapped particle and attenuated galactic cosmic ray (GCR) components. The two primary limitations 

in the LEO trapped environmental models AP8MIN and AP8MAX as discussed by Wilson et al. [2003b] is the 

assumption that the trapped particles are isotropic (resulting from the omnidirectional fluence description) and poor 

representation of the dynamic behavior. These omnidirectional models have been relatively successful in describing 

the radiation environment aboard the highly maneuverable shuttle spacecraft wherein anisotropies tend to be 

averaged. Such models have been found to be less accurate in the formation flying of the ISS, mainly oriented in the 

local horizontal plane along the velocity vector as was demonstrated elsewhere [Hugger et al. 2003; Wilson et al. 

2005; Nealy et al. 2006; Slaba et al. 2011a, 2013]. A dynamic and anisotropic trapped proton environmental model

was developed for future use in LEO shield design and operations [Wilson et al. 2003b]. More recently, an updated 

3
trapped proton model, AP9, has been made available and preliminary validation comparisons have shown reduced 

uncertainties [Badavi et al. 2014].

These environmental models are placed in a suitable form for evaluation of the incident radiation on the 

bounding surface of the six-degree of freedom motion of an orbiting spacecraft for shield evaluation [Wilson et al. 

2005]. To test the dynamic behavior, shuttle thermoluminescent dosimeter (TLD) data from 1983 to 2000 were 

used, giving good coverage for nearly two solar cycles. A sample of shuttle TLD measurements and the 

computational model results is given in Table 1. With the use of a normalization procedure intended to represent 

environmental model uncertainty, Badhwar et al. [1996, 2001] showed that the root mean square error of both 

observed and calculated dose rate and dose equivalent rate is within 15%. The normalization procedure is based on 

scaling the model trapped proton contribution to the TLD measurements. Recent validation studies by Slaba et al. 

[2011a, 2013] did not utilize normalization procedures and focused entirely on the GCR component of the LEO 

environment. Those studies revealed larger differences exceeding 40% if pion production and the associated 

electromagnetic cascade were not included, as was the case in most of the previous validation efforts. These 

comparisons highlight the usefulness of HZETRN within the straight-ahead approximation in many cases, but that 

further improvement is still necessary.

Table 1. Comparison of straight-ahead HZETRN transport model with shuttle flight data. [Wilson et al. 2003b].

Flight Date DRNM* DLOC TLD† (PGy/d) Calc. (PGy/d)

STS-41A 11/83 6421 3 64.6 59.6

STS-51D 4/85 6661 4 917.4 889.3

STS-31 4/90 5701 1 2141 2290

STS-43 8/91 5894 4 20.7 18.6

STS-62 3/94 6771 1 94.3 89.2

STS-65 7/94 6822 2 28.3 25.1

STS-67 3/95 6925 3 250.8 238.1

STS-80 11/96 6973 4 264.4 256.5

STS-82 2/97 7074 1 2978 3080

STS-91 6/98 6894 1 89.1 83.2

STS-101 5/00 6460 2 140.8 131.1

STS-92 10/00 6417 2 165.9 153.4

* Deep River Neutron Monitor count rate, † GCR corrected TLD-100 data

Most of the previous validation calculations were performed in a homogeneous media (aluminum alloy 

2219), but radiation protection requires, at a minimum, the calculation of the environment within the human body, 

yielding an inherently inhomogeneous problem.

Applying radiation exposure constraints to the space radiation shield design process requires not only the 

estimation of the internal spacecraft environment but how that environment is transported in human tissues of a 

characteristic human shape having location and orientation within the spacecraft [Wilson et al. 1993]. To examine 

this problem experimentally, a human skeleton (head and torso) was covered and filled with a tissue equivalent 

plastic and sectioned to allow instrumentation for in-flight measurements (see Fig. 2). The torso phantom was 

Computerized Tomography (CT) scanned to define the geometry for use as a computational model, especially in 

regard to organ location. The instrumentation used TLD, nuclear track foils and Tissue Equivalent Proportional 

Counters (TEPCs) to allow the mapping of the dose and dose equivalent within the phantom torso model. The 

instrumented torso phantom was flown on STS-91 followed by data analysis and dosimetric analysis using the 

HZETRN transport code. The results of the measurements expressed herein as the average quality factor are given in 

Table 2. If the environmental model is again normalized to measurements, the calculated dose equivalent was within 

12% of the measurements except for the stomach which was 19%. The cause of the differences was in large part 

due to the need to correct the galactic cosmic ray under-response of the TLDs [Cucinotta 2004], and the corrected 

data were all within 12% of the computational models.

The usefulness of the TEPC in evaluating dose equivalent in the space environment being established, 

Badhwar provided a test rig of five TEPCs (Fig. 3) to study dose equivalent attenuation in the LEO environment 

within the shuttle spacecraft. The bare detector allows characterization of the ambient field within the shuttle while 

the shielded detectors (spheres in the figure) sample the modified environment due to the shielding material. A 

study of the LEO radiation attenuation in polyethylene was conducted in the shuttle bay on STS-81 and are 

compared with a similar study of aluminum shielding on STS-89. 

4
Fig. 2. The JSC phantom torso as flown on STS-91 [Cucinotta 2004].

Fig. 3. Apparatus for shield attenuation studies during shuttle flight [Cucinotta 2004].

Table 2. Calculated and “measured” (STS-91) quality factors. All dose equivalents were within 12% except for 

stomach which was -19.7%.

Organ Qcalc Qmeas

Brain 2.10 1.8 + 0.1

Colon 2.17 2.0 + 0.2

Heart 2.15 2.0 + 0.4

Stomach 2.20 1.7 + 0.3

Thyroid 1.89 1.7 + 0.2

Skinbreast 2.10 1.8 + 0.2

Skinabdomen 2.11 1.9 + 0.2

Note that a direct comparison of aluminum shielding effectiveness to that of polyethylene cannot be made 

as the two flights are in radically different environments. STS-81 is in a high inclination orbit (51.6o

) with large 

contributions from galactic cosmic radiation and STS-89 is in a low inclination orbit (28.5o

) where the environment 

is dominated by trapped protons. The results are compared with the HZETRN model calculations in Fig. 4. 

One TEPC failed on STS-89, but it is still clear that good agreement was found between the HZETRN and 

the measured attenuation curves. This is also true for the STS-81 results for polyethylene shields as shown in the 

figure. It is clear that polyethylene has superior attenuation characteristics even in environments with large 

contributions from galactic cosmic radiation (STS-81). As a result, polyethylene is the material of choice for 

augmenting the ISS shield in spite of its need to be encapsulated to guard against fire hazards. Most important is the 

ability of HZETRN, even in the straight-ahead approximation, to evaluate correctly the dose equivalent from all 

components important to human rated systems and provide estimates of ISS augmentation requirements and the 

value of the benefits. Consequently, the layout of the augmentation in the ISS was evaluated using the HZETRN 

code within the straight-ahead approximation [Shavers et al. 2004].

5
Fig. 4. Experimentally measured and calculated attenuation curves in aluminum and polyethylene in shuttle flight

[Cucinotta 2004].

It is clear that the straight-ahead approximation provides accurate dosimetric results and a good 

approximation to the particle fluence in many circumstances (see Badhwar et al. [1995]). One limitation of the 

straight-ahead approximation is near the front boundary where the calculated neutron fluence vanishes (unless 

neutrons are part of the external radiation environment). This results from the straight-ahead approximation and can 

be improved by employing 3D corrections containing backward leakage to the front boundary. While leakage is an 

important process of neutron transport in space systems, it is less important to the charged particle transport due to 

range/energy relations. 

In the early transport studies of Wilson and Lamkin [1974, 1975] and Lamkin [1974], it was demonstrated 

that once neutrons are produced by a proton beam, that the recoupling of neutrons back to the proton field is limited, 

but the neutron fields tend to build-up with increasing penetration depth. This mainly occurs since many protons 

produced by neutron-induced reactions are of lower energy than the forward components and are quickly removed 

by atomic interactions (especially the isotropically produced protons). The probability of even a 100 MeV proton 

undergoing a nuclear reaction is only several percent, and the probability decreases substantially at lower energies.

The low-energy neutrons, on the other hand, have no effective atomic interactions and propagate until a nuclear 

interaction occurs or they escape from the media. Yet, these low-energy protons and other light ions produced in 

tissue through neutron collisions contribute significantly to biological injury and must be accounted. This was the 

basis of the studies of Heinbockel et al. [2003] starting with Clowdsley et al. [2000, 2001] using the NUCFRG2 

database [Wilson et al. 1995, 1998], with recent improvements provided by Slaba et al. [2010a]. Herein, we will 

concentrate on extending the transport formalism to full 3D using the NUCFRG3 [Adamczyk et al. 2012] database. 

Future work will be a direct evaluation of equation (3) quantifying the propagated error of the 3D neutron code.

Extension of the transport formalism begins by writing the nuclear reactive and scattering differential cross 

sections in the following form 

, ,

1 ( , ', , ') ( , ') ( , ', , ') 4

T TT jk jk iso jk for EE EE EE

Q

88 88   , (9)

where the first term is isotropic and associated with lower-energy particles produced, including target fragments, 

defined by 

, 2 ( , ') 2 ( , ', , ') ' jk iso jk B

EE EE d Q T T   ̈ 88 8 , (10)

where 2SB represents the backward hemisphere, and the second term in equation (9) is highly peaked in the forward 

direction and is associated mainly with direct quasi-elastic events and projectile fragmentation products [Wilson 

1977, Wilson et al. 1988]. In this context, quasi-elastic refers to collisions in which the energy transfer is small 

compared to the incident energy, and the secondary products are highly peaked in the forward direction. 

6
Surprisingly, even nucleon-induced reactions follow the simple form of equation (9), and the isotropic term 

extends to relatively high energies (see Fig. 5). Further details on the nature of the elastic and reactive cross sections 

used in HZETRN are given in Wilson et al [1991] and Cucinotta et al. [1994, 1998].

Fig. 5. Isotropic (iso) and forward (for) neutron spectra produced by 500 MeV protons in aluminum.

The first step towards a more general solution uses a bi-directional neutron transport approximation as first 

investigated by Clowdsley et al. [2000, 2001] using multigroup methods and recoupling to the light ions using an 

analytic solution as a quadrature over the induced light ion source term. These types of approximations were further 

studied and improved by Slaba et al. [2010a] and Heinbockel et al. [2011a, 2011b].

In the most recent bi-directional model described by Slaba et al. [2010a], the charged particle components 

were first solved in the straight-ahead approximation with the neutron coupling represented only by the forward 

component. The isotropically produced neutrons are then solved using a coupled bi-directional transport formalism 

to complete the neutron solution [Slaba et al. 2010a]. These isotropically produced neutrons further collide with the 

media, giving rise to an isotropic source of ions that are transported assuming only atomic interactions as an 

approximation to the complete solution, justified by the limited range of these ions. The governing equation is given 

by equation (1) which is solved using equation (8) for charged components and assuming the separation of the 

neutron field into a forward and isotropic component as in equation (9). The charged components are solved using 

the straight-ahead approximation as

, , ( , , )] ( , ') ( , , ') ' j for jk k for E k

G TG E E E E dE d#< 8 8   ̈ x x , (11)

for j z n, and for j = n, only the forward produced neutrons, given by equation (9), are represented as

, ,, [ ( )] ( , , ) ( , ') ( , , ') ' n n for nk for k for E k

TG T G E E E E E dE d   ̈ 88 8 < x x . (12)

Equations (11) and (12) are solved simultaneously using methods previously developed [Wilson et al. 2006; Slaba et 

al. 2010b]. Note that the ions generated by the isotropically produced neutrons are neglected in the above but will be 

treated later. Whereas the full straight-ahead cross section of equation (8) appears in the ion transport of equation 

(11), only the forward components of the cross sections appear in the transport of neutrons in equation (12) given by

, , 4 ( , ') ( , ', , ') ' nk for nk for EE EE d Q T T   ̈ 88 8 . (13)

The isotropic source of neutrons are treated as a perturbation given by

7
, , 4

, , 4

[ ( )] ( , , ) ( , ', , ') ( , ', ') ' '

( , ', , ') ( , ', ') ' ' .

n n iso nn n iso E

nk iso k for E k

E E E E E d dE

E E E d dE

Q

Q

TG T G

T G

d

d

 



 ̈  ̈

 ̈  ̈

< x x

x

8 8 88 8 8

88 8 8

(14)

Couplings between the isotropic neutron field and the ion fields of equation (11) are not explicitly included in 

equation (14) but are treated once the isotropically produced neutrons are evaluated. Clearly, equations (11) and (12)

can be solved independent of equation (14), but the forward solution of equations (11) and (12) enters into equation 

(14) generating the isotropic neutron source. Whereas equations (11) and (12) take on value only if an external 

source is present at the boundary, equation (14) has no external sources, and the resulting field develops from the 

internal isotropic neutron source. Development of an algorithm for solving equation (14) is expedited by 

approximating the double differential cross sections as

, , ( , ', , ') ( , ') ( ') ( , ') ( ') T TE TE nn nn f nn b EE EE EE 88 88 88   , (15)

where 

, ,, 2 ( , ') ( , ', , ') ' ( , ') 2 ( , ') / 4 nn f nn nn for nn iso F

EE EE d EE EE Q T T T QT Q    ̈ 88 8 , (16)

with 2SF as integration over the forward hemisphere and 

, , 2 ( , ') ( , ', , ') ' 2 ( , ') / 4 nn b nn nn iso B

EE EE d EE Q T T QT Q    ̈ 88 8 , (17)

with 2SB as integration over the backward hemisphere. Similar relations for nk iso , T in equation (14) are given by

, , ( ) , () ( , ', , ') ( , ') ( ') ( , ') ( ') T T ET E nk iso nk iso f nk iso b EE EE EE 88 88 88   , (18)

where

, () , , 2 ( , ') ( , ', , ') ' 2 ( , ') / 4 nk iso f nk iso nk iso F

EE EE d EE Q T T QT Q    ̈ 88 8 , (19)

with 2SF as integration over the forward hemisphere and 

, () , , 2 ( , ') ( , ', , ') ' 2 ( , ') / 4 nk iso b nk iso nk iso B

EE EE d EE Q T T QT Q    ̈ 88 8 , (20)

with 2SB as integration over the backward hemisphere. Using equations (15) to (20) in equation (14) results in

, ,,

, ,

, () ,

, () ,

[ ( )] ( , , ) ( , ') ( , , ') '

( , ') ( , , ') '

( , ') ( , , ') '

( , ') ( , , ') ' ,

n n iso nn f n iso E

nn b n iso E

nk iso f k for E k

nk iso b k for E k

E E E E E dE

E E E dE

E E E dE

E E E dE

TG T G

T G

T G

T G

d

d

d

d

 

 



 

 ̈

 ̈

 ̈

 ̈

88 8 <

8

8

8

x x

x

x

x

(21)

8
where the backward moving isotropic induced neutron field , (, , ) Gn iso x 8 E satisfies a similar equation obtained by 

replacing : by -: in equation (21) as

, ,,

, ,

, () ,

, () ,

[ ( )] ( , , ) ( , ') ( , , ') '

( , ') ( , , ') '

( , ') ( , , ') '

( , ') ( , , ') ' .

n n iso nn f n iso E

nn b n iso E

nk iso f k for E k

nk iso b k for E k

E E E E E dE

E E E dE

E E E dE

E E E dE

TG T G

T G

T G

T G

d

d

d

d

    



 



 ̈

 ̈

 ̈

 ̈

88 8 <

8

8

8

x x

x

x

x

(22)

Note that k for , G in equations (21) and (22) represents the forward directed flux along 8 , and consequently, 

, ( , , ') 0 Gk for x   8 E . The isotropically produced neutrons induce a source of isotropic light ions in collisions with 

the media resulting in a light ion fluence given by

, ,

1 ( , , ) ( ) ( , , ) ( , , ) ( , ') ( , , ') ' j j j j iso jn n iso E

j

E S E E E E E E dE

A E G G G TG

   ̄ s d ¡ °      ̄ ¢ ± s ¢ ±  ̈ 88 8 8 8 < x xx x , (23)

where j z n, and , ( , , ') Gn iso x 8 E is the isotropically produced neutron fluence. Equations (21) and (22) are solved 

simultaneously for no external fluence incident on the media directed along : (see Fig. 1). Note that the implicit 

assumption of equations (21) and (22) is that lateral diffusion from the ray : is fully compensated by diffusion from 

adjacent rays. This is strictly true for a semi-infinite slab and represents the use of periodic boundary conditions as 

an importance sampling method common to MC simulation practice. Note that the structure of these equations is the 

same as described by Wilson et al. [2014] for uniform media in simple geometry, and herein, the continuity of 

solution across the interface of the two media is enforced in the numerical solution. The formalism described by 

equations (11) - (23) is hereafter referred to interchangeably as the bi-directional solution, or the N = 2 solution. The 

N = 2 terminology was discussed by Wilson et al. [2014] and indicates that the isotropic neutron field is described 

by two ray directions, forward and backward.

The solution is now examined near the interfaces of a 30 g/cm2 semi-infinite slab of ICRU tissue (mass 

composition of 76.2% O, 11.1% C, 10.1% H, 2.6% N) [ICRU 1993] shielded and backed by 20 g/cm2 of aluminum 

exposed to the Webber solar particle event (SPE) model [Webber 1966] given as 

\ ^

9 10 ( 938) ( ) exp 239.1 ( ) / 100 200 ( ) p

E E pE p E G       ̄ ¢ ± , (24)

with pE EE ( ) ( 1876)   . Transport procedures were carried out using the bi-directional (N = 2) solution with 

results shown in Fig. 6.

Results are shown for the usual biologically significant depths assigned by the ICRU [1993] as 0.007, 0.3, 

and 1 g/cm2 but include fluence in the aluminum shield adjacent to the forward (near) and the backward (distal) 

interface. One first observes that the neutron fluence varies little across both regions. The proton and alpha fluences 

are nearly constant in the aluminum shield before the near interface but vary rapidly within the tissue after passing 

the interface and approach a new equilibrium spectrum by 1 g/cm2

. The main change is in greater thermalization in 

tissue. The proton and alpha fluences in tissue before the distal interface are nearly constant but vary rapidly in the 

aluminum past the interface. 

The solutions for charged components at low energies are proportional to the inverse stopping power of the 

media, and the transition of the charged particle spectra in the first several hundred microns of the interface is clear 

from the figure. Note that no transition occurs in the aluminum before the near interface and in tissue before the

distal interface resulting from the use of the straight-ahead approximation for charged components.

9
Fig. 6. Particle fluence induced by the Webber SPE in 30 g/cm2 ICRU tissue slab shielded and backed by 20 g/cm2

aluminum slabs in the vicinity of the forward (left pane) and backward (right pane) interface. The depths of 0 and 30 

correspond to the front and back interfaces, respectively. Fluences are shown at the biologically significant depths 

(solid) in tissue defined by ICRU 51 [1993] and similar depths (dashed) in the aluminum shield near the interfaces.

An expanded view of the alpha fluence on a linear scale is given in Fig. 7 for several depths in the 

neighborhood of both the near and distal interfaces. It is clear that a large transition occurs as a result of the change 

in both stopping power and alpha production within the two media. Improvements in the transition description will 

require, at a minimum, a bi-directional treatment of the charged particle fluence.

This bi-directional approximation (N = 2) in slab geometry is compared with the MC codes, Geant4, 

PHITS, and FLUKA, in Fig. 8. All Monte Carlo fluence results presented herein are plotted as mean values with 

estimated statistical uncertainty represented by one standard deviation and shown as vertical error bars placed at the 

energy bin midpoint. Details of the MC simulation setups are given in Appendix A for Figs. 8 and 17-23. Although 

the bi-directional approximation makes significant improvements to the neutron fluence in slab geometry, additional 

improvements with further 3D corrections are expected, especially for the charged particle fluences near interfaces. 

The success of the straight-ahead approximation under validation in human rated systems (shuttle, ISS) clearly 

meets the requirements for a useful solution of the transport problem, although the bi-directional solution of the 

neutron fluence makes marked improvements. Considered in the subsequent section is the transport in finite 

inhomogeneous objects in which lateral leakage can play an important role.

Fig. 7. Alpha fluences from Fig. 6 with the vertical axis on a linear scale and the horizontal axis restricted below 10 

MeV/n.

10
Fig. 8. Comparison of fluence spectra at 0 g/cm2 (left pane) and 30 g/cm2 (right pane) in a 30 g/cm2 ICRU tissue slab 

shielded and backed by 20 g/cm2 aluminum slabs exposed to the Webber SPE spectrum.

3D Marching Procedures

A sequence of approximations that span the formalism from the simple straight-ahead approximation to 

more complex propagation algorithms is described in this section. The governing equations in the CSDA ignoring 

Coulomb scattering are given by equations (1) and (3), and the cross sections include nuclear elastic and reaction 

processes. The double differential cross sections are found to be closely approximated by a forward component and 

an isotropic component given in equation (9), where the isotropic cross section is given by equation (10). The 

corresponding forward component is found from 

, ,

1 ( , ', , ') ( , ', , ') ( , ') 4

TTT jk for jk jk iso EE EE EE

Q

88 88   . (25)

As was done in the development of the bi-directional transport formalism, the solution is divided into forward and 

isotropic components where the isotropic neutron term is treated as a perturbation given by the following equations 

(26) - (29)

, ,, ( , , )] ( , ') ( , , ') ' j for jk for k for E k

G TG E E E E dE d#< 8 8   ̈ x x . (26)

On the scale of a vehicle or habitat, the space environment is nearly uniform in position and varies slowly in angle 

(if not isotropic) and each directional component :0 (see Fig. 9) is propagated according to equation (26) or 

according to the corresponding Volterra equation as

0

, 0 00

, ,0 '

( )( ) ( , , ) , ), , () ()

( ') ( ', '') ( , , '') '' ', () ()

j j

j for j

j j

E j

j jk for k for E E

j j k

S E PE

E E

S EP E

A

P E E E E dE dE

S EP E

H

H H

H G G

T G d

    ̄ ¡ ° ¢ ±

  ̈  ̈ 8

8 (	8 8

8

x x

x

(27)

where 0 0 [ ( ) ( ')] j j x8   x RE RE 8 . The means of handling the zero neutron charge is discussed elsewhere 

[Wilson et al. 1991]. 

11
Fig. 9. Geometry of 3D marching procedure.

Whereas the full straight-ahead cross section of equation (8) appears in the ion transport of equation (26),

only the forward propagating components appear in the transport of neutrons given by

4 ,

, 4

( , ', , ') ' , ( , ') ( , ', , ') ' ,

jk

jk for

jk for

EE d j n E E

EE d j n Q

Q

T T T

£¦

¦ v  ¦

¤

¦

¦  ¦¥

 ̈

 ̈

88 8

88 8 . (28)

In equation (26), the use of the straight-ahead approximation allows : to be viewed as a parameter. The equation 

may subsequently be solved according to the boundary condition set as the incident fluence from the direction : = 

:0 as shown in Fig. 9. The coupling to the isotropically produced neutron field for each :0 is given by

, , 4

, 0, 0

[ ( )] ( , , ) ( , ', , ') ( , ', ') ' '

( , ', , ) ( , , ') '

n n iso nn n iso E

nk iso k for E k

E E E E E d dE

E E E dE

Q TG T G

T G

d

d

 



 ̈  ̈

 ̈

8 8 88 8 8 <

88 8

x x

x

. (29)

Once the forward fluence is evaluated for direction :0 over the domain of x and E, one must yet solve equation (29)

for the induced isotropic neutron fluence within the shield configuration. Note that the last term in equation (29) is 

an isotropic source of neutrons from collisions of the forward propagating fluence given as 

, 0 , 0, 0 ( , , ) ( , ', , ) ( , , ') ' n iso nk iso k for E k

Y TG E E E E dE d8
 8 8 8 8   ̈ x x , (30)

for which the forward fluence , 0 ( , , ') Gk for x 8 E is a function of the penetration depth z (Fig. 9) along the direction 

:0. This is denoted by introducing , 0, 0 [ ( ), , ] ( , , ) n iso n iso D Y zE E x x 8
 8 8
 8  where z(x) is the penetration depth 

along :0 to the point x. Equation (29) is rewritten as

, , 4

, 0

[ ( )] ( , , ) ( , ', , ') ( , ', ') ' '

[ ( ), , ]

n n iso nn n iso E

n iso

E E E E E d dE

z E

Q TG T G

D

d

 



 ̈  ̈ 8 8 88 8 8 <

8
 8

x x

x

. (31)

Note that equation (31) also provides a basis for introducing anisotropic source terms as well. Whereas the 

forward propagating solution has its source fixed by the boundary condition, the isotropic solution (having no 

inbound fluence) is driven by internal sources generated by the collisions of the inbound fluence as modified by the 

forward propagator in penetrating to x along direction :0. In evaluating the penetration depth, one must account for 

12
differences in the material composition that determines the attenuation in reaching x from the boundary in addition 

to the strength of the neutron source term at location x.

The direction : may be viewed as a parameter in equation (31), allowing the equation to be solved along an 

arbitrary number of directions used to construct a numerical representation of , (, , ) Gn iso x 8 E . The : dependence of 

, (, , ) Gn iso x 8 E comes from the distance to the boundary t(:) and the penetration depth z along :0. The isotropic 

nature of the neutron source of equation (31) effectively decouples the incident direction :0 from the neutron 

direction of propagation :. To solve equation (31), one must know the neutron fluence, , (, , ) Gn iso x 8 E , at every x

from all directions :' approaching x from adjacent locations, as was the case of equation (14). Although the use of 

the approximation in (15), based on the assumption that losses along the ray : are compensated by diffusion from 

adjacent rays, is a demonstrated useful assumption in a semi-infinite slab (see Fig. 8 and Wilson et al. [2014], Slaba 

et al. [2010a, 2011b]), it is generally an overestimate in a sphere (or any finite object). Using equation (15) in 

evaluation of equation (31) seems a reasonable approximation in most human rated vehicles, leading to an extension 

of transport along an arbitrary ray : as a solution to the bi-directional equations

, ,,

, ,

, 0

( )] ( , , ) ( , ') ( , , ') '

( , ') ( , , ') '

[ ( ), , , ] ,

n n iso nn f n iso E

nn b n iso E

n iso

E E E E E dE

E E E dE

z E

TG T G

T G

D

d

d

 

 



 ̈

 ̈

<8 8 8 <

8

8 8

x x

x

x

(32)

where the backward moving field satisfies a similar equation obtained by replacing : by -: in equation (32) as

, ,,

, ,

, 0

( )] ( , , ) ( , ') ( , , ') '

( , ') ( , , ') '

[ ( ), , , ] .

n n iso nn f n iso E

nn b n iso E

n iso

E E E E E dE

E E E dE

z E

TG T G

T G

D

d

d

   



 

 ̈

 ̈

<8 8 8 <

8

8 8

x x

x

x

(33)

These isotropically produced neutrons produce isotropic light ions in collisions with the media resulting in a light 

ion fluence approximately given by

, ,

1 ( , , ) ( ) ( , , ) ( , ') ( , , ') ' j j j iso jn n iso E

j

E S E E E E E dE

A E G G TG

   ̄ s d ¡ °  

s ¢ ±  ̈ 88 8 8 < xx x , (34)

where , ( , , ') Gn iso x 8 E is the isotropically produced neutron fluence. Equations (32) and (33) are solved 

simultaneously for no additional external fluence incident on the boundary of the media directed along :0 (see Fig. 

9). Note that the implicit assumption of equations (32) and (33) is that lateral diffusion from the ray : is fully 

compensated by diffusion from adjacent rays. This was demonstrated for slab geometry as a reasonable result as 

shown in Figs. 6-8. The lateral leakage for incident particles along :0 in this case is found to a first order 

approximation by solving equations (32) and (33) along each :i between the ray limits {-t(-:i), t(:i)} for a given 

set of N directions :i. A set of directions from N = 1 to N = 22 was introduced by Wilson et al. [2014], and given in 

Appendix B. Convergence tests are applied in the next section to verify the ability of discrete ray numbers in 

describing the isotropic radiation fields. 

A Simple Inhomogeneous Geometry

Consider now a 15 g/cm2 radius tissue equivalent sphere (ICRU) shielded by a 20 g/cm2 spherical shell of 

aluminum, as shown in Fig. 10. The penetration depth now generally includes aluminum and tissue thicknesses. The 

solution along any direction :i must generally include propagation through both media from the outer limits of the 

object along the ray r:i to the evaluation point. In this example, the solution within the ICRU sphere serves as 

13
proxy for the astronaut (ICRU phantom, see ICRU [1993]). In particular, the D*

(td, :0) and H*(td, :0) as defined by 

ICRU are used. With verification of the methods in this simplified geometry, these methods will be ready to apply to 

more realistic astronaut and spacecraft geometry in subsequent studies. Simple geometry is chosen to allow efficient 

MC simulation for comparison with the calculation of quantities using the above formalism. In this way, the solution 

methodology is verified and experience gained in the formalism in preparation for connecting to the general 

geometric descriptions generated by the engineering design process for which MC simulation is difficult to 

implement. Target points (detectors) are placed at depths along the z–axis (x = 0, y = 0). Most results will be shown 

at or near detectors located at the top of the tissue sphere (0 g/cm2

) and bottom of the tissue sphere (30 g/cm2) as 

shown in the right pane of Fig. 10. These depths are being shown to highlight transition effects that occur near 

material interfaces, as previously shown in Figs. 6 and 7.

The external radiation environment (external source boundary condition) is assumed to be anti-parallel with 

the z-axis, uniform in the x-y plane, positioned above the inhomogeneous sphere, and directed down onto the top of 

the sphere as shown in Fig. 10. The added complexity in the calculation is the separate evaluation of the tissue 

thickness and shield thickness along each :, along with the evaluation of the penetration through the shield and 

tissue along :0 for each : for use in the transport procedure. Note that the penetration depths for sources in the 

lower hemisphere will generally have double penetration through the shield, once at the top and again at the bottom 

after passing through the tissue equivalent sphere. 

Fig. 10. Geometric relations for a sphere of tissue shielded by a spherical shell of aluminum. The left pane shows the 

geometry with variables related to ray-tracing description. The right pane shows dimensions of geometry (in g/cm2

)

and locations of detectors where results will be shown and compared. Both geometries drawn in cm scale.

Field values are to be evaluated at a point x where the thickness of material t is known over an array of 

directions, :i (i.e. ti = t(:i)). Note that each :i for i odd is backed by :i+1 = -:i. To complete the solution 

methodology, the directions of propagation, :i, must be selected, and the penetration depth, zp, for a given :0 must 

be evaluated. The N directions, :i, are chosen such that each ray accounts for the same solid angle, ': = 4S/N, on 

the unit sphere. The evaluation will be made for a single direction of incidence, :0, but this could include an array 

of incident directions for which the solutions would be integrated after all transport procedures have been 

completed. In that low-energy neutron transport is dominated by diffusive processes, it is anticipated that the 

solution will rapidly converge for low values of N.

The source of neutrons induced within the inhomogeneous sphere by the Webber SPE event is shown in 

Fig. 11. Note that the excess attenuation in tissue relative to aluminum for the primary protons is seen as the slight 

shadow cast by the tissue sphere, and the weaker neutron source strength within the tissue shifts the neutron 

production to lesser values easily seen at the top of the tissue sphere. The source of neutrons induced by four 

components of the 1977 solar minimum GCR ions of hydrogen, helium, carbon, and iron is shown in Fig. 12. The 

greater penetrating power of the GCR is apparent in the figures resulting from the much greater energy of the GCR 

ions. Also, note that the fluence of the forward propagating light ions increases with penetration depth, forming the 

most intense neutron production near the exit. While the GCR hydrogen and helium result in a shadow behind the 

tissue sphere, the more complex ions (carbon and iron in the present example) are more effectively fragmented on 

the hydrogen contained in the tissue, forming a more intense source of neutrons behind the tissue sphere. 

14
Fig. 11. Integral (E > 1 MeV) isotropic neutron source (particles/(g-event)) induced by the Webber SPE incident on 

an ICRU sphere and aluminum spherical shell.

Fig. 12. Integral (E > 1 MeV) isotropic neutron source (particles/(g-day)) induced by 1977 solar minimum ions of 

hydrogen, helium, carbon, and iron incident on an ICRU sphere and aluminum spherical shell.

The resulting neutron fluence depends not only on the distance, t(:), to the boundary but the intensity of 

the source along the ray in reaching the boundary. Clearly for solar particle events, the most intense regions must be 

adequately represented in the domain of the solution. The convergence rate for increasing number of rays N will be 

tested in the following sections. The number of rays, :i, for the current example is now discussed.

15
N = 1, Straight-ahead approximation

For each :0, only one value of :i is allowed and is taken to be that of the incident direction :0. The transport 

distance x along :1 and penetration depth z are the same, reducing complication in the formulation for N = 1. In this 

case, leakage from the media is only through the distal face. The structure of this code is quite different than the 

original HZETRN code, although computational results are identical, and the original HZETRN code provides a 

verification test-case of the new 3DHZETRN code. For N = 1, the comparisons in the prior sections describing the 

verification and validation of HZETRN are fully applicable. If one restricts interest to N = 1, then the HZETRN code 

is preferred in that computational efficiency is higher, but the added overhead of 3DHZETRN for N = 1 is minimal. 

In the following sections, comparisons between 3DHZETRN and MC simulations will also include results from the 

HZETRN code (N =1) to highlight shortcomings of the straight-ahead approximation and improvements of the 3D 

corrections presented herein.

N = 2, Bi-directional approximation

The first step towards a full 3D code is the treatment of backward propagating particles. This was first 

accomplished by Clowdsley et al. [2000, 2001] and later improved by Slaba et al. [2010a] with demonstrated 

improved results [Slaba et al. 2010a, 2011b; Heinbockel et al. 2011a, 2011b]. In this case, propagation of the 

isotropic neutrons is represented in two directions by :1 = :0 and :2 = -:0. Therefore, the penetration depth z is 

equal to x along :1 and equal to t(:1) + t(:2) - x along :2. Note that lateral leakage is not represented in this 

approximation, although leakage at the forward and backward boundaries is represented. The neutron transport in 

semi-infinite slab geometry has extensive MC results for verification for which this bi-directional approximation, 

ignoring lateral leakage, is most appropriate as demonstrated in Fig. 8. This approximation is appropriate for slab 

geometry and works well in applications in which the shield geometry is a collection of thin plates. The next step for 

treatment of lateral leakage is now discussed.

N > 2, N-stream approximation

For general N, the transport formalism requires the penetration depth, zp, to be known as a function of 

transport depth, xT, along the propagation direction : for each :0. In order to simplify algebraic manipulation and 

expressions, the center of the sphere is assumed to be positioned at the origin with the positive z-axis anti-parallel to 

the source direction :0, as shown in the left pane of Fig. 10. Due to the simplified geometry and source orientation, 

an analytic expression for the penetration depth through the aluminum shield and tissue target may be determined. 

The point x0 may be written as 

(1) (2)

0 ( )( ) st T T tt x x        ̄ ¡ ° ¢ ± x x 8 , (35)

where ts and tt are the known shield and target thicknesses along direction -: from x to the outer boundary, 

respectively. The symbols (1)

Tx and (2)Tx represent the transport depth through the shield and target along :,

respectively. The points (1)

( x and (2)( x may be written in terms of x0 as 

(1) (2) (1)

0 p z x x ( (   8 , (36)

(2) (2)

0 0 p z x x (   8 , (37)

where (1)

p z and (2)p z are the penetration thicknesses through the aluminum shield and tissue target, respectively. The 

expression for (2)

p z is obtained first by writing equation (37) in component form and substituting into the equation of 

a tissue sphere with radius 15 g/cm2

. This yields

16
2 (2) 2 2

00 0 00 15 ( ) p z   x xx < < 8 8 . (38)

The expression for (1)

p z is obtained by writing equation (36) in component form and substituting into the equation of 

an aluminum shell with thickness 20 g/cm2 surrounding a tissue sphere with radius 15 g/cm2

. This yields

2 (1) 2 2 (2) (2) (2)

0 0 (15 20) ( ) p z     ( (( x xx < < 8 8 . (39)

For preliminary testing and benchmark comparisons against MC codes, the solutions for N = 6, 10, 14, 18, 

and 22 are considered separately, in addition to the simple N = 1 and N = 2 cases discussed above. The distributions 

for all N are shown in Fig. 13. Note that the N = 6, 14, and 22 cases include rays perpendicular to the z-axis (i.e. 

tangent to the sphere surface at the top and bottom detector locations). On the surface of the sphere, these rays 

correspond to shielding by aluminum only with no intervening tissue. This is corrected in the N = 10 and N = 18 

distributions. The directional elements for each N are given in the Appendix B, along with a brief description of how 

the elements were determined. For each N, the omni-directional flux is obtained by integrating the fluxes along each 

ray direction after the transport procedure. 

Fig. 13. Directional components for N = 1, 2, 6, 10, 14, 18, and 22.

Preliminary testing

One very large solar particle event occurred on February 23, 1956 and has served as an important event to 

not only space operations but for high altitude aircraft as well [Foelsche et al. 1974]. As a result, this event has 

become a test environment used over many decades and has been (at times) approximated by the Webber spectrum 

given by equation (24). Of interest has been the dose and dose equivalent within the tissue sphere for this Webber 

spectrum. Dosimetric quantities are evaluated at depths in the tissue sphere at td = 0, 0.007, 0.3, 5, 25, 29, 29.7, 

29.993, 30 g/cm2 with results in Tables 3 and 4 for N = 1 through N = 22. The dose even on the surface of the tissue 

sphere is seen to converge rapidly with increasing N. The local fields appear dominated by transport through the 

aluminum shell giving improved convergence behavior over the case concerning the surface of the aluminum shield 

[Wilson et al. 2014]. The rate of convergence of dose equivalent is somewhat slower as the details on the penetration 

of charged particles produced by neutrons are more dependent on angular factors. 

Dose and dose equivalent convergence rates are reasonably rapid for all locations within the tissue sphere 

except near the tissue/aluminum interface as shown in Tables 3 and 4. At each depth, results appear to 

asymptotically approach a common value, as will be further demonstrated for transported fluence spectra. The 

induced particle fluence is shown in Fig. 14 at the near and distal interfaces. Also shown are the lowest order N = 1 

and N = 2 solutions of the HZETRN code. The N = 1 solution has no contributions from the backward leaking 

neutrons that have been moderated within the tissue sphere, as evidenced by the incorrect spectral shape at the 

lowest energies. The incorrect spectral shape is largely corrected by the N = 2 bi-directional approximation as seen 

in Fig. 14. 

17
Table 3. Webber SPE dose as a function of number of rays, N, at various depths including those defined by ICRU.

d (g/cm2

) Dose (cGy/event)

N = 1 N = 2 N = 6 N = 10 N = 14 N = 18 N = 22

0 7.01 7.01 7.03 7.03 7.02 7.03 7.03

0.007 7.00 7.00 7.08 7.03 7.05 7.03 7.04

0.3 6.64 6.70 6.73 6.71 6.72 6.72 6.72

1 6.07 6.12 6.14 6.13 6.13 6.13 6.13

5 3.83 3.87 3.86 3.86 3.86 3.86 3.86

25 0.72 0.73 0.71 0.71 0.71 0.71 0.71

29 0.56 0.56 0.54 0.55 0.55 0.55 0.55

29.7 0.53 0.54 0.52 0.52 0.52 0.52 0.52

29.993 0.52 0.53 0.51 0.51 0.51 0.51 0.51

30 0.52 0.52 0.50 0.51 0.51 0.51 0.51

Table 4. Webber SPE dose equivalent as a function of number of rays, N, at various depths including those defined 

by ICRU.

d (g/cm2

) Dose equivalent (cSv/event)

N = 1 N = 2 N = 6 N = 10 N = 14 N = 18 N = 22

0 11.33 11.15 11.61 11.35 11.41 11.45 11.42

0.007 11.18 11.01 12.10 11.31 11.58 11.42 11.51

0.3 10.05 10.45 10.92 10.62 10.71 10.69 10.70

1 9.25 9.63 9.82 9.68 9.73 9.73 9.73

5 6.11 6.36 6.20 6.17 6.20 6.20 6.19

25 1.46 1.50 1.29 1.35 1.34 1.33 1.34

29 1.19 1.23 1.06 1.11 1.11 1.10 1.10

29.7 1.15 1.19 1.05 1.08 1.09 1.08 1.07

29.993 1.13 1.16 1.09 1.06 1.09 1.06 1.07

30 1.13 1.12 0.98 1.01 1.01 1.01 1.00

Fig. 14. Webber SPE convergence test in 30 g/cm2 diameter tissue sphere within a 20 g/cm2 spherical shell of 

aluminum at the two interfaces.

The solution methodology appears to be converged for N t 10. This is supported in Tables 3 and 4 by 

comparing exposure values for different N at a fixed depth. The variation is found to be less than 1% for dose and 

less than 2% for dose equivalent for N > 10 at all depths. Convergence for N > 10 is also supported in Fig. 14 in that 

spectra for neutrons, protons, and alphas are nearly indistinguishable. The fluence spectra in the neighborhood of the 

near and distal interface are shown in Fig. 15. A rapid transition occurs in passing from aluminum to tissue in the 

18
near interface and in passing from tissue to aluminum in the distal interface. The alpha fluence spectra are shown in 

the neighborhood of the near and distal interfaces on a linear scale in Fig. 16. Qualitatively, the spectral 

characteristics are similar to those near the slab interfaces in Fig. 7. For example, the alphas produced in the 

aluminum shield rapidly form a new equilibrium spectrum within the first few millimeters of tissue as observed in 

the near interface. Similarly, at the distal interface, the leakage of alphas from the tissue into the aluminum shield 

quickly forms a new equilibrium spectrum within the first few millimeters of aluminum. 

Fig. 15. Webber SPE induced fluence within the ICRU sphere/aluminum spherical shell neighboring the near (left 

panel) and distal (right panel) tissue/shield interface with N = 22.

Fig. 16. Webber induced alpha fluence within the ICRU sphere/aluminum spherical shell neighboring the near (left 

panel) and distal (right panel) tissue/shield interface with N = 22.

MC benchmarks will now be used to verify as much of the solution as possible. It is found that the details 

of the secondary alpha spectra for the Webber SPE boundary condition (for example, see Fig. 16) are difficult to 

verify with MC codes in the spherical geometry without excessively large run-times. For the same reason, defining 

the transition regions near the interfaces for neutrons and protons is also difficult. Consequently, the larger features 

are used to verify the solution of the neutron and proton fluences at the interfaces and results deeper within the tissue 

sphere will be the focus. If the trends in the solution over larger distances are verified, it is expected that the details 

revealed by 3DHZETRN are reasonably correct. The N = 22 will be used in the MC benchmarks discussed in the 

next section.

19
Monte Carlo Benchmarks

In 2005, a systematic benchmark exercise was administered as a means of verifying various computational 

procedures resulting from the first phase of code development (1996-2004). One component of transport code 

verification uses established benchmarks using standard environments: the Webber SPE spectrum given in equation 

(24) and the 1977 solar minimum GCR environment. The standard configuration of 30 cm semi-infinite slab of

water shielded in front by 20 g/cm2 of either aluminum or iron utilizing only dose and dose equivalent (ICRP 26 

quality factor) has been used in the past [Wilson et al. 1991]. It was generally concluded that HZETRN agreed with 

the MC codes as well as the various MC codes agreed among themselves [Wilson et al. 2006; Heinbockel et al. 

2011a, b; Slaba et al. 2010a]. There are currently no compelling reasons to attempt to use MC methods in design and 

operational applications. However, as will be further demonstrated below, the MC codes provide an indispensable 

guide in the development of the 3DHEZTRN code, which does provide a useful basis for space radiation protection. 

Attention is now directed towards benchmarking 3DHZETRN in which the added complications of

finite/inhomogeneous geometry are addressed.

Webber SPE Benchmark

The benchmark SPE spectrum is the Webber [1966] approximation of the February 23, 1956 SPE (not an 

accurate representation but a historically useful test spectrum, see [Scott and Alsmiller 1967]). For the present study, 

the induced neutron and proton fluence spectra were evaluated in the 70 g/cm2 diameter aluminum/tissue sphere by 

three MC codes (Geant4, PHITS, and FLUKA) and 3DHZETRN (N = 22) with the results at 0 and 30 g/cm2 depths 

in the tissue sphere shown in Fig. 17. Monte Carlo results for protons below 1 MeV have been removed since the 

statistical fluctuations were exceedingly large.

The proton fluence predictions among the three codes are in reasonable agreement above 1 MeV with very 

large statistical uncertainty below 1 MeV (not shown). The larger discrepancy lies with the neutron spectra where 

disagreement among the MC is on the same order as for 3DHZETRN. The FLUKA and 3DHZETRN results exhibit 

a similar inflection in the neutron spectrum near 10 to 60 MeV that may be related to the use of an early FLUKA 

approximation to the secondary neutron spectrum used in HZETRN [Ranft 1980, Wilson et al. 1988] and carried 

over into 3DHZETRN. Some of the differences in proton spectra on the back surface may be associated with the 

simplified assumptions of the straight-ahead approximation in light ion transport.

Fig. 17. Particle spectra at detector locations 0 g/cm2 (left pane) and 30 g/cm2 from top of tissue sphere (radius 15 

g/cm2

) surrounded by aluminum shell (thickness 20 g/cm2) exposed to the Webber SPE spectrum.

20
A principal difference in the three codes is the computational efficiency. The time required by the codes to 

run this benchmark is given in CPU seconds used and listed in the first column of Table 5. The 3DHZETRN (N =

22) time of ~1 minute allows for the development of affordable means of design optimization software on a 

commercial laptop while Geant4, PHITS, and FLUKA codes require a multiprocessing cluster environment of a 

thousand or more processors to do a single design evaluation and would still be impractical for design optimization. 

It should be noted that source biasing, physics biasing, or other optimization strategies could significantly reduce the 

Monte Carlo runtimes shown in Table 5. However, the fully optimized run times would still be orders of magnitude 

larger than the 3DHZETRN run times and still require a multi-processor cluster computing environment. 

Table 5. Total CPU seconds required for benchmark calculations. Total number of histories ran in Monte Carlo 

simulations are shown in parentheses.

Webber SPE GCR hydrogen GCR helium GCR carbon GCR iron

3DHZETRN (N = 22) 69 952 955 953 954

Geant4 2×108 (1×1011) 1×108 (5×1010) 2×108 (3×109) 2×108 (1×109) 3×108 (3×108)

FLUKA 1×108 (4×1010) 4×107 (1×1010) 2×107 (3×109

) 6×107 (1×109) 2×107 (1×108)

PHITS 3×107 (2×1010) 3×107 (4×109

) 8×107 (1×109) 7×107 (4×108) 1×108 (1×108)

1977 Solar Minimum GCR benchmark

The GCR benchmarks utilized the Badhwar-O’Neill 2010 model [O’Neill 2010] to represent the 1977 solar 

minimum conditions (solar modulation parameter set to 475 MV). The hydrogen, helium, carbon and iron spectra 

were considered separately as external source boundary conditions. The results for the neutron and proton fluence 

spectra for an incident 1977 solar minimum hydrogen spectrum are shown in Fig. 18 at the interface depths of 0 and 

30 g/cm2

.

Fig. 18. Same as Fig. 17 but for incident solar minimum GCR hydrogen spectrum.

The corresponding MC results for secondary alpha particles had large statistical uncertainty across the energy 

domain and are not shown. The GCR hydrogen results are qualitatively similar to results obtained for the Webber 

SPE, only now 3DHZETRN agrees more closely with PHITS over the high energy range. The 3DHZETRN proton 

fluence below 100 MeV shows better agreement with Geant4. The MC computation time used for transporting GCR 

hydrogen is slightly shorter than that for the Webber SPE transport in the same geometry mainly due to the slower 

convergence of the MC solution at lower energies. The 3DHZETRN (N = 22) run time of ~15 minutes is very 

manageable whereas the MC codes take several orders of magnitude longer.

21
The neutron, proton and alpha spectra generated by GCR helium ions are shown in Fig. 19. A large 

contribution to the alpha fluence spectrum is from the incident helium ions transmitted without suffering a nuclear 

interaction. There is generally better agreement with Geant4 results than those generated by GCR hydrogen. PHITS 

and FLUKA have a tendency toward greater production of neutrons and protons below 1 GeV, whereas Geant4 is in 

better agreement with 3DHZETRN. 

Fig. 19. Same as Fig. 17 but for incident solar minimum GCR helium spectrum.

Fig. 20. Same as Fig. 17 but for incident solar minimum GCR carbon spectrum.

22
Fig. 21. Same as Fig. 17 but for incident solar minimum GCR iron spectrum.

Similar results are shown for GCR carbon and iron ions in Figs. 20 and 21. The comparisons are more 

mixed in these cases. The carbon ion results in Fig. 20 exhibit closest agreement with FLUKA and Geant4 below 10 

MeV for neutrons and with PHITS and FLUKA above 1 GeV. A similar result is found for proton spectra. This is 

likely due to differences in carbon breakup cross sections. The alpha results for 3DHZETRN agree more closely 

with PHITS after penetrating the top aluminum shield and after passing through the tissue sphere. The 3DHZETRN 

results for incident iron tend to agree more closely with the FLUKA results. It is clear that the breakup of HZE into 

light fragments requires additional study.

To gain further understanding of the light ion flux problem, the 'Z = 0, 1, 2 flux induced by GCR carbon and 

iron environments is evaluated as shown in Figs. 22 and 23. Note that 'Z refers to the charge difference between the 

primary and secondary particles. The Z = 6 flux contains the surviving carbon primaries while Z = 5, 4 are purely 

secondary fragments. It is especially clear that there are large differences in nuclear fragment cross sections between 

the codes, and that this is in part the cause of differences in the light ion flux in Fig. 20. Similar conclusions are 

drawn from the GCR iron ion induced flux shown in Fig. 23, resolving these differences is of great importance to 

shielding from these heavy ions.

Fig. 22. Particle spectra for 'Z = 0, 1, 2 at detector locations 0 g/cm2 (left pane) and 30 g/cm2 from top of tissue 

sphere (radius 15 g/cm2

) surrounded by aluminum shell (thickness 20 g/cm2) exposed to the 1977 solar minimum 

GCR carbon spectrum.

23
Fig. 23. Same as Fig. 22 but for incident solar minimum GCR iron spectrum.

Although it is difficult to favor one code over another on the basis of the veracity of the solutions obtained as 

significant differences in the light ion spectra exist, there are even greater differences in computational time required 

to evaluate fluence spectra even in this simple geometry (inhomogeneous spherical aluminum shell containing a 

tissue sphere). Still, MC methods provide an important test point for developing codes capable of meeting 

operational and especially design requirements. In addition, computational speed has proven to be important in 

spaceflight validation in LEO where the time structure of the environment can be used to test various environmental 

components as was done using ISS and the Liulin instrument [Wilson et al. 2007, Slaba et al. 2011a, 2013]. The 

main limitation on such studies remains the uncertainty in the environmental models (especially for the trapped 

environment), uncertainty in nuclear cross sections, and vehicle geometry and composition.

Conclusions

The HZETRN and 3DHZETRN radiation transport codes for space design are still incomplete in both the 

physical description and corresponding implementation into computational procedures. Even so, the current 

development appears as state-of-the-art in computational procedures for spacecraft design. The present version of 

3DHZETRN has some inherent limitations (aside from the simplified geometry) that carry over from the available 

bi-directional methods in which the main share of light ions are still treated in the straight-ahead approximation as 

coupled to the forward propagating neutrons. Only the isotropically produced neutrons are first transported in 3D 

without the coupling to the isotropically produced light ion components. Once the isotropically produced neutrons 

are evaluated, the final coupling to the light ions is implemented (without further nuclear reactions) in a full 3D 

manner. Hence, only a fraction of the light ions have a 3D treatment of transport effects of unquantified magnitude. 

Future work will require extension to complex geometries and treatment of the cross-coupling of 3D neutron/light 

ion transport. Even so, the present expanded treatment of 3D effects shows great promise in the path to a fully 3D 

transport algorithm.

Further limitations result from the current HZE nuclear database. The current database used in 3DHZETRN 

only accounts for the breakup of the projectile ions without regard to the struck nucleus breakup. Furthermore, the 

projectile fragments are assumed to have the same velocity as the initiating projectile. Adding these neglected target 

fragments and correcting the projectile fragment spectra will result in a more reliable simulation.

Acknowledgements

This work was supported by the Human Research Program under the Human Exploration and Operations 

Mission Directorate of NASA and by NASA Grant number NNX09AR20A.

24
Appendix A: Monte Carlo simulations

In this section, a brief description is given of the MC codes Geant4, FLUKA, and PHITS used to simulate 

particle spectra shown in Figs. 8 and 16 - 23. A discussion of the source, geometry, and simulation setup used for the 

codes is also given.

Geant4 version 9.4.6 [Agostinelli et al. 2003] was used in the present calculations with the physics list 

QGSP_BERT_HP. This physics list utilizes the quark gluon string model (QGS) for high energy interactions with 

nucleons, pions, and nuclei. Post-interaction nuclear de-excitation is handled by the precompound (P) model. The 

Bertini cascade (BERT) model is used for interactions below 10 GeV. A high-precision (HP) tracking model is used 

for neutrons below 20 MeV. Nucleus-nucleus collisions are represented by the native Quantum Molecular Dynamics 

(QMD) model [Koi 2008]. More information about the chosen physics list can be found at the Geant4 website 

[Geant4 2012a]. The QGSP_BERT_HP physics list has been recommended by the Geant4 collaboration for high 

energy physics applications [Geant4 2012b] and has been used, with some slight variation, in other space radiation 

related studies [Bernabeu and Casanova 2007, Hayatsu et al. 2008, Martinez and Kingston 2012, Slaba et al. 2013].

FLUKA is a MC program that contains fully integrated physics and performs the transport of elementary 

particles and ions through materials. It has the ability to transport and calculate interactions with all elementary 

hadrons, light and heavy ions, and electrons and photons over an energy range which extends up to 104 TeV for all 

particles, and down to thermal energies for neutrons [Battistoni et al. 2007]. The nuclear models are integrated into 

the software for all particles and energies, with the exception of neutrons below 20 MeV, where standard 

international evaluated data files are used. For nucleus-nucleus interactions above 5 GeV/n, FLUKA utilizes 

DPMJET-III with a modified initialization procedure, and for energies between 0.1 and 5 GeV/n, a modified RQMD 

model is used [Roesler et al. 2001, Sorge et al. 1989a, Sorge et al. 1989b]. For low energy nucleus interactions 

below 0.1 GeV/n, the Boltzmann Master Equation is used [Cerutti et al. 2006]. All models are followed by a 

generalized intra-nuclear cascade process followed by a pre-equilibrium stage and then an evaporationfragmentation-fission stage. This code has been extensively benchmarked against available accelerator and cosmic 

ray experimental data, at beam energies as low as a few MeV and as large as cosmic ray energies [Aarnio et al. 

1993; Andersen et al. 2004].

PHITS has been used extensively in space radiation transport. It simulates ionization through transport 

processes and uses the mean free path to determine when collisions occur [Niita et al. 2006, Sihver et al. 2007]. 

Ionization transport can include angle and energy straggling, showing good agreement for Bragg peaks and 

fragmentation tails [Niita et al. 2006]. Tabulated nuclear cross section data are used for neutron interactions below 

20 MeV, and evaluated data are also used for photon transport below 100 GeV and electron/positron transport below 

100 MeV [Sato et al. 2013]. Medium energy neutron and light ion induced nuclear interactions are simulated 

through the Intra-Nuclear Cascade model of Liege [Boudard et al. 2013], while high energy (>3 GeV) nuclear 

interactions involving neutrons and other hadrons are simulated through the Jet AA Microscopic transport model 

(JAM) [Sato et al. 2013]. The Japan Atomic Energy Research Institute (JAERI) Quantum Molecular Dynamics 

model (JQMD) addresses nuclear interactions induced by deuterons, tritons, helions, and alphas at energies above 3 

GeV/n and heavier charged particle induced nuclear interactions from 10 MeV/n to 100 GeV/n [Sato et al. 2013]. 

Once the nuclear interaction dynamic portion is complete, the General Evaporation Model is used to address nuclear 

de-excitation through particle emission in a statistical manner [Niita et al. 2006, Sihver et al. 2007, Sato et al. 2013]. 

The default total nucleus-nucleus cross sections used in PHITS are the same as those that have been used in versions 

of HZETRN [Niita et al. 2006, Sato et al. 2013]. 

In Fig. 8, the geometry setup was for a 30 g/cm2 ICRU tissue slab shielded in front and back by 20 g/cm2 of 

aluminum. All slabs are assumed to have infinite lateral dimensions. Particle spectra were tallied through surfaces at 

the depths indicated in the figure. In Figs. 17 - 21, the geometry setup was for an ICRU tissue sphere with radius 15 

g/cm2 surrounded by an aluminum shell with thickness 20 g/cm2

. The external source was applied uniformly onto 

the top of the sphere (covering the entire sphere) parallel with the z-axis as shown in Fig. 10. Particle spectra were 

tallied at spherical tissue detectors (radius 2.5 mm) placed along the z-axis down through the sphere. Convergence 

tests on the detector radius were performed to arrive at the final choice. In general, the chosen detector radius is 

large enough to allow reasonable tally statistics, but small enough to effectively prevent fragmentation events from 

occurring within the detector. 

For all of the GCR simulations, the source energy spectrum was sampled directly from the input spectrum 

(i.e. no source energy biasing). For the SPE simulations, the source energy spectrum was biased so that source 

sampling was performed from a uniform distribution covering 0.01 MeV to 2500 MeV. This biasing technique 

places greater emphasis on the higher energy portion of the SPE spectrum improving tally statistics for secondary 

neutrons. 

25
In Figs. 22 and 23, the simulation setup is identical to the other GCR simulations, except for the following 

changes to improve efficiency and statistical uncertainty. First, all particles produced in the simulation with 'Z > 2 

(i.e. anything lighter than Be for the GCR carbon benchmark and anything lighter than Cr for GCR iron benchmark) 

were terminated. Second, instead of covering the entire geometry, the source irradiation plane was restricted to a 

radius of 5 cm. Both modifications greatly reduced simulation run time allowing more histories to be simulated, 

thereby reducing statistical uncertainties. The first modification is justified since HZE particles produced within 

shielding receive little or no contributions from lighter mass particles [Wilson and Badavi 1986, Wilson et al. 1991]. 

The exception to this occurs at low energies, less than a few MeV/n, which receive contributions from nucleoninduced target fragmentation. These energies are not shown in Figs. 22 and 23, and MC codes generally have 

difficulty resolving this part of the spectrum due to the short residual ranges. The second modification is justified 

since HZE particles are highly peaked in the forward direction [Wilson and Badavi 1986, Wilson et al. 1991].

Convergence testing with Geant4 revealed that no particles were scored in a detector (out of 1 billion primaries) if 

the primary ion struck the geometry with a radius greater than 2 cm. The larger source radius of 5 cm was selected to

ensure even the largest production angles in HZE collisions were accounted.

Appendix B: Evaluation of :i for N = 6, 10, 14, 18, 22

The :i are herein evaluated with :0 as the top hemispheric pole so that :1 = :0, while the bottom pole is at 

:2 = -:0 and are the central rays of the two polar regions. The remaining directions denoted by :i for i = 3...N are 

symmetrically distributed about the mid-latitudes. Generally in application, the :i’s must be rotated to the body 

coordinates of the transport media. This is accomplished by a rotation matrix R(T0, I0) determined from the polar 

vector :0 relative to the body coordinate system (note, T denotes co-latitude with values on 0 to S, and I denotes 

longitude running from 0 to 2S). The regional boundaries are chosen such that

4 / i

i d N Q % ̈  8

8 , (40)

ZKHUHǻ:i denotes the i

th region, N is the number of regions and is related to the number of NT segments of 'cosTi

and the number of NI segments 'Ii plus the two polar regions such that

N N N 2   G R , (41)

and defined such that

1 1

cos

N N

i ii

i i

G R

 

  %  %% 8 . (42)

The polar caps are bound by 

1 cos 1 2 / R   N , (43)

2 cos 1 2 / R   N , (44)

and the corresponding central rays :1 and :2 lie along the rz axis. The middle latitudes are divided into NINT

regions extending from the polar cap to the equator in each hemisphere. The corresponding longitudinal boundaries 

are set by the NI such that the 

2 / i NG %  G Q . (45)

Note that NI is generally greater than or equal to 4. The remaining boundaries in latitude are given as 

26
1 2 cos (cos cos ) / i NR %  R RR . (46)

For N = 6, there is only one midpoint in latitude taken as

1 2 cos (cos cos ) / 2 0 i R RR   , (47)

for i t 2. The four longitudinal points are given as

G Q QQ  \ ^ 0, / 2, , 3 / 2 . (48)

The directional elements for each solid angle region are given by

\ ^ cos sin , sin sin , cos i il il l 8  GR GR R , (49)

and are given for N = 6 in Table B.1. Similar results for other N used in this paper are given in Tables B.2 - B.5.

Table B.1 Directional elements of the :i for N = 6.

i 123456

ex 0 0 1 -1 0 0

ey 0 0 0 0 1 -1

ez 1 -1 0 0 0 0

Table B.2 Directional elements of the :i for N = 10.

i 1 2 3 4 5 6 7 8 9 10

ex 0 0 0.9165 -0.9165 0 0 -0.9165 0.9165 0 0

ey 0 0 0 0 0.9165 -0.9165 0 0 -0.9165 0.9165

ez 1 -1 0.4 -0.4 0.4 -0.4 0.4 -0.4 0.4 -0.4

Table B.3 Directional elements of the :i for N = 14.

i 1 3 5 7 9 11 13

ex 0 0.8207 0 -0.8207 0 1 0

ey 0 0 0.8207 0 -0.8207 0 1

ez 1 0.5714 0.5714 0.5714 0.5714 0 0

i 2 4 6 8 10 12 14

ex 0 -0.8207 0 0.8207 0 -1 0

ey 0 0 -0.8207 0 0.8207 0 -1

ez -1 -0.5714 -0.5714 -0.5714 -0.5714 0 0

Table B.4 Directional elements of the :i for N = 18.

i 1 3 5 7 9 11 13 15 17

ex 0 0.7454 0 -0.7454 0 0.9750 0 -0.9750 0

ey 0 0 0.7454 0 -0.7454 0 0.9750 0 -0.9750

ez 1 0.6667 0.6667 0.6667 0.6667 0.2222 0.2222 0.2222 0.2222

i 2 4 6 8 10 12 14 16 18

ex 0 -0.7454 0 0.7454 0 -0.9750 0 0.9750 0

ey 0 0 -0.7454 0 0.7454 0 -0.9750 0 0.9750

ez -1 -0.6667 -0.6667 -0.6667 -0.6667 -0.2222 -0.2222 -0.2222 -0.2222

27
Table B.5 Directional elements of the :i for N = 22.

i 1 3 5 7 9 11 13 15 17 19 21

ex 0 0.6863 0 -0.6863 0 0.9315 0 -0.9315 0 1 0

ey 0 0 0.6863 0 -0.6863 0 0.9315 0 -0.9315 0 1

ez 1 0.7273 0.7273 0.7273 0.7273 0.3636 0.3636 0.3636 0.3636 0 0

i 2 4 6 8 10 12 14 16 18 20 22

ex 0 -0.6863 0 0.6863 0 -0.9315 0 0.9315 0 -1 0

ey 0 0 -0.6863 0 0.6863 0 -0.9315 0 0.9315 0 -1

ez -1 -0.7273 -0.7273 -0.7273 -0.7273 -0.3636 -0.3636 -0.3636 -0.3636 0 0

References

Aarnio, P.A., Fasso, A. Ferrari, A., Moehring, H.J., Ranft, J., Sala, P.R., Stevenson, G.R., Zazula, J.M., FLUKA: 

Hadronic benchmarks and applications. Proceedings of the MC93 International Conference on Monte Carlo 

Simulation in High Energy and Nuclear Physics, 1993.

Adamczyk, A.M., Norman, R.B., Sriprisan, S.I., Townsend, L.W., Norbury, J.W., Blattnig, S.R., Slaba, T.C., 

NUCFRG3: Light ion improvements to the nuclear fragmentation model. Nucl. Instrum. & Methods A 678: 21-32; 

2012.

Agostinelli, S. et al., Geant4--a simulation toolkit. Nucl. Instrum. & Methods A 506: 250-303; 2003.

Alsmiller, R.G., Irving, D.C., Kinney, W.E., Moran, H.S., The validity of the straightahead approximation in space 

vehicle shielding studies. In Second Symposium on Protection Against Radiations in Space, Arthur Reetz, ed. NASA 

SP 71: 177-181, 1965.

Alsmiller, R.G., High-energy nucleon transport and space vehicle shielding. Nucl. Sci. & Eng. 27: 158–189; 1967.

Alsmiller, R.G., Santoro, R.T., Barish, J., Claiborne, H.C., Shielding of manned space vehicles against protons and 

alpha particles. ORNL-RSIC-35, 1972.

Andersen, V., Ballarini, F., Battistoni, G., Campanella, M., Carboni, M., Cerutti, F., Empl, A., Fasso, A., Ferrari, A., 

Gadioli, E., Garzelli, M.V., Lee, K., Ottolenghi, A., Pelliccioni, M., Pinsky, L.S., Ranft, J., Roesler, S., Sala, P.R., 

Wilson, T.L., The FLUKA code for space applications: recent developments. Adv. Space Res. 34: 1338–1346; 2004.

Armstrong, T.W. and Colburn, B.L., Predictions of secondary neutrons and their importance to radiation effects 

inside the International Space Station. Rad. Meas. 33: 229-234; 2001.

Badavi, F.F., Validation of the new trapped environment AE9/AP9/SPM at low Earth orbit. Adv. Space Res. 54:

917-928; 2014.

Badhwar, G.D., Patel, J.U., Cucinotta, F.A., Wilson, J.W., Measurements of the secondary particle energy spectra in 

the space shuttle, Rad. Meas. 24: 129–138; 1995.

Badhwar, G.D., O’Neill, P.M., Cucinotta, F.A., In-flight radiation measurements on STS-60. Rad. Meas. 26: 17-34; 

1996.

Badhwar, G.D., Atwell, W., Badavi, F.F., Yang, T.C., Space radiation absorbed dose distribution in human 

phantom. Rad. Meas. 33: 235-241; 2001.

Battistoni, G., Muraro, S., Sala, P.R., Cerutti, F., Ferrari, A., Roesler, S.,Fasso, A., Ranft, J., The FLUKA code: 

Description and benchmarking. Proceedings of the Hadronic Shower Simulation Workshop 2006, 896: 31–49; 2007.

Bernabeu, J. and Casanova, I., Geant4-based radiation hazard assessment for human exploration missions. Adv. 

Space Res. 40: 1368-1380; 2007.

28
Boudard, A., Cugnon, J., David, J.C., Leray, S., Mancusi, D., New potentialities of the Liege intranuclear cascade 

model for reactions induced by nucleons and light charged particles. Phys Rev C. 87: 014606; 2013.

Cerutti, F., Battistoni, G., Capezzali, G., Colleoni, P., Ferrari, A., Gadioli, E., Mairani, A., Pepe, A., Low energy 

nucleus-nucleus reactions: the BME approach and its interface with FLUKA. Proceedings of the 11th International 

Conference on Nuclear Reactions Mechanisms, 2006.

Clowdsley, M.S., Heinbockel, J.H., Kaneko, H., Wilson, J.W., Singleterry, R.C., Shinn, J.L., A comparison of the 

multigroup and collocation methods for solving the low-energy neutron Boltzmann equation. Can. J. Phys. 78: 45-

56; 2000.

Clowdsley, M.S., Wilson, J.W., Shinn, J.L., Badavi, F.F., Heinbockel, J.H., Atwell, W., Neutron environment 

calculations for low Earth orbit. SAE ICES 01ICES2327, 2001.

Cucinotta, F.A. and Dubey, R.R., Alpha-cluster description of excitation energies in 12C(12&Į; DW $ *H9

Phys. Rev. C 50: 979–984; 1994.

Cucinotta, F.A., Wilson, J.W., Shinn, J.L., Tripathi, R.K., Assessment and requirements of nuclear reaction data 

bases for GCR transport in the atmosphere and structures, Adv. Space Res., 21(12): 1753-1762; 1998. 

Cucinotta, F.A., Personal communication, 2004.

Fasso, A., Ferrari, A., Ranft, J., Sala, P.R., FLUKA: A multi-particle transport code, CERN-2005-10, INFN/TC 

05/11, SLAC-R-773, 2005.

Foelsche, T., Mendell, R.B., Wilson, J.W., Adams, R.R., Measured and calculated neutron spectra and dose 

equivalent rates at high altitudes: Relevance to SST operations and space research. NASA TN D-7715, 1974.

Geant4 collaboration. Geant4 Reference Physics Lists. Website: 

http://www.geant4.org/geant4/support/physicsLists/referencePL/referencePL.shtml; 2012a.

Geant4 collaboration. Geant4 Recommended Physics Lists. Website: 

http://www.geant4.org/geant4/support/physicsLists/referencePL/useCases.shtml; 2012b.

Hayatsu, K., Hareyama, M., Kobayashi, S., Yamashita, N., Miyajima, M., Sakurai, K., Hasebe, N., Radiation doses 

for humans exposed to galactic cosmic rays and their secondary products on the lunar surface. Bio. Sci. in Space 22:

59-66; 2008.

Heinbockel, J.H., Feldman, G.A., Wilson, J.W., Singleterry, R.C., Leakeas, C.L., Clowdsley, M.S., Solutions to the 

low energy neutron Boltzmann equation for space applications. SAE ICES 03ICES339, 2003.

Heinbockel, J.H., Slaba, T.C., Blattnig, S.R., Tripathi, R.K., Townsend, L.W., Handler, T., Gabriel, T.A., Pinsky, 

L.S., Reddell, B., Clowdsley, M.S., Singleterry, R.C., Norbury, J.W., Badavi, F.F., Aghara, S.K., Comparison of the 

transport codes HZETRN, HETC and FLUKA for a solar particle event. Adv. Space Res. 47: 1079-1088; 2011a.

Heinbockel, J.H., Slaba, T.C., Tripathi, R.K., Blattnig, S.R., Norbury, J.W., Badavi, F.F., Townsend, L.W., Handler, 

T., Gabriel, T.A., Pinsky, L.S., Reddell, B., Aumann, A.R., Comparison of the transport codes HZETRN, HETC and 

FLUKA for galactic cosmic rays. Adv. Space Res. 47:1089-1105; 2011b.

Hugger, C.P., Qualls, G.D., Wilson, J.W., Cucinotta, F.A., Shavers, M.R., Zapp, N., ISS as a platform for 

environmental evaluation. AIP Conf. Proc. 654: 1011-1017; 2003.

ICRU, Internal Commission on Radiation Units and Measurements, Quantities and units in radiation protection 

dosimetry, ICRU Report 51, 1993.

29
Koi, T., New native QMD code in Geant4. IEEE Nuclear Science Symposium Conference Record, 2008.

Lambiotte, J.J., Wilson, J.W., Filippas, T.A., PROPER3C: A nucleon-pion transport code. NASA TM 2158, 1971.

Lamkin, S.L., A theory for high-energy nucleon transport in one-dimension. Masters Thesis, Old Dominion 

University, June 1974.

Martinez, L.M. and Kingston, J., Space radiation analysis: Radiation effects and particle interaction outside the 

earth's magnetosphere using GRAS and GEANT4. Acta Astronaut. 72: 156-164; 2012.

Nealy, J.E., Cucinotta, F.A., Wilson, J.W., Badavi, F.F., Dachev, Ts.P., Tomov, B.T., Walker, S.A., de Angelis, G., 

Blattnig, S.R., Atwell, W., Pre-engineering spaceflight validation of environmental models and the 2005 HZETRN 

simulation code. Adv. Space Res. 40(11): 1592-1610; 2006.

Niita, K., Sato, T., Iwase, H., Nose, H., Nakashima, H., Sihver, L., 2006 PHITS - A particle and heavy ion transport 

code system, Rad. Meas. 41: 1080 – 1090; 2006.

O'Neill, P.M., Badhwar-O'Neill galactic cosmic ray flux model – Revised. IEEE Trans. Nuc. Sci., 57: 3148–3153; 

2010.

Pinsky, L., Carminati, F., Ferrari, A., Simulation of space shuttle neutron measurements with FLUKA. Rad. Meas.

33: 335-339; 2001.

Qualls, G.D., Wilson, J.W., Sandridge, C.A., Cucinotta, F.A., Nealy, J.E., Heinbockel, J.H., Hugger, C., Verhage, J., 

Anderson, B.M., Atwell, W., Zapp, N., Barber, R., International space station shielding model development. SAE 

ICES 2001-01-2370, 2001.

Ranft, J., The FLUKA and KASPRO hadronic cascade codes. Computer Techniques in Radiation Transport and 

Dosimetry, W. R. Nelson and T.M. Jenkins, eds, Plenum Press, pp. 339-371, 1980

Roesler, S., Engel, R., Ranft, J., The Monte Carlo event generator DPMJET-III. Proceedings of the Monte Carlo 200 

Conference: 1033-1038, 2001.

Sato, T., Niita, K., Iwase, H., Nakashima, H., Yamaguchi, Y., Sihver, L., 2006 Applicability of particle and heavy 

ion transport code PHITS to the shielding design of spacecrafts, Rad. Meas. 41: 1142 – 1146; 2006.

Sato, T., Niita, K., Matsuda, N., Hashimoto, S., Iwamoto, Y., Noda, S., Ogawa, T., Iwase, H.,Nakashima, H., 

Fukahori, T., Okumura, K., Kai, T., Chiba, S., Furuta T., Sihver, L., Particle and heavy ion transport code system 

PHITS, version 2.52, J. Nucl. Sci. Technol. 50: 913-923; 2013.

Scott, W.W. and Alsmiller, R.G., Comparisons of results obtained with several proton penetration codes. ORNLRSIC-17. 1967.

Shavers, M.R., Zapp, N., Barber, R.E., Wilson, J.W., Qualls, G., Toupes, L., Ramsey, S., Vinci, V., Smith, G., 

Cucinotta, F.A., Implementation of ALARA radiation protection on the ISS through polyethylene shielding 

augmentation of the Service Module Crew Quarters, Adv. Space Res. 34: 1333-1337; 2004.

Shinn, J.L., Cucinotta, F.A., Simonsen, L.C., Wilson, J.W., Badavi, F.F., Badhwar, G.D., Miller, J., Zeitlin, C., 

Heilbronn, L., Tripathi, R.K., Clowdsley, M.S., Heinbockel, J.H., Xapsos, M.A., Validation of a comprehensive 

space radiation transport code, IEEE Trans. Nuc. Sci., 45: 2711-2719; 1998.

Sihver, L., Mancusi, D., Sato, T., Niita, K., Iwase, H., Iwamoto, Y., Matsuda, N., Nakashima, H., Sakamoto, Y., 

2007 Recent developments and benchmarking of the PHITS Code. Adv. Space Res., 40: 1320 – 1331; 2007.

Slaba, T.C., Blattnig, S.R., Aghara, S.K., Townsend, L.W., Handler, T., Gabriel, T.A., Pinsky, L.S., Reddell, B., 

Coupled neutron transport for HZETRN. Rad. Meas. 45: 173-182; 2010a.

30
Slaba, T.C., Blattnig, S.R., Badavi, F.F., Faster and more accurate transport procedures for HZETRN. Comp. Phys.

229: 9397-9417; 2010b.

Slaba, T.C., Blattnig, S.R., Badavi, F.F., Stoffle, N.N., Rutledge, R.D., Lee, K.T., Zapp, E.N., Dachev, Ts. P., 

Tomov, B.T., Statistical validation of HZETRN as a function of vertical cutoff rigidity using ISS measurements. 

Adv. Space Res. 47: 600-610; 2011a.

Slaba, T.C., Blattnig, S.R., Clowdsley, M.S., Variations in lunar neutron dose estimates. Rad. Res. 176: 827-841; 

2011b.

Slaba, T.C., Blattnig, S.R., Reddell, B., Bahadori, A., Norman, R.B., Badavi, F.F., Pion and electromagnetic 

contribution to dose: Comparisons of HZETRN to Monte Carlo results and ISS data. Adv. Space Res. 52: 62-78; 

2013.

Sorge, H., Stoecker, H., Greiner, W., Poincare invariant Hamiltonian dynamics: Modelling multi-hadronic 

interactions in a phase space approach. Ann. Phys. 192: 266-306; 1989a.

Sorge, H., Stoecker, H., Greiner, W., Relativistic quantum molecular dynamics approach to nuclear collisions at 

ultrarelativistic energies. Nucl. Phys. A498, 567c-576c; 1989b.

Tweed, J., Walker, S.A., Wilson, J.W., Cucinotta, F.A., Tripathi, R.K., Blattnig, S., Mertens, C.J., Computational 

methods for the HZETRN code. Adv. Space Res. 35: 194-201; 2006a.

Tweed, J., Walker, S.A., Wilson, J.W., Tripathi, R.K., Cucinotta, F.A., Badavi, F.F., An improved Green’s function 

code for HZE transport. SAE ICES 2006-01-2147, 2006b.

Webber, W.R., An evaluation of solar-cosmic-ray events during solar minimum. D2-84274-1, Boeing Co. 1966.

Wilson, J.W. and Khandelwal, G.S., Proton dose approximation in convex geometry. Nucl. Tech. 23: 298-305; 

1974.

Wilson, J.W. and Lamkin, S.L., Perturbation approximation to charged particle transport. Trans. Am. Soc., 19: 443; 

1974.

Wilson, J.W. and Lamkin, S.L, Perturbation theory for charged-particle transport in one dimension. Nucl. Sci. & 

Eng. 57(4): 292-299; 1975.

Wilson, J.W., Analysis of the theory of high-energy ion transport. NASA TN D-8381, 1977.

Wilson, J.W., Townsend, L.W., Bidasaria, H.B., Schimmerling, W., Wong, M. Howard, J., 20-Ne depth-dose 

relations in water. Health Phys. 46: 1101-1111; 1984.

Wilson, J.W. and Badavi, F.F., Methods of galactic heavy ion transport. Rad. Res. 108: 231-237; 1986.

Wilson, J.W., Townsend, L.W., Badavi, F.F., A semiempirical nuclear fragmentation model. Nucl. Instrum. &

Methods B 18: 225-231; 1987a.

Wilson, J.W., Townsend, L.W., Badavi, F.F., Galactic cosmic ray propagation in earth's atmosphere. Rad. Res. 109:

173-183; 1987b.

Wilson, J.W., Chun, S. Y., Buck, W.W., Townsend, L.W., High energy nucleon data bases. Health Phys., 55: 817-

819; 1988.

31
Wilson, J.W., Townsend, L.W., Schimmerling, W., Khandelwal, G.S., Khan, F., Nealy, J.E., Cucinotta, F.A., 

Simonsen, L.C., Shinn, J.L., Norbury, J.W., Transport methods and interactions for space radiations, NASA RP 

1257, 1991.

Wilson, J.W., Nealy, J.E., Wood, J.S., Qualls, G.D., Atwell, W., Shinn, J.L., Simonsen, L.C., Exposure fluctuations 

of astronauts exposure due to orientation. NASA TP 3364, 1993.

Wilson, J. W., Shavers, M.R., Badavi, F.F., Miller, J., Shinn, J.L., Costen, R.C., Nonperturbative methods in HZE 

propagation. Rad. Res. 140: 241-244; 1994.

Wilson, J.W., Tripathi, R. K., Cucinotta , F. A., Shinn, J. L., Badavi, F. F., Chun, S. Y., Norbury, J. W., Zeitlin, 

C. J., Hielbronn, L. H., Miller, J., NUCFRG2, an evaluation of the semiempirical nuclear fragmentation database, 

NASA TP-3533, 1995.

Wilson, J.W., Cucinotta, F.A., Tai, H., Shinn, J.L., Chun, S.Y., Tripathi, R.K., Sihver, L., Transport of light ions in 

matter. Adv. Space Res. 21(12): 1763-1771; 1998.

Wilson, J.W., Tripathi, R.K., Qualls, G.D., Cucinotta, F.A., Prael, R.E, Norbury, J.W., Heinbockel, J.H., Tweed, J., 

de Angelis, G., Advances in space radiation shielding codes, Rad. Res. 43: Suppl., S87-S91; 2002.

Wilson, J.W., Korte, J.J., Sobieski, J., Badavi, F.F., Chokshi, S.M., Martinovic, Z.N., J. Cerro, J., Qualls, G.D., 

Radiation shielding, MDO processes, and single stage to orbit design, AIAA Space 2003, Sept 22, 2003, Long 

Beach, California, AIAA 2003-6259, 2003a.

Wilson, J.W., Nealy, J.E., De Angelis, G., Badavi, F.F., Hugger, C.P., Cucinotta, F.A., Kim, M.Y., 

Dynamic/anisotropic low earth orbit environmental models. AIAA Space Conference paper 2003-6221, 2003b.

Wilson, J.W., Cucinotta, F.C., Schimmerling, W.S., Emerging radiation health-risk mitigation technologies. AIP 

Conf. Proc. 699, 913-924; 2004a.

Wilson, J.W., Tripathi, R.K., Qualls, G.D., Cucinotta, F.A., Prael, R.E., Norbury, J.W., Heinbockel, J.H., Tweed, J., 

Space radiation transport methods development, Adv. Space Res. 34(6): 1319-27; 2004b.

Wilson, J.W., Tripathi, R.K., Mertens, C.J., Blattnig, S.R., Clowdsley, M.S., Verification and validation: High 

charge and energy (HZE) transport codes and future development. NASA TP-2005-213784, 2005.

Wilson, J.W., Tripathi, R.K., Badavi, F.F., Cucinotta, F.A., Standardized radiation shield design method: 2005 

HZETRN, SAE ICES 2006-01-2109, 2006.

Wilson, J.W., Nealy, J.E., Dachev, Ts.P., Tomov, B.T., Cucinotta, F.A., Badavi, F.F., De Angelis, G., Atwell, W., 

Leutke, N., Time serial analysis of the induced LEO environment within the ISS 6A. Adv. Space Res. 40(11): 1562-

1570; 2007. 

Wilson, J.W., Slaba, T.C., Badavi, F.F., Reddell, B.D., Bahadori, A.A., Advances in NASA space radiation 

research: 3DHZETRN, Life Sci. Space Res. 2: 6-22; 2014.

32
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Publication

 4. TITLE AND SUBTITLE

3D Space Radiation Transport in a Shielded ICRU Tissue Sphere

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Wilson, John W.; Slaba, Tony C.; Badavi, Francis F.; Reddell, Brandon; 

Bahadori, Amir A.

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-20473

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 93

Availability: NASA CASI (443) 757-5802

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

A computationally efficient 3DHZETRN code capable of simulating High Charge (Z) and Energy (HZE) and light ions 

(including neutrons) under space-like boundary conditions with enhanced neutron and light ion propagation was recently 

developed for a simple homogeneous shield object. Monte Carlo benchmarks were used to verify the methodology in slab and 

spherical geometry, and the 3D corrections were shown to provide significant improvement over the straight-ahead 

approximation in some cases. In the present report, the new algorithms with well-defined convergence criteria are extended to 

inhomogeneous media within a shielded tissue slab and a shielded tissue sphere and tested against Monte Carlo simulation to 

verify the solution methods. The 3D corrections are again found to more accurately describe the neutron and light ion fluence 

spectra as compared to the straight-ahead approximation. These computationally efficient methods provide a basis for software 

capable of space shield analysis and optimization. 

15. SUBJECT TERMS

3DHZETRN; HZETRN; Monte carlo; Radiation; Transport

18. NUMBER

 OF 

 PAGES

40

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

651549.02.07.10

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/TP-2014-218530

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 09 - 2014