Cosmic radio noise absorption events associated with equatorward drifting arcs during a substorm growth phase

Abstract. Cosmic radio noise absorption (CNA) events associated with equatorward drifting arcs during a substorm growth phase are studied by using simultaneous optical auroral, IRIS imaging riometer and EISCAT incoherent scatter radar measurements. The CNA is generally attributed to energetic particle precipitation in the D-region. However, it has been argued that plasma irregularities or enhanced electron temperature (Te) in the E-region could also produce CNA. Both of the latter mechanisms are related to intense electric fields in the ionosphere. We present two events which occur during a substorm growth phase in the evening MLT sector. In both of the events, an auroral arc is drifting equatorward, together with a region of CNA (auroral absorption bay) located on the equatorward side and outside of the arc. Both of the events are associated with enhanced D-region electron density on the equatorward side of the auroral arc, but in the second event, a region of intense electric field and enhanced electron temperature in the E-region is also located on the equatorward side of the arc. We show that in the studied events neither plasma instabilities nor enhanced Te play a significant role in producing the measured CNA, but the CNA in the vicinity of the equatorward drifting arcs is produced by D-region energetic electron precipitation. Key words. Ionosphere (auroral ionosphere; particle precipitation; electric fields and currents)


Introduction
The height-integrated ionospheric attenuation of high frequency cosmic radio waves, i.e. cosmic radio noise absorption (CNA) is measured by riometers operating in the frequency range of 20-50 MHz. Under disturbed conditions the CNA is generally attributed to the enhancement of Dregion ionization produced by particle precipitation. The Dregion can be ionized by energetic electrons (≥30 keV) or Correspondence to: J. R. T. Jussila (jouni.jussila@oulu.fi) protons of MeV range. It has also been suggested that CNA could be produced in the E-region by plasma irregularities (D'Angelo, 1976) or by an enhanced electron temperature (Stauning, 1984). Both of these E-region mechanisms are related to strong ionospheric electric fields.
The aim of this paper is to study the relative importance of the conventional D-region precipitation and the proposed E-region processes producing CNA in the vicinity of quiet auroral arcs occurring during a substorm growth phase.
It is well known that auroral arcs are drifting equatorward during the growth phase of a substorm. The development of auroral absorption substorms has been studied in several papers (Hargreaves et al., 1975;Pytte et al., 1976;Ranta et al., 1981;Pellinen et al., 1982;Ranta and Ranta, 1984;Ranta et al., 1999). Studies show a pre-onset absorption bay that often precedes an onset of a substorm. The absorption bay is extended in the east-west direction, moves slowly equatorward and is typically 50-100 km wide latitudinally. Optical data has very seldom been available in absorption studies. The study by Pellinen et al. (1982) showed that the maximum of CNA was located on the equatorward side of the growth phase auroral arc. Also, many studies have utilized wide-beam riometer antennas, and, therefore, it has not been possible to determine the extent of the actual absorption area with a high accuracy. Stoker et al. (1997) studied the relationship of CNA and discrete auroral arcs or pulsating auroral patches by an imaging riometer during disturbed conditions and they concluded that the CNA occurred adjacent to the auroral forms. Sergeev et al. (1990) proposed that the observed CNA during a substorm growth phase was due to pitch-angle scattering of energetic electrons in the equatorial current sheet and the equatorward movement of the CNA region was due to configuration changes in the magnetotail. On the other hand, Williams et al. (1990), Opgenoorth et al. (1990), Aikio et al. (1993) and Aikio et al. (2002) have shown that an enhanced meridional electric field and an accompanying enhanced Eregion electron temperature are often observed outside of an auroral arc, on the equatorward side of an eveningside arc. In addition, there is evidence that CNA can occur outside of an arc and that the enhanced electric field associated with drifting quiet arc is also located outside of the arc. Our aim is to study whether these two phenomena (CNA and enhanced electric field) are related causally.
We have looked for events of equatorward drifting quiet growth phase arcs with simultaneous measurements of 1) auroras by the Kilpisjärvi (KIL) all-sky camera (ASC), 2) CNA by the IRIS imaging riometer at Kilpisjärvi and 3) plasma parameters (N e , T e , T i and E) by the EISCAT incoherent scatter facility in Tromsø within the field-of-view (fov) of both the KIL ASC and IRIS. We selected two good events: 7 December 1999 and 8 April 1998. In the latter event the auroral arc is associated with a region of enhanced electric field and increased T e in the E-region on the equatorward side of the arc, as discussed by Aikio et al. (2002).
For these events, the CNA measured by IRIS will be compared to the CNA calculated from the EISCAT data using the equation of the conventional nondeviative absorption. This equation is formulated mainly for D-region absorption phenomena and if E-region instabilities play a major role, the calculated CNA should be significantly less than the measured one.

