HYDROFLASH : A 2-D Nuclear EMP Code Founded on Finite Volume Techniques

The basic mechanisms that govern the generation of an electromagnetic pulse (EMP) following a nuclear detonation in the atmosphere, including heights of burst (HOB) relevant to surface bursts (0 km), near surface bursts (0-2 km), air bursts (2-20 km) and high-altitude bursts (> 20 km), are reviewed. Previous computational codes developed to treat the source region and predict the EMP are discussed. A new 2-D hydrodynamic code (HYDROFLASH) that solves the fluid equations for electron and ion transport in the atmosphere and the coupled Maxwell equations using algorithms extracted from the Conservation Law (CLAW) package for solving multi-dimensional hyperbolic equations with finite volume techniques has been formulated. Simulations include the ground, atmospheric gradient, and an azimuthal applied magnetic field as a first approximation to the geomagnetic field. HYDROFLASH takes advantage of multiprocessor systems by using domain decomposition together with the Message Passing Interface (MPI) protocol for parallel processing. A detailed description of the model is presented along with computational results for a generic 10 kiloton (kT) burst detonated at 0 and 10 km altitude.


Nuclear EMP
The detonation of a nuclear weapon is accompanied by the release of tremendous amounts of energy (several to thousands of tera-Joules, depending on yield) in a very short period of time (less than a microsecond) and in the form initially of X-rays, γ -rays, neutrons, and radioactive debris.Approximately 70-80% of the weapon yield is emitted as radiation (primarily X-rays but also neutrons and gammas) while 20-30% comes out as debris kinetic energy.For a detonation in sea-level air the X-rays are absorbed within a few meters of the burst and the air temperature in this volume is instantaneously raised to millions of degrees.This deposited energy then evolves over milliseconds and seconds into a hot, expanding fireball that consists of blast and shock waves, that emits intense optical and UV radiation, and that, for near surface bursts, vaporizes or entrains copious amounts of solid materials.While most of the energy produced by a nuclear detonation is dissipated during this radiative and hydrodynamic phase, very little coherent electromagnetic radiation is produced.In fact, it is only during the very early prompt phase (t < a few tens of microseconds) of the detonation, when weapon leakage γrays and neutrons are radiated and absorbed in the air, that a strong coherent electromagnetic pulse (EMP) is generated.The γ-rays comprise a small fraction of the total yield ranging from 0.1% -0.5% depending on the weapon design and have absorption scale lengths in air at sea level in the range of hundreds of meters depending on their energy.Compton scattering of the γ-rays in air leads to the production of a Compton current (mainly relativistic electrons) that in turn radiates an electromagnetic pulse.
Conversion of the prompt, unscattered γ-ray energy to EMP occurs within one microsecond of the detonation.A second phase of EMP production occurs in the time frame from 1 µs to roughly 10 ms after the detonation.During this time the neutrons produced by the detonation (roughly one percent of the weapon yield) generate γ-rays through neutron inelastic scattering and neutron capture.These secondary or delayed gammas along with the scattered prompt gammas produce EMP by means of a Compton current as described previously.A third phase of EMP production occurs for high-altitude bursts above 100 km and involves a completely different mechanism characterized by formation of a magnetic bubble in the Earth's field through magnetohydrodynamic (MHD) expansion of the conducting debris plasma created by the detonation and by distortion of the geomagnetic field resulting from heave of the ionized and heated air that originates in the x-ray deposition region at lower altitudes (~ 70-100 km).The motion of the geomagnetic field induces electric fields in the ground approximately 1 second after the burst and for as long as 100 s.Only the first two phases are addressed in this paper.A spherically symmetric radial current system produces zero net electromagnetic radiation therefore a break in the current symmetry must occur.The Earth's magnetic field, differential absorption of the γ-rays due to the presence of the ground or the decrease in air density and water vapor content with height, and weapon design asymmetries can all contribute to symmetry breaking.The turning of Compton electrons in the geomagnetic field (synchrotron radiation) produces the dominant 'fast' component of the EMP while for surface bursts the ground-air asymmetry yields the strongest 'slow' component.The peak of the latter component can be more than a factor of ten larger than that of the former.The turning of Compton electrons in the geomagnetic field and the self-consistent magnetic fields produced by electron currents in the source region together generate synchrotron radiation.The strongest fields however are generated by the vertical currents near the ground-air interface.A strong coherent pulse is emitted along the ground by these currents and this radiation is generally confined to within a few degrees of the ground.The contribution of various asymmetries to the net EMP depends on altitude.For surface bursts (HOB = 0 km) the ground asymmetry dominates for an observer near the ground as noted previously while synchrotron radiation still contributes to the high-frequency component near the ground and to the total radiated field at angles greater than roughly 5-10 degrees above the surface.For near surface bursts (HOB = 0-2 km), the ground asymmetry is a weaker component and the air asymmetry begins to take over the slow component of the EMP.Synchrotron radiation is an important contributor to the fast component of EMP throughout the altitude range of interest.For air bursts (H0B = 2-20 km) the air-asymmetry component dominates the lower frequency EMP while synchrotron radiation increases in amplitude dramatically with height.The EMP from highaltitude bursts (HOB > 20 km) is dominated by synchrotron radiation.As altitude increases it is also necessary to consider the increase in mean free paths for both gammas and neutrons (increases the volume of the source region) and the fact that inertial effects in the production and motion of secondary (lower energy) electrons become important.Xrays can also play a role in absorbing and generating EMP for high altitude bursts.Our present understanding of the physics inherent to generation of nuclear EMP was first elucidated in the initial publications of Karzas and Latter [1]- [4] and the lecture notes of Longmire [5]- [8].The early US computational codes (cf., [9]- [14]) addressed all aspects of the problem but incorporated many approximations (cf., [15]- [20]) that in part depended on the altitude regime addressed by the code.Since the 1970s additional codes have been developed with improvements made to the air chemistry model as new data has become available, to the treatment of X-rays, to the treatment of secondary electrons (including inertial effects) at high altitude, and to the numerical algorithms as computing power has increased; however, the majority of approximations still persist in available codes into this decade, see [21] for a recent review of high-altitude EMP (HEMP) and [22] for an earlier review.More recently there has been an effort to develop a 2D nuclear EMP code that uses parallel processing and modern Finite Difference Time Domain (FDTD) techniques to solve Maxwell's equations [23].We are also aware of additional work underway at various institutions in the US, Russia, and Great Britain to develop 3-D Monte Carlo and quasi-kinetic codes coupled to FDTD solutions of Maxwell's equations; however, much of this work is unpublished and cannot be addressed in this manuscript.While a self-consistent kinetic treatment is desirable in order to fully address all aspects of nuclear EMP generation at all altitudes, the 2-D hydrodynamic code, HYDROFLASH, described in this paper represents a first step towards this goal and provides an advance over previous work by establishing a robust numerical foundation rooted in finite volume techniques that more easily incorporate sharp boundaries (such as the air-ground interface and material interfaces in an urban environment) and that allow for straightforward expansion to a kinetic formulation.In addition, HYDROFLASH includes the ground and air asymmetries, an azimuthally applied magnetic field (as a rough approximation to the geomagnetic field), runs at all altitudes, and takes advantage of modern multiprocessor systems by using domain decomposition together with the MPI protocol for parallel processing.Further details regarding the physical approximations and numerical methods inherent to HYDROFLASH are provided in Section 2.

