2-3-6 Ionosphere Simulation
SHINAGAWA Hiroyuki
Space weather forecasts are about to enter a stage incorporating numerical forecasts basedon full-scale numerical simulation, in addition to conventional methods used by forecasters tomake predictions based on observational data and experience. NICT is currently running a realtime integrated simulator for space weather developed under the leadership of the Space Environment Group. This system consists of three real-time simulators used for the sun and solar
wind, the magnetosphere, and ionosphere, respectively. In the system, the simulation data is
made available to the public immediately. This paper gives an overview of the ionospheric model
used in this system and presents an example of real-time ionospheric simulation results.

Keywords
Space weather, Ionosphere, Numerical forecast, Simulation, Magnetic storm, Real time

1 Introduction
The atmospheric region at an altitude ofmore than 80 km is called the upper atmosphere and is a mixture of neutral gas and ionized gas consisting of ions and electrons. The
thermosphere is an extremely thin atmosphericregion with density less than 1/100,000 of theatmosphere at the ground surface. The ionosphere has even lower density by several digits
and is so thin that it can almost be called avacuum. However, it is an important region
for human activity. For example, communications and radio broadcasts using short waves
are conducted through reflection in the ionosphere. If the ionosphere becomes unusually
ionized through X-rays resulting from a solar
flare or falling energetic particles from magnetic storms, it will disrupt communications
and broadcasts, perhaps even affecting ourlives. Moreover, when a severe magnetic
storm occurs, significant energy flowing into
the ionospheric region in the form of auroral
particles causes strong heating and expansionof the thermosphere, resulting in a dramaticincrease of the atmospheric density. In such an
event, a low-altitude earth orbiter at about 400

km or lower will experience difficulty in controlling its orbit because of exposure to strong
atmospheric drag, and may even fall into theearth in the worst-case scenario. By understanding the current conditions of the upper
atmosphere and predicting the future state,appropriate measures could then be taken toprevent or alleviate damage. At the moment,
certain research institutes in the world aredeveloping numerical prediction models forvarious ionospheric disturbances.

2 Ionospheric simulation
2.1 Development of the ionosphericmodels
Theoretical studies on the ionospherebegan in the 1960s. With the development ofcomputers, realistic simulation models werefirst built around 1970. The first model—thevertical one-dimensional diffusion model—
was developed by Schunk et al. at Utah State
University, who successfully reproduced thebasic structure and daily changes of the midlatitude ionosphere for the first time［1］. Laterin the 1980s, global three-dimensional models
were developed by several groups in the

SHINAGAWA Hiroyuki

199

USA［2］. This made it possible to examine theoretically how polar region disturbances associated with magnetic storms affect the global
ionospheric environment［3］–［5］.
The ionosphere is characterized by the factthat it is heavily affected by the magnetospheric electric field, high-energy particles(auroral particles) and waves propagatingfrom the upper and lower atmospheres. In thefirst model, the effects of the thermosphereand magnetosphere were given from an empirical model. The 1990s witnessed the development of both a magnetospheric model and athermospheric model. The first magnetosphere-ionosphere-thermosphere coupledmodel was developed around 2000［6］. In thepresent century, the models in each region
were further improved in terms of precision,
while efforts were initiated to develop a modelof solar surfaces and solar winds, and a “solarterrestrial environmental integrated model”
that includes all regions from the lower atmosphere to the thermosphere［7］.
However, numerical models based on thefluid equations cannot, in principle, represent
the effects of microscopic phenomena such as
turbulent flows and the particle-level effects
such as auroral particles. For that reason,
empirical input data, empirical equations, andother parameters were introduced into themodel with their effects being taken in. Such atechnique (called parameterization) is commonly used in many global models. To
increase the precision of numerical models, it
is essential to compare observational data andmodels, and enhance the technique of parameterization.
2.2 Method of ionospheric modeling
Many ionospheric models have beendeveloped by research institutes in various
parts of the world, along with various types ofequation systems and methods of modeling,ranging from very simple ones to those highly
rigorous. Here, the most widely used and relatively simple method is presented. The NICT
model is also based on the method presented
below.

