A Dynamic Fluid-Kinetic (DyFK) Model for Ionosphere-Magnetosphere Plasma Transport: Effects of Ionization and Thermal Electron Heating by Soft Electron Precipitation

In order to efficiently and self-consistently describe the transport of high-latitude ionospheric plasma into the magnetosphere, we have devel­ oped a Dynamic coupled Fluid-Kinetic (DyFK) model for describing plasma flow along a magnetic flux tube. The collision-dominated ionospheric plasma is treated with a low�speed fluid approach for altitudes between 120-1100 km� while a generalized semikinetic approach is used for the topside and higher altitudes, starting at 800 km. This paper presents a description of the new DyFK model, along with illustrative results obtained from model­ ing the effects of soft (100 e V) auroral electron precipitation using the new model. The results demonstrate the accuracy and efficiency of the new DyFK model and illustrate how soft electron precipitation provides a mecha­ nism to drive observed upflows of high-latitude ionospheric plasma along geomagnetic field lines.


INTRODUCTION
More than thirty years ago, Dessler and Michel (1966) predicted a "gentle" loss of H+ ions from the p-0lar ionosphere due to an evaporation of lighter ions which they called the "polar breeze". Axford (1968) later considered the effect of ions being accelerated by an ambipolar electric field associated with the dominant Q+ ions and the thermal electrons, and suggested that the flow would become supersonic, in analogy with solar wind models of that period.
Thus, Axford suggested calling these light ion outflows the "polar wind". Holzer (1968, 1969) used nonlinear continuity and momentum equations to calculate this polar wind. In 1970, Hoffman (1970 used Explorer 31 data to confirm the existence of the polar wind. H+ ion flow velocities parallel to the magnetic field lines of up to 10 km/sec and fluxes 1Center for Space Plasma and Aeronomic Research, The University of Alabama in Huntsville, USA 2Mission Research Corporation, Nashua, New Hampshlre, USA 492 TAO, Vol. 10, No. 3, September 1999 of about 108 cm-2 sec-1 were observed above 2500 km altitude. Dynamics Explorer I satellite measurements in the 1980's verified the supersonic nature of the polar wind (e.g., Gurgiolo and Burch, 1982). Both supersonic H+ flows (Nagai et al., 1984) and Q+ ion flows (Waite et al., 1985) were observed.
From the above early calculations to the present, simulations of the polar wind were basi cally cast in terms of whether the ion description was either moment-based or kinetic-based.
Earlier hydrodynamic computer models solved the continuity and momentum equations for ions (e.g., Holzer, 1968, 1969), and later included the effects of energy, heat, momentum transfer processes (e.g., collisions), and ionospheric photoelectrons (e.g., Holzer et al., 1971). In the mid 1980's, time-dependent hydrodynamic (Gombosi et al., 1985) and 16 moment (Ganguli et al., 1987) models were developed to simulate the dynamic polar plasma outflow. Early kinetic models, such as that by Lemaire and Scherer (1971), utilized a guiding center approximation, Liou ville mapping and quasi neutrality conditions to construct steady state models of the high-altitude collisionless polar wind. Later kinetic models considered two-dimensional features associated with restricted sources, EX B transverse velocities (Horwitz and Lockwood, 1985), EXB-driven centrifugal acceleration (Horwitz et al., 1994;Ho et al., 1997), and effects of photoelectrons (e.g., Tam et al., 1995).
The idea of constructing a model which treats the lower altitude collision-dominated part of the ionosphere by a moment-based approach while using a kinetic treatment in the higher altitude collisionless regime has been considered for many years. The earliest coupled ap proach, by Lemaire (1972) used a fluid approach for low altitudes and a kinetic approach for high altitudes. In their three-dimensional time-dependent model of the polar wind, Sojka (1989, 1997) combined a low-altitude (90-800 km) low-moment ionosphere-atmo sphere and a high-altitude (500-9000 km) higher-moment hydrodynamic model. Their model also includes a region of overlap between the two different treatments, necessary to assure reliable boundary conditions between them. Lie-Svendsen and Rees (1996) describe a model in which they join a kinetic approach (including a Fokker-Planck collision operator) for the collisionless regime with standard moment equations in the lower part.
There have been a few partial comparisons of hydrodynamic or generalized transport model results with those of semikinetic treatments, although these have been limited to the baseline polar wind and limited deviations from it. Holzer et al. (1971 ), Demars and Schunk (1992), and Leer et al. (1996) compared steady-state kinetic approaches with various moment based approaches. The only available comparison between hydrodynamic and semikinetic models which included time-dependent effects was that of Ho et al. (1993), which compared both steady-state and time-dependent versions of transport and semikinetic collisionless mod els of the evolving polar wind expansion into a vacuum. These limited comparisons suggested that for steady-state conditions, moment-based generalized transport models adequately ap proximate kinetic results. For time-dependent phenomena involving sudden bursts, the con sequent velocity-dispersion effects are inadequately described by moment-based treatments.
Obviously, the need for a kinetic description of the high-altitude regions becomes acute when such effects as perpendicular ion heating, creating conics, are to be considered. Winske and Omidi (1996), in their tutorial article on kinetic models, pointed out that one way to model the transition from collision dominated to a collisionless plasma is to link a fluid 493 solution at low altitude with an essentially collisionless kinetic solution at high altitude. They pointed out that the assumptions used to close t he higher order moment equations are not valid for most collisionless plasmas. Linking two treatments this way will involve using the results of the fluid solution as boundary conditions for the kinetic equations. Winske and Omidi also discussed that this method would conserve energy and momentum locally.
The first such developed ·coupled fluid-semikinetic model was introduced by Su et al., (1998) who used a new steady-state coupled fluid-generalized semikinetic model to investi gate the photoelectron-driven polar wind. This new coupled model used a fluid transport approach in the low-altitude portion of the simulation space (120-800 km), a generalized semikinetic treatment from 800 km to 2 RE, and a steady-state collisionless semikinetic ap proach from 2 � to 9�. Using the new model, Su et al., investigated ionospheric reaction to the photoelectrons for solar minimum and maximum conditions. The new modeling approach incorporated lower-altitude sources of thermal plasma in a self-consistent manner, as opposed to imposing artificial boundary conditions at the base of the simulation space.
However, the new model does not allow for any dynamic changes in the lower-altitude sources to be communicated to the upper regions. It allows for simulations at solar maximum or minimum (or any discrete condition between), but does not provide information on how changes are effected over time. This needed development is the subject of this paper. How ever, we first describe the salient features of this steady-state version.

THE STEADY-STATE FLUID-SEMIKINETIC MODEL
Here we briefly review the Su et al., (1998) effort to construct a steady-state coupled fluid-semikinetic model; the next section describes the· extension to time-dependence. The moment-based treatment for the low altitude portion of the model is based on the Field Line Interhemispheric Plasma (FLIP) Model, a one-dimensional time-dependent chemical-physi cal model of the ionosphere-plasmasphere which has been developed over the past fifteen years Young et al., 1980). Solving time-dependent energy, momentum, continuity, and photoelectron transport equations, the moment-based model provides ion den sities, ion and electron temperatures, and ion and electron velocities over the altitude range of the model. It should be noted that the fluid model uses a low-speed approximation which neglects stress and nonlinear acceleration terms as well as density and temperature gradients perpendicular to the geomagnetic field lines. The set of transport equations is closed by as suming isotropic electron and ion temperatures. This assumption is valid for the lower altitude collision-dominated part of the high-latitude ionosphere, as long as there are no convection driven anisotropies and toroidal or otherwise highly non-Maxwellian velocity distribution func tions (e.g., Wilson, 1994).
The generalized semikinetic section of the coupled model treats ions as gyrocenters mov ing parallel to the geomagnetic field line, while electrons are treated as a massless neutralizing fluid. This portion of the coupled model is based on a one-dimensional time-dependent treat ment first developed by Wilson et al., (1990) and Brown et al., (1991), but also including effects of various collisions (e.g., Wilson, 1998;Ho et al., 1997). It does not solve the Boltzmann equation for ions directly, but uses a particle in cell (PIC) technique in which each simulation particle represents a large number of real particles, and the movements of these simulation particles are followed. Each simulation particle is acted on by the Lorentz force, the gravita tional force, and the magnetic mirror force. The ambipolar electric field is determined through the differentiation of the electric potential, which is determined through the assumption of quasineutrality and the use of either a Boltzmann electron distribution (e.g., Wilson et al., 1990) or a more elaborate electron fluid treatment (e.g., Brown et al., 1991).
Chemical reactions and collisions between different ions and neutral species are included in the semikinetic part of the model. Reactions of ions with 0, N2, and 02, collisions with 0, polarization collisions with N 2 , self coulomb collisions, and o•-H+ coulomb collisions are handled using a Monte Carlo technique. The probability of a reaction between two particles is determined by using the thermal velocity of the simulation particles with the reaction rate coefficient or the collision cross section. The Monte Carlo technique compares randomly generated numbers with the probability of a reaction or collision occurring. If the random number falls within the range of the probability of a reaction, the number of simulation par ticles is adjusted accordingly. When a collision occurs, the simulation particle is given an impulse at a randomly chosen angle.
The coupled fluid-generalized semikinetic treatment is illustrated by Figure 1 from Su et al., (1998). Boundary conditions in the steady-state version of the model are specified for both the fluid and the semikinetic treatment at 800 km altitude ( Figure 2). A steady-state is initially obtained using the fluid treatment for a full flux tube. Bulk parameters from the top cell of the fluid zone are used to determine the distribution of particles in the bottom cell of the semikinetic   Su et al., (1995). A low speed fluid approach is used between 120-800 km, while a generalized semikinetic approach is used for the topside and higher altitudes. zone.