HYDROFLASH Model
As noted previously, the detonation of a nuclear device in the atmosphere leads to the immediate release of gamma rays, X-rays, neutrons, and radioactive debris.The earlytime portion of nuclear EMP is dominated by the interaction of the gammas rays with the air and the ground.Neutrons play an important role at later times but are not included in our formulation at present.X-rays are confined initially to a small volume (meters to tens of meters in extent) near the detonation and do not play a significant role at low altitudes but can be important at high-altitudes (> 20 km) in generating and absorbing nuclear EMP.X-rays are not included presently in HYDROFLASH.The radioactive debris has a small contribution during intermediate times and very late times but has not yet been incorporated in our model.Gamma ray transport in the air including reflections off the ground is modeled based on the treatment advanced by Wyatt [24].In this formulation the gamma photons are assumed to undergo full attenuation (absorption and scattering) along the ray path and then a build-up factor is applied to account for those photons that are scattered into the line-of sight from neighboring rays and off the ground.Wyatt provides analytic forms for the impulse responses for gamma ray deposition in air based on Monte Carlo calculations for mono-energetic photons over a range of energies.The impulse responses are folded against an input gamma ray pulse and weighted according to a specified energy spectrum at every mesh point in the simulation volume.HYDROFLASH uses the 2D Clawpack version 4.3 collection of routines [25]- [27] written in Fortran to solve the hydrodynamic equations in spherical coordinates centered on the burst point for four charged species and the Maxwell equations simultaneously.The four relevant constituents in the source region include Compton electrons with energies > 10 keV, secondary electrons with energy < 100 eV, positive ions, and negative ions.The appropriate hydrodynamic equations for each species can be written in general as: where n is the fluid density, u is the fluid radial velocity, v is the fluid velocity in the theta direction, f r = nmu is the fluid momentum flux in the radial direction, f θ = nmv is the fluid momentum flux in the theta direction, F Lr and F Lθ are the Lorentz forces acting on the charged fluid in the r and θ direction, respectively, m is the mass of the fluid particle, S and L are source and loss terms respectively, and S fr , L fr , S fθ , and L fθ are source and loss terms for the radial and theta fluid fluxes, respectively.Multiplying Eqs. ( 1) -( 3) by r 2 sinθ and performing a coordinate transformation into the retarded time frame that follows the prompt gamma pulse out radially from the burst point at the speed of light i.e., τ = t-r/c, r′ = r, and θ' = θ, where τ, r′, and θ' are the new time, radial, and theta coordinates, we find: where we have dropped the primes on the new coordinates, the bar designates a variable multiplied by r 2 sinθ, and β r is the fluid velocity in the radial direction divided by the speed of light (β r = u/c).In the retarded time frame it is possible to adopt a spatial resolution that is adequate for the absorption scale length (a hundred meters or more in sea level air) of the gammas while at the same time maintaining resolution in retarded time to accommodate the large radial gradients associated with the fast rise of the gamma pulse.
After moving the Lorentz force terms to the right hand side (RHS) of Eqs. ( 5) and ( 6), Eqs. ( 4) -( 6) take on a form that can be solved explicitly using standard finite-volume methods for hyperbolic equations [27], see pp. 170-171 and pp.460-463.Dimensional splitting allows one to extend one dimensional high-resolution methods to more dimensions by applying these methods along each coordinate direction (r and θ in our case).In HYDROFLASH high-resolution Godunov splitting with Van Leer limiting is invoked to solve Eqs. ( 4) - (6).Each individual fluid is modeled differently with additional approximations and with appropriate source and loss terms as follows:

Compton Electrons
The high energy Compton electron density and fluxes are evolved in time by solving Eqs. ( 4) -( 6) with the Lorentz force given as: where e is the electron charge, E r and E θ are the selfconsistent electric fields in the source region, and B φ is the self-consistent magnetic field plus applied field in the φdirection.The mass m corresponds to the relativistic mass of the electron which is equal to the rest mass times the Lorentz factor i.e., m = γm 0 with  = 1/�(1 −  2 ) , m 0 = rest mass of the electron, and β = the speed of the electron divided by the speed of light.For the Compton electron fluid we assume that the mean energy per particle is constant and equal to half the mean energy of the gamma photons emitted by the detonation.In this model, forces (both frictional and accelerating) that act on the fluid serve only to reduce the number density and/or change the direction of motion.Hence there is no explicit energy equation for the Compton electrons.This basic model is somewhat crude but sufficiently accurate because the gamma ray production of Compton electrons at a given energy is rapid and is followed by the loss of those particles to the production of low energy secondary electrons.Our intention is to replace this single fluid model with a fully kinetic treatment for all species including the gamma rays in future versions of HYDROFLASH.The source and loss terms for the Compton electrons are given by: where n c is the Compton electron density, f γ is the gamma flux in number of photons /m 2 /s, η is the fractional gamma efficiency, Y is the total yield of the device in kT, ε γ is the gamma energy in MeV, S γ is the gamma output rate normalized to one (i.e. the integral of S γ over time is equal to 1, units of S γ are MeV/s), ν E is the energy loss rate of for Compton electrons, �� is the Compton electron total speed, F D is the dynamical friction force for energetic electrons in air [28], ε c is the Compton electron energy, u 0 and v 0 are the radial and theta speeds with which Compton electrons are born, ν s is the scattering rate of Compton electrons in air, e is the charge of an electron, z (= 7.2) is the mean atomic charge of air and D (= 8.078x10 19 ) is a normalization constant to accommodate MKS units.MKS is used throughout the formulation of HYDROFLASH.

