- Email: [email protected]
A numerical solution of the time-varying diffusion equation for the F2-layer
Planet.
Space
Sci.
1967, Vol.
15. pp. 701 to 713.
A NUMERICAL DIFFUSION
Pergamon
Press
Ltd.
Printed
in Northern
Ineland
SOLUTION OF THE TIME-VARYING EQUATION FOR THE F2-LAYER
ROGER G. BAXTER Department of Applied Mathematics, University of Sheffield (Received 14 November 1966) Abstract-The continuity equation for the electron density in the R-region of the ionosphere is considered. This is a linear second order non-homogeneous partial differential equation. A general method of solution of the equation when diffusion, production, loss and timedependence are included is given. This is achieved by making two substitutions and integrating along a field line, using the Crank Nicolson method to advance the integration timewards. A world wide distribution of the maximum electron density is obtained and the results are discussed in connection with results obtained by other authors. 1. INTRODUCTION In 1960, Rishbeth and Barron”) solved the relevant equation for conditions of equilibrium such as may be approached in the daytime F-region in moderate latitudes. The diffusion equation was treated numerically in a number of cases using Chapman-type production functions, and exponential formulae for the diffusion and loss coefficients including electrodynamic drift. They deduced that the maximum electron density and the height at which it occurs are such that the magnitude of the production, loss and diffusion processes are approximately equal. The modifications produced by vertical electrodynamic drifts were also discussed. Martyn,t2) Dungeyc3) and Duncant4) discussed the behaviour of the FZ-layer at night, when it is under the influence of loss processes, diffusion and vertical movement. Yonezawa(5) considered the daytime model which includes the production of ionization by solar radiation and solved the relevant equations for a particular model atmosphere which led to solutions in terms of analytic functions. Also, Gliddon@) obtained an analytic time dependent solution in the case of moderate latitude, for straight infkritely long field lines. Ferraro and t)zdogan(‘) discussed the effect on the vertical distribution of ionization in an isothermal Chapman region. They assumed that electrons disappeared according to an attachment type law of recombination, the attachment coefficient K being a constant and confined their investigation to the equator at the equinox. Using the straight field line model, Gliddon and Kendall@) calculated world curves for two cases in order to give a baseline from which the true extent of the equatorial anomaly could be determined. They detected the apparent influence of diffusion at moderate latitudes, but the model could be of little help near the magnetic equator. Later, Gliddon and Kendall@) generalized the solution to include uniform upward or downward drift velocity. The electron density was expressed in the form of a definite integral, as a function of three parameters representing attachment-type recombination, diffusion and electrodynamic drift velocity. However, once again the solution was only valid at moderate latitudes. Rishbeth,o”) using an analogue computer, solved the continuity equation for electron density in the F-region and demonstrated the respective roles of production, loss and diffusion in determining the time-varying electron density. In 1964, Rishbeth”” considered a model quiet F2-layer at mid-latitudes which combined time dependent solutions of the F2-layer continuity equation in a static atmospheredo) with 701
702
ROGER
G. BAXTER
the time-varying atmospheric models of Harris and Priester.02) Rishbeth, Lyon and Pearto3) calculated the equilibrium distribution of ionization in the equatorial Fregion on the assumption that the predominant processes were diffusion along the Earth’s magnetic field together with the usual production by photoionization and loss by recombination. Moffett and Hanson, Bramley and Peart(1s*16)and Hanson and Moffettu7) obtain solutions for the steady-state continuity equation jn the equatorial F-region, incorporating an electrodynamic drift in an upward direction perpendicular to the magnetic field lines. They suggest that the results show that it would be possible to develop a quantitive explanation of the Appleton anomaly. Thomas(1e) gave results of calculations of the electron density distribution in the Fregion when thermal non-equilibrium exists between the neutral gas, electrons and ions. The equation included photoionization, loss and diffusion for a model in which the neutral gas and electron temperatures were independent of height and the ion temperature increased from that of the neutral gas at low heights and approached that of the electrons at great heights. The investigation indicates how the height distribution of ionization at any time is determined by the neutral gas and electron temperatures and by the variation of the ion temperature with height. An atmospheric model similar to that taken by these previous authors which includes production, loss, diffusion and time dependence (but not electrodynamic drift) is considered here. The curvature of the field lines is taken into account and a numerical solution for the maximum electron density is found, giving world wide distributions valid near the equator as well as at moderate latitudes. 2. EQUATIONS The
OF THE PROBLEM
continuity equation for the distribution of electron density in the F2-layer is given by aN
-=q-L+
(1)
D.92N-91N
at
where q is the production rate, L the loss, gS the diffusion operator, D the coefficient of diffusion and .& the electrodynamic drift operator. Clearly, full time dependent calculations are desirable, including both electrodynamic drift and the curvature of the field lines. However, in the present paper, the effects of drift are omitted, i.e. we omit the QIN term. Otherwise the model is physically consistent. The d@sion operator used is that given by Kendall(20) and Lyonf21) for diffusion along the field lines of a dipole approximation to the geomagnetic field and can be written as Nsin2 Z .92N=2Ha+2HsinZ
sinlgf-+(9+
+ &
aN COS ZaN sinI% +-r ae
3
cos z a sin I aN Y&+-r ae )(
15.u2)sinZ 2rA2
-l---
COS 1
aN
r
ae
, cosZaN r ae
b2(9 + 15~~) - sina e]
12)
SOLUTION
where I = H = 0= /J = A= =
OF THE TIME-VARYING
DIFFUSION
EQUATION
703
magnetic dip angle scale height of the neutral ionizable gas which is assumed to be atomic oxygen colatitude cos 0 1 + 3~~.
In deriving equation (2) it has been assumed that the atmosphere is isothermal,, with equal electron, ion and neutral particle temperatures. The production term q is taken to be the modified expression given by Chapman(22*23) and is given by
4=
1
(3)
qO 0 exp ( 1 - ‘$
- secXexp(-‘+))
5O
where z is the height above some fixed reference, q,, is the maximum production rate for overhead sun and occurs at the height z, (z,, is taken to be 180 km) and x denotes the angular zenith distance of the Sun from the point of observation. The values 0 < 1x1< 5 correspond to daylight 7l
5
< 1x1< n correspond to night.
Also, we only consider equinoctial conditions, in which case cos x = sin 0 sin +
(4)
where 4 is the time measured in radians. The loss term is of an attachment typef”) and is given by L=KN
(5)
where K is the coefficient of attachment and given by K=$exp(-17)
(6)
where /I and 2 are constants and T is the number of seconds in a radian. 3. SUBSlTI’UTIONS
We now make the following substitutions throughout Y=exp(-_)
(7)
N=yu
(8)
and where r = a sin2 8
(9)
704
ROGER G. BAXTER
is the equation of a typical field line and r,, is the base of the Fregion (r,-,= 6370 + 180 km: the mean radius of the Earth is taken as 6370 km). We also use t = T+
(IO)
where 4 is time measured in radians and T is the number of seconds in a radian. Making the above substitutions and using equations (2), (3), (5) and (6), equation (1) becomes (I - 10~” - 15~4) av aA4H sin2 0 3 + 7qOeyexp (-7” cosec
e cmec +) - /?yaAv
(11)
Finally, we make the assumption that the coefficient of diffusion D is inversely proportional to n, the density of the neutral atoms, i.e. DCwhere b is a constant.
b n
Also, for an isothermal atmosphere n = n, exp (-z/H)
(13)
where n, is a constant, and thus we have
We define a dimensionless parameter y by &!!
(15)
Y=HP rDy2
(16)
and using equation (12)
Hence, equation (11) finally becomes
av
-= a+
sin2 z asv + -4~ aY2
(1 - 10~2- 15jP)H 2 ya sin2 0A4y ay
+ 7qOeyexp (-y2 cosec e cosec 4) - &%
(17)
which is to be solved for v. Note that this is an equation with v in terms of y and + for each field line. 4. METHOD OF SOLUTION
For solving equation (17) the method used is the Crank and Nicolson’26) which is known to be a stable method [see also Moden Computing Methods p. 116’28)]. This states that
(18)
SOLUTION OF THE TIME-VARYING
DIFFUSION
EQUATION
705
where C, is the difference correction for 4. (v)+ means v evaluated at time 4. It is standard practice on a fast computer to take S+, the step length in time, small enough so that C, may be neglected. A check on the truncation error may be easily carried out by halving the step length 84, keeping the other parameters fixed, repeating the computations, and examining the difference between the two sets of results. The equation is integrated over the semi-infinite rectangle 0 =ZZ y Q 1, 0 Q +. The range of y is split into. (n - 1) equal intervals each interval, being h, apart, and the time interval r$ is split into equal intervals of S+. Thus, (u,).+ means the value of u at y = y,, + rh, at time 4 where y,, is one of the boundaries of y to be described later. 5. BOUNDARY
CONDlTIONS
Three boundary conditions are required (I) at y = y0 for all 4, (II) at y = yn_l for all 4 and (III) at rj = #,,, i.e. a starting condition is required. (I) Boundary condition at the equator About the equator 0 = 7r/2 or a = 0, the electron density N is assumed to be symmetrical, where o! is the latitude (a = 742 - 0) and y=y,=exp(-9) This gives aN
~=oa~=o
av
(20)
Now consider the variable y, given by equation (7) which can be written as y = exp
-
a co9 a - r, 2H
(21)
This is symmetrical about a = 0, and so it would appear that a boundary condition has been lost. However, using central differences, from equation (20) (03” = (t)-,)
(22)
where the subscript 0 represents the value at the equator and where the points corresponding to (v$ and (u_$ are each separated by a distance h, from the equatorial value. Now at the equator
(~)~=o=~ou(~),,+Y(~)
(23)
which, from equation (21) and the fact that
asvo
( ) -= aa2
2(v,” - vo=) h,2
(24)
becomes aaN
(4 aa 8
a=o
(25)
706
ROGER G. BAXTER
Also, since N=yu we have at the equator
If N = N,” = yOv,,is the equatorial value, then if N,” coincides with yrui, equations (25) and (26) may be equated to give 2W1 - %) ay&Z
a% -= aY
(27)
To make equation (27) valid, we have at a = 0 vi at h, = u1 at h, Now sinh,=
a - r, + 2Hlog yl a
cm
where (29)
YI= YOf h, But since h, is small, the condition that equation (27) is valid is hma=
a - r, + 2Hlog yl a
(30)
This is the boundary condition taken at the equator until very high field lines are considered (i.e. a is very large) in which case y. --c 0. In this case the boundary condition taken is 0, = 0,
(31)
a sina 8 = r,
(32)
Y= 1
(33)
(II) Bounabry condition at low height At low height which implies The atmosphere there is considered to be in photoequilibrium,
i.e.
q=L
(34)
which gives, from equations (3), (4), (5) and (6) N=
1
j{~qeexpC
- cosec 8 cosec 4)))
(35)
(III) Starting condition We start the computations at noon where we assume equilibrium conditions, i.e. aN z=O’x-
au
av
-“=-%=o
(36)
SOLUTION OF THE TIME-VARYING
DIFFUSION
EQUATION
707
With this condition equation (17) can be solved immediately to give initial values. 6. PROCEDURE
For each set of integrations, a particular field line is taken by choosing (I = constant and integrating from the equator 0, = J+,) along the field line until it sinks below a certain height 0, = 1) which is below the F24ayer. Thus, the curved field line has been mapped into a straight line from y = y. toy = 1 by the mapping in equation (7). The integration, if carried out for the period of one complete day, gives a world wide distribution of N. Since the required solution of N is known to be periodic, with a period of 24 hr, the computations are carried on until this so called quasi-steady state is reached. The results along the field line are printed out on the hour throughout the period of one day. At each time, N attains a maximum value Nmsx at a certain latitude. Thus, repeating the calculations for a wide distribution of field lines (i.e. various values of a), a complete distribution of N max with latitude may be obtained. From these distributions of Nmax, world contours of Nmax can be drawn. VALUES AND DEXUSSION Two sets of results (a) and (b) will be considered. (a) The@ set of results was obtained by using typical F2-layer parameters. The scale height H = 50 km, fi = 10, y = 10, 1= 1, the number of equations = 90 and the steplength in time = 6 min. As stated in the previous section, the computations were carried out until a quasisteady state was reached. For this to happen, using the equilibrium solution at noon as the starting condition, it was found that the computations had to be continued for three or four times round the Earth (i.e. 3 or 4 days). On investigation it was seen that the equilibrium approximation was about 25 per cent higher in equatorial regions than the periodic solution finally obtained. Figure 1 shows N maxplotted against latitude for noon day values. The dotted curve I is the initial equilibrium solution and the dotted curve IV is the periodic solution obtained 3 or 4 days later. However, in order to speed up the computations, a linear multiple of the equilibrium solution was taken as the actual starting value, which 7. NUhlEIUCAL
FIG. 1. Nmax
AGAINST
LATITUDE
AT
Comparison of starting values.
NOON.
708
ROGER G. BAXTER
gave a periodic solution after one revolution round the Earth, thus halving computer time. Curve II shows the linear multiple of the equilibrium solution actually used for starting the computations and curve III shows the solution after one revolution. It must be noted that owing to the nature of the equation, it does not matter what starting values are taken, the same result will be obtained after computing long enough. However, these results would appear to indicate that equilibrium solutions tend to be too high. World contours of N,,, Figure 2 shows world contours of NmsX. The numbers on the contours which are isolines of ALax are N x 102/7q,e where N is the electron density.
Fm. 2. N,,,,, WORLD ~~R~Es-I~~LINES OF Nmax. The numbers on the contours are N x 101/rq0 e where N is the actual electron density, /?=10,y=10,1=1,N=50km.
This Figure is compared with the world curves derived for the equinox by Gliddon and Kendall@) where the maximum electron density was calculated for a few latitudes and where the results were known to be valid only for higher latitudes. The two sets of curves have a general resemblance with a drawing out of the electron density in evening and nighttime and a burst of concentration at sunrise. It will be noticed that around the equator at certain times there are small equatorial troughs. Since they are very shallow (i.e. the ratio of the maximum to the equatorial value is small ~1.07:1), and also occur within 10”N and S latitude, they are not particularly significant. The behaviour of the troughs may be seen by examining Fig. 2. Diurnal variation We next examine the diurnal variation of the equatorial ALax. Figure 3 shows ALaX at the equator throughout a day. The special features to note here are the minimum at about
SOLUTION OF THE TIME-VARYING
DIFFUSION
EQUATION
709
6 hours, and maximum at 16 hours, which is in agreement with previous results. The diurnal variation is considered further with the second set of results. (b) A second set of results was also obtained, for more realistic FZ-layer parameters. A scale height H of 80 km was chosen which gives an electron and ion temperature T of about 1300°K which is for a medium sun-spot number. The value of ;I was taken to be l-75 which
007-
N_ TV
OQ5-
004
Oa-
0
I
I
I
I
I
Ia
6
12
16
20
2&
4
FIG.3. DIURNAL VARLOION OF Nmar AT THE EQUATOR. p = 10, y = 10, I = 1, H = 50 km.
would be so if nitrogen N, were the molecular gas participating subsequent dissociative recombination of the molecular ions, i.e.
in the formation
O++N,+NO++N At a temperature is of the order
and (37)
of 1300°K the coefficient of diffusion at noon at an altitude of 300 km DW w 4.3 x lOlo cme set-l
and at the same altitude the coefficient of attachment K [equation (5)] taken to be KW = 4 x 1O-4se+ (see Rishbethuu p. 683, Fig. 14 where results from various models are plotted together with collected data from various sources). This gives the values of /I and y to be approximately 76 and O-5respectively. These vary greatly with sunspot number and in fact the values /? = 35 and y = 1 were used in the calculations. The number of equations and the step length in time were taken as in (a) (90, 6 mm respectively). World contours of Nmax
Figure 4 shows world contours for Nmax. The numbers on the contours, which are isolines of I&,, are N x 102/7qoe where N is the actual electron density. Once again it will be noticed that there are small equatorial troughs but since the ratio of the equatorial value to the peak value is at most 1: l-06 occurring at 09 hours, and also the peaks are within IO“ latitude of the equator they are not particularly significant. However, it should be noted the difference in behaviour of the troughs between these results and the results (a) Fig. 2. The burst of electron density at sunrise and the tailing off after
710
ROGER
G. BAXTER
SunSet iS Once again evident and the peak value of the electron density occurs just before sunset. The height of the maximum Figure 5 shows the variation with latitude a of the height of the maximum h, above z = zo for noon and midnight. On the noon curve, h, initially rises slightly with increasing
60-
OP Nm.=. FIG. 4. Nmax WORLD ~~R~JB-I~~LINW N x lo*/7 q. e where N is the actual electron density. 1.75,H = 80 km. Taking q. = 1000 cm-* WC-‘,(sq,,e)lO-’H 3.7 x lo6 cm-* e.g. for contour number 10, N N 3.7 x 10‘ cm-*.
The numbers on the contours are
/? = 35, y = 1,1=
latitude away from the equator and then falls to an approximately constant level until about 70” and then rises again. The rise again at 70” may be due to the production function used here which is known to become less valid as the poles are approached. At midnight, there is no “equatorial trough” in h,, which in fact falls from a maximum at the equator to a constant level at about 45” latitude and remains constant out to the poles. Also, the height of the maximum at midnight is higher than at noon. A possible explanation of this is that at noon, production, loss and diffusion are active but at night there is no production with the result the loss term becomes more effective and “eats” into the electron density at lower heights causing the peak to rise during the night. Diurnal variation Figure 6 shows the diurnal variation of Nmax at four different latitudes (0”, 30”, 50”, 80”). It is seen that on the equator, the maximum occurs at about 17 hours and as latitude increases, it occurs earlier until at 80” the maximum is at about 134 hours. A possible explanation of this is that at noon on the equator the loss coefficient at the peak is smaller than for higher latitudes. The difference between the equatorial and 30” curves is much greater than the difference between the 30” and 50” curves. In equatorial regions h, is
SOLUTION
100
OF THE TTME-VARYING
I 0
10
DIFFUSION
EQUATION
I
I
I
I
I
I
I
P
30
40
50
60
m
M)
G_
711
lalllU&
Ro. 5. LATITUDE VARIATIONOF THE HEIGHTh,
Ttn
8”
OF THE FUCIRON PEAK.
hO”lr
FIG. 6. DIIJRNAI.VAIUTION OF Nmax. /3= 35, y = 1. A = 1.75, H = 80 km.
higher than for larger latitudes (see Fig. 5), thus the loss coefficient is smaller and so the time of decay of ionization takes longer. Between 30” and 50”, h, is approximately constant and so the loss coefficient is nearly constant, thus the time of decay will be approximately the same. Also it will be noticed that the peak electron density decreases steadily from the equator towards the poles. These results are in general agreement with the results obtained by Gliddon and Kendall.@)
712
ROGER 8. GENERAL
REMARKS
G. BAXTER ON THE TRAWFORMATIONS
Although the electrodynamic drift has been omitted from the above calculations, they do serve as a more consistent base from which the effect of time dependent electrodynamic drifts will cause departure. They are also of interest because of the fruitful transformation y = exp (-2/2H) which is used along a field line, and the subsequent transformation N = yv Owing to the nature of the equation being solved, the transformations as a stabilizing effect, for at great height
(38) and (39) act
v = constant. The transformation to y simplifies the algebra considerably and also y is a solution of the diffusion operator (2). Also, a solution in y gives a much better distribution of points which are relevant to the calculation, especially for very high field lines (i.e. for very large values of a). Further, the transformation maps the field lines into straight lines giving straight boundaries to the problem. (See Section 4.) Finally, using the transformations (38) and (39), we have been able to solve the timedependent diffusion equation, giving a world wide distribution of the electron density from the equator out to the poles. Acknowledgenzents-I should like to thank Professor P. C. Kendall who su ted this investigation. The results were carried out on the Ferranti “Mercury” computer here at SheffieY d. This research was partly sponsored by the Air Force Cambridge Research Laboratories under Grant AF EOAR 64-74 through the European Office of Aerospace Research (OAR), United States Air Force. REFERENCES 1. H. Rrstrn~nr and D. W. BARRON,J. afmm. few. Phys. 18,234 (1960). 2. D. F. MARTYN,Ausr. J. Phys. 9, 161 (1956). 3. J. W. DUNGEY.J. atmos. terr. Phvs. 9.90 (1956). 4. R. A. DUNCAN, Aust. J. Phys. 9,236 (1958). ’ 5. T. J. YONEZAWA,Radio Res. Labs. 5, 165 (1958). 6. J. E. C. GLIDDON,Q. J. mech. appl. Math. 12, 340 (1959). 7. V. C. A. FERRAROand I. ~ZD~GAN, J. atmos. terr. Phys. 12,140 (1958). 8. J. E. C. GL~DD~Nand P. C. KENDALL,J. armos. ten. Phys. 18,48 (1960). 9. J. E. C. GLIDDONand P. C. KENDALL,J. armos. terr. Phys. 24, 1073 (1962). 10. H. RISHBETH, Proc. phys. Sot. 81,65 (1963). 11. H. RISHBETH,J. atrues. ferr. Phys. 26,657 (1964). 12. I. HARRISand W. PRDSTER,NASA Technical Note D-1444 (1962). 13. H. Rrsrn~s.m, A. J. LYON and M. PEART,J. geophys. Res. 68,2559 (1963). 14. R. J. MOFFETT and W. B. HANSON,Nutwe, L0m.f. 206,705 (1965). 15. E. N. BRAMLEYand M. PEART,Nature, Lord. 206,1245 (1965). 16. E. N. BRAMLEYand M. PEART,J. atmos. ferr. Phys. 27, 1201 (1965). 17. W. B. HANSONand R. J. Mom, Private communication (1966). 18. E. V. APPLETON,Nature, Land. 157.691 (1946). 19. L. THOMAS,J. geophys. Res. 71, 1357 (1966). 20. P. C. KENDALL.J. atmos. terr. Phvs. 24.805 (1962). 21. A. J. LYON, J.geophys. Res. 68, i531 (i963). 22. S. CHAPMAN,Proc. Phys. Sot. 43,26 (1931). 23. S. CHAPMAN,Proc. Phys. Sot. 43,483 (1931). G, C. S. G. K. SEITY, and J. 0. THOMAS,Phil. 7kuns. R. Sot. A248, 24. J. A. RATCLIPFE, E. R. SC-
621 (1956). 25. J. CRANKand P. NKXXSON, Proc. Camb. Phil. Sot. 43,50 (1947). 26. National Physical Laboratory, Modem Computing Methods H.M.S.O. (1961).
SOLUTION
OF THE TIME-VARYING
Pf33IOMhPaCCMaTpEiBaeTClI B
CJIOe
F,i?
,qm_$@epeKqtfanbKoe CJlyYae, aT0
KOrAa
4aIOT U
IIyTeM
cnoco6
JJl@I#y3HR, AByX
pe3yJIbTaTbI .
061qnti
MeTOA
UpOABHFKeHHH
paClIpeAeneHKe CBHBIl
C
3JIeKTpOHHOi
¶aCTEl¶HO
ypaBHeHHR
Ii 3aBHCKMOCTb
BAOJIb
pe3yJIbTaTaMU,
IIJIOTHOCTH
HeOAHOpOAHOe,
II0
IIOJl~eAHbIMEi
TOM
IlpKMeHfWI
BpeMeHK.
3JleKTpOHHOt
B
OT BpeMeHK.
JIEIHHH IIOi’lH,
IiHTerpaqHH
lIpeAeJIbHOlf
713
EQUATION
peIJ.IeHHH
IlOTHpR
H KHTerpHpOBaHHeM
AJIU B
AJIfi
IIOpJIJJKa,
IIpOH3BOACTB0,
3aMeH
HUKOJIbCOHa
06CyxAaIOTCH
BTOpOl’O
AaeTCFI
IllEpOKO-paCIlpOCTpaHeHHOe
aBTOpaMEl
Hepa3pblBHOCTEi
aTO-JlHHettHOe,
ypaBKeme. BKJ’fIOUeHbI
AOCTHPaeTCR
KOJleH=iaTbl&i
ypaBHeHHe
EIOHOCtjjepbI.
DIFFUSION
nOJIy-
IIJIOTHOCTH, ApyrHMH
Space
Sci.
1967, Vol.
15. pp. 701 to 713.
A NUMERICAL DIFFUSION
Pergamon
Press
Ltd.
Printed
in Northern
Ineland
SOLUTION OF THE TIME-VARYING EQUATION FOR THE F2-LAYER
ROGER G. BAXTER Department of Applied Mathematics, University of Sheffield (Received 14 November 1966) Abstract-The continuity equation for the electron density in the R-region of the ionosphere is considered. This is a linear second order non-homogeneous partial differential equation. A general method of solution of the equation when diffusion, production, loss and timedependence are included is given. This is achieved by making two substitutions and integrating along a field line, using the Crank Nicolson method to advance the integration timewards. A world wide distribution of the maximum electron density is obtained and the results are discussed in connection with results obtained by other authors. 1. INTRODUCTION In 1960, Rishbeth and Barron”) solved the relevant equation for conditions of equilibrium such as may be approached in the daytime F-region in moderate latitudes. The diffusion equation was treated numerically in a number of cases using Chapman-type production functions, and exponential formulae for the diffusion and loss coefficients including electrodynamic drift. They deduced that the maximum electron density and the height at which it occurs are such that the magnitude of the production, loss and diffusion processes are approximately equal. The modifications produced by vertical electrodynamic drifts were also discussed. Martyn,t2) Dungeyc3) and Duncant4) discussed the behaviour of the FZ-layer at night, when it is under the influence of loss processes, diffusion and vertical movement. Yonezawa(5) considered the daytime model which includes the production of ionization by solar radiation and solved the relevant equations for a particular model atmosphere which led to solutions in terms of analytic functions. Also, Gliddon@) obtained an analytic time dependent solution in the case of moderate latitude, for straight infkritely long field lines. Ferraro and t)zdogan(‘) discussed the effect on the vertical distribution of ionization in an isothermal Chapman region. They assumed that electrons disappeared according to an attachment type law of recombination, the attachment coefficient K being a constant and confined their investigation to the equator at the equinox. Using the straight field line model, Gliddon and Kendall@) calculated world curves for two cases in order to give a baseline from which the true extent of the equatorial anomaly could be determined. They detected the apparent influence of diffusion at moderate latitudes, but the model could be of little help near the magnetic equator. Later, Gliddon and Kendall@) generalized the solution to include uniform upward or downward drift velocity. The electron density was expressed in the form of a definite integral, as a function of three parameters representing attachment-type recombination, diffusion and electrodynamic drift velocity. However, once again the solution was only valid at moderate latitudes. Rishbeth,o”) using an analogue computer, solved the continuity equation for electron density in the F-region and demonstrated the respective roles of production, loss and diffusion in determining the time-varying electron density. In 1964, Rishbeth”” considered a model quiet F2-layer at mid-latitudes which combined time dependent solutions of the F2-layer continuity equation in a static atmospheredo) with 701
702
ROGER
G. BAXTER
the time-varying atmospheric models of Harris and Priester.02) Rishbeth, Lyon and Pearto3) calculated the equilibrium distribution of ionization in the equatorial Fregion on the assumption that the predominant processes were diffusion along the Earth’s magnetic field together with the usual production by photoionization and loss by recombination. Moffett and Hanson, Bramley and Peart(1s*16)and Hanson and Moffettu7) obtain solutions for the steady-state continuity equation jn the equatorial F-region, incorporating an electrodynamic drift in an upward direction perpendicular to the magnetic field lines. They suggest that the results show that it would be possible to develop a quantitive explanation of the Appleton anomaly. Thomas(1e) gave results of calculations of the electron density distribution in the Fregion when thermal non-equilibrium exists between the neutral gas, electrons and ions. The equation included photoionization, loss and diffusion for a model in which the neutral gas and electron temperatures were independent of height and the ion temperature increased from that of the neutral gas at low heights and approached that of the electrons at great heights. The investigation indicates how the height distribution of ionization at any time is determined by the neutral gas and electron temperatures and by the variation of the ion temperature with height. An atmospheric model similar to that taken by these previous authors which includes production, loss, diffusion and time dependence (but not electrodynamic drift) is considered here. The curvature of the field lines is taken into account and a numerical solution for the maximum electron density is found, giving world wide distributions valid near the equator as well as at moderate latitudes. 2. EQUATIONS The
OF THE PROBLEM
continuity equation for the distribution of electron density in the F2-layer is given by aN
-=q-L+
(1)
D.92N-91N
at
where q is the production rate, L the loss, gS the diffusion operator, D the coefficient of diffusion and .& the electrodynamic drift operator. Clearly, full time dependent calculations are desirable, including both electrodynamic drift and the curvature of the field lines. However, in the present paper, the effects of drift are omitted, i.e. we omit the QIN term. Otherwise the model is physically consistent. The d@sion operator used is that given by Kendall(20) and Lyonf21) for diffusion along the field lines of a dipole approximation to the geomagnetic field and can be written as Nsin2 Z .92N=2Ha+2HsinZ
sinlgf-+(9+
+ &
aN COS ZaN sinI% +-r ae
3
cos z a sin I aN Y&+-r ae )(
15.u2)sinZ 2rA2
-l---
COS 1
aN
r
ae
, cosZaN r ae
b2(9 + 15~~) - sina e]
12)
SOLUTION
where I = H = 0= /J = A= =
OF THE TIME-VARYING
DIFFUSION
EQUATION
703
magnetic dip angle scale height of the neutral ionizable gas which is assumed to be atomic oxygen colatitude cos 0 1 + 3~~.
In deriving equation (2) it has been assumed that the atmosphere is isothermal,, with equal electron, ion and neutral particle temperatures. The production term q is taken to be the modified expression given by Chapman(22*23) and is given by
4=
1
(3)
qO 0 exp ( 1 - ‘$
- secXexp(-‘+))
5O
where z is the height above some fixed reference, q,, is the maximum production rate for overhead sun and occurs at the height z, (z,, is taken to be 180 km) and x denotes the angular zenith distance of the Sun from the point of observation. The values 0 < 1x1< 5 correspond to daylight 7l
5
< 1x1< n correspond to night.
Also, we only consider equinoctial conditions, in which case cos x = sin 0 sin +
(4)
where 4 is the time measured in radians. The loss term is of an attachment typef”) and is given by L=KN
(5)
where K is the coefficient of attachment and given by K=$exp(-17)
(6)
where /I and 2 are constants and T is the number of seconds in a radian. 3. SUBSlTI’UTIONS
We now make the following substitutions throughout Y=exp(-_)
(7)
N=yu
(8)
and where r = a sin2 8
(9)
704
ROGER G. BAXTER
is the equation of a typical field line and r,, is the base of the Fregion (r,-,= 6370 + 180 km: the mean radius of the Earth is taken as 6370 km). We also use t = T+
(IO)
where 4 is time measured in radians and T is the number of seconds in a radian. Making the above substitutions and using equations (2), (3), (5) and (6), equation (1) becomes (I - 10~” - 15~4) av aA4H sin2 0 3 + 7qOeyexp (-7” cosec
e cmec +) - /?yaAv
(11)
Finally, we make the assumption that the coefficient of diffusion D is inversely proportional to n, the density of the neutral atoms, i.e. DCwhere b is a constant.
b n
Also, for an isothermal atmosphere n = n, exp (-z/H)
(13)
where n, is a constant, and thus we have
We define a dimensionless parameter y by &!!
(15)
Y=HP rDy2
(16)
and using equation (12)
Hence, equation (11) finally becomes
av
-= a+
sin2 z asv + -4~ aY2
(1 - 10~2- 15jP)H 2 ya sin2 0A4y ay
+ 7qOeyexp (-y2 cosec e cosec 4) - &%
(17)
which is to be solved for v. Note that this is an equation with v in terms of y and + for each field line. 4. METHOD OF SOLUTION
For solving equation (17) the method used is the Crank and Nicolson’26) which is known to be a stable method [see also Moden Computing Methods p. 116’28)]. This states that
(18)
SOLUTION OF THE TIME-VARYING
DIFFUSION
EQUATION
705
where C, is the difference correction for 4. (v)+ means v evaluated at time 4. It is standard practice on a fast computer to take S+, the step length in time, small enough so that C, may be neglected. A check on the truncation error may be easily carried out by halving the step length 84, keeping the other parameters fixed, repeating the computations, and examining the difference between the two sets of results. The equation is integrated over the semi-infinite rectangle 0 =ZZ y Q 1, 0 Q +. The range of y is split into. (n - 1) equal intervals each interval, being h, apart, and the time interval r$ is split into equal intervals of S+. Thus, (u,).+ means the value of u at y = y,, + rh, at time 4 where y,, is one of the boundaries of y to be described later. 5. BOUNDARY
CONDlTIONS
Three boundary conditions are required (I) at y = y0 for all 4, (II) at y = yn_l for all 4 and (III) at rj = #,,, i.e. a starting condition is required. (I) Boundary condition at the equator About the equator 0 = 7r/2 or a = 0, the electron density N is assumed to be symmetrical, where o! is the latitude (a = 742 - 0) and y=y,=exp(-9) This gives aN
~=oa~=o
av
(20)
Now consider the variable y, given by equation (7) which can be written as y = exp
-
a co9 a - r, 2H
(21)
This is symmetrical about a = 0, and so it would appear that a boundary condition has been lost. However, using central differences, from equation (20) (03” = (t)-,)
(22)
where the subscript 0 represents the value at the equator and where the points corresponding to (v$ and (u_$ are each separated by a distance h, from the equatorial value. Now at the equator
(~)~=o=~ou(~),,+Y(~)
(23)
which, from equation (21) and the fact that
asvo
( ) -= aa2
2(v,” - vo=) h,2
(24)
becomes aaN
(4 aa 8
a=o
(25)
706
ROGER G. BAXTER
Also, since N=yu we have at the equator
If N = N,” = yOv,,is the equatorial value, then if N,” coincides with yrui, equations (25) and (26) may be equated to give 2W1 - %) ay&Z
a% -= aY
(27)
To make equation (27) valid, we have at a = 0 vi at h, = u1 at h, Now sinh,=
a - r, + 2Hlog yl a
cm
where (29)
YI= YOf h, But since h, is small, the condition that equation (27) is valid is hma=
a - r, + 2Hlog yl a
(30)
This is the boundary condition taken at the equator until very high field lines are considered (i.e. a is very large) in which case y. --c 0. In this case the boundary condition taken is 0, = 0,
(31)
a sina 8 = r,
(32)
Y= 1
(33)
(II) Bounabry condition at low height At low height which implies The atmosphere there is considered to be in photoequilibrium,
i.e.
q=L
(34)
which gives, from equations (3), (4), (5) and (6) N=
1
j{~qeexpC
- cosec 8 cosec 4)))
(35)
(III) Starting condition We start the computations at noon where we assume equilibrium conditions, i.e. aN z=O’x-
au
av
-“=-%=o
(36)
SOLUTION OF THE TIME-VARYING
DIFFUSION
EQUATION
707
With this condition equation (17) can be solved immediately to give initial values. 6. PROCEDURE
For each set of integrations, a particular field line is taken by choosing (I = constant and integrating from the equator 0, = J+,) along the field line until it sinks below a certain height 0, = 1) which is below the F24ayer. Thus, the curved field line has been mapped into a straight line from y = y. toy = 1 by the mapping in equation (7). The integration, if carried out for the period of one complete day, gives a world wide distribution of N. Since the required solution of N is known to be periodic, with a period of 24 hr, the computations are carried on until this so called quasi-steady state is reached. The results along the field line are printed out on the hour throughout the period of one day. At each time, N attains a maximum value Nmsx at a certain latitude. Thus, repeating the calculations for a wide distribution of field lines (i.e. various values of a), a complete distribution of N max with latitude may be obtained. From these distributions of Nmax, world contours of Nmax can be drawn. VALUES AND DEXUSSION Two sets of results (a) and (b) will be considered. (a) The@ set of results was obtained by using typical F2-layer parameters. The scale height H = 50 km, fi = 10, y = 10, 1= 1, the number of equations = 90 and the steplength in time = 6 min. As stated in the previous section, the computations were carried out until a quasisteady state was reached. For this to happen, using the equilibrium solution at noon as the starting condition, it was found that the computations had to be continued for three or four times round the Earth (i.e. 3 or 4 days). On investigation it was seen that the equilibrium approximation was about 25 per cent higher in equatorial regions than the periodic solution finally obtained. Figure 1 shows N maxplotted against latitude for noon day values. The dotted curve I is the initial equilibrium solution and the dotted curve IV is the periodic solution obtained 3 or 4 days later. However, in order to speed up the computations, a linear multiple of the equilibrium solution was taken as the actual starting value, which 7. NUhlEIUCAL
FIG. 1. Nmax
AGAINST
LATITUDE
AT
Comparison of starting values.
NOON.
708
ROGER G. BAXTER
gave a periodic solution after one revolution round the Earth, thus halving computer time. Curve II shows the linear multiple of the equilibrium solution actually used for starting the computations and curve III shows the solution after one revolution. It must be noted that owing to the nature of the equation, it does not matter what starting values are taken, the same result will be obtained after computing long enough. However, these results would appear to indicate that equilibrium solutions tend to be too high. World contours of N,,, Figure 2 shows world contours of NmsX. The numbers on the contours which are isolines of ALax are N x 102/7q,e where N is the electron density.
Fm. 2. N,,,,, WORLD ~~R~Es-I~~LINES OF Nmax. The numbers on the contours are N x 101/rq0 e where N is the actual electron density, /?=10,y=10,1=1,N=50km.
This Figure is compared with the world curves derived for the equinox by Gliddon and Kendall@) where the maximum electron density was calculated for a few latitudes and where the results were known to be valid only for higher latitudes. The two sets of curves have a general resemblance with a drawing out of the electron density in evening and nighttime and a burst of concentration at sunrise. It will be noticed that around the equator at certain times there are small equatorial troughs. Since they are very shallow (i.e. the ratio of the maximum to the equatorial value is small ~1.07:1), and also occur within 10”N and S latitude, they are not particularly significant. The behaviour of the troughs may be seen by examining Fig. 2. Diurnal variation We next examine the diurnal variation of the equatorial ALax. Figure 3 shows ALaX at the equator throughout a day. The special features to note here are the minimum at about
SOLUTION OF THE TIME-VARYING
DIFFUSION
EQUATION
709
6 hours, and maximum at 16 hours, which is in agreement with previous results. The diurnal variation is considered further with the second set of results. (b) A second set of results was also obtained, for more realistic FZ-layer parameters. A scale height H of 80 km was chosen which gives an electron and ion temperature T of about 1300°K which is for a medium sun-spot number. The value of ;I was taken to be l-75 which
007-
N_ TV
OQ5-
004
Oa-
0
I
I
I
I
I
Ia
6
12
16
20
2&
4
FIG.3. DIURNAL VARLOION OF Nmar AT THE EQUATOR. p = 10, y = 10, I = 1, H = 50 km.
would be so if nitrogen N, were the molecular gas participating subsequent dissociative recombination of the molecular ions, i.e.
in the formation
O++N,+NO++N At a temperature is of the order
and (37)
of 1300°K the coefficient of diffusion at noon at an altitude of 300 km DW w 4.3 x lOlo cme set-l
and at the same altitude the coefficient of attachment K [equation (5)] taken to be KW = 4 x 1O-4se+ (see Rishbethuu p. 683, Fig. 14 where results from various models are plotted together with collected data from various sources). This gives the values of /I and y to be approximately 76 and O-5respectively. These vary greatly with sunspot number and in fact the values /? = 35 and y = 1 were used in the calculations. The number of equations and the step length in time were taken as in (a) (90, 6 mm respectively). World contours of Nmax
Figure 4 shows world contours for Nmax. The numbers on the contours, which are isolines of I&,, are N x 102/7qoe where N is the actual electron density. Once again it will be noticed that there are small equatorial troughs but since the ratio of the equatorial value to the peak value is at most 1: l-06 occurring at 09 hours, and also the peaks are within IO“ latitude of the equator they are not particularly significant. However, it should be noted the difference in behaviour of the troughs between these results and the results (a) Fig. 2. The burst of electron density at sunrise and the tailing off after
710
ROGER
G. BAXTER
SunSet iS Once again evident and the peak value of the electron density occurs just before sunset. The height of the maximum Figure 5 shows the variation with latitude a of the height of the maximum h, above z = zo for noon and midnight. On the noon curve, h, initially rises slightly with increasing
60-
OP Nm.=. FIG. 4. Nmax WORLD ~~R~JB-I~~LINW N x lo*/7 q. e where N is the actual electron density. 1.75,H = 80 km. Taking q. = 1000 cm-* WC-‘,(sq,,e)lO-’H 3.7 x lo6 cm-* e.g. for contour number 10, N N 3.7 x 10‘ cm-*.
The numbers on the contours are
/? = 35, y = 1,1=
latitude away from the equator and then falls to an approximately constant level until about 70” and then rises again. The rise again at 70” may be due to the production function used here which is known to become less valid as the poles are approached. At midnight, there is no “equatorial trough” in h,, which in fact falls from a maximum at the equator to a constant level at about 45” latitude and remains constant out to the poles. Also, the height of the maximum at midnight is higher than at noon. A possible explanation of this is that at noon, production, loss and diffusion are active but at night there is no production with the result the loss term becomes more effective and “eats” into the electron density at lower heights causing the peak to rise during the night. Diurnal variation Figure 6 shows the diurnal variation of Nmax at four different latitudes (0”, 30”, 50”, 80”). It is seen that on the equator, the maximum occurs at about 17 hours and as latitude increases, it occurs earlier until at 80” the maximum is at about 134 hours. A possible explanation of this is that at noon on the equator the loss coefficient at the peak is smaller than for higher latitudes. The difference between the equatorial and 30” curves is much greater than the difference between the 30” and 50” curves. In equatorial regions h, is
SOLUTION
100
OF THE TTME-VARYING
I 0
10
DIFFUSION
EQUATION
I
I
I
I
I
I
I
P
30
40
50
60
m
M)
G_
711
lalllU&
Ro. 5. LATITUDE VARIATIONOF THE HEIGHTh,
Ttn
8”
OF THE FUCIRON PEAK.
hO”lr
FIG. 6. DIIJRNAI.VAIUTION OF Nmax. /3= 35, y = 1. A = 1.75, H = 80 km.
higher than for larger latitudes (see Fig. 5), thus the loss coefficient is smaller and so the time of decay of ionization takes longer. Between 30” and 50”, h, is approximately constant and so the loss coefficient is nearly constant, thus the time of decay will be approximately the same. Also it will be noticed that the peak electron density decreases steadily from the equator towards the poles. These results are in general agreement with the results obtained by Gliddon and Kendall.@)
712
ROGER 8. GENERAL
REMARKS
G. BAXTER ON THE TRAWFORMATIONS
Although the electrodynamic drift has been omitted from the above calculations, they do serve as a more consistent base from which the effect of time dependent electrodynamic drifts will cause departure. They are also of interest because of the fruitful transformation y = exp (-2/2H) which is used along a field line, and the subsequent transformation N = yv Owing to the nature of the equation being solved, the transformations as a stabilizing effect, for at great height
(38) and (39) act
v = constant. The transformation to y simplifies the algebra considerably and also y is a solution of the diffusion operator (2). Also, a solution in y gives a much better distribution of points which are relevant to the calculation, especially for very high field lines (i.e. for very large values of a). Further, the transformation maps the field lines into straight lines giving straight boundaries to the problem. (See Section 4.) Finally, using the transformations (38) and (39), we have been able to solve the timedependent diffusion equation, giving a world wide distribution of the electron density from the equator out to the poles. Acknowledgenzents-I should like to thank Professor P. C. Kendall who su ted this investigation. The results were carried out on the Ferranti “Mercury” computer here at SheffieY d. This research was partly sponsored by the Air Force Cambridge Research Laboratories under Grant AF EOAR 64-74 through the European Office of Aerospace Research (OAR), United States Air Force. REFERENCES 1. H. Rrstrn~nr and D. W. BARRON,J. afmm. few. Phys. 18,234 (1960). 2. D. F. MARTYN,Ausr. J. Phys. 9, 161 (1956). 3. J. W. DUNGEY.J. atmos. terr. Phvs. 9.90 (1956). 4. R. A. DUNCAN, Aust. J. Phys. 9,236 (1958). ’ 5. T. J. YONEZAWA,Radio Res. Labs. 5, 165 (1958). 6. J. E. C. GLIDDON,Q. J. mech. appl. Math. 12, 340 (1959). 7. V. C. A. FERRAROand I. ~ZD~GAN, J. atmos. terr. Phys. 12,140 (1958). 8. J. E. C. GL~DD~Nand P. C. KENDALL,J. armos. ten. Phys. 18,48 (1960). 9. J. E. C. GLIDDONand P. C. KENDALL,J. armos. terr. Phys. 24, 1073 (1962). 10. H. RISHBETH, Proc. phys. Sot. 81,65 (1963). 11. H. RISHBETH,J. atrues. ferr. Phys. 26,657 (1964). 12. I. HARRISand W. PRDSTER,NASA Technical Note D-1444 (1962). 13. H. Rrsrn~s.m, A. J. LYON and M. PEART,J. geophys. Res. 68,2559 (1963). 14. R. J. MOFFETT and W. B. HANSON,Nutwe, L0m.f. 206,705 (1965). 15. E. N. BRAMLEYand M. PEART,Nature, Lord. 206,1245 (1965). 16. E. N. BRAMLEYand M. PEART,J. atmos. ferr. Phys. 27, 1201 (1965). 17. W. B. HANSONand R. J. Mom, Private communication (1966). 18. E. V. APPLETON,Nature, Land. 157.691 (1946). 19. L. THOMAS,J. geophys. Res. 71, 1357 (1966). 20. P. C. KENDALL.J. atmos. terr. Phvs. 24.805 (1962). 21. A. J. LYON, J.geophys. Res. 68, i531 (i963). 22. S. CHAPMAN,Proc. Phys. Sot. 43,26 (1931). 23. S. CHAPMAN,Proc. Phys. Sot. 43,483 (1931). G, C. S. G. K. SEITY, and J. 0. THOMAS,Phil. 7kuns. R. Sot. A248, 24. J. A. RATCLIPFE, E. R. SC-
621 (1956). 25. J. CRANKand P. NKXXSON, Proc. Camb. Phil. Sot. 43,50 (1947). 26. National Physical Laboratory, Modem Computing Methods H.M.S.O. (1961).
SOLUTION
OF THE TIME-VARYING
Pf33IOMhPaCCMaTpEiBaeTClI B
CJIOe
F,i?
,qm_$@epeKqtfanbKoe CJlyYae, aT0
KOrAa
4aIOT U
IIyTeM
cnoco6
JJl@I#y3HR, AByX
pe3yJIbTaTbI .
061qnti
MeTOA
UpOABHFKeHHH
paClIpeAeneHKe CBHBIl
C
3JIeKTpOHHOi
¶aCTEl¶HO
ypaBHeHHR
Ii 3aBHCKMOCTb
BAOJIb
pe3yJIbTaTaMU,
IIJIOTHOCTH
HeOAHOpOAHOe,
II0
IIOJl~eAHbIMEi
TOM
IlpKMeHfWI
BpeMeHK.
3JleKTpOHHOt
B
OT BpeMeHK.
JIEIHHH IIOi’lH,
IiHTerpaqHH
lIpeAeJIbHOlf
713
EQUATION
peIJ.IeHHH
IlOTHpR
H KHTerpHpOBaHHeM
AJIU B
AJIfi
IIOpJIJJKa,
IIpOH3BOACTB0,
3aMeH
HUKOJIbCOHa
06CyxAaIOTCH
BTOpOl’O
AaeTCFI
IllEpOKO-paCIlpOCTpaHeHHOe
aBTOpaMEl
Hepa3pblBHOCTEi
aTO-JlHHettHOe,
ypaBKeme. BKJ’fIOUeHbI
AOCTHPaeTCR
KOJleH=iaTbl&i
ypaBHeHHe
EIOHOCtjjepbI.
DIFFUSION
nOJIy-
IIJIOTHOCTH, ApyrHMH