495
Simulation particles are injected upward into the semikinetic zone of the model according to a distribution constructed from bulk parameters obtained from the fluid zone.
(1) profile in the semikinetic zone is found by assuming a thermal conductivity proportional to T51 2 (Spitzer, 1956) and the profile is determined by heat conduction, assuming a constant electron heat flow above 800 km. Heating due to Coulomb collisions with electrons is in cluded by calculating the heating rate ( r) due to electrons. During each time step the simula tion ions are given a random impulse in perpendicular and parallel energy from the Gaussian where Q is heat flux, n is density, vis flow velocity, C. and C. are the collisional heating and er m cooling rates for the ions with electrons and neutrals respectively, land q are production and loss terms, ll s is the cross sectional area of the flux tube, Tis temperature, and KT 512 is the I I ion thermal conductivity, with (Banks and Kocharts, 1973) K; = l.2X10 4 [n(o+)+2n(He+)+4n(H+)] n electron ( 3 ) in units of eV cm-1s-1K-712• The moment-based part of the model is then run to a steady-state using a truncated flux tube, one which extends only to 800 km altitude.
This process of passing boundary conditions back and forth between the. fluid and semikinetic parts of the model is repeated until the bulk parameters show good agreement at the interface between the two models.

THE DYNAMIC FLUIDwSEMIKINETIC (DyFK) MODEL
In order to generalize the steady-state version of the coupled model into a time-dependent model or Dynamic Fluid-Semikinetic (Dy PK) Model, a suitable scheme for advancing each treatment component must be developed. Initially the simulation will be run to asymptotic steady-state for the given input geophysical conditions, but once a time-dependent process is initiated, one must not allow a steady state to be achieved separately in either the moment based or the semikinetic zone. Consideration needs to be given to how the boundary condi tions at the altitude where the two different treatments are coupled together will be treated in a Initially, one might consider allowing the interface between the two treatments to remain fixed at 800 km altitude but not pinned at one specific value. Then, by iterating back and forth between the two treatments, allowing each iteration to last only a single time step, the bulk parameters at 800 km altitude will have changed, and these changes can then be passed be tween the two treatments. This method is the same as that used by the steady-state version of the model, except one does not allow either treatment to achieve a steady-state condition.
Keeping in mind that the goal is to communicate c;lynamic changes between the two treat ments, one does not wish to restrict the values at the interface. However, by allowing the parameters at the interface to change, the connection between the two different treatments is lost ( Figure 3). A method is needed which will allow for dynamic changes in the parameters to happen while maintaining a smooth connection between the two different model treatments.
It might be argued that by decreasing the magnitude of the time step one will minimize the differences in parameters between two iterations, resulting in a smoother connection between the two model approaches. However, there is a limitation on the amount that the time step can be decreased. Experimenting with various time steps it was found that in the low speed fluid part of the model there is a minimum time step for which the velocities and temperatures remain stable" (in other words, they do not need multiple iterations to achieve a stable value).
In the semikinetic regime, the time step must be smaller than the collision times which vary * When the time step becomes too small the model does not iterate to a solution within a reasonable amount of time. This is an artifact of the numerical method used to solve the fluid equations (see A useful solution is to allow for an "overlap" region between the two different model treatments ( Figure 4). By overlap, one should understand that the fluid part of the model will extend to 1100 km altitude, but boundary conditions will be passed to the semikinetic part of the model at the 800 km interface. Similarly, the semikinetic part of our approach will extend down to 780 km altitude while passing values at 1100 km altitude. The boundary values will be kept constant while the remainder of the overlap area is allowed to evolve dynamically.
This maintains a connection between the two different modeling regions without limiting the dynamic behavior of either part of the model. semikinetic part of the model, using a four second time step.

SOFT-ELECTRON PRECIPITATION EFFECTS ON IONOSPHERIC UPFLOWS
In this initial DyFK-based report, we will use the code to illustrate effects of soft-electron precipitation in driving ionospheric upflows.
Satellite spectrometers have detected Q+ fluxes of up to 109 cm·2 s·1 flowing along geo magnetic field lines from the high-latitude ionosphere into the magnetosphere (e.g., Shelley and Collin, 1991). Seo et al., (1997), in analyzing DE�2 observations for the topside iono sphere, found a relationship between field-aligned ion upflow velocities and fluxes and the characteristic energy of precipitating soft electrons (less than 1 keV). The observations of Seo et al., (1997) generally supported the notion that soft-electron precipitation effects are major drivers for high-latitude ionospheric upflows.
Thermal velocities of Q+ are normally low enough that the ions are gravitationally bound, resulting in fluxes of only about 105 cm2 s·1 (e.g., Horwitz et al., 1994). Soft (less than 150 eV) electron precipitation results in increased ionizatioI). and electron heating at higher altitudes in the ionosphere. The resulting increased ion densities and enhanced ambipolar electric field associated with the thermal electron heating increases the amount of ions flowing into the magnetosphere from the high-latitude ionosphere. The increase in o+ density could also en hance the accidentally resonant charge exchange between H-0+, resulting in an enhanced H+ density and flux.
Modeling efforts have provided further evidence that soft electron precipitation is a major contributor to the observed increased outflows of Q+ in the ionosphere. Brown et al., (1995) used a time-dependent semi-kinetic treatment to model the effects of 100 eV precipitating electrons. Their results showed a pulse in the o• drift velocity which flowed upward, as well as an enhanced H+ density. Richards (1995) used a fluid approach to model the response of the high-latitude ionosphere to 300 eV precipitating electrons. The electron precipitation re sulted in immediate increases in H+ and Q+ fluxes. Richards explained the flow as responding to a sudden increase in temperature, noting that the H+ upflow resulted from increased ion production. Liu et al., (1995) also used a fluid approach and was able to reproduce observed (DE-2) o+ fluxes of 3x109 cm·2 s·1 by introducing soft electron precipitation of about 100 eV energy and 3 ergs cm·2 s·1 energy flux in the high-latitude ionosphere. Caton et al., (1996) in modeling European Incoherent Scatter Radar (EISCAT) observations produced results indi cating that soft electron preCipitation play a dominant role in driving F-region and topside ion outflows.
Most recently, systematic fluid-based modeling of soft-electron precipitation effects was carried out by Su et al., (1999). Using a transport model and a two-stream auroral electron model they distilled the results of numerous simulations, varying the electron energy fluxes and characteristic energy levels. Their results indicated that the largest o+ upflow velocities and increases in electron temperature in the topside are to be found when the characteristic precipitating electron energy is lowest (30 e V) and when the energy flux was held constant. This results, at least partially, from the fact that more electrons are available for constant energy flux to enhance the electron temperature through heating, and thus the ambipolar elec-tric field. It was also found that an increase in the electron energy flux (while keeping the characteristic energy constant at 125 e V) resulted in a linear increase in o+ upflow at the topside. In the DyFK model described here, precipitating electrons are modeled using a two-stream energy loss code developed by Richards and Torr (1990). We designate a characteristic elec tron energy and a total energy flux. The model accounts for effects of ionization, backscatter, secondary electron production, and thermal electron heating, returning the o+ production rate due to the precipitating electrons as a function of altitude.