Instrumentation
The CNA was measured by the IRIS (Imaging Riometer for Ionospheric Studies) system, which is located at Kilpisjärvi (69.05 • N, 20.79 • E). IRIS consists of 49 (7×7) narrow beams with widths between 13 • and 16 • , and it operates at a frequency of 38.2 MHz (Browne et al., 1995). The basic time resolution is 1 s, but in this study we have used integration times of 10 s. The calculated absorption was mapped to the altitude of 80 km in both events, and the data was interpolated to produce a grid comprising 67. 8-70.2 • N and 17.75-23.75 • E, with a resolution of 0.1 • N×0.25 • E.
Auroras were recorded by an all-sky camera of the MIR-ACLE network (Syrjäsuo et al., 1998) at Kilpisjärvi. The IRIS fov is located inside that of the ASC. The wavelength of 557.7 nm originating from the O 1 S emission was used for the interference filter of the ASC. The time resolution of the ASC during both of the events was 20 s. The exposure time of each image was 1 s. ASC images were mapped to the geographical grid by assuming the lower border of auroras at 100 km (Event 1) or 120 km (Event 2).
The EISCAT (European Incoherent SCATter radar) UHFradar (Risbeth and van Eyken, 1993) was measuring during both of the events in the CP-1 magnetic field-aligned mode. The basic time resolution of the data was 5 s, but the data were post-integrated to 60 or 120 s resolution. The altitude resolution of the data used was 3 km. The CP-1 experiment contained several pulse codes: long pulses, alternating codes and power profiles. The use of EISCAT data is explained in more detail in the following section.
2.2 Method to calculate the CNA from EISCAT measurements The nondeviative absorption of a vertically incident radio wave is given, for example, by Hargreaves (1992) A(db) = 4.6 · 10 −5 h N e ν en ν 2 en + (ω − e ) 2 dh , where N e is the electron density in m −3 , ν en the electronneutral collision frequency, ω the angular frequency of the radio wave and e the electron gyrofrequency.
To calculate the absorption with Eq. (1), the electron density was obtained from the EISCAT measurement. The alternating code pulse scheme extended from 86 to 270 km (Event 1) or to 191 km (Event 2), and it provided true electron densities. At lower heights, only the power profile data was available. The relationship between the true electron density N e and the so-called raw electron density N e raw provided by the power profiles is where λ d =(ε 0 k B T e /N e e 2 ) 1/2 is the Debye length and k is the wave number of the EISCAT UHF-radar (λ=0.322 m −1 ).
It can be seen that if the Debye term can be omitted and the electron and ion temperatures are equal, then the raw and true electron densities are also equal. In order to make a correction for the raw electron density, we have estimated the Debye length by using the raw electron density as N e and by assuming T e =T i =T n , which is usually a good approximation in the D-and lower E-region. The T n profiles were obtained from the MSIS model (Hedin, 1991). Then the power profile data, together with Eq. (2), was used to estimate N e between 70 and 86 km. The alternating code provided N e between 86 and 140 km.
To calculate the absorption, the electron-neutral collision frequency profile was also needed and it was calculated by ν en = ν eO + ν eO 2 + ν eN 2 (3) (Doe et al., 1995), where the individual collision frequencies for atomic and molecular oxygen, and molecular nitrogen are given by where the number densities of |O], [O 2 ] and [N 2 ] are taken from the MSIS model. The electron temperature for the height region 100-140 km was taken from the EISCAT measurement and below this from MSIS, assuming T e =T n . The electron gyrofrequency was calculated by using the IGRF magnetic field model for Tromsø. To estimate how sensitive the results were to the electron-neutral collision frequency profile selected, we also calculated the absorption by using an equation ν en =5.8·10 −11 NT 5/6 e from Gurevich (1978), where N is the neutral density from the MSIS model, but the difference was negligible.