The ionosphere can be represented by afluid equation consisting of ions and electrons.
The general equation system used in ionospheric models is shown below.
Continuity equation…………………（1）

Ion momentum equation…………………（2）Electron momentum equation………………（3）Ion energy equation………………（4）Electron energy equation……………………（5）where, n denotes density, v velocity, P rateof particles generated, L loss rate of particles,m mass, vab momentum transfer collision frequency of “a” particles with regard to “b” particles, p pressure, g gravity acceleration, qheat flow velocity, and Q net heating rate. Suffixes i, n and e denote i-type ions, n-type neutral particles, and electrons, respectively, andsuffix t also represents all particles.
In the continuity equation (1), the production of ions involves photoionization where
the solar extreme ultraviolet radiation (EUV)
directly ionizes neutral particles, together withthe contribution of secondary ionization where
ionization-induced electrons (photoelectrons)
ionize the surrounding neutral atmosphere.
Since a rigorous calculation of ionizationinduced by photoelectrons is rather complex, asimple equation of approximation is common-

200 Journal of the National Institute of Information and Communications Technology Vol.56 Nos.1-4 2009

ly used. Moreover, there are many chemical
reactions between ions and neutral particles,
so that a certain ion produces another kind of
ion also included in the generation. Ions disappear through an ion chemical reaction, but dissociative recombination eventually combines
ions with electrons, thereby turning them backinto neutral particles and disappearing. The
photochemical process is predominant below
200 km in altitude, while in the F layer of theionosphere above 200 km, the term of transportation in the first term on the right side inequation (1) is important. In that region, thediffusion process of ions and the effects of
convection by the electrical field are predominant.
In the ion momentum equation (2), temporal changes in the first term on the left sideand the inertia term in the first term on theright side are relatively small and therefore
ignored at times, but such changes are not necessarily small in the upper ionosphere. Moreover, in strict terms, a viscosity term for ions
may be added. In the electron energy equation(3), temporal changes and the term of inertia
can normally be ignored.
In the equations for ion and electron energy (4) and (5), local temporal changes in thefirst term on the left side are small and may beignored; these changes can be expressed in theform of balance between the term of heating
and cooling, and the term of thermal conduction. As an uncertain factor, there is a heat flux
given as an upper boundary condition.
Neglecting to add this condition is known toresult in unrealistic ion or electron temperatures, so that an appropriate value is given to
match the observations.
In ionospheric simulations, the basic procedure involves solving these equations, but inthe case of the earth’s ionosphere, the intrinsic
magnetic field is sufficiently strong at ionospheric altitudes, and one can safely assume
that the form of magnetic lines remains
unchanged. In solving the momentum equations, therefore, the method described below isused in many cases for the sake of convenience.