MODEL RESULTS
Here, the new DyFK model is used to investigate the response of a high-latitude flux tube to time dependent changes in the fluid zone induced by ionization and heating by soft electron precipitation. With the DyFK model one is able to produce the time evolution of o+ and H+ velocity distributions as well as the time evolution of bulk parameters.
A magnetic flux tube at 80.4 degrees latitude and 1830 local time is allowed to achieve a steady state. The local time is kept constant to eliminate the diurnal effects of changing solar radiation. Once a steady state has been achieved, soft electron precipitation is turned on with an average energy of 100 e V and energy flux at 800 km of 3 ergs cm-2 s-1• The precipitation is allowed to continue for one hour geophysical time. The electron precipitation produces an increase in the heating of the electrons in the F-region, resulting in an increase in the ambipolar electric field, increasing the o+ upflow velocities. Localized density enhancements of the ions, due to increased production in the F-region, also occur and are expected to propagate upwards along the magnetic field lines.   upwards, those above 1000 km all beginning at the start of the precipitation event. These enhancements can be attributed to both increased ionization in the F-region due to the precipi tation, and to local production as the Q+ enhancement moves upward. Between 600 km and 1000 km it appears that density enhancements begin later as one moves down in altitude. For lower altitudes, below 600 km, the density decreases as the precipitation event continues. One can also see somewhat localized velocity enhancements which propagate upwards at around 0.3 km/sec. Figure 8 shows H+ phase space distributions at 3300 km altitude for 160 sec., 600 sec., 1528 sec., 2080 sec., 2720 sec., and 3360 sec. The formation of a downward tail on the H+ distribution can be seen 600 seconds into the precipitation event. Figure 9 shows H+ phase space distributions at 4900 km altitude at the same times. A downward tail can be seen form-

