- Email: [email protected]
Numerical modelling of scintillation variations from interplanetary disturbances
Planet. Space Sci., Vol. 35, No. 3, pp. 271-283, 1987 Printed in Great Britair~.
0032~)633/87 $3.00+0.00 © 1987 Pergamon Journals Ltd.
NUMERICAL MODELLING OF SCINTILLATION VARIATIONS F R O M I N T E R P L A N E T A R Y DISTURBANCES S. J . T A P P I N *
Mullard Radio Astronomy Observatory, Cavendish Laboratory, Madingley Road, Cambridge CB3 0HE, U.K.
(Received 14 August 1986) A~tract--The interpretation of interplanetary scintillation observations requires the use of models in order to deduce 3-D structures from observations of a line-of-sightintegral of the scattering. In this paper a system of computer programs which can be used to calculate the enhancements or decrements of interplanetary scintillation caused by arbitrary interplanetary disturbances is described. A number of the resulting calculations are shown and an example of the application of the procedure to a real disturbance is given.
I. I N T R O D U C T I O N
2. T H E P R I N C I P L E S
Observations made over many years have shown interplanetary scintillation (IPS) to be able to monitor large-scale variations in the level of turbulence in the solar wind. These are seen as regions of enhanced or decreased scintillation on the sky (e.g. Houminer and Hewish, 1972; Erskine et al., 1978; Vlasov, 1981; Gapper et al., 1982; Tappin et al. 1983, 1984). These have been found to correlate well with large-scale density structure in the solar wind, with the scintillation enhancement factor (#) and the plasma density (N) obeying the relation g oc N °5 (Tappin, 1986). However, IPS measures a line-of-sight integral of the scattering power and so it is not possible to calculate the 3-D structure of the disturbances directly from observations made from one point (the Earth). This difficulty can be overcome by fitting models to the observations and using our knowledge of what is physically plausible. In order to do this effectively it is essential to be able to calculate the expected appearance of the models and quickly and efficiently so as to allow the iterative process of refinement to proceed as quickly as possible, In this paper a system of computer programs which has been developed to perform the calculations necessary for the interpretation of the Cambridge 81.5 MHz IPS observations is described. A number of example calculations from a set of idealized disturbances which can be used as starting points in the interpretation of real disturbances are presented. In the final section of the paper the application of the modelling to a real disturbance is illustrated,
OF THE MODELS
The model calculations were based on the solar wind model of Readhead et al. (1978) in the region of weak scattering (that is to say, where the r.m.s, phase variations are less than 1 rad, which at 81.5 MHz is anywhere at elongations > 35°). This approximation was chosen because it gives a good theoretical fit to the observations over its region of applicability. No attempt was made to extend the model calculations into the strong scattering region for three reasons. (i) There is no adequate theory to handle the calculations in this region. (ii) The computational cost of including strong scattering calculations, with whatever theory is available, in the models would have greatly outweighed the extra information gained, as will be explained later. (iii) Observations have shown that there is much less variation of scintillation at small elongations than occurs in the weak scattering region (Tappin, 1984). In the calculations the expression for the scintillation index (m) was used in the form given by Kemp (1979) (which is equivalent to that of Readhead et al., 1978),with an extra factor ( f ) to specify the structure of the disturbance being calculated, i.e. ioo (4Z)2 m2 = ~ v f ( r , O, (a)~(r) (1 + h2) - i (1 + h2) 2+ (4Z) 2 dz
*Present address: Space Environment Laboratory, 325 Broadway, Boulder, CO 80303, U.S.A. 271
(1)
272
S.J. TAPPIN
where h = source size term (h 2 = 2z20~/a 2) Z = distance term (Fresnel filter) = z/xa 2 =
z2/2zca2 a = scale size of the medium x = wavenumber of the radio waves z = distance from Earth (not from the point
00 =
fl(r) = q~0 =
f(r, 0, 4) =
of closest approach to the Sun as used by Kemp, 1979) angular radius of a Gaussian source scattering power = d~a~/dz r.m.s, phase variations on the incoming wave enhancement factor of scattering power at position (r, 0, ~b) in the solar wind.
Kemp's (1979) values were used for the scale size and scattering power of the solar wind, namely : a(r) = 175r 1/2 km for0.2 < r < 0.6a.u.
a(r) = 250rS/4kmforr > 0.6a.u. fl(r) = 2.73 × 10-3;t2r4 rad 2 a . u . - ' (2 in m, r in a.u.), The expression (1) assumes that the scale size in the medium does not change significantly in disturbances, which was found to be the case by Houminer (1971). 3. COMPUTATIONALMETHOD In order to be able to calculate the variations of scintillation caused by arbitrarily shaped disturbances in the solar wind, it was necessary to store the integrand of (1) in a grid and integrate through the grid, rather than integrate analytic forms numerically as Readhead et al. (1978) and Kemp (1979) were able to do when calculating the variation of scintillation with elongation for the quiet solar wind. The size and shape of the grid through which the integration was to be performed was determined by the need to satisfy certain constraints. (i) It should be easy to make structures evolve in time in realistic ways. This dictated a heliocentric system of spherical coordinates to facilitate both radial expansion and corotation, (ii)A major disturbance reaching from the outer edge of the grid to oo should not cause a significant enhancement of scintillation. For the parameters appropriate to the observations and models being used this required that the grid should reach out to about 3 a.u. from the Sun. (iii) The whole grid and the programs needed to perform the calculations had to fit into the
available computer memory (500 kbyte~ = 125 kwords), thus limiting the total size of the grid. (iv) Additionally it was desirable that the cells of the grid should be approximately cubic at I a.u. in order to allow the maximum accuracy in representing corotating regions, which typically have a spiral angle of about 45 ° at Earth. The most convenient grid which satisfied all these criteria simultaneously was one with a cell size of 0.1 a.u. x 6 ° × 6 °, reaching out to 3 a.u. from the Sun. This gives a size of 30 x 30 x 60 cells, but of these the innermost 5 shells could be omitted as the calculations were not carried out at elongations < 30 ° so none of the lines of sight integrated went within 0.5 a.u. of the Sun, thus reducing the size required to 25 × 30 × 60 = 45,000 cells. The integrand in (1) is linear in the scattering power (ffl) and so the full integrand can be stored in each cell rather than just ffl, which saves a considerable a m o u n t of computing time. Had strong scattering calculations been included this could not have been done. In addition, the slowing of computation would have been exacerbated by the large numbers of cells traversed by lines of sight at small etongations as they pass through the small cells close to the Sun. The source size enters (1) in a non-linear way and it is therefore necessary to choose a source size for all the calculations performed on any given model. Most o f the calculations presented in this paper use a size of 0':6 [that is, 00 = 0'.'3 in (1)], this being a "typical" source size. The value used does not make a major difference to the appearance of the model, as is shown later. In order to avoid recalulating the integrand of (1) for every model, the values for f = 1 were stored on disc and for any given calculation this was read by the programs and then multiplied by the values o f f The integration was performed by evaluating the s u m m 2 = Zldzo~, over all the cells on a line of sight, where dzce, is the distance travelled by the line of sight in the cell and I is the integrand of (1). The details of the calculations of dz~, are given in the Appendix to this paper. In order to speed up the calculations use was made of the symmetry of the grid, which means that for any line of sight there are three others with the same values of dz~, and simply-related cell numbers, thus reducing by a factor of four the n u m b e r of trigonometrical calculations needed, which is important as these dominate the computational costs. The integration was carried out for lines of sight at 5 ° intervals in latitude and longitude, and elongations > 30 °. The values of m obtained for f = !, which will subsequently be referred to as m0, agree with those of
Modelling scintillation disturbances Readhead et al. (1978) to about 1%. The final values output from the program are of the scintillation enhancement factor g (defined as m/mo), which is what was plotted in the observational maps of Gapper et al. (1982) and of Tappin et al. (1983, 1984). In order to facilitate plotting the maps g was assigned the value 1 at elongations < 30 °.
4. T H E STANDARD MODELS
A large number of calculations were made for idealized disturbances involving density enhancements only. These can then be used as a starting point in the more detailed analysis of individual disturbances, as described by Tappin et al. (1983, 1984) and also in the next section of this paper. In addition they can be used to help make rough estimates of the sizes, speeds, etc. of disturbances which are being considered in less detail, These so-called standard models were all calculated for a source size of 0'.'6. The models calculated were spherical shells, with a wide range of sizes and directions to the Sun Earth line, and corotating sectors with a variety of extents. F o r the streams two different radial variations are used: (a) f = const, and (b)
a)
b)
273
f - 1 oc r 2 for r < 1.2 a.u. and f = const for r > 1.2. In all cases f = 9 at 1 a.u. Some shell models are displayed in Fig. 1, where they are displayed as contour plots o f g in ecliptic coordinates with the Sun at the centre, using Mollweide's equal-area projection. Each frame represents the whole sky at one particular time. F o r a disturbance speed of 430 km s i the frames are 1 day apart. Some of the streams are plotted in Fig. 2 using the same format, but here a speed of 350 km s ~ was used. In order to confirm that the source size in unimportant in determining the main features of the region of enhanced scintillation, a number of models have been calculated for other sizes as well as the standard 0':6. The similarity of the observed structure is confirmed by the plots shown in Fig. 3 where the same model is shown for three different source sizes. There is some difference in the low level contours after the disturbance had passed Earth (i.e. where the Fresnel filter is most important) but nothing of any significance at times important for the determination of structure. The distance to which we are able to see disturbances depends on the strength of the disturbance ; however, for a disturbance of average size or more
c)
d)
@@@ <2> <2><2> FIG. I. (a~t). Caption overleaf.
274
S . J . TAPPIN
el
fl
g)
h)
@ @ @ @
FIG. 1. CONTOUR PLOTSOF SONEOF THE SI-~LL-TYPESTANDARDMODELS. All the models are simple shells centred on the Sun with specified semi-vertical angles, directions to the Sun-Earth line and thicknesses. All lie in the ecliptic plane and all have speeds o f 430 km s - i if the interval between frames is taken to be 24 h. In all the maps the contours are at g = 1.1, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 3.0, 3.5. Below are the parameters of the shells listed in the order: semi-vertical angle, direction, thickness.
(a) (b) (c) (d) (e) (f) (g) (h)
18° 18° 24 ° 24 ° 36 ° 36 ° 48 ° 48 °
with a tripling o f d e n s i t y ( n o t a 9-fold increase as in t h e " s t a n d a r d m o d e l s " ) the e n h a n c e m e n t will be d e t e c t a b l e if t h e n e a r e s t p a r t c o m e s n e a r e r to the S u n E a r t h line t h a n a b o u t 45 °. T h e m a i n d i s t i n g u i s h i n g feature o f d i s t u r b a n c e s lying so far f r o m t h e S u n E a r t h line is t h a t they d o n o t pass b e y o n d a b o u t 90 ° e l o n g a t i o n , b u t r a t h e r t h e y a p p e a r to slow d o w n a n d fade o u t at o r n e a r 90°; d i s t u r b a n c e s n e a r to the S u n - E a r t h line, o n the o t h e r h a n d , are seen at large e l o n g a t i o n s after they have p a s s e d b e y o n d 1 a.u. f r o m the Sun.
18° 30 ° 0° 48 ° 18° 54 ° 48 ° 96 °
0.1 0.1 0.1 0.2 0.2 O. 1 0.2 0.2
5. APPLICATION TO A REAL DISTURBANCE In this section a real d i s t u r b a n c e is c o n s i d e r e d a n d a m o d e l p r o d u c e d for it. T h e d i s t u r b a n c e , for w h i c h t h e m a p s are s h o w n in Fig. 4, was first seen o n 18 S e p t e m b e r 1980. T h e m a p s used here are similar in principle to t h o s e o f T a p p i n etal.(1983, 1984) b u t here a grid o f s o m e 2500 sources is used r a t h e r t h a n the 900 o f the earlier p a p e r s . T h e m a p s s h o w the w h o l e r e g i o n o f the sky visible using the 3.6 h a A r r a y at C a m b r i d g e . T h e y are p l o t t e d in
@@@@@@@@@@@
r
276
S . J . TAPPIN
el
f)
FIG. 2. CONTOUR MAPS OF SOME OF THE COROTATING SECTOR STANDARD MODELS.
The format is similar to that of Fig. 1. All the models have a speed of 350 km s- ~. The parameters are : the longitude range, the latitude range (half angle) and whether there is any radial variation. If there is radial variation then it takes the form ( f - 1) oc r 2. (a) (b) (c) (d) (e) (f)
12° 12° 24 ° 24 ° 48 ° 48 °
30° 30° 30° 30° 30° 30°
no yes no yes no yes
The observations Mollweide's equal-area projection. Each pixel corresponds to a 5 ° square at the ecliptic a n d to equivalent areas at higher latitudes. Because the maps are displayed in ecliptic coordinates a n d time tracks in celestial coordinates it is necessary to m a k e a decision o n where to cut the maps. If this is done in such a way as to prevent the occurrence o f any lines o n the m a p across which time is discontinuous, t h e n it is necessary t h a t some o f the d a t a should lie outside the region defined by the ellipse o f 180 ° longitude. In the maps used here we have used the simplest solution, which is to m a k e the cut along the m i d n i g h t line a n d then follow a great circle from the celestial pole to the ecliptic pole. This scheme has the a d v a n t a g e t h a t each m a p refers to exactly one calendar day.
T h e m a i n features to note in the m a p s s h o w n in Fig. 4 are as follows. 18 S e p t e m b e r : the m a i n feature was a loop o f e n h a n c e d scintillation with its apex at a b o u t 60 ° E, 60 ° N. 19 September : the loop h a d e x p a n d e d to form a n a r r o w b a n d o f e n h a n c e d scintillation a r o u n d 90 ° elongation, b u t curving in towards the Sun at the tails. 20 September: the e n h a n c e d scintillation was still visible, b u t it h a d been followed o u t by a region of decreased scintillation, which filled the inside of the loop. 21 September: the decreased scintillation h a d n o w completely replaced the enhancement.
Modelling scintillation disturbances
a)
b)
277
c)
FIG. 3. CONTOUR PLOTS OF THE CALCULATED SCINTILLATION ENHANCEMt!NTS FOIl T I ~ DISTURBANCE BUT DIFFERENT SOURCE SIZES.
SAME MODEL
The disturbance was a shell 0.2 a.u. thick with a semi-vertical angle of 48 °, just grazing Earth. The source sizes are (a) 0'.'0, (b) 0';6 and (c) 1';0. The contour levels are the same as in Figs 1 and 2.
Direct deductions Certain deductions as to the shape of this disturbance can be made from the observations prior to
Earth it is necessary that the centre of the disturbance should be near 30 ° N, 20 ° E. The details of the models used were :
any attempt to fit a model, It is clear that the disturbance lay to the N E of the S u n - E a r t h line, but that the disturbance hit the Earth sometime on the 19th or 20th. It is also very clear that the disturbance does not fit into the pattern of a spherical shell centred on the Sun, as such a shape cannot give rise to so narrow a band of enhancement as was seen on the 19th, nor to so strongly curved a region as was seen on the 18th. Both of these require that the leading edge should be more strongly curved than such a model would imply. F r o m the levels of scintillation it is apparent that the turbulence was approximately doubled in the shell (at 1 a.u.), and almost halved in the low density region. These correspond to densities of about 4 times and about 0.3 times the normal values, respectively,
outer shell 0.1 a.u. thick at leading edge; scattering
The model F r o m the deductions above it is necessary to find a form which can be reasonably readily calculated in order to provide a starting point for the model. Such a form is that of a spherical bubble, such as is shown in section in Fig. 5. In order to achieve consistency of position angle and the timing of the contact with the
power increased 6-fold at 1 a.u. (should be 4) and proportional to (1 + r 2) to 1.2 a.u. ; symmetry axis at 30 ° N, 20 ° E ; interior scattering power reduced 5-fold at all radii. The maps computed from this model are shown in Fig. 6. The best identification of days is that the third frame in the model corresponds to the 18th, and that the sixth frame corresponds to the 19th. This gives a speed of about 0.3 a.u. d a y - ~ or about 520 km s- ~. What did the real disturbance look like? The model gives a reasonable fit to the observations (apart from the overestimate of the levels). The general character of a shell of enhanced density, followed by a region of low density, is similar to a great many other IPS disturbances, which we believe to be caused by the eruption of short-lived solar wind steams from coronal holes (Hewish et al., 1985 ; Hewish and BravoNunez, 1986). Therefore a good estimate of the real shape is as illustrated in Fig. 7, with the curvature resulting from a fast flow with a smooth velocity profile, varying from the ambient speed at the edges to a high value at the centre.
278
S . J . TAPPIN
FIG. 5. A SECTION THROUGH THE MERIDIONAL PLANE OF THE MODEL FITTED TO THE DISTURBANCE FIRST SEEN ON 18 SEPTEMBER. The shells are s e p a r a t e d 0.1 a.u. at their leading edges a n d for each frame in Fig. 6 the region between a pair o f adjacent shells h a d e n h a n c e d scattering a n d the interior h a d the scattering reduced.
FastdrMngflow
Ecliptic Plane
_----~~Sun FIG. 7. A MORE REALISTIC PICTURE OF THE SHAPE OF THE 18 SEPTEMBER DISTURBANCE, AGAIN IN MERIDIONAL SECTION.
M o d e l l i n g scintillation d i s t u r b a n c e s
FIG. 4. MAPS OF THE DISTRIBUTION OF g ACROSS THE SKY FOR 17 21 SEPTEMBER 1980. The m a p s cover the r a n g e g = 1.5 (red) to y = 0.67 (blue). E a c h pixel is 5 ° s q u a r e at the ecliptic and o f c o r r e s p o n d i n g area at higher latitudes : typically there are 3 5 sources in each pixel.
279
280
S . J . TAPPIN
FIG. 6. MAPS OF THE y-VALUESCALCULATEDFROMTHE MODELSHOWNIN FIG. 5. The format is similar to Fig. 4, but with a greyscale representation ; however, each pixel corresponds to a single integrated line of sight (i.e. 5° intervals in latitude and longitude in all positions).
Modelling scintillation disturbances 6. CONCLUDING REMARKS I n c o n c l u s i o n t h e following p o i n t s o u g h t to be emphasized, (i)It is s h o w n t h a t a simple s y s t e m o f empirical m o d e l l i n g is able to a c c o u n t for the m a j o r features o f the i n t e r p l a n e t a r y d i s t u r b a n c e s seen by IPS. (ii) T h e l e a d i n g edge o f a spherical shell c e n t r e d o n t h e S u n will a p p e a r to follow a line o f c o n s t a n t e l o n g a t i o n ; m o r e s t r o n g l y c u r v e d sky p r o j e c t i o n s r e q u i r e m o r e s t r o n g l y c u r v e d shells, ( i i i ) D i s t u r b a n c e s travelling at large angles to the S u n - E a r t h line d o n o t p a s s 90 ° e l o n g a t i o n , b u t r a t h e r they slow d o w n a n d fade o u t n e a r 90 ° . (iv) M o s t d i s t u r b a n c e s w h i c h c o m e w i t h i n 45 ° o f the S u n - E a r t h line will be seen by IPS. (v) C o m p a r i s o n o f IPS o b s e r v a t i o n s w i t h t h e set o f standard models enables the gross features of m a n y large-scale h e l i o s p h e r i c d i s t u r b a n c e s to be d e t e r m i n e d in a s e m i - q u a n t i t a t i v e way.
Acknowledgements--I thank Professor A. Hewish for suggesting this work and Dr A. Purvis for his contribution to reducing the observational data. I acknowledge the support of an Isaac Newton Studentship at the time this work was
281
a major interplanetary disturbance. Planet. Space Sci. 31, 1171. Tappin, S. J., Hewish, A. and Gapper, G. R. (1984) Tracking a high latitude corotating stream for more than half a solar rotation. Planet. Space Sci. 32, 1273. Vlasov, V..I, (1981)Interplanetary shock waves from observations of scintillation of radio sources. Geomagn. Aeronomy 21, 686. APPENDIX : THE STEPLENGTHS 1. The radial step From Fig. 8 and a simple application of Pythagoras' theorem it is obvious that for zs > 0 then dzr = [(r+dr)2+p~]l/2-zs and for zs < 0 then dzr = --[(r--dr)2-p2]l/2-zs and in the casep > r - d r then dzr = 2(r 2 _p2)~/2
2. Azimuthal step Here it is clear from Fig. 9 and the sine rule that dz, r' r sin dq~ sin q~ sin (q +d~b) for a line of sight lying in the plane of the figure, so for a line of sight an an elevation ~ we have
carriedassociateshipOUt, andatNOAA/ERL/SEL.the current support of a N R C research
dz,
r~ sin dq~ 1 sin (r/+dqb) cos ~k
where re is the distance from the Sun projected into the plane of the ecliptic as shown in Fig. 10 such that REFERENCES
r e = z tan 0 sin ~O.
Erskine, F. T., Cronyn, W. M., Shawhan, S. D., Roelof, E. C. and Gotwols, B L. (1978) Inter:planetary scintillations at large elongation angles : response to solar wind density structure. J. yeophys. Res. 83, 4153. Gapper, G. R., Hewish, A., Purvis, A. and Duffett-Smith, P. J. (1982) Observing interplanetary disturbances from the ground. Nature, Lond. 296, 633. Hewish, A. and Bravo-Nunez, S. S. (1986) The sources of large-scale heliospheric disturbances. Solar Phys. 106, 185. Hewish, A., Tappin, S. J. and Gapper, G. R. (1985) The origin of strong interplanetary shocks. Nature, Lond. 314, 137. Houminer, Z. (1971) Corotating plasma streams revealed by interplanetary scintillation. Nature Phys. Sci. 231, 165. Houminer, Z. and Hewish, A. (1972) Long-lived sectors of enhanced density irregularities in the solar wind. Planet. Space ScL 20, 1703. Kemp, M. C. (1979) Density fluctuations in the solar wind and the interpretation of IPS observations. PhD thesis, University of Cambridge. Readhead, A. C. S., Kemp, M. C. and Hewish, A. (1978) The spectrum of small scale density fluctuations in the solar wind. Mon. Not. R. astr. Soe. 185, 207. Tappin, S. J. (1984) Transient disturbances in the solar wind. PhD Thesis, University of Cambridge. Tappin, S. J. (1986) Interplanetary scintillation and plasma density. Planet. Space Sei. 34. 93. Tappin, S. J., Hewish, A. and Gapper, G. R. (1983) Tracking
3. Latitudinal step The best way to obtain the step in the latitudinal direction is to calculate the distance from the Earth to each of the two relevant cones of constant latitude and then calculate the difference ; now : tan 0 - (1 + z2 cos: ~O- 2zcos qJcos ~) I/2 zsin~k using the notation of Fig. I0, which can be rearranged to give.
z2(cos 2~b- sin 2q, tan 2 0 ) - 2z cos ~bcos ~ + 1 = 0. Hence cos e + [cos2 e + (sin 2 qJtan z 0 - cos 2 if)] ~/2 (cos2 ~ksin 2 ~, tan 2 0) where cos e = cos ~ cos ~k.The appropriate root is the positive or the negative one according to whether Iq~lis greater or less than ~/2. Then dz = z(O2)-z(G). z-
4. Combining the step lengths The grid was set up so that the Earth was at an intersection of all three types of boundary. So determining the first step was simply a matter of calculating the three types of step away and then choosing the shortest. Then for subsequent steps, recalculate the last type of step used, subtract the previous step from the other two, select the smallest and repeat.
282
S . J . TAPPIN
~
dz~
Zs2
r-dr
r
dz~
~
~
Zsl
FIG. 8.
CONSTRUCTION ILLUSTRATING THE DETERMINATION OF THE RADIAL STEP IN THE MODEL CALCULATIONS.
J r FIG. 9. CONSTRUCTION ILLUSTRATING THE DETERMINATION OF THE ECLIPTIC PROJECTION OF THE LONGITUDINAL STEP IN THE MODEL CALCULATIONS.
Modelling scintillation disturbances
283
~
°
i'i~iil~ i~i!iii!if!i!i if!!iiiii i ........ o~o~'~~ ~
~
!
~
~
| S
E FIG. 10.
CONSTRUCTION ILLUSTRATING THE CALCULATION OF THE LONGITUDINAL STEP FROM ITS ECLIPTIC PROJECTION, AND THE CALCULATION OF THE LATITUDINAL STEP.
0032~)633/87 $3.00+0.00 © 1987 Pergamon Journals Ltd.
NUMERICAL MODELLING OF SCINTILLATION VARIATIONS F R O M I N T E R P L A N E T A R Y DISTURBANCES S. J . T A P P I N *
Mullard Radio Astronomy Observatory, Cavendish Laboratory, Madingley Road, Cambridge CB3 0HE, U.K.
(Received 14 August 1986) A~tract--The interpretation of interplanetary scintillation observations requires the use of models in order to deduce 3-D structures from observations of a line-of-sightintegral of the scattering. In this paper a system of computer programs which can be used to calculate the enhancements or decrements of interplanetary scintillation caused by arbitrary interplanetary disturbances is described. A number of the resulting calculations are shown and an example of the application of the procedure to a real disturbance is given.
I. I N T R O D U C T I O N
2. T H E P R I N C I P L E S
Observations made over many years have shown interplanetary scintillation (IPS) to be able to monitor large-scale variations in the level of turbulence in the solar wind. These are seen as regions of enhanced or decreased scintillation on the sky (e.g. Houminer and Hewish, 1972; Erskine et al., 1978; Vlasov, 1981; Gapper et al., 1982; Tappin et al. 1983, 1984). These have been found to correlate well with large-scale density structure in the solar wind, with the scintillation enhancement factor (#) and the plasma density (N) obeying the relation g oc N °5 (Tappin, 1986). However, IPS measures a line-of-sight integral of the scattering power and so it is not possible to calculate the 3-D structure of the disturbances directly from observations made from one point (the Earth). This difficulty can be overcome by fitting models to the observations and using our knowledge of what is physically plausible. In order to do this effectively it is essential to be able to calculate the expected appearance of the models and quickly and efficiently so as to allow the iterative process of refinement to proceed as quickly as possible, In this paper a system of computer programs which has been developed to perform the calculations necessary for the interpretation of the Cambridge 81.5 MHz IPS observations is described. A number of example calculations from a set of idealized disturbances which can be used as starting points in the interpretation of real disturbances are presented. In the final section of the paper the application of the modelling to a real disturbance is illustrated,
OF THE MODELS
The model calculations were based on the solar wind model of Readhead et al. (1978) in the region of weak scattering (that is to say, where the r.m.s, phase variations are less than 1 rad, which at 81.5 MHz is anywhere at elongations > 35°). This approximation was chosen because it gives a good theoretical fit to the observations over its region of applicability. No attempt was made to extend the model calculations into the strong scattering region for three reasons. (i) There is no adequate theory to handle the calculations in this region. (ii) The computational cost of including strong scattering calculations, with whatever theory is available, in the models would have greatly outweighed the extra information gained, as will be explained later. (iii) Observations have shown that there is much less variation of scintillation at small elongations than occurs in the weak scattering region (Tappin, 1984). In the calculations the expression for the scintillation index (m) was used in the form given by Kemp (1979) (which is equivalent to that of Readhead et al., 1978),with an extra factor ( f ) to specify the structure of the disturbance being calculated, i.e. ioo (4Z)2 m2 = ~ v f ( r , O, (a)~(r) (1 + h2) - i (1 + h2) 2+ (4Z) 2 dz
*Present address: Space Environment Laboratory, 325 Broadway, Boulder, CO 80303, U.S.A. 271
(1)
272
S.J. TAPPIN
where h = source size term (h 2 = 2z20~/a 2) Z = distance term (Fresnel filter) = z/xa 2 =
z2/2zca2 a = scale size of the medium x = wavenumber of the radio waves z = distance from Earth (not from the point
00 =
fl(r) = q~0 =
f(r, 0, 4) =
of closest approach to the Sun as used by Kemp, 1979) angular radius of a Gaussian source scattering power = d~a~/dz r.m.s, phase variations on the incoming wave enhancement factor of scattering power at position (r, 0, ~b) in the solar wind.
Kemp's (1979) values were used for the scale size and scattering power of the solar wind, namely : a(r) = 175r 1/2 km for0.2 < r < 0.6a.u.
a(r) = 250rS/4kmforr > 0.6a.u. fl(r) = 2.73 × 10-3;t2r4 rad 2 a . u . - ' (2 in m, r in a.u.), The expression (1) assumes that the scale size in the medium does not change significantly in disturbances, which was found to be the case by Houminer (1971). 3. COMPUTATIONALMETHOD In order to be able to calculate the variations of scintillation caused by arbitrarily shaped disturbances in the solar wind, it was necessary to store the integrand of (1) in a grid and integrate through the grid, rather than integrate analytic forms numerically as Readhead et al. (1978) and Kemp (1979) were able to do when calculating the variation of scintillation with elongation for the quiet solar wind. The size and shape of the grid through which the integration was to be performed was determined by the need to satisfy certain constraints. (i) It should be easy to make structures evolve in time in realistic ways. This dictated a heliocentric system of spherical coordinates to facilitate both radial expansion and corotation, (ii)A major disturbance reaching from the outer edge of the grid to oo should not cause a significant enhancement of scintillation. For the parameters appropriate to the observations and models being used this required that the grid should reach out to about 3 a.u. from the Sun. (iii) The whole grid and the programs needed to perform the calculations had to fit into the
available computer memory (500 kbyte~ = 125 kwords), thus limiting the total size of the grid. (iv) Additionally it was desirable that the cells of the grid should be approximately cubic at I a.u. in order to allow the maximum accuracy in representing corotating regions, which typically have a spiral angle of about 45 ° at Earth. The most convenient grid which satisfied all these criteria simultaneously was one with a cell size of 0.1 a.u. x 6 ° × 6 °, reaching out to 3 a.u. from the Sun. This gives a size of 30 x 30 x 60 cells, but of these the innermost 5 shells could be omitted as the calculations were not carried out at elongations < 30 ° so none of the lines of sight integrated went within 0.5 a.u. of the Sun, thus reducing the size required to 25 × 30 × 60 = 45,000 cells. The integrand in (1) is linear in the scattering power (ffl) and so the full integrand can be stored in each cell rather than just ffl, which saves a considerable a m o u n t of computing time. Had strong scattering calculations been included this could not have been done. In addition, the slowing of computation would have been exacerbated by the large numbers of cells traversed by lines of sight at small etongations as they pass through the small cells close to the Sun. The source size enters (1) in a non-linear way and it is therefore necessary to choose a source size for all the calculations performed on any given model. Most o f the calculations presented in this paper use a size of 0':6 [that is, 00 = 0'.'3 in (1)], this being a "typical" source size. The value used does not make a major difference to the appearance of the model, as is shown later. In order to avoid recalulating the integrand of (1) for every model, the values for f = 1 were stored on disc and for any given calculation this was read by the programs and then multiplied by the values o f f The integration was performed by evaluating the s u m m 2 = Zldzo~, over all the cells on a line of sight, where dzce, is the distance travelled by the line of sight in the cell and I is the integrand of (1). The details of the calculations of dz~, are given in the Appendix to this paper. In order to speed up the calculations use was made of the symmetry of the grid, which means that for any line of sight there are three others with the same values of dz~, and simply-related cell numbers, thus reducing by a factor of four the n u m b e r of trigonometrical calculations needed, which is important as these dominate the computational costs. The integration was carried out for lines of sight at 5 ° intervals in latitude and longitude, and elongations > 30 °. The values of m obtained for f = !, which will subsequently be referred to as m0, agree with those of
Modelling scintillation disturbances Readhead et al. (1978) to about 1%. The final values output from the program are of the scintillation enhancement factor g (defined as m/mo), which is what was plotted in the observational maps of Gapper et al. (1982) and of Tappin et al. (1983, 1984). In order to facilitate plotting the maps g was assigned the value 1 at elongations < 30 °.
4. T H E STANDARD MODELS
A large number of calculations were made for idealized disturbances involving density enhancements only. These can then be used as a starting point in the more detailed analysis of individual disturbances, as described by Tappin et al. (1983, 1984) and also in the next section of this paper. In addition they can be used to help make rough estimates of the sizes, speeds, etc. of disturbances which are being considered in less detail, These so-called standard models were all calculated for a source size of 0'.'6. The models calculated were spherical shells, with a wide range of sizes and directions to the Sun Earth line, and corotating sectors with a variety of extents. F o r the streams two different radial variations are used: (a) f = const, and (b)
a)
b)
273
f - 1 oc r 2 for r < 1.2 a.u. and f = const for r > 1.2. In all cases f = 9 at 1 a.u. Some shell models are displayed in Fig. 1, where they are displayed as contour plots o f g in ecliptic coordinates with the Sun at the centre, using Mollweide's equal-area projection. Each frame represents the whole sky at one particular time. F o r a disturbance speed of 430 km s i the frames are 1 day apart. Some of the streams are plotted in Fig. 2 using the same format, but here a speed of 350 km s ~ was used. In order to confirm that the source size in unimportant in determining the main features of the region of enhanced scintillation, a number of models have been calculated for other sizes as well as the standard 0':6. The similarity of the observed structure is confirmed by the plots shown in Fig. 3 where the same model is shown for three different source sizes. There is some difference in the low level contours after the disturbance had passed Earth (i.e. where the Fresnel filter is most important) but nothing of any significance at times important for the determination of structure. The distance to which we are able to see disturbances depends on the strength of the disturbance ; however, for a disturbance of average size or more
c)
d)
@@@ <2> <2><2> FIG. I. (a~t). Caption overleaf.
274
S . J . TAPPIN
el
fl
g)
h)
@ @ @ @
FIG. 1. CONTOUR PLOTSOF SONEOF THE SI-~LL-TYPESTANDARDMODELS. All the models are simple shells centred on the Sun with specified semi-vertical angles, directions to the Sun-Earth line and thicknesses. All lie in the ecliptic plane and all have speeds o f 430 km s - i if the interval between frames is taken to be 24 h. In all the maps the contours are at g = 1.1, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 3.0, 3.5. Below are the parameters of the shells listed in the order: semi-vertical angle, direction, thickness.
(a) (b) (c) (d) (e) (f) (g) (h)
18° 18° 24 ° 24 ° 36 ° 36 ° 48 ° 48 °
with a tripling o f d e n s i t y ( n o t a 9-fold increase as in t h e " s t a n d a r d m o d e l s " ) the e n h a n c e m e n t will be d e t e c t a b l e if t h e n e a r e s t p a r t c o m e s n e a r e r to the S u n E a r t h line t h a n a b o u t 45 °. T h e m a i n d i s t i n g u i s h i n g feature o f d i s t u r b a n c e s lying so far f r o m t h e S u n E a r t h line is t h a t they d o n o t pass b e y o n d a b o u t 90 ° e l o n g a t i o n , b u t r a t h e r t h e y a p p e a r to slow d o w n a n d fade o u t at o r n e a r 90°; d i s t u r b a n c e s n e a r to the S u n - E a r t h line, o n the o t h e r h a n d , are seen at large e l o n g a t i o n s after they have p a s s e d b e y o n d 1 a.u. f r o m the Sun.
18° 30 ° 0° 48 ° 18° 54 ° 48 ° 96 °
0.1 0.1 0.1 0.2 0.2 O. 1 0.2 0.2
5. APPLICATION TO A REAL DISTURBANCE In this section a real d i s t u r b a n c e is c o n s i d e r e d a n d a m o d e l p r o d u c e d for it. T h e d i s t u r b a n c e , for w h i c h t h e m a p s are s h o w n in Fig. 4, was first seen o n 18 S e p t e m b e r 1980. T h e m a p s used here are similar in principle to t h o s e o f T a p p i n etal.(1983, 1984) b u t here a grid o f s o m e 2500 sources is used r a t h e r t h a n the 900 o f the earlier p a p e r s . T h e m a p s s h o w the w h o l e r e g i o n o f the sky visible using the 3.6 h a A r r a y at C a m b r i d g e . T h e y are p l o t t e d in
@@@@@@@@@@@
r
276
S . J . TAPPIN
el
f)
FIG. 2. CONTOUR MAPS OF SOME OF THE COROTATING SECTOR STANDARD MODELS.
The format is similar to that of Fig. 1. All the models have a speed of 350 km s- ~. The parameters are : the longitude range, the latitude range (half angle) and whether there is any radial variation. If there is radial variation then it takes the form ( f - 1) oc r 2. (a) (b) (c) (d) (e) (f)
12° 12° 24 ° 24 ° 48 ° 48 °
30° 30° 30° 30° 30° 30°
no yes no yes no yes
The observations Mollweide's equal-area projection. Each pixel corresponds to a 5 ° square at the ecliptic a n d to equivalent areas at higher latitudes. Because the maps are displayed in ecliptic coordinates a n d time tracks in celestial coordinates it is necessary to m a k e a decision o n where to cut the maps. If this is done in such a way as to prevent the occurrence o f any lines o n the m a p across which time is discontinuous, t h e n it is necessary t h a t some o f the d a t a should lie outside the region defined by the ellipse o f 180 ° longitude. In the maps used here we have used the simplest solution, which is to m a k e the cut along the m i d n i g h t line a n d then follow a great circle from the celestial pole to the ecliptic pole. This scheme has the a d v a n t a g e t h a t each m a p refers to exactly one calendar day.
T h e m a i n features to note in the m a p s s h o w n in Fig. 4 are as follows. 18 S e p t e m b e r : the m a i n feature was a loop o f e n h a n c e d scintillation with its apex at a b o u t 60 ° E, 60 ° N. 19 September : the loop h a d e x p a n d e d to form a n a r r o w b a n d o f e n h a n c e d scintillation a r o u n d 90 ° elongation, b u t curving in towards the Sun at the tails. 20 September: the e n h a n c e d scintillation was still visible, b u t it h a d been followed o u t by a region of decreased scintillation, which filled the inside of the loop. 21 September: the decreased scintillation h a d n o w completely replaced the enhancement.
Modelling scintillation disturbances
a)
b)
277
c)
FIG. 3. CONTOUR PLOTS OF THE CALCULATED SCINTILLATION ENHANCEMt!NTS FOIl T I ~ DISTURBANCE BUT DIFFERENT SOURCE SIZES.
SAME MODEL
The disturbance was a shell 0.2 a.u. thick with a semi-vertical angle of 48 °, just grazing Earth. The source sizes are (a) 0'.'0, (b) 0';6 and (c) 1';0. The contour levels are the same as in Figs 1 and 2.
Direct deductions Certain deductions as to the shape of this disturbance can be made from the observations prior to
Earth it is necessary that the centre of the disturbance should be near 30 ° N, 20 ° E. The details of the models used were :
any attempt to fit a model, It is clear that the disturbance lay to the N E of the S u n - E a r t h line, but that the disturbance hit the Earth sometime on the 19th or 20th. It is also very clear that the disturbance does not fit into the pattern of a spherical shell centred on the Sun, as such a shape cannot give rise to so narrow a band of enhancement as was seen on the 19th, nor to so strongly curved a region as was seen on the 18th. Both of these require that the leading edge should be more strongly curved than such a model would imply. F r o m the levels of scintillation it is apparent that the turbulence was approximately doubled in the shell (at 1 a.u.), and almost halved in the low density region. These correspond to densities of about 4 times and about 0.3 times the normal values, respectively,
outer shell 0.1 a.u. thick at leading edge; scattering
The model F r o m the deductions above it is necessary to find a form which can be reasonably readily calculated in order to provide a starting point for the model. Such a form is that of a spherical bubble, such as is shown in section in Fig. 5. In order to achieve consistency of position angle and the timing of the contact with the
power increased 6-fold at 1 a.u. (should be 4) and proportional to (1 + r 2) to 1.2 a.u. ; symmetry axis at 30 ° N, 20 ° E ; interior scattering power reduced 5-fold at all radii. The maps computed from this model are shown in Fig. 6. The best identification of days is that the third frame in the model corresponds to the 18th, and that the sixth frame corresponds to the 19th. This gives a speed of about 0.3 a.u. d a y - ~ or about 520 km s- ~. What did the real disturbance look like? The model gives a reasonable fit to the observations (apart from the overestimate of the levels). The general character of a shell of enhanced density, followed by a region of low density, is similar to a great many other IPS disturbances, which we believe to be caused by the eruption of short-lived solar wind steams from coronal holes (Hewish et al., 1985 ; Hewish and BravoNunez, 1986). Therefore a good estimate of the real shape is as illustrated in Fig. 7, with the curvature resulting from a fast flow with a smooth velocity profile, varying from the ambient speed at the edges to a high value at the centre.
278
S . J . TAPPIN
FIG. 5. A SECTION THROUGH THE MERIDIONAL PLANE OF THE MODEL FITTED TO THE DISTURBANCE FIRST SEEN ON 18 SEPTEMBER. The shells are s e p a r a t e d 0.1 a.u. at their leading edges a n d for each frame in Fig. 6 the region between a pair o f adjacent shells h a d e n h a n c e d scattering a n d the interior h a d the scattering reduced.
FastdrMngflow
Ecliptic Plane
_----~~Sun FIG. 7. A MORE REALISTIC PICTURE OF THE SHAPE OF THE 18 SEPTEMBER DISTURBANCE, AGAIN IN MERIDIONAL SECTION.
M o d e l l i n g scintillation d i s t u r b a n c e s
FIG. 4. MAPS OF THE DISTRIBUTION OF g ACROSS THE SKY FOR 17 21 SEPTEMBER 1980. The m a p s cover the r a n g e g = 1.5 (red) to y = 0.67 (blue). E a c h pixel is 5 ° s q u a r e at the ecliptic and o f c o r r e s p o n d i n g area at higher latitudes : typically there are 3 5 sources in each pixel.
279
280
S . J . TAPPIN
FIG. 6. MAPS OF THE y-VALUESCALCULATEDFROMTHE MODELSHOWNIN FIG. 5. The format is similar to Fig. 4, but with a greyscale representation ; however, each pixel corresponds to a single integrated line of sight (i.e. 5° intervals in latitude and longitude in all positions).
Modelling scintillation disturbances 6. CONCLUDING REMARKS I n c o n c l u s i o n t h e following p o i n t s o u g h t to be emphasized, (i)It is s h o w n t h a t a simple s y s t e m o f empirical m o d e l l i n g is able to a c c o u n t for the m a j o r features o f the i n t e r p l a n e t a r y d i s t u r b a n c e s seen by IPS. (ii) T h e l e a d i n g edge o f a spherical shell c e n t r e d o n t h e S u n will a p p e a r to follow a line o f c o n s t a n t e l o n g a t i o n ; m o r e s t r o n g l y c u r v e d sky p r o j e c t i o n s r e q u i r e m o r e s t r o n g l y c u r v e d shells, ( i i i ) D i s t u r b a n c e s travelling at large angles to the S u n - E a r t h line d o n o t p a s s 90 ° e l o n g a t i o n , b u t r a t h e r they slow d o w n a n d fade o u t n e a r 90 ° . (iv) M o s t d i s t u r b a n c e s w h i c h c o m e w i t h i n 45 ° o f the S u n - E a r t h line will be seen by IPS. (v) C o m p a r i s o n o f IPS o b s e r v a t i o n s w i t h t h e set o f standard models enables the gross features of m a n y large-scale h e l i o s p h e r i c d i s t u r b a n c e s to be d e t e r m i n e d in a s e m i - q u a n t i t a t i v e way.
Acknowledgements--I thank Professor A. Hewish for suggesting this work and Dr A. Purvis for his contribution to reducing the observational data. I acknowledge the support of an Isaac Newton Studentship at the time this work was
281
a major interplanetary disturbance. Planet. Space Sci. 31, 1171. Tappin, S. J., Hewish, A. and Gapper, G. R. (1984) Tracking a high latitude corotating stream for more than half a solar rotation. Planet. Space Sci. 32, 1273. Vlasov, V..I, (1981)Interplanetary shock waves from observations of scintillation of radio sources. Geomagn. Aeronomy 21, 686. APPENDIX : THE STEPLENGTHS 1. The radial step From Fig. 8 and a simple application of Pythagoras' theorem it is obvious that for zs > 0 then dzr = [(r+dr)2+p~]l/2-zs and for zs < 0 then dzr = --[(r--dr)2-p2]l/2-zs and in the casep > r - d r then dzr = 2(r 2 _p2)~/2
2. Azimuthal step Here it is clear from Fig. 9 and the sine rule that dz, r' r sin dq~ sin q~ sin (q +d~b) for a line of sight lying in the plane of the figure, so for a line of sight an an elevation ~ we have
carriedassociateshipOUt, andatNOAA/ERL/SEL.the current support of a N R C research
dz,
r~ sin dq~ 1 sin (r/+dqb) cos ~k
where re is the distance from the Sun projected into the plane of the ecliptic as shown in Fig. 10 such that REFERENCES
r e = z tan 0 sin ~O.
Erskine, F. T., Cronyn, W. M., Shawhan, S. D., Roelof, E. C. and Gotwols, B L. (1978) Inter:planetary scintillations at large elongation angles : response to solar wind density structure. J. yeophys. Res. 83, 4153. Gapper, G. R., Hewish, A., Purvis, A. and Duffett-Smith, P. J. (1982) Observing interplanetary disturbances from the ground. Nature, Lond. 296, 633. Hewish, A. and Bravo-Nunez, S. S. (1986) The sources of large-scale heliospheric disturbances. Solar Phys. 106, 185. Hewish, A., Tappin, S. J. and Gapper, G. R. (1985) The origin of strong interplanetary shocks. Nature, Lond. 314, 137. Houminer, Z. (1971) Corotating plasma streams revealed by interplanetary scintillation. Nature Phys. Sci. 231, 165. Houminer, Z. and Hewish, A. (1972) Long-lived sectors of enhanced density irregularities in the solar wind. Planet. Space ScL 20, 1703. Kemp, M. C. (1979) Density fluctuations in the solar wind and the interpretation of IPS observations. PhD thesis, University of Cambridge. Readhead, A. C. S., Kemp, M. C. and Hewish, A. (1978) The spectrum of small scale density fluctuations in the solar wind. Mon. Not. R. astr. Soe. 185, 207. Tappin, S. J. (1984) Transient disturbances in the solar wind. PhD Thesis, University of Cambridge. Tappin, S. J. (1986) Interplanetary scintillation and plasma density. Planet. Space Sei. 34. 93. Tappin, S. J., Hewish, A. and Gapper, G. R. (1983) Tracking
3. Latitudinal step The best way to obtain the step in the latitudinal direction is to calculate the distance from the Earth to each of the two relevant cones of constant latitude and then calculate the difference ; now : tan 0 - (1 + z2 cos: ~O- 2zcos qJcos ~) I/2 zsin~k using the notation of Fig. I0, which can be rearranged to give.
z2(cos 2~b- sin 2q, tan 2 0 ) - 2z cos ~bcos ~ + 1 = 0. Hence cos e + [cos2 e + (sin 2 qJtan z 0 - cos 2 if)] ~/2 (cos2 ~ksin 2 ~, tan 2 0) where cos e = cos ~ cos ~k.The appropriate root is the positive or the negative one according to whether Iq~lis greater or less than ~/2. Then dz = z(O2)-z(G). z-
4. Combining the step lengths The grid was set up so that the Earth was at an intersection of all three types of boundary. So determining the first step was simply a matter of calculating the three types of step away and then choosing the shortest. Then for subsequent steps, recalculate the last type of step used, subtract the previous step from the other two, select the smallest and repeat.
282
S . J . TAPPIN
~
dz~
Zs2
r-dr
r
dz~
~
~
Zsl
FIG. 8.
CONSTRUCTION ILLUSTRATING THE DETERMINATION OF THE RADIAL STEP IN THE MODEL CALCULATIONS.
J r FIG. 9. CONSTRUCTION ILLUSTRATING THE DETERMINATION OF THE ECLIPTIC PROJECTION OF THE LONGITUDINAL STEP IN THE MODEL CALCULATIONS.
Modelling scintillation disturbances
283
~
°
i'i~iil~ i~i!iii!if!i!i if!!iiiii i ........ o~o~'~~ ~
~
!
~
~
| S
E FIG. 10.
CONSTRUCTION ILLUSTRATING THE CALCULATION OF THE LONGITUDINAL STEP FROM ITS ECLIPTIC PROJECTION, AND THE CALCULATION OF THE LATITUDINAL STEP.