Secondary Electrons
The secondary electrons are produced primarily by impact ionization of the air by Compton electrons.These lowenergy (< 100 eV) electrons are highly collisional and are assumed to equilibrate to drift speeds and energies defined by the instantaneous electric field as described in Section 2.4.As a result only Eq. ( 4) is solved for the secondary electron density.In addition we note that the drift speed of secondary electrons is small; so that, the transport of electrons across a grid cell in the times of interest is negligible.Therefore only the time derivative on the left hand side of Eq. ( 4) is retained.The drift speed, characteristic energy, ionization rate (secondary electrons producing additional secondary electrons), and attachment rate as a function of electric field strength are derived from "swarm" measurements [28]- [29].The source and loss terms in Eq. ( 4) are given in this case by: =    +    +   (14) where ν i is the ionization rate, ε c is the Compton electron energy in eV in this case, ν E is the energy loss rate of Compton electrons as defined in Eq. ( 10), n c is the Compton electron density, α is the sum of the 2-body and 3-body attachment rates for electrons in air, and   is the electronion recombination rate [30].The factor, 34, in Eq. ( 13) corresponds to the energy loss (in eV) per ion pair produced by Compton electrons.Note that this formulation does not allow for formative lag as described in Section 2.3.

Negative Ions
The negative ions (primarily small ions, i.e., 2 − ,  − ,  3 − ,  − ,  2 − ,  3 − , and  2  − ) are assumed to equilibrate to the instantaneous electric field with a fixed mobility, defined as µ − = 2.27x10 -4 (m 2 /V/s) at sea level, that is characterized by an inverse scaling with atmospheric density [31]- [32].The ion temperature is assumed constant and equal to the air temperature.In the same way as for secondary electrons, only Eq. ( 4) is solved for the negative ion density and only the time derivative on the left hand side of this equation is retained.The source and loss terms in this case are given as: where   is the total ion-ion recombination rate that in turn is the sum of a 3-body rate ( 3 = 210 −37  6 /) and a 2body rate (  2 = 210 −13  3 / ),  + is the positive ion density,  − is the negative ion density, and all other quantities have been defined previously.The recombination rates are derived from [33].