First, ignore the term of temporal changes
in the momentum equation and divide thevelocity into the component parallel to themagnetic lines (suffix //) and the vertical component (suffix⊥).
………………………………（6）The motion of ions perpendicular to the magnetic lines can be expressed as follows:……………（7）This can be solved with regard to vi⊥ to determine the velocity.
On the other hand, the motion parallel tothe magnetic lines can be expressed withregard to electrons and ions as follows:
…（8）

…（9）

Delete E// from that and one will obtain thevelocity. Now, think of a simple case andassume that ion i is the main ion. Using theapproximations ni = nc, vin >> ven, and mi >>me, one can obtain the result as follows:……（10）The motion parallel to the magnetic lines
can be expressed in the form of a diffusion
equation in a magnetic field flux tube, and the
direction perpendicular to the magnetic lines
is expressed as motion E×B to which an electrical field is added. In that case, E becomesthe convection electrical field of the magnetosphere at high altitudes and the dynamo electrical field at mid-to-low altitudes.

3 Real-time simulation model ofthe ionosphere at NICT
3.1 Real-time ionospheric simulation
NICT has developed an ionospheric model
on its own and is now conducting ionosphericsimulations on a real-time basis. This model is

SHINAGAWA Hiroyuki

201designed to solve the three-dimensional basic
equations of the ionosphere mentioned above
in real time, in order to determine ionosphericquantities. Background neutral atmosphericstructure is given by an empirical model. The
ionospheric model includes a thermosphericmodel, which is designed to determine perturbation quantities of the neutral atmosphere,
and can express thermospheric changes. In
addition, the effects of the magnetosphere areincluded by using the electric conductivity andpotential of the polar regions as obtained fromthe real-time magnetospheric model operating
at the same time. The real-time magnetospheric model uses real-time solar wind dataobtained from the ACE solar observation
satellite.
As such, real-time solar wind data is used
as input, so that the system can reproduce thepresent conditions of the ionosphere and thermosphere. This system has been in operation
since July 2008, and information obtained by
the system is available to the public on theWeb. The following describes an overview ofthis ionospheric model.
3.2 Ionospheric model
The coordinate system employs a sphericalcoordinate system (vertical, latitude, longitude) for both the thermosphere and the ionosphere. The vertical extent of the neutral
atmosphere ranges from 0 km (ground surface) to 600 km with the 10-km grid interval.
The ionospheric model ranges from 0 km
(ground surface) to 3000 km, including the
10-km grid interval up to an altitude of
600 km, and above that, the interval gradually
increases until becoming 100 km near theupper boundary. The grid points in the directions of latitude and longitude are set at intervals of 1 and 5 degrees, respectively. A spherical coordinate system used on global simulations is generally disadvantageous in that theinterval between grid points close to the poles
narrows rapidly, but it is not necessarily a disadvantage in ionospheric simulation because
upper atmospheric phenomena are often
accompanied by small-scale phenomena suchas auroras at high altitudes. However, the grid
points narrow down more than necessary near
the poles, so that regions within 3 degrees
from the poles are deprived of grid points andthe values are calculated by interpolation.
The ionosphere involves five types of ions+(N2 , O2+, NO+, O+, H+). The method described
is used regarding the equation of motion, andthe density ratio of O+ and H+ is empiricallygiven and regarded as a one-fluid ion divided
into two velocity components (diffused in themagnetic line direction and drifted by E×B in
the direction perpendicular to the magneticline), and the velocity is determined based onequations (7) and (10). Moreover, photochemical equilibrium is assumed to determine thedensity of molecular ions.
Regarding the equation of energy, the ionsare designed not to be within the term of temporal changes on the left side of equation (3)with thermal conductivity excluded. For electrons as well, the term of temporal changes on
the left side of equation (5) are ignored and, inthe lower ionosphere (at about 150 km or
lower), an exchange of energy occurs throughcollisions and photoelectron-induced heating.
In the upper ionosphere (at about 250 km or
higher), thermal conduction alone is set to beeffective and given by an analytical solution,and the solutions for the upper and lower
ionospheres are smoothly connected. The
shape of the magnetic field is approximatedwith the dipole field. Moreover, a dynamo
electric field generated by the neutral wind isincluded using a separately built dynamo
model［8］.
Model input utilizes a solar EUV spectral
model EUVAC［9］ based on F10.7 using solar
radio intensity as a proxy. The ionization rate,heating rate, and electric potential in the polar
regions are given by an empirical formulabased on the output data of the magnetospheric model.
3.3 Thermospheric model
The ionospheric simulator incorporates aneutral atmospheric model using an empiricalmodel as a basic state. This part of the model

202 Journal of the National Institute of Information and Communications Technology Vol.56 Nos.1-4 2009can be replaced with another atmosphericmodel. This model is characterized by the fact
that, in order to solve small-scale phenomena
temporally and spatially with high precision,
the compressible fluid equation system for anonhydrostatic equilibrium is used instead of amodel of hydrostatic approximation, which isused in many global atmospheric models.
The general equation system used in neutral atmospheric models is shown below.
Continuity equation…………………………（11）Neutral momentum equation……（12）

Neutral energy equation……………（13）Equation of state……………………………………（14）where, ρ denotes mass density, vn velocity,p pressure, R gas constant, T temperature, Cv
specific heat at constant volume, μ molecularvelocity, κ thermal conductivity coefficient, Q
net heating rate, and Ω angular velocity ofrotation.
The atmosphere is designed as single fluidwith its composition consisting of three elements: O, N2, and O2. For the basic field, theatmospheric density and temperature are givenby empirical model NRL-MSISE00［10］, while
the ratio of composition is calculated based on
the composition ratio of the basic field［11］.
The basic field for velocity is designed asgiven by empirical model HWM93［12］. In
actual calculations, this basic field is assumed
to be in a steady-state equilibrium, and theamount of change (such as ρ´= ρ−ρ0 ) fromthat state is solved. The method of numericalsolution is the finite difference method, while
the advection terms are solved using the CIP
method.
3.4 Coupling to the magnetosphere
The response of the ionosphere to magnetospheric inputs and feedback processes fromionospheric variations to the magnetosphere(i.e., the process of interaction between themagnetosphere and ionosphere) remain
unsolved and are a major challenge. Numerical simulation is an effective means ofaddressing this issue. However, in terms ofmodel connection, the issue is to find an effective way to exchange mutual information. Of
particular importance is how to give electric
conductivities of the polar ionosphere corresponding to ionization caused by auroral particles falling from the magnetosphere. Our magnetospheric model is based on a sphericalshell-shaped ionosphere and projects the parameters of the magnetosphere’s inner boundaries (r = 3.5 Re) onto the ionosphere. The
electric conductivity of the ionosphere isgiven as a height-integrated value, and is supposed to be determined by three components:(1) an electric conductivity component corresponding to solar extreme ultraviolet radiation(EUV), (2) a component corresponding to diffuse aurora (and set to a function of plasmapressure and mass density on the plasma sheet
in the magnetospheric model), and (3) a component corresponding to discrete aurora (andset to a function of field-aligned currents inthe polar regions of the magnetosphericmodel). This function is empirically given to
reproduce observations as effectively as possible.
Using the height-integrated electric conductivity, altitude distribution of the electric
conductivity is determined and used as threedimensional data in the ionospheric model.
This electric conductivity is then used to estimate the ion production rate, the plasma heating rate, and the neutral particle heating rate.
This series of procedures is called parameterization—a method necessary for entering

SHINAGAWA Hiroyuki

203microscopic effects into a global model.
The electric potential of the polar regions
is another important parameter for the ionosphere. Since it is directly given by the magnetospheric model, this parameter is used as
input for the polar electric field of the ionospheric model. For ionized components associated with EUV on the dayside, input parameters are calculated in the ionospheric model.
With this method, using the solar wind as an
input enables the execution of a model that
reproduces ionospheric and thermosphericconditions on a real-time basis.
3.5 Example of real-time simulation
results
Here we introduce an example of calculations obtained from the real-time model.
Figure 1 shows the actual results obtained
from a real-time simulation on July 23, 2008.
It shows the electron density and plasma
velocity at an altitude of 120 km. With a risein the southward component of the solar windmagnetic field, disturbances occurring in themagnetosphere result in disturbances in theionosphere as well. The night side (bottom) of
the figure shows a crescent-shaped region
with high electron density; this is an “aurora”
reproduced on a computer. The bright region

Fig.1 Electron density distribution andplasma velocity at an altitude of
120 km as obtained in real-time
ionospheric simulation
The crescent-shaped bright region on thenighttime side (bottom) represents the aurora
in the simulation.

at the top is due to the effects of ionization
stemming from solar ultraviolet radiation on
the dayside. Given the insufficient resolutionof the model, the aurora microstructure cannot
be reproduced, but its overall structure agrees
well with the representative observation
results shown in Fig. 2. The convection velocity of ionospheric plasma excited in the electric field of the magnetosphere exceeds
1 km/s, thereby indicating the occurrence offast convection.
Figure 3 shows increases in temperature
and wind speed in the thermosphere at an altitude of 300 km at the same time of day as inFig. 1. Both figures show the amount of deviations from the atmosphere in a steady state.
When a magnetic storm occurs over the polar
regions, the energy of the solar wind flows
into the upper atmosphere, and dramaticallyincreases thermospheric and ionospheric temperatures. This heating enhances the thermospheric pressure in the polar regions, thereby
driving strong winds from the poles to flow
toward the equator. These winds elevate theionosphere, resulting in increased electron
density at mid-latitudes. Such a disturbanceFig.2 Auroral image taken in outer spacefrom an artificial satellite (courtesy:
Space Sciences Laboratory, University of California, Berkeley)This is an image taken on a day other thanwhen the image in Fig. 1 was taken, but does
show the characteristic aurora structure.

204 Journal of the National Institute of Information and Communications Technology Vol.56 Nos.1-4 2009

involving increased electron density is called apositive storm.
As the thermospheric atmosphere becomesheated, it also expands, thereby increasing theneutral atmospheric density in high-altitude

Fig.3 Increases in temperature (color con-tours) and wind velocity (arrow) in thethermosphere at an altitude of 300 km
in the Northern Hemisphere
Both represent the amount of changes from theatmosphere in a steady state. The red regions
represent the heated portions, along with risingatmospheric density. Thermospheric winds flow
toward the equator from these heated regions.
The center represents the North Pole; the top
represents the daytime side.

regions. In the event of extreme disturbances,
atmospheric drag will alter artificial satellite
orbits, and in the worst cases cause a satellite
to fall from orbit. Such drag may even change
the distribution of space debris, which is now
regarded as a serious problem. Expansion of
the thermospheric atmosphere is an extremely
important phenomenon in terms of space
weather as well as in other respects.
Figure 4 shows an example of positive
storms actually reproduced in the ionospherethrough real-time ionospheric simulation. This
disturbance occurred on June 25, 2008. This
event is an example of heating that occurs in
the polar regions, resulting in thermosphericwinds directed toward the equator, and causing an increase in the total electron content(TEC) of the ionosphere. The disturbance isrelatively small in scale, but when comparedwith the trend of changes in TEC over
Kokubunji as observed with GPS satellites,
the model reproduces a trend where its TEC
increased around the end of the 14th andbecame higher than the average on the 15th,
thus showing that the simulation agrees wellwith the observations.

Fig.4 Changes in the total electron content (in 10 16m−2) in the ionosphere on days of ionospher-

ic disturbance (June 14 and 15, 2008) (red line) near Tokyo(Top) Calculation results of the real-time space weather simulator (bottom) total electron content of the ionosphere obtained in observation. The dotted line represents the variation during quiet days in the ionosphere. The
model reproduces the enhancement in electron density compared with that on quiet days from the end of June
14 to June 15.

SHINAGAWA Hiroyuki

205

4 Conclusion
This paper gave an overview of the ionospheric model and related techniques. At present, continual operation is conducted on thereal-time space weather integrated simulator
of NICT, while comparing ionosphericchanges with various observational data on theionosphere, and thereby validating the magnetosphere-ionosphere coupling. Although thissimulation needs to be further improved inquantitative terms, ionospheric changes in thesimulation agree fairly well with observations
regarding many disturbances stemming fromthe magnetosphere. This leads us to consider
the possibility of developing a numerical prediction model for the ionosphere covering aperiod of about several hours. Moreover,ionospheric disturbances stemming from thelower atmosphere are also important at mid to
low latitudes. To adopt such disturbances, it isnecessary to use an extended atmosphericmodel including meteorological phenomena as

well. To that end, efforts are now under way todevelop a model［13］ and compare it with［15］.
observational data［14］Coping sufficiently with space environment disturbances in practical applications
often requires predictions half a day to about
several days in advance. To that end, it isindispensable to perform numerical predictions not only on the magnetosphere andionosphere but also on the solar wind on the
upstream side. Since 2008 we have been operating a real-time simulator for the sun and thesolar wind as well. The data obtained is now
being compared with solar wind data obtained
from the ACE satellite in low-earth orbit toverify the model. Since it takes a few days forthe solar wind to reach the earth’s orbit fromthe solar surface, it seems theoretically possible to predict space weather a few days inadvance with a degree of practical precision
by enhancing the precision of the simulator forthe sun and the solar wind, and integrating itwith a model of the magnetosphere, ionosphere and atmosphere.

References
01 Schunk, R. W. and J. C. G. Walker, “Theoretical ion densities in the lower ionosphere,” Planet. Space
Sci., Vol. 21, pp. 1875–1896, 1973.

02 Sojka, J. J. and R. W. Schunk, “A theoretical study of the global F-region for June solstice, solar maximum, and low magnetic activity,” J. Geophys. Res., Vol. 90, pp. 5285–5298, 1985.

03 Fuller-Rowell, T. J. et al., “Interactions between neutral thermospheric composition and the polar ionosphere using a coupled ionosphere-thermosphere model,” J. Geophys. Res., Vol. 92, pp. 7744–7748,
1987.

04 Roble., R. G., E. C. Ridley, A. D. Richmond, and R. E. Dickinson, “A coupled thermosphere/ionospheregeneral circulation model,” Geophys. Res. lett., Vol. 15, No. 12, pp. 325–1328, 1988.

05 Roble., R. G. and E. C. Ridley, “A thermosphere-ionosphere-mesosphere-electrodynamics general circulation model (TIME-GCM): Equinox solar cycle minimum simulations (30–500 km),” Geophys. Res.
lett., Vol. 21, No. 6, pp. 417–420, 1994.

06 Raeder, J., Y. L. Wang, and T. J. Fuller-Rowell, “Geomagnetic storm simulation with a coupled magnetosphere-ionosphere-thermosphere model,” in Space Weather: Progress and Challenges in Researchand Applications, pp. 377–384, P. Song, H.J. Singer, and G. Siscoe, editors, Geophysical Monograph,No. 125, AGU, Washington, D.C., 2001.

07 Tóth, G., D. L. De Zeeuw, T. I. Gombosi, W. B. Manchester, A. J. Ridley, Igor V. Sokolov, and I. I. Roussev, “Sun-to-thermosphere simulation of the 28–30 October 2003 storm with the Space Weather Modeling Framework,” Space Weather, Vol. 5, S06003, doi: 10.1029/2006SW000272, 2007.

206 Journal of the National Institute of Information and Communications Technology Vol.56 Nos.1-4 2009

08 Jin, H., Y. Miyoshi, H. Fujiwara, and H. Shinagawa, “Electrodynamics of the formation of ionosphericwave number 4 longitudinal structure,” J. Geophys. Res., Vol. 113, No. A09307, doi:
10.1029/2008JA013301, 2008.

09 Richards, P. G., J. A. Fennelly, and D. G. Torr, “EUVAC: A Solar EUV Flux Model for Aeronomic Calculations,” J. Geophys. Res., Vol. 99, No. A5, pp. 8981–8992, 1994.

10 Hedin, A. E., “Extension of the MSIS thermosphere model into the middle and lower atmosphere,”
J. Geophys. Res., Vol. 96, pp. 1159–1172, 1991.

11 Fuller-Rowell, T. J. and D. Rees, “Derivation of a conservation equation for mean molecular weight for a
two constituent gas within a three-dimensional, time-dependent, model of the thermosphere,” Planet.
Space Sci., 31, 1209–1222, 1983.

12 Hedin, A. E., E. L. Fleming, A. H. Manson, F. J. Schmidlin, S. K. Avery, R. R. Clark, S. J. Franke, G. J. Fraser,
T. Tsuda, F. Vial, and R. A. Vincent, “Empirical wind model for the upper, middle and lower at mosphere,” J. Atmos. Terr. Phys., Vol. 58, pp. 1421–1447, 1996.

13 Miyoshi, Y. and H. Fujiwara, “Day-to-day variations of migrating diurnal tide simulated by a GCM fromthe ground surface to the exobase,” Geophys. Res. Lett., Vol. 30, No. 5, pp. 1789–1792, doi:
10.1029/2003GL017695, 2003.

14 Liu, H., H. Lühr, V. Henize, and W. Köhler, “Global distribution of the thermospheric total mass density
derived from CHAMP,” J. Geophys. Res., Vol. 110, No. A04301, doi: 10.1029/2004JA010741, 2005.

15 Liu, H., S. Watanabe, and T. Kondo, “Fast thermospheric wind jet at the Earth’s dip equator,” Geophys.
Res. lett., Vol. 36, No. L08103, doi: 10.1029/2009GL037377, 2009.

SHINAGAWA Hiroyuki, Ph.D.
Senior Researcher, Space Environment
Group, Applied Electromagnetic
Research Center
Ionospheric Physics

SHINAGAWA Hiroyuki

207