Studied events
3.1 Event 1, 7 December 1999 The first event occurred between 19:40 and 20:05 UT on 7 December 1999 during a growth phase of a small substorm. Selected images of ASC and IRIS data from this event are shown in Figs. 1 and 2. Figure 1 shows the 557.7 nm emis-sion from the KIL ASC in geographic coordinates by assuming an altitude of 100 km as the lower border of optical auroras. The altitude was selected on the basis of EISCAT electron density profiles. The location of EISCAT is marked with a magenta dot. The optical data shows a wide region of diffuse aurora with a weak discrete arc at the poleward edge, which becomes visible at 19:48 UT. At the poleward horizon, more intense aurora exist. The weak arc and the diffuse aurora drift slowly equatorward. Also, scanning photometer data at 557.7 nm from Karesuvanto (68.47 • N, 22.44 • E) show this equatorward drift clearly (data not shown). New arcs form on the poleward side of the studied arc, as seen at 20:00 UT. The onset of the substorm occurs at about 20:08 UT.
The top panel of Fig. 2 shows the CNA measured by IRIS in a geographic grid, assuming an altitude of 80 km for absorption. The data shows a longitudinally stretched absorption region, which drifts equatorward. The total width of this region is all the time below 100 km, and the half-width is of the order of 50 km. The absorption region contains structures changing with time, in agreement with the observations of pulsating absorption patches by Ranta et al. (1999). Furthermore, absorption shows intensity variations and it is the most intense with a maximum value of 1 dB between 19:55 and 19:59 UT, when the maximum of the CNA region is no longer within the EISCAT beam (red circle in the figure).
The middle panel of Fig. 2 shows the ASC data in geographic coordinates from the same time intervals as the IRIS data in top panel. Again, the structures are mapped to the altitude of 100 km. All auroral structures drift equatorward and the velocity of the studied arc is about 170 m s −1 . The 20:02 UT image shows the intensification of the arcs to the east of EISCAT and IRIS.
The bottom panel of Fig. 2 combines the IRIS and optical images from the top and middle panels. Here the location of the arc at the poleward boundary of the diffuse aurora is shown by a red line. Since the mapping height of the arc was 100 km and the IRIS measurement 80 km, an accurate comparison would require mapping of the arc location 20 km further down along magnetic field lines. Since the inclination of the magnetic field lines is only 12 • at KIL, this correction would account for a shift of the arc of only about 4 km to the poleward direction, which would not be visible in the figure, so that correction is omitted. The figure clearly demonstrates that the absorption region is located on the equatorward side of the weak arc, within the same region as the weak diffuse aurora. Figure 3 shows the EISCAT measurement and the measured and calculated CNA. The top panel shows the electron temperature from 100 to 140 km altitude. In the calculation the values of T e were used from the MSIS model for the height range of 70-100 km. The electron temperature remains rather small throughout the event, except for 19:59-20:00 UT, on the poleward side of the studied arc at an altitude of about 110 km. We argue that the T e enhancement is not related to the studied arc but to the arc on the poleward side of it, since the analysis of the electric field shows that at this location the northward electric field reaches a value of 50 mV m −1 and thus that arc obeys the typical electrodynamic pattern of some eveningside arcs (e.g. Aikio et al., 2002).
The second panel shows the raw electron density between 70 and 86 km, and the electron density from the alternating code between 86 and 140 km, as used in the calculation of the CNA. As explained in the previous section, the raw electron density was corrected for the Debye effect, which accounted for a correction factor of the order of 1.1 for the density. The diffuse aurora is visible as the electron density enhancement mainly between the 110 and 130 km altitude. The studied arc is measured between 19:54-19:57 UT by EISCAT, and it is visible both as an electron density enhancement and as a decrease of the lower edge down to 100 km. After 20:00 UT there are small electron density enhancements within a broad range of altitudes above 110 km, which are related to other auroral arcs drifting to the EISCAT beam. Between 19:45 and 19:55 UT, on the equatorward side of the studied arc, there are enhancements of electron density at lower altitudes, with a maximum at about 83 km. They are obviously produced by energetic electron precipitation.
The third panel shows the total intensity of the electric field measured by EISCAT. The electric field is enhanced from 19:58 UT onward in the region of decreased electron density in the E-region poleward of the arc.
The bottom panel shows the calculated total absorption from the EISCAT data (red curve) and that measured by IRIS (blue curve) at the location of the EISCAT beam. The absorption from EISCAT shows two peaks corresponding to the two strong enhancements of electron densities between 80 and 90 km. The IRIS data shows the same two peaks and an additional middle peak at about 19:49 UT. These peaks reflect the temporal and spatial variation in absorption. When analysing EISCAT data with a 20-s integration time, D-region electron density values showed significant variations within one minute time, e.g. 19:49:00-19:49:20 UT, where a clear maximum was obtained. Obviously, the D-region precipitation varied in time and space and this, together with the narrow beam of EISCAT (0.6 • ) compared to the wide beam of IRIS (13 • ), resulted in the absence of a clear middle peak from the calculated CNA. During 19:50-19:53 UT, CNA was small throughout the IRIS fov and after 19:56 UT absorption reached its maximum value, but since the region of CNA had drifted away from the EISCAT beam, it is not visible in Fig. 3. It can be further noticed that during the time interval of high E-region Te, 19:59-20:00 UT, the CNA produced between 100 and 120 km is only 0.006 dB. The two curves show clearly the same characteristics with similar background absorption values. The measured CNA values are larger than the calculated ones, but the difference is rather small, with a mean value of 0.1 dB (average over 19:44-20:00 UT). The possible reasons for this difference are discussed in Sect. 4. Figure 4 shows altitude profiles of electron density, electron temperature, electron-neutral collision frequency and incremental absorption (absorption per unit height), calculated from the EISCAT data at the first peak of CNA in Fig. 3, 19:46-19:47 UT. Absorption is produced below 90 km with a maximum close to 81 km, close to the maximum of the Dregion electron density, which is located at 87 km altitude.