Positive Ions
The positive ions (primarily small ions, i.e.,  2 + ,  2 + ,  + ,  + ,  + ,  2 + , and  2  + ) are assumed to equilibrate to the instantaneous electric field with a fixed mobility µ + = 1.6x10 -4 (m 2 /V/s) at sea level and with the same inverse scaling with atmospheric density [31]- [32].The ion temperature is assumed constant and equal to the air temperature.In the same way as for negative ions, only Eq. ( 4) is solved for the positive ion density and only the time derivative on the left hand side of this equation is retained.The source and loss terms in this case are given as: =    +  − +    +   (18) where all quantities have been defined previously.
Maxwell's equations in spherical coordinates and when transformed to the retarded time frame can be written: where   and   are the radial and theta electric field components,   is the φ-component of the magnetic field,   and   are the radial and θ-components of the Compton current,  is the secondary electron plus ion conductivity, ε is the electric permittivity, µ is the magnetic permeability, c is the vacuum speed of light, and we have only considered the transverse magnetic mode for the EMP problem.Eqs. ( 19) -( 21) can be further simplified with the substitution: or and the corresponding equations: F and G represent the strength (in Tesla, MKS units of the magnetic field are used) of the electromagnetic waves that flow inward and outward, respectively, times the radial distance (see Eq. 24).The magnetic permeability is assumed to be a constant both spatially and temporally while the electric permittivity changes spatially and takes a jump in moving from the air to the ground.As a result the speed of electromagnetic waves changes in the ground relative to that of the air and free space.The conductivity σ is derived from solutions of the hydrodynamic equations and the constitutive relation, j s = σE = − n s eu s + n + µ + eE − n − µ − eE, based on the assumption that the secondary electrons and ions are collision dominated.We note that the high-frequency approximation described in Section 2 is not invoked; so that, the solution to Eqs. ( 26)-( 28) applies for all time provided the source model is valid and complete.Eqs. ( 26)-( 28) can be solved using standard finite volume techniques for hyperbolic equations as described in [28], see pp. 424-434 for the acoustics problem in two dimensions.A thorough presentation of multi-dimensional finite volume techniques for solving Maxwell's equations is contained in the thesis written by Sankaran [34].A perfectly matched layer (PML) is applied at the outer radial boundary to minimize unwanted reflections back into the simulation volume.The basic method is described in [34].
The simulation volume is represented as a uniform mesh in both the radial and polar coordinates.HYDROFLASH does have the capability to accept a non-uniform grid in either or both coordinates while Clawpack also has an integrated adaptive mesh capability.Neither of these options has been exercised or thoroughly vetted.The hydrodynamic variables are all assumed to be continuous at both the radial and polar boundaries.For Maxwell's equations the boundary conditions differ from those chosen for finite difference time domain calculations.Because FVTD solutions with dimensional splitting are focused on specifying fluxes associated with waves propagating in and out of cells, the radial boundary values are chosen so that waves do not flow along or in and out of the boundary.This condition is satisfied by setting the magnetic field at the radial boundaries to zero (i.e., F = -G at the two radial boundaries).
The radial electric field is chosen to be continuous across the radial boundary.In the polar direction the theta component of the electric field is set to zero (i.e.F = G) at both boundaries as dictated by symmetry and in order to prevent an unphysical buildup of charge along the axis.The radial electric field is continuous across the polar boundaries.
A variable time step is used in the calculations and chosen so as to satisfy the Courant condition and to be a fraction (= 0.2) of the inverse of the dominant rate for production and loss of particles and momentum for all cells in the simulation volume and for all time.Parallel processing using the MPI protocol and domain decomposition is used to minimize run time on multiprocessor systems.The results of two representative cases are presented in the following section.

Results
The parameters chosen for the two HYDROFLASH runs discussed in this section are summarized in Table 1.The primary differences between the two runs are the height of burst and the physical scale of the simulation.At higher altitudes it is necessary to increase the radial extent in order to accommodate larger scale lengths that accompany the decrease in atmospheric density.The spatial resolution is decreased for the same reason; however, the number of radial cells was increased (see Table 1) in order to place an observer on the ground at a reasonable horizontal range from the burst point.Note that the water vapor content within the source region was set to 2% for Case 1 and 0 for the higher altitude run.The simulations were run out to 30 µs (retarded time) with a variable time step (< 1.7x10 -9 s).Each case was run on a 32 core server with 64 GB of random access memory (RAM) and 16 central processing units (CPUs) running at 2.9 GHz.The CPU time for Case 1 was 4009 s and 4872 s for case 2.