20
0t·· · · · · · ----· · · · · .. ·· ·-j   Note the tail which forms at 1528 sec into the precipitation event (later than at 3300 km altitude). The enhanced distribution near the core of velocity space results from H•-o• collisional drag associated with the enhancement of o• density at these altitudes.
ing 1528 seconds into the precipitation event, becoming much more apparent later. The popu lation of the "core" region of velocity space is due to H+-Q+ collisional drag in the presence of enhanced Q+ during the event.

SUMMARY AND CONCLUSIONS
In this report we have presented results from DyFK simulations of the effects of soft electron precipitation on high-latitude ionospheric plasma transport. The DyFK model results suggest that soft ( 100 e V) electron precipitation provides a mechanism, through enhanced ionization and electron temperatures in the F-region, which can drive the upflow of high lati tude ionospheric plasma along the geomagnetic field lines. Model results indicate Q+ fluxes which are well above the normal values of 105 cm-2 s-1 (e.g., Horwitz et al., 1994), and are in accord with satellite observations of Q+ ion fluxes of magnitudes up to 109 cm-2 s-1 (Shelley and Collin, 1991). Good agreement in the behavior of the o+ bulk .parameters is found with the results of other modeling studies which used a strictly fluid approach (e.g., Richards, 1995;Liu et al., 1995;Caton et al., 1996). H• ions were found to respond to the increase in o+ through the accidentally resonant charge exchange reaction Q+ + H ¢:::? H+ + 0, resulting in H+ density enhancements. H+ is also found to respond to drag due to o+-H+ collisions, as evi denced by the formation of a tail on the H+ velocity distributions.
The current preliminary simulations demonstrate the newly developed DyFK model, which will provide an exciting new tool for efficiently modeling the dynamic transport of high-lati tude plasma in the ionosphere-magnetosphere. It allows for fully self-consistent plasma and energy coupling with the ionosphere but with full descriptions of non-Maxwellian distribu tions at topside and higher altitudes. The ability to generate time series of bulk profiles and distribution functions will allow for the modeling of various types of geophysical situations of interest. These include modeling the polar wind as driven by time-dependent convection and photo-electron processes, auroral plasma outflows influenced by soft-electron precipitation and wave-driven topside ion heating, and plasmasphere refilling processes on closed flux tubes.