Event 2, 8 April 1998
The second event occurred between 19:45 and 20:20 UT on 8 April 1998, also during a growth phase of a small substorm. Figure 5 shows IRIS and optical data in the same format as Fig. 2 but now the auroral arc has been mapped to 120 km, as deduced by the EISCAT data. The top panel shows again a region of east-west elongated and patchy CNA drifting equatorward. The absorption area is slightly wider than in Event 1, the total width exceeds 120 km and the half-width is about 60 km. In this case the CNA is very weak, with a maximum value of about 0.3 dB. The all-sky camera was switched on at 20:11 UT, just prior to a faint discrete arc which drifted to the beam of EISCAT. The arc continued to drift equatorward until the substorm onset at 20:20 UT, when the arc intensified and a fold developed in the arc. The average drift speed of the arc was 120 m s −1 , and the width of the arc was 15 km (Aikio et al., 2002). Again, the region of CNA was located on the equatorward side of the arc, which can be seen from the bottom panel. Figure 6 shows EISCAT and IRIS data in the same format as Fig. 3. Electron temperature (top panel) shows a significant enhancement centered at 110 km just on the equatorward side of the arc, with values reaching 1700 K during 20:10:45-20:11:45 UT. Within the same region, the mainly northward electric field magnitude reaches 51 mV/m, as shown in the third panel. The discrete auroral arc is visible as an increased electron density above 120 km during 20:11:45-20:12:45 UT (second panel), after which there is a short data gap. On the equatorward side of the arc there is some electron precipitation to the E-as well as the D-region. The D-region electron density has a maximum value at about 90 km, between 19:51 and 20:11 UT. The bottom panel of Fig. 6 shows the values of total absorption from IRIS (blue line) and EISCAT (red line). The values of CNA measured with IRIS were small, with a maximum value of about 0.3 dB. At the location of the EISCAT beam the maximum was 0.25 dB. The enhanced absorption is co-located with the enhanced electron density in the D-region from 19:51 to 20:14 UT. During this time interval, the measured and calculated absorption are almost equal. However, the calculated background values are about 0.1 dB larger than the measured ones.
Within the narrow region of the enhanced electric field and electron temperature on the equatorward side of the arc, there is no increase in the measured CNA. Also, the calculated CNA is at the same level as inside the arc. Figure 7 shows the altitude profiles from this time interval. The T e enhancement centered at the altitude of 110 km produces a small enhancement to the incremental absorption profile. The enhancement amounts to about 0.01 dB between altitudes of 100 and 120 km, so its effect to the total absorption is small, and the main contribution comes below the altitude of 90 km.