Case 1 -Ground Burst
The results are more readily understood in the rest frame and we have developed an interpolation algorithm for converting from the retarded frame.Because only a limited number of snapshots are output over the entire duration of the run and because the grid resolution is large compared to the spatial scale of the gamma pulse itself, both the temporal and spatial resolution are coarse in a rest frame presentation.Nevertheless it is possible to resolve the main features of the simulation results.
The spatial distributions of the Compton electron density at four times are shown in Figs.1(a)-(d).Note that the ground in all spatial plots for this case lies between polar angles θ = 90 0 and θ = 180 0 .The problem is not modeled for ranges r < R min (= 135 m in this case) in order to avoid numerical issues associated with simulating a point source and this region in the plots is blank (white).The Compton scattering produced by the gamma pulse is seen as an expansion with time of the Compton electron density as the pulse moves outward.
Because of the relatively small gamma absorption scale length, the density in the ground is limited to a very small volume defined in the simulation by the first theta angle in the ground (θ=90.5 0 ).Near the burst point these cells are extremely small and hard to see (an artifact of the spherical geometry of the grid).The focusing of the density towards higher altitudes observed at late times is associated with the air density gradient (loss rate decreases with height), the background applied magnetic field (in the φ-direction), and the theta electric fields above the ground (see Fig. 3).The peak density, N c = 1x10 14 /m 3 , occurs earlier in time at approximately 500 ns (in the fixed frame) into the simulation.The density enhancements that occur in the ground at roughly 1 and 3 km radial distance from the burst are associated with strong theta electric fields that drive electrons into the ground (see Fig. 4).
The secondary electrons are generated by the Compton electrons and indeed their density follows the same general spatial distribution.The attachment rates are also lower at higher altitudes and this fact contributes to relatively higher densities with increasing height at late times.The maximum secondary electron density also occurs at ~ 500 ns into the simulation with a value of N s ~ 3x10 18

Summary
The 2-D hydrodynamic code, HYDROFLASH, described in this paper represents a first step towards development of a fully kinetic model that can address all aspects of nuclear EMP.HYDROFLASH uses algorithms extracted from the Conservation Law (CLAW) package for solving multidimensional hyperbolic equations with finite volume techniques and provides an advance over previous work by establishing a robust numerical foundation that more easily incorporates sharp boundaries (such as the air-ground interface and material interfaces in an urban environment) and that allows for straightforward expansion to a kinetic formulation.As illustrated by the results presented for a surface burst and a detonation at 10 km altitude, HYDROFLASH easily accommodates sharp boundaries such as the ground-air interface, can include externally applied fields such as an azimuthally applied magnetic field, runs at all altitudes, and takes advantage of modern multiprocessor systems by using domain decomposition together with the MPI protocol for parallel processing.The transformation to the retarded time frame reduces the computational time by eliminating the need to carry out the simulation over the entire time for radiation to reach the spatial limits of the grid as required in a fixed frame calculation.The transformation back to the frame can be performed rapidly by simple interpolation of the retarded time results as illustrated in Figs.1-3 and Figs.5-7.
Future plans for enhancing HYDROFLASH include: 1) incorporating the neutron output of the device and neutron transport, 2) adding the energetic X-ray output, 3) solving the full Boltzmann equation for electron transport with finite volume techniques as described in [35], [28]- [29], 4) solving the time-dependent radiation transport equations for gamma rays [36], 5) incorporating additional materials relevant to an urban environment, and 6) developing non-uniform meshes to accommodate sharp boundaries throughout the grid.

3 . 2
Fig.3(c).This field drives Compton electrons into the ground.Note also that the θ-electric field develops initially near the ground and then propagates along the radial and theta directions.The fields near the ground are associated in

Table 1 :
Run parameters for two cases.