Discussion
We have studied two stable drifting arcs, which were associated with a region of CNA on the equatorward side of them. The CNA close to magnetic midnight hours is generally attributed to the enhancement in the D-region electron densities, which, together with the high electron-neutral collision frequency in the D-region, is capable of absorbing a part of the cosmic radio noise (Collis et al., 1996;Stauning, 1996).
Previous studies have also shown that drifting auroral arcs during a substorm growth phase are often associated with an enhanced electric field and enhanced electron temperature in the E-region on the equatorward side of them, indicative of a plasma instability Aikio et al., 1993). The frequency range generally in use by riometers corresponds to the few-meter-length Eregion irregularities produced by the Farley-Buneman (FB) instability. These waves are known to be generated in the electrojet layer during periods of strong electric fields (Farley, 1963;Bunemann, 1963). Initially, the mechanism has been proposed to be operated in the equatorial electrojet (D'Angelo, 1976) and then the conclusion has been extended to the auroral electrojet (D'Angelo, 1978;D'Angelo, 1981;Mehta and D'Angelo, 1980). However, Farley et al. (1981) pointed out that the irregularities are highly field-aligned, and since the magnetic field is nearly vertical in the auroral zone, even very strong irregularities will be unable to scatter back into space cosmic noise coming from more or less overhead. Analysis of riometer observations by Siren (1982) supported the idea that scattering of cosmic radio waves is unable to quantitatively account for the observed absorption intensities at high latitudes.   6. The EISCAT data for the Event 2, with an integration time of 60 s for T e and N e , and 120 s for the electric field. The data format is the same as in Fig. 3. The panels from top to bottom are electron temperature, electron density, total intensity of electric field and absorption, both from IRIS and EISCAT. The E-region absorption has also been attributed to the enhanced electron temperature directly. Stauning (1984) explained the absorption event with values below 1.0 dB in the summer cusp by an enhanced classical collision frequency in the strongly heated E-region, in association with strong ionospheric electric fields exceeding 50 mV/m. This value of the electric field magnitude corresponds quite well to the experimentally determined "threshold" condition which should be met in order for strong electron heating to start operating Haldoupis et al., 1993;Aikio et al., 2002). It was also pointed out by Stauning (1984) that an event of this type requires a relatively dense E-region and substantially elevated electron temperatures to cause visible effects in the measured absorption. Due to this, daytime profiles were selected in his work from the polar cusp, and during the selected events, the density at 105 km height stayed close to 2·10 11 m −3 . During similar conditions, Olesen et al. (1986) have shown that the frequency selective absorption (lacuna) at the E-region critical frequency might well be explained by the enhanced deviative absorption caused by the increased electron-neutral collision frequency in the heated E-region. High electron temperature, therefore, seems a possible explanation for lacuna and associated cosmic noise absorption events.
The analysis of EISCAT data for these two events showed that within the regions of CNA, enhanced electron densities in the D-region existed. However, on the equatorward side of an arc an elevated electron temperature in the E-region and enhanced electric field was found in Event 2. The nondeviative absorption was calculated between 70 and 140 km, and it was compared to the measured absorption at the point of EISCAT observation. In Event 1, CNA varied with time, and this temporal variation was similar, both in the IRIS data and the calculated CNA. However, the calculated CNA was on the average 0.1 dB smaller than the measured one in Event 1, and there were some differences in the behaviour. There are several possible reasons for these differences. First, CNA was spatially structured, and this may have some effect, since the beam width of EISCAT is only 0.6 • as compared to the IRIS beam of about 13 • . Second, the EISCAT beam is not located exactly at any single IRIS beam position but the IRIS measurements have been interpolated to the location of EIS-CAT. Third, there are uncertainties in the calculated CNA. We used the MSIS model to estimate the electron temperature for the lowest altitudes, and they were used in calculating the electron-neutral collision frequencies and the Debyecorrection for the electron density. The number densities of neutral molecules are from the MSIS model and thus, the υ en is probably the largest source of error in the calculation. Also, the measured parameters have certain confidence limits. Despite these uncertainties, the calculated and measured CNA are in good agreement in both of the events.
The maximum absorption was produced at an altitude of about 80 km in both events. The maximum altitude of the D-region electron density enhancement took place at 84 km (Event 1) and at 90 km (Event 2). Both the altitude of the density enhancement and the peak absorption remained rather constant throughout the events. The electrons causing density enhancement in D-region may originate from the radiation belts but the discussion of the precipitation mechanisms is beyond the scope of this paper and we refer to the paper by Sergeev et al. (1990).
In Event 2, a region of elevated T e , centered at 110 km, was located on the equatorward side of the arc, with the value reaching 1700 K, about 1300 K above the background level. The width of the region of enhanced T e was of the order of the arc width, 15 km. The altitude range of the temperature enhancement was about 10 km (height resolution was 3 km), so it was very narrow. The incremental absorption profile showed a clear enhancement between 100 and 120 km, completely accounting for the absorption of the order of 0.01 dB. However, the main part of the absorption comes from altitudes below 90 km. The small effect of the enhanced T e on the absorption is partly due to the very low N e within the same region: the enhanced electric field and T e tend to occur in the auroral cavity (Aikio et al., 2002).
In Event 2, the enhancement of T e , with a value of 1300 K over the background in the narrow region just equatorward of the arc, was well correlated with the electric field value exceeding 50 mV/m, in accordance with FB instability. We can argue that during this time the electric field is responsible for the T e enhancement in relatively low density plasma. However, in this narrow region no measurable enhancement in CNA due to plasma instability can be detected. In Event 1, the magnitude of the electric field remains small (≤20 mV/m) within the region of CNA. Poleward of this region the arc electric field intensifies to a level of about 50 mV/m and somewhat above. Also, the T e around 110 km is increased within this region, and the maximum (500 K above background) is obtained at a point where the electric field is about 50 mV/m but the variance is very large. No significant increase in CNA is observed within this region.
If the wave scattering mechanism would have played a significant role in these events, the calculated absorption would have been clearly smaller than the measured one. In that case it had been necessary to replace the classical collision frequency with an effective collision frequency caused by the collective interaction of electrons with waves.

Conclusions
We have presented simultaneous data from optical auroral, cosmic radio noise absorption and plasma parameter measurements for two events, 7 December 1999 (Event 1) and 8 April 1998 (Event 2). Our aim was to study whether the conventional D-region or the alternative E-region processes produce the measured absorption bay during the substorm growth phase. The proposed mechanisms to cause absorption are (1) the enhancement of D-region electron density by precipitated energetic electron or protons, (2) the wave scattering of cosmic radio noise due to instabilities in the Eregion and/or (3) the effect of the highly enhanced electron temperature in the E-region.
Both of the events studied occurred during the growth phase of an auroral substorm in the evening MLT sector, when a quiet auroral arc and a region of CNA in the vicinity of the arc drifted equatorward. It was shown that the region of the most intense absorption in the absorption bay was located on the equatorward side and outside of the optical auroral arc. Both of the events studied were associated with enhanced N e in the D-region, with a maximum between the altitudes of 80 to 90 km. Event 2 was also associated with an intense electric field and an enhanced electron temperature centered at an altitude of 110 km, on the equatorward side of the arc. The analyzed data show clearly that the effect of both the electric field exceeding 50 mV/m, and the enhanced T e at 110 km altitude with values up to 1700 K, for the total absorption, is small.
The measured value of the electric field exceeded the threshold value of Farley-Bunemann instability. However, IRIS measured no increase in CNA within the latitudinally narrow region of the high electric field and increased T e . This is in accordance with the theoretical predictions that field-aligned irregularities are unable to scatter back a significant amount of cosmic noise (Farley et al., 1981). Our calculations also showed that even though the value of T e was tripled, it didn't produce an observable effect on the total absorption by using the equation of nondeviative absorption. Due to the location of the T e enhancement in the E-region, the effect on the electron-neutral collision frequency remains small. In addition, the enhancement is located in an auroral density cavity, which further decreases its effect. The situation may be different in the dayside, as indicated by observations of Stauning (1984) and Olesen et al. (1986).
As a conclusion, we argue that the absorption observed in the vicinity of auroral arcs during a substorm growth phase, i.e. the absorption bay can entirely be explained by precipitation of electrons with energies over 30 keV from the plasma sheet, increases the electron density in the D-region.
all-sky cameras. The Ph. D. work of the first author, J. R. T. Jussila, is funded by Vilho, Yrjö and Kalle Väisälä Foundation of Finnish Academy of Science.
Topical Editor M. Lester thanks T. Kikuchi and another referee for their help in evaluating this paper.