LWPC ANALYSIS OF LIGHTNINGINDUCED SFERICSâ€™ ELF PROPAGATION
VELOCITY
By
SANDEEP RANJAN SARKER
BSEE, University of Alabama, 2015
A thesis submitted to the Faculty of the Graduate School of the University of Colorado in partial fulfillment of the requirements for the degree of Master of Science Electrical Engineering Program
2017
ii
This thesis for the Master of Science degree by Sandeep Ranjan Sarker has been approved for the Electrical Engineering Program by
Stephen Gedney, Chair Marek Golkowski Miloje Radenkovic
Sarker, Sandeep Ranjan (MS, Electrical Engineering)
LWPC Analysis of Lightninginduced Sfericsâ€™ ELF Propagation Velocity Thesis directed by Assistant Professor Marek Golkowski
ABSTRACT
The Long Wave Propagation Capability (LWPC) software package is a comprehensive simulation tool with a longstanding track record of simulating the propagation of very low frequency (VLF: 33 kHz) electromagnetic waves in the Earthionosphere waveguide. The code is able to simulate waves launched at any geographic location and takes into account effects of parameters that vary spatially around the world, such as the Earthâ€™s magnetic field, the ground conductivity, and the zenith angle of the sun.
In this work, we apply LWPC to the extremely low frequency (ELF : 0 3kHz) band and focus on the group velocity of electromagnetic impulses generated by lightning discharges, known as â€œsfericsâ€. The group velocity shows variation as a function of the electron density profile of the lowest region of the ionosphere known as the Dregion. The Dregionâ€™s electron density profile can generally be fitted as an exponential function quantified by two key parameters, 1) the scaling factor for the reflection height of the Earthionosphere waveguide 2) the sharpness of the exponential profile. We analyze the change of the group velocity for profiles with reflection heights ranging from 70  90 km as well as sharpness magnitudes ranging from 0.4  0.6 km1. These parameters correspond to day and night conditions. The variation of the propagation velocity in an Earthionosphere waveguide containing a nonexponential electron density profile is also analyzed.
IV
Other factors that are analyzed in this work are the path lengths of the propagation, the conductivities of the ground, different directions of propagation, geographic locations of the transmitter and receiver, propagation across the daynight terminator, as well as the trend of propagation velocity change under the event of a sudden ionospheric disturbance (SID). A SID occurs when the sun emits a short wavelength Xray carrying a highlevel of energy towards the Earth. These Xrays carry enough energy to penetrate into the Dregion and induce photoionization. This photoionization process lowers the reflection height and increases the sharpness of the profile. Two models for predicting the change of the electron density profile under a SID event have been developed. The first developed model simply changes linearly with the logarithm of the strength of the Xray flux. This model is relatively accurate in predicting the change of the amplitude and phase, but does not capture the relaxation of the ionosphere after the SID event took place. The relaxation time depends on the number of electrons that were freed during the event and the duration of time needed for these electrons to reattach to the ions in the Dregion. The second model captures the relaxation of the ionospheric profile after the event. The results show that ELF propagation velocity increases during a SID event.
The form and content of this abstract are approved. I recommend its publication.
Approved: Marek Golkowski
v
To my parents, Krishna and Sanjoy, for putting up with all of my tantmms, giving me guidance, and always believing in me. To my brother, Anindo, for all of his guidance throughout my life. To my cousin, Amrit Bandyopadhai, for inspiring me to become an electrical engineer.
VI
TABLE OF CONTENTS
I. Introduction..............................................................1
Lightning Source Waveform and Radio Atmospherics...........................2
EarthIonosphere Waveguide.................................................4
Ionospheric Profdes..............................................4
Mode Theory......................................................6
Long Wave Propagation Capability Software (LWPC)..........................10
PRESEG..........................................................10
MODEFNDR........................................................11
FASTMC..........................................................12
II. Modeling of Lightning Induced Sferics Using LWPC............................13
Lightninginduced Sferics.................................................13
Modeling the Lightning Source in LWPC.....................................15
ELF Propagation Theory....................................................16
Obtaining the ELF Propagation Velocity from LWPC..........................23
III. LWPC Results.................................................................25
Initial Results...........................................................2
EW vs NS.................................................................22
Land vs Sea..............................................................27
IV. LWPC vs Experimental Data.....................................................29
Variations of ELF Propagation Velocity Over DayNight Terminator..........29
LWPC vs Greifinger Theory.................................................31
Diagnosing A Sudden Ionospheric Disturbance Using ELF Propagation Velocity Variation..35
V. Conclusion...................................................................39
Bibliography........................................................................40
vii
viii
I.
Introduction
Lightning from thunderstorms involves a powerful discharge current over a long distance and therefore generates wideband electromagnetic radiation. The electromagnetic radiation from lightning is known as radio atmospherics often shortened to â€˜sfericsâ€™. These signals have a broad frequency range, from a few Hertz to MHz. Because of the timesignatures of the lightning current source as well as the attenuation of propagation, most of the energy radiated that can observed at long distances is within the extremely low frequency (ELF) and very low frequency (VLF) bands. Here we take the EFF frequency band to range from 03kHz and the VFF band to range from 330 kHz. Over long distances, the sferic energy in the EFF/VFF bands is guided along the surface of the Earth and propagates as if the signals were traveling in a very large parallel plate waveguide. The bottom of the waveguide is the ground. The â€œroofâ€™ of the waveguide is the upper part of the atmosphere, what is specifically referred to as the ionosphere. This configuration is also known as the Earthionosphere waveguide. As discussed below, at ELF/VLF frequencies both the ground and the ionosphere act as conductors.
In addition to neutral atmospheric particles, the ionosphere contains ions and electrons with an overall net charge that is equal to zero. Such a collection of ions and electrons is also known as plasma, which is regarded as the fourth state of matter. Plasma is created when a neutral gas molecule exceeds a threshold energy. The ionosphere is dynamic and is constantly varying based on day vs. night, solar cycles, seasonal cycles, and numerous other space weather effects. The varying distribution of the electron density in the ionosphere causes the ELF/VLF reflection height of the Earthionosphere waveguide to be constantly changing. The changing reflection height affects the propagation of ELF/VLF waves including sferics.
1
Lightning Source Waveform and Radio Atmospherics
Lightning is one of the most frequently occurring natural processes in the Earth. During any moment, there are about 2000 thunderstorms taking place. The amount of the Earth, which is covered by thunderstorms, is about ten percent of the Earthâ€™s surface. The average number of lightning strikes per a second globally is about 45. Certain geographical areas experience thunderstorms at higher rate than other areas. Typical regions that experience a heavy number of thunderstorms are the Carribean, southeast Asia, and central Africa [Uman, 3739]. Most of this activity occurs during the summer time. This work will mainly cover the effects of lightning strikes generated in the east coast of the United States.
Several different lightning generation processes occur in nature, many of which are still poorly understood. For the purpose of this work the most important process to understand is the cloudtoground lightning return stroke waveform. By creating a vertical current path, the charge moving as a current from cloud to ground acts as the source for all of the ELF and VLF sferics that will be covered in this work. At close to 8ps in rise time, the current waveform of a typical cloud to ground discharge suddenly increases in a step sequence, and then decays exponentially with a â€œslowtail.â€ The overall current relaxes back to the initial value around 500ps.
2
Figure 1.1 Lightning returnstroke
We briefly comment on three major classifications of transient luminous events, which occur after thunderstorms, and are known as sprites, elves, and jets. Sprites are large flashes of light above a thunderstorm, which occur immediately after a cloudtoground (CG) lightning strike. The altitude region at which a sprite occurs is generally around 5090 km. Elves are expanding light rings that occur in the Dregion of the ionosphere. The general altitude range is 8090 km. These rings are centered above a CG lightning stroke, and are emitting red photons at intensities up to tens of megarayleighs. Rather than being a direct consequence of lightning like sprites and elves, jets are discharges occurring due to the quasielectrostatic field created by the cloud in the ionosphere. The two types of jets are blue and gigantic. Blue jets are relatively frequent discharges ranging from the cloud top to 40 km altitude region. Gigantic jets are much more infrequent occurrences. The altitude range extends from 80 km to the cloud top. The rings exhibit streamerlike properties [Blaes, 810]. Although transient luminous events and their associated electromagnetic effects have received considerable attention in recent years, in this work we focus on the propagation of sferics in the Earthionosphere waveguide.
3
The EarthIonosphere Waveguide
Ionospheric Profiles
The ionosphere is defined as the upper layer of the atmosphere where various ionizing processes occur for there to be a significant number of free electrons and ions. Due to the amount of free electrons and ions, the ionosphere acts as a good conductor to reflect ELF and VLF waves toward the ground in the Earthionosphere waveguide configuration. The ionosphere is generally divided into three (sometimes four) regions: D, E, and F. The Dregion is the lowest altitude region of the ionosphere. The altitude range above the ground of the Earth that the Dregion encompasses is about 50150 kilometers. The configuration of the number of electrons in the Dregion greatly affects the reflection height for the ELF/VLF frequency band in the Earthionosphere waveguide. For this reason, Dregion variations will be solely analyzed.
The electron density profile can be characterized by a two parameter exponential formula, which gives a general exponential fit of the electron density vs altitude in the Dregion of the ionosphere:
JVe(/i) = iV0exp(â€”.15/T)exp[(/? â€” .15)(/i â€” h')] . (1.1)
In this expression, h is the altitude, lr is the reference height of saturation, and /? is the magnitude of saturation for the ionospheric profile, it is also referred to as the sharpness. In realtime, h â€™ and /? are constantly varying due to daynight variations, seasonal cycles, solar cycles, and numerous other space weather effects [Schmitter, 2013]. The electron density profile of the Dregion can almost always be estimated with this exponential model. General daytime values of the reference height are between 6875 km, and nighttime values are around 7890 km.
4
General daytime values of the sharpness are around 0.350.45 inverse kilometers, and nighttime values are around 0.50.6 inverse kilometers.
Ionospheric profiles(daytime vs nighttime)
Figure 1.2 Daytime and Nighttime Electron Density Profiles of the Dregion
Figure 1.2 displays two different electron density profiles versus altitude. Altitude is the yaxis in kilometers, and electron density is the xaxis. The blue line is for the daytime profile of h â€™ = 70 and /? = .4. The orange line is for the nighttime profile of h â€™ = 84 and /? = .6. This work will also investigate the variations of the ionosphere under the occurrence of a sudden ionospheric disturbance. This phenomenon occurs when the sun emits a solar flare towards the Earth adding energy to the ionosphere. This solar flare is composed of Xrays, which add to the density of electrons and the ions in the ionosphere. Under the occurrence of a sudden ionospheric disturbance, the sharpness increases and the reference height decreases.
5
Mode Theory
The widely accepted theory of the Earthionosphere waveguide configuration begins with a parallelplate waveguide (rather than a spherical waveguide), and then introduces corrections for the spherical geometry of the configuration. James Waitâ€™s formulation of the Earthionosphere waveguide provides an accurate model for radiation in the Earthionosphere configuration (Wait, 1970). In this work, he shows that the fields can be reproduced in terms of a sum of modes. Rather than using impedance matching conditions to solve for the eigenvalues of each mode in the waveguide, Wait focuses on the boundary conditions of the permittivity, permeability, and conductivity profile of the ionosphere. By matching tangential magnetic field and normal electric field conditions, Wait formulates a Legendre polynomial series solution for the fields in the ground, freespace, and ionosphere.
The source of the radiation across this waveguide is modeled as a vertical dipole source of strength II with a frequency of a) oriented in the zdirection. K.G. Budden shows that this source can be equivalently modeled as the sum of an infinite number of line quadrapole sources parallel to the yaxis [Budden, 1962], The Hertz vector U, which is much like the magnetic vector potential A, is derived as follows:
kll
Uz = â€” jc exp[â€”jk(xcos9 + zsin9]cos9d9 (1.2)
In this formulation, 0 is the angle in the xz axis. A general model for a freespace region in the model of two reflecting layers shows that these reflecting layers can be stratified since they are only defined by their reflection coefficients. In the freespace region, the operator ^/'q is
equivalent to zero for all cases, so the individual modes can be classified as either transverse electric (TE) or transverse magnetic (TM) propagation [Cummer, 1997]. The ionosphere is a
6
magnetized plasma. This means the fields are coupled at the roof of the waveguide configuration, therefore the incident wave produces both TE and TM components so that there are no purely transverse magnetic or transverse electric fields propagating through the Earthionosphere waveguide. This coupling means that the reflection coefficients of the upper boundary are not a single scalar quantity, but rather a 2x2 matrix of the different types of polarizations. The ground is isotropic in terms permittivity and conductivity values and can also be treated as a 2x2 diagonal matrix [Budden, 1962], The reflection coefficients of the ionosphere is shown be:
II Rtf') 1R(.6)
\R(9) iR(9)
(1.3)
, =rF(0)
L 0 {R'(9)
The top and bottom subscripts of the matrixâ€™s indices denote the incident and reflected polarization of the wave, respectively. The matrix Rv denotes the downward reflection coefficients of the ionosphere. The matrix, RL, denotes the upward reflected wave from the lower region of the waveguide. A modified version for the Hertz vector potential is derived below:
Uz = â€” fc exp[â€”jk(xcos9 + zsin9] â– [1 0](/ + RV)W(I + RL)[^]cos9d9 (1.4)
7
The multiple reflections occurring in a single path of propagation are modeled using a series of
equivalent image sources [Cummer, 1997].
Figure 1.3 Imagesource equivalent of propagation through the Earthionosphere waveguide
W is equal to (7  RVRL)"1. W is singular or the inverse is equal to zero at 0n, which is where the poles of the source are. This means:
\I  Ru(en)RL(0n)\ = o (1.5)
Using this relation, one can solve for the eigen angles of each mode of a wave propagating through the Earthionosphere waveguide [Kraus 1992, pg 691]. Using a series of asymptotic expansions, the transverse horizontal magnetic field at z=0 is expressed as:
By(x) = Jk^U exp^/q) AnA/singnexp()/cxsingn) (1.6)
8
This asymptotic assumption comes from redefining 0n as the angle of incidence relative to the roof of the waveguide (cos 0 sin 0). The mode excitation factor is given by
An = [l
U+Ru(en)]x[i+RL(en)] rin
UJ 3A/ , LnJ
= 9n U
dAj
/<
(1.7)
An is defined as an excitation factor for the corresponding mode. A(0) is expressed as det(PV_1). X is expressed as WA(0) 0=0n. Using multiple matrix manipulations, and exploiting the fact that A = 0 at 0n, the excitation factor of the mode can be rewritten as:
An â€”
(l+R02(ljR'iR) II d dA  llR,S sin9'9=9n
(1.8)
The physical properties of the magnetic field can be described in more detail. The constant is frequency dependent, and it also depends on the currentmoment strength II of the dipole source. The excitation factor, An, is widely known as the excitation factor of the mode. This quantity describes the efficiency of excitation from the vertical dipole source [Cummer, 1997]. The summations described above are over an infinite number of modes, but in practice are used for only the modes with significant contribution to the fields.
The x1/2 dependence of By(x) in equation (1.6) describes the attenuation due to the spreading
of energy over the entire waveguide. In order to make spherical corrections, this term must
become [REsin(â€”)]1/2, where RE is the radius of the Earth. The corresponding factor that the Re I
I ~lR
field must be multiplied by for this specific spherical correction is Jjn(x/ER y If the propagation
distance, x, is much smaller than the radius, this factor is close to unity. It is only for large propagation distances, where this impact is tmly quantifiable. Condition (1.5) must also be
9
modified for the geometry of the Earth. The modified refractive index nmod becomes
7l + 2 z/Re.
This derivation accounts for a vertical dipole source on the ground oriented in the zdirection. In order to account for various sources of an extended altitude range with different orientations, new parameters must be introduced. Papert and Ferguson [1986] introduce new functions known as heightgain factors. In the aforementioned work, these functions fex(z), fez(z), and fty(z) are normalized. These functions are the factors for their respective field components, which is denoted by the subscript. However, in this work, these functions will stay unnormalized for the ease of computation. Once this correction is included, a corrected expression for the output field at the receiver is derived.
F = C(F)i7^exP(;T)SrlAtnArnexp()/cxsin(0n) (1.9)
The software package MODEFNDR, which solves the for the eigenangles of a various modes of any single propagation path, does not list heights in reference to the ground as this derivation does. Instead, the software choses an arbitrary reference height at which the modal solutions are much easier to compute. With this change, the excitation factor from the transmitter can be represented as:
Atn = â€”Asin(0n) cos(y) fby(zt) + flsin(y) cos(0)/ez(zt) + Asin(y) sin(0) fex(zt) (1.10)
The coefficients A, B, and the heightgain factors are all explicitly defined in Papert and Ferguson [1986]
Long Wave Propagation Capability Software (LWPC)
10
The Naval Oceans Systems Center successfully implemented this formulation in a simulation software package known as the Long Wave Propagation Capability code (LWPC). This simulation software is comprised of three separate software routines, which are written in FORTRAN: PRESEG, MODEFNDR, and FASTMC [Cummer, 1997]
PRESEG
The program PRESEG.f is the first phase of the LWPC. The modified LWPC, which was adapted by Robert Moore and Neal Dupree at the University of Florida to simulate ELF paths rather than the original navy version that simulates VLF paths, takes in the user input parameters and finds the necessary waveguide parameters as inputs for the MODEFNDR and the FASTMC. In this version of the LWPC, the user inputs are the transmitter location, the receiver location, the transmitterâ€™s power and frequency, the electron density profile for each segment of the path, and the collisional frequency profile for each segment of the path. From these inputs, PRESEG determines the conductivity, and permittivity along each section of the Earth. The program also determines the magnitude and direction of the Earthâ€™s magnetic field. If the user inputs an inhomogeneous ionosphere for many segments of the path, the program segments the waveguide into many sections, where the boundaries are determined by when the ground conductivity or the ambient field changes.
MODEFNDR
MODEFNDR is the most computational extensive program of the LWPC. This can also be considered the pillar of the LWPC propagation theory and computation. Using the inputs determined by PRESEG, MODEFNDR searches for the eigen angles that satisfy condition 1.5. The reflection coefficients required to solve this equation are computed based on the electron
11
density profile, the ion densities, and the collisional frequency profile for separate angles of incidence. Once these angles are determined, the separate modal solutions and overall output field can be computed based on the derivation. This program also computes the excitation factor of each mode, which is needed to calculate the overall magnitude of the field.
FASTMC
The FASTMC computes the modal solutions of the field, and then sums these parameters to calculate an overall resultant field. If the user inputs in homogeneities in the ionosphere, this program computes the propagation based on the discontinuities in the waveguide. The output of the FASTMC is field strength, which is in decibels for one microvolt per meter based on one kilowatt of source power, and phase, which is outputted in degrees, of the vertical electric field.
12
II.
Modeling of Lightning Induced Sferics Using LWPC
Since the closure of the US Navy Project ELF facility in Clam Lake, Wisconsin, there are no ELF transmitters in existence. The primary source of ELF waves is that induced by lightning strikes. The electromagnetic radiation from a lightning strike is called a â€˜sfericâ€™. Much of the focus of this work will be comparing the propagation velocities that were found experimentally with simulation models developed using the LWPC or other theoretical models. To properly model these sources, an understanding of the lightninginduced sferic waveform must be developed.
Lightninginduced Sferics
Atmospheric lightning can be modeled as a currentmoment dipole source. The current is defined as the amount of charge travelling from cloudtoground through a certain channel cross section per second. The current moment is defined as charge multiplied by the length of the channel. The process of lightning discharges is complex, with many separate processes occurring on various timescales. Because of this, lightning emits energy across multiple frequency bands. There is a significant amount of electromagnetic energy released in the ELF and VLF bands.
Although the Earthionosphere waveguide has many similarities to an ideal parallel plate waveguide, the curvature of the Earth and the finite and anisotropic conductivity of the ionosphere cause departures from the ideal model. For propagation of the Earthionosphere waveguide, there are no pure transverse electric or transverse magnetic modes to the xz plane,
13
rather the modes are referred to as quasitransverse electric or quasitransverse magnetic modes. In the ELF propagation, there is only a single mode contributing to the propagation. This mode is referred to as the quasitransverse electromagnetic (QTEM) mode. For this reason, the part of the lightningspectra that must be modeled for this work is the portion contributing to this most basic mode of propagation in the waveguide.
The current flowing from cloud to ground is impulsive in nature and starting at initial time t = 0. This percentage of the currentmoment frequency band that is analyzed in this work ranges from 01 kHz. The amount of charge flowing from cloudtoground that will be analyzed in this work is estimated to have the waveform of:
Igit) = Ig0(e~at  e~bt). (2.1)
Reasonable values are Ig0 = 20kA, = 2 x 104 s1, and b = 2 x 105 s1. Assuming that this negative cloudtoground lightning strike acts as a vertical dipole source on the ground oriented in the zaxis, the current moment can also be estimated.
The pulse of the returnstroke is assumed to be propagating up the lightning channel. The propagation velocity can be expressed as vUghtning(t) = v0jiightningexp(yt). By integrating this velocity function with respect to time, the length of the channel can be obtained. After this integration has taken place, the function for the dipole moment is:
i(t) â– l(t) = ^ [e~at  e~bt][1  e~Yt], (2.2)
The amount of charge that is transferred from cloudtoground can also be estimated by integration of the current waveform:
q = C Ig(t)dt = igo(i  Â£). (2.3)
14
What is relevant to this work is the amount of ELF energy radiated from one of these lightninginduced sferics [Cummer, 1997]. Generally, a lower propagation velocity up the lightning channel yields a larger portion of ELF energy radiated from the sferic. Modeling this source in LWPC is not a straightforward matter as the LWPC is not a timedomain solver simulation software. The LWPC code was designed to model the radiation from US Navy VLF transmitters. Hence the inputs are the source is the location, the power, and the frequency assuming the sourceâ€™s waveform is a single frequency sinusoidal signal. Details of how the sferic waveform is input into LWPC are provided below.
Modeling the Lightning Source in LWPC
As mentioned above, the LWPC is not a timedomain solver. It is mode theory solver. Since there is not a timedomain metric in the LWPC, to obtain the ELF propagation velocity, the analysis must be done in the frequency domain. The main inputs of the LWPC are the transmitterâ€™s power, frequency, source and receiver locations (in longitude and latitude), as well as the ionospheric profile along the path. The LWPC only assumes a single sinusoidal source, so any other type of source must be modeled as a sum or sequence of sinusoidal sources. In order to simulate lightning source in the timedomain as figure 1.1 shows, a Fourier transform is performed on the output field electric field from the lightning current. Taking the square of the field, a power profile versus frequency is obtained.
Using equation 2.1 and 2.2, an output field from any observation point can be derived. The expression for an electric field radiated from a lightning return stroke is:
E = li o
hesini; dl 2nR dt
(2.4)
15
where he is the height of the initial lightning charge. is the angle of the electric field in relation to the local zenith angle, and R is the distance from observation to source. In the case of obtaining the amount of charge flowing from cloudtoground that is radiating the electric field, both the expression for current and the dipolemoment from the lightning return stroke can be used in place of the term I. The farfield spectral density that is obtained from the current and dipolemoment, respectively, is:
SiO) â€” K0
co2(ab)2 (a>+a)(a>+by
(2.5)
c ( \ â€” is fvÂ°)2 (S2+4a>2)(aycobya>)2
W ~ AÂ°Vy/ (a2 +co2)(b2 +co2)(a2 +co2)(P2 +co2)'
(2.6)
In these expressions, the term K0 includes all of the constants from the electric field and current equations mentioned.
K0
(2.7)
a = a + y, (3 = b+ y,8 = a + b+ Y
The units of this spectral density is watts times meters inversesquared, a) is the radial frequency with the units of radians per second. The units of the intrinsic impedance, Z0, is ohms, and the units of the permeability of free space, /i0, is Henries per meter [Cotts et ah, 2011]. Through the spectral density, a power profile versus frequency can be obtain to input into the LWPC in order to model the propagation of lightning sferic.
ELF Propagation Theory
Carl and Phyllis Greifinger developed a formulation to determine ELF eigenvalues in the Earthionosphere waveguide configuration [Greifinger and Greifinger, 1978]. This theory used
16
the magnetic vector potential created by a lightning source to derive the fields propagating in the waveguide. The most relevant part of their formulation relating to this work is that the fields have a proportionality relation to the propagation distance is shown below.
E,B
qJIlSqP 1 P2
(2.8)
where p is the propagation distance. The phase constant, S0, is a complex quantity where the imaginary component multiplied by the wave number of free space is the attenuation of the wave propagating in the waveguide. The real component multiplied by the frequency is the effective wave number.
ke/f â€” kS0 (2.9)
In this equation, k is the wave number of free space. Using these two basic electromagnetic relations,
vp =
keff
(2.10)
_ l
va ~ dkeff/
'do
(2.11)
it is possible to derive the ELF Propagation velocity for a single frequency sinusoidallyvarying wave in space and time propagating in the Earthionosphere waveguide.
According to Greifingerâ€™s theory, the square of the complex phase constant is equal to
5
2
0
h^
. n7'
hÂ°~J^
(2.12)
17
From these relations, both the group and phase velocity of ELF waves in the waveguide configuration can be obtained. First, boundary conditions of the ionosphere for ELF waves must be analyzed to determine the quantities h0, h1> and Â£ where Â£ is the conductivity scale height, and is defined to be
< =
a
do/dz
(2.13)
In this case, a is the conductivity of the ionosphere. The conductivity profile analyzed in this work is the DC conductivity profile of the ionosphere,
W(z)q2 mv{z) â€˜
(2.14)
For this conductivity profile, q and m are the charge and mass of individual electrons in the ionosphere, respectively. The term, N(z), is the electron density profile of the ionosphere with respect to altitude where N is the number of electrons, z is the altitude, and v(z) is the collisional frequency profile of the ionosphere with respect to altitude. Equation 1.1 is the widelyaccepted 2 parameter electron density profile. The collisional frequency profile accepted in most literature is
u(z) = v0e 15z (2.15)
In this equation the initial collisional frequency, v0, is 1.86ell collisions per second. Inserting Eqaution 1.1 and 2.15 into 2.14, the conductivity profile of the ionosphere is found to be
a(z) =M!eP(zft') . (2.16)
mv0
Plugging 2.16 into 2.14, the conductivity scale height is obtained.
18
This is the single quantity in the phase constant where both the reference height and the sharpness do not vary. Since the sharpness is the direct inverse of the conductivity scale height, it can be inferred that varying the sharpness varies the ELF propagation velocity to a greater extent than varying the reference height. After determining the conductivity scale height of the ionosphere, it is now possible to determine both the reflection heights of the electric field and magnetic field, which are h0 and hÂ±, respectively. In order to determine the reflection height of the electric field, the conduction boundary condition must be understood. The conduction boundary condition is
where a0 is the permittivity of free space. Inserting 2.16 into this boundary condition, the reflection height of the electric field in a standard conductivity profile is obtained as
What is found is that the reflection height of the magnetic field is equivalent to the reference height of the ionosphere plus a certain spatial step. This spatial step is dependent of the frequency and the sharpness. To obtain the reflection height of the magnetic field, the diffusion boundary must be analyzed. The condition for the diffusion boundary condition is listed below.
^Oo) =
(2.18)
(2.19)
(2.20)
From plugging the obtained conductivity scale height into this boundary condition, the reflection height of the magnetic field is found.
(2.21)
hi
1
P
mv0p2
4N0q0q20)
Now that the reflection heights of the electric and magnetic fields have been obtained, along with the conductivity scale height, the phase constant can be determined. This work decomposes the phase constant into numerators and denominators in terms of both magnitude and phase as depicted in this equation below.
S0 = Â£ exp [/ i tan1 Q] (2.22a)
The terms A, B, d, and e are listed below.
A = (
(Iln + h,)( i,n( ) + _ (p)T + +
\P V N0q2 ) ) \ p \4N0q0q2coJ ) \2(3/ \2(3/ YB \4Nlq*q0) V
(2.22b)
B =
( LUJ!^E.) + hâ€™)2 + (JL)
V B \4N0il0q2a)J ) \2B)
1/2
(2.22c)
v2 B)yB V4 N2q*qJ J
(2.22d)
e = a In (SMf) + kâ€™) pin ( mvÂ°P ) + h')  (iV
\B V N0q2 ) ) \B \4N0q0q2coJ ) \2/?/
(2.22e)
Plugging in this expression to 2.9 and 2.10, the expression for the resultant phase velocity of an ELF wave travelling in the Earthionosphere waveguide is obtained. It should be noted that Greifinger and Greifinger estimated the proportionality of the ELF propagation velocity to the speed of light to be the inverse of the real component of the phase constant.
20
After obtaining the phase constant, in order to determine the group velocity of the ELF wave, a differential of the wave number with respect to frequency must be done as Equation 2.11 shows. The expression for the differential of the wave number with respect to frequency is
d^e// Â£o __ ()so
dco c dco â€˜
(2.23)
The differential of phase constant can be written as
dSo
dco
3A
dco
â– Bâ€”A
SB
dco
cos [. 5tan1 (J)]  f f sinlftan^f)]
dd
dco
ed
de
dco
'+(?)
(2.24a)
To determine the differential of the phase constant, the differential of the magnitude and phase coefficients with respect frequency must be computed.
â€” â€” r ~dco ~ ^
[Â© Â©Â©+',o (> Â©sÂ©+ Â© ] + Â© Â© &+ Â©> [Â© Â©Â©+'*') (> Â©Â©)+hâ€™)  Â©1 (Â© Â© (Â©)] <224b>
dB
dco
211/2
(Â©Â©Â©Â©') +Â©J
(2.24c)
(2.24d)
de _ 1 . ,/?2c2
dco ~ (32com^Aco22
(2.24e)
21
Inserting 2.23 and 2.24 into 2.11, an expression for the ELF group velocity is obtained. The two figures below show the frequency versus the group velocity for two separate profiles.
Figure 2,1a Daytime frequency vs group velocity (left), Figure 2,1b Nighttime frequency vs
group velocity (right)
Figure 2.1(a) shows the group velocity versus frequency for a daytime exponential profile where the reference height is equal to 70 kilometers, and the sharpness is equal to 0.4 inverse kilometers. Figure 2.1(b) is the same plot for a nighttime ionospheric profile where the reference height is equal to 84 kilometers, and the sharpness is equal to 0.6 inverse kilometers. The units for the frequency along the vaxis is Hertz. The propagation velocity on the yaxis is normalized by the speed of light. As shown below, the daytime profile yields a group velocity of ~ 70% of the speed of light at 1000 Hz. The nighttime profile yields a group velocity of 83% of the speed of light at 1000 Hz. The conclusions are that nighttime ionospheric profile yields a higher velocity than the daytime profile because the sharpness of the nighttime profile is higher, and both profiles in the ELF band yield lower velocities than the speed of light.
22
Obtaining ELF Propagation Velocity from the LWPC
In order to model the lightning sferic in the LWPC, each power to the corresponding frequency in the power profile derived in Section 2.2 must be entered into a simulation run in the LWPC. Generally, this work is composed of simulations from 301000 Hz in 10 Hz resolution or 30330 Hz in 3 Hz resolution. Once each simulation has been executed, an inverse fast Fourier transform is performed on the LWPCoutput amplitudes and phases of the receiver to the corresponding operating frequencies of the transmitter to retrieve the output waveform in the time domain. Once that is done, the maximum is found and the corresponding time is defined as the timedelay. Dividing the distance between transmitter and receiver by the timedelay gives the output propagation velocity. Figure 2.2 is the FFT of the source, the source electric field in the frequency domain obtained from the source, the output signals in the frequency domain, and output signals at the receivers in the timedomain computed from the power profile.
Time (ms)
Figure 2,2 a) Current waveform of sferic vs frequency (topleft), b) source Efield vs frequency (topright), c) output data in frequency domain for three different distances, (middle), d) output waveform at the receiver in timedomain (bottom)
23
24
III. LWPC Results
Initial Results
To develop an understanding of how the ELF propagation varies with ionospheric parameters, an initial set of simulations were run were over a wide range of electron density profiles. There were exactly 421 profiles run in this set of simulations for one path. The profiles were over the range of reference height being set from 70 to 90 kilometers in 1 kilometer resolution and the sharpness being set from 0.4 to 0.6 inverse kilometers in 0.01 resolution. The source coordinates are 22.1781Â°N, 297.2852Â°E, and the receiver coordinates correspond to the ELF receiver in Hugo, Colorado (103.4078Â°N, 38.8918Â°W). These plots are made using the same power profiles and process for retrieving the velocity discussed in Section 2.4.
90
88
86
84
82
78
76
74
72
70
relative velocity (v/c)
1
0.99
0.98
0.97
0.96
0.95
0.94
0.93
0.92
0.91
0.9
0.4 0.45 0.5 0.55 0.6
beta
Figure 3.1 Relative Velocity for 2parameter profiles inputted into the LWPC to Hugo receiver
Figure 3.1 shows the results in a 2 dimensional color plot where the vaxis is the sharpness, the yaxis is the reference height, and the intensity of the color is the propagation velocity. The color
25
scale on the right of the plot is the corresponding velocity value listed as a percentage of the speed of light. Analyzing the plot in Figure 3.1, there is a linear relation along the vaxis between the sharpness and the propagation velocity. It should be noted that the dark space on the bottomright corner represent ionospheric profiles that the LWPC did not yield a valid propagation solution. Along the yaxis, the propagation velocity shares more of an oscillatory relationship with the reference height. This is in good agreement with the formulations developed in section 2.3.
relative velocity (v/c)
.99 .98 57 .96 35 .94 53 .92 51 5
0.4 0.42 0.44 0.46 0.48 0.S 0.S2 0.54 0.56 0.58 0.6
beta
Figure 3.2 LWPCgenerated Relative Velocity for 2parameter profiles to Poland receiver
Figure 3.2 shows the same set of simulations as Figure 3.1 with the only difference being the receiver coordinates. The receiver is located in Hylaty, Poland at the coordinates 49.209Â°N, 22.553Â°E. For this path, the LWPC yielded very little difference between ionospheric profiles. Every transmitter in this chapter uses the same power profile; the frequency range of this power
26
profile is 301000 Hz in 10 Hz resolution. This means that the frequency range and resolution must be lowered in order to notice a difference for certain paths simulated through the LWPC
East/West vs North/South Propagation
While Greifinger and Greifingerâ€™s theory [Greifinger and Greifinger, 1978] is a very good analysis of the propagation in the Earthionosphere waveguide, the LWPC is a more complete simulation tool. All of the conductivity values of the ground, as well as the orientations of the geomagnetic field on the propagation path, are geographically mapped in the simulation software. Moreover, the ionospheric profile does not have to be constant along the path. In the following we analyze the effect of the orientation of the Earthâ€™s magnetic field by comparing propagation paths in the cardinal directions over land and sea. What varies in each of the propagation paths is the orientation of the geomagnetic field of the Earth relative to the
27
propagation path. This varies the conductivity tensor of the ionosphere.
WORLD Latitude and Longitude 7? nr ns* tar iiw w ir nr Â«â€¢ w GREENLAND 1 0N M  is cf ir *r i? w ;r v inv tar 135* isr w W* Mir ? Â«v*4*. ARCTIC 0 C E A l\ 3 \ \ l_ . W V m,
<*:33â€˜ Alaska IU.SJU Bap Arctic Circle ^ i 41 A ft A Hilliwu KtLAN > , â€™ I U Â£ g i ^ RUSSIAN FEDERATION
^ NORTH CANADA /*Â« iv UNITED KINGDOM THAN tenmark Laima POLAND BELARUS GERMANY . ***"**Â£Â£ WAmMâ€˜N MONGOLIA SERBIA UZBEKISTAN r*\ ITALY â€˜ KYRGYZSTAN gSg? NORTH \ GflEECC Tllfl*F' AFGHANISTAN nxciitr â€¢niM^A Sâ„¢A lftAN ru.MA SOUTH MPAN 1 4 C 11 U ' tomsia jyyjf, mAo CHINA kbma OCEAN .kr
P AC in C OCEAN 30Â° UNITED STATES NORTH OF AMERICA ATLANTIC OCT AN PORTUGAL y,A MOROCCO
Zyjh' Tropic of Can HÂ»wi Idandt I3> IDS) f â– *â€œ qiT&U. â€œHj " â€œ UBV* SAUDI U.A.E INDIA , â„¢ v â€¢Â» r pjwllco MAURITANIA MALI CHAD ARABIA OMAN aANGLACFSM GUATFMA' A MONhuflAS 1151 SFNftiAL i ^CrfR SKlAN vsutw THAILAND VIETNAM IP
OP Equator MCMACU* TWOHOftTOMCd COtOlfK GUINEA g LIBERIA ' M** COAT NtGERW ETHIOPIA IaMBCO* PULL***S CAMP KtWJflUlL SSUDAN SOMALIA 5IUAHKA ^ v A, Y S j, BOON & U5AAO* JGMOA MALDV15 ^ ^ 4 Equator P
NflBATI If fCUAD3l Pt4l) B R A 7 < L GABON $ WANDA Â«*ya VCYolBJLB INOONESIA Mu nfo _ SLmjHOI ** TM â€™ WUÂ«MÂ» , CONGO Z(UWA mczakwx INDIAN OCEAN east ' '
2J*aar r ><>/> 30* * . scxma; f of Capricorn Paraguay SOUTH NAMBIA ZIMBAEWE MAURT1US g  â– 2T2* MADAGASCAR AUSTRALIA 3ir
b U 1 r â€ cmc; URUGUAY J â€¢ PACIFIC MutmiM Off 1 \ LAN DCF. A TIC N SOOTH AftKA if / 4S*
mr SOUTHERN OCEAN ji
Copyright* 201? wwwmapsotworktc Antarctic 0 â‚¬00 1800 3000 4200 Kma. W nr I2T 105Â° V 75* Â«â€¢ 45 *JrrJ< 30r 3 Â£ 15* Antarctic Circle j  ttfVS ANTARCTICA ir ar 4r wr 75 Â«cr ics* 12^ iav inr thsÂ° ^Rir
Figure 3.3 Propagation paths analyzed in this work along the east/west direction
To isolate the effect of the geomagnetic field, we attempt to hold all other factors constant. The first factor is the propagation distance. The other factor being the transmitter location. East and west propagation is tested over both land and sea. This is significant because the conductivity values of the ground vary between land and sea. For east and west propagation over land, the transmitterâ€™s location was at 90Â°E, 60Â°N. The westward receiver was located at 30Â°E, 60Â°N, and the eastward receiver was located at 150Â°E, 60Â°N. For propagation paths over sea, receivers and transmitters were held at the conjugate point so the propagation distances were not varied. The transmitter was held at 90Â°E, 60Â°S. The westward receiver was held 30Â°E, 60Â°S, and the eastward receiver was held at 150Â°E, 60Â°S. For a visual representation of the propagation paths,
28
refer to figure 3.3. A larger contributing factor than propagation over land and sea is most likely the propagation over the two hemispheres.
relative velocity (v/c)
0.4 0.45 0.5 0.55 0.6
beta
relative velocity (v/c)
0.4 0.45 0.5 0.55 0.6
beta
0.4 0.45 0.5 0.55 0.6
beta
90
88
86
84
82
78
76
74
72
70
relative velocity (v/c)
0.4 0.45 0.5 0.55 0.6
beta
0.98
0.97
0.96
0.95
0.94
0.93
0.92
0.91
0.9
Figure 3.4 2D color plots for the propagation paths corresponding to Figure 3.3
Figure 3.4 displays 2D color plots similar to figure 3.1 and 3.2 corresponding to the propagation paths in Figure 3.3. For all of these plots, the vaxis is the sharpness ranging from 0.4 to 0.6 inverse kilometers in 0.01 resolution. The yaxis is the reference height from 70 to 90 kilometers in one kilometer resolution. The intensity of the brightness denotes the propagation velocity with respect to the speed of light in a vacuum. The plots corresponding to west and east propagation over land are topleft and topright plots, respectively. The plots corresponding to west and east propagation over sea are bottomleft and bottomright, respectively. There is
29
again an oscillatory relationship between the reference height and the velocity in all four panels. There is also a linear relationship between the sharpness and the propagation velocity, which agrees with the formulation in Chapter 2. It should be noted that the eastward propagation path over land possesses a magnetic dip angle of 296.5651Â°. The westward propagation path over land possesses a magnetic dip angle of 63.4349Â°. The westward propagation path over sea possesses a magnetic dip angle of 116.5651Â°, which is 180Â° out of phase with the eastward propagation over land. The eastward propagation over sea is 180Â° out of phase with the westward propagation over land at 243.4349Â°.
The next set of simulations were held in the north/south direction varying the lattitude rather than longitude.
Figure 3.5 Propagation paths for the north/south propagation simulated in the LWPC
For comparing land and sea propagation, the transmitters are located on the equator at 30Â°E and 30Â°W longitude, respectively. For propagation over sea, the receiver for northward propagation is located at 30Â°E, 30Â°N, and the southward receiver is held at 30Â°E, 30Â°S. For propagation over land, the receiver for northward propagation is held at 30Â°W, 30Â°N, and the
30
receiver for southward propagation is at 30Â°W, 30Â°S. The magnetic dip angle for both northward paths is 0Â°. The southward paths are, of course, 180Â° out of phase with the northward paths.
Figure 3.6 Propagation paths analyzed in this work along the north/south direction
Again, the plots in Figure 3.6 correspond to the paths denoted in Figure 3.5. The topleft plot is for propagation in the northward direction over sea. The bottomleft plot is for propagation in the southward direction over sea. The topright plot is for northward propagation over land, and the bottomright plot is for southward propagation over land. The vaxis and yaxis are the same for all of the 2D color plots. The relationships between the sharpness and the propagation velocity are the same for all of the past figures in this chapter. Southward
31
propagation seems to share more of a constant propagation velocity across reference heights holding the sharpness constant. All of the other figures, the same oscillatory relationship between the reference height and the propagation velocity exists.
Land vs Sea
Relative Velocities Over Land Propagation
East
West
North
South
1 1.5 2 2.5 3 3.5 4
1h =77,beta=.4 2h =77,beta=.6 3h =85,beta=.4 4h =85,beta.6
Figure 3.7 Scatter plot for propagation over land. Each point on the xaxis denotes a separate profile. (1 â€” h' = 77km, (3 = .4/rm_1 2 â€” h'77km,fi = .6km1 3 â€” h' = QSkm.fi =
.4/cm_1 4 â€” h' = 85 km.fi = .6 km1)
Figure 3.7 shows what each point on the vaxis denotes, the yaxis is the propagation velocity with respect to the speed of light in a freespace medium. This makes sense since seawater possesses a higher conductivity value. Between points 1 and 2 and points 3 and 4, the reference height is held constant, and the sharpness is increased. The LWPC yields a higher
32
propagation velocity when the sharpness is increased for all paths. Between points 1 and 3 and points 2 and 4, the reference height is increased, and the sharpness is held constant. The propagation velocity decreases or is held constant between these points of reference. The reference height shares a less direct relationship with the magnitude of the velocity. It is important note that the northward propagation velocity is the lowest due to the fact that the propagation path is completely in phase with the Earthâ€™s geomagnetic field. Southward propagation generally possesses the highest propagation velocity being 180Â° out of phase with the geomagnetic field of the Earth.
Relative Velocities Over Sea Propagation
u.yo ' ' â€™ o T â€™ ^ 1 1
0.95 o 
0.94 1 + o 
'Tj 0 * â€¢
5 0.93
0
0.92 o + East 
0 West
0.91 O North 
â€¢ South
0.9 ( 1 1 1 1 1 1 1 1 1
) 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
1h =77.beta=.4 2h =77,beta=.6 3h =85,beta=.4 4h =85,beta=.6
p P P P
Figure 3.8 Scatter plot for propagation over sea. Each point on the xaxis denotes a separate profile. (1 â€” h' = 77km,(1 = .4km_1 2 â€” h'77km,fi = .6km1 3 â€” h' = 85km,(3 =
.4/cm_1 4 â€” h' = QSkm.fi = .6 km1)
Figure 3.8 is a scatter plot in all of the cardinal directions over sea. It should be noted that the average propagation velocity is greater over sea than over land. This makes sense because the ground conductivity of seawater is greater than the conductivity values of any land terrain. The
33
relationships for when either the sharpness is increased and the reference height is held constant or viceversa are very similar both over land and over sea. The magnetic dip angle variation does not seem to vary these relationships very much, either.
34
IV. LWPC vs Experimental Data
Variations of ELF Propagation Velocity Over DayNight Terminator
In this chapter, the ELF propagation velocity variation of sferics is employed as a diagnostic for daynight transitions. The variation of this quantity is also analyzed as a diagnostic for the occurrence of a solar induced sudden ionospheric disturbance (SID). As mentioned in Section 2.3, because the sharpness (and not the reference height) is much higher for nighttime ionospheric profiles rather than daytime profiles, it can be assumed that the magnitude of the propagation velocity is much higher at nighttime than at daytime. To validate the simulations, data from ELF receivers in Poland (49.209N, 22.553E) and Hugo, Colorado (103.4078N, 38.8918W) was used. These receivers continuously record in the ELF band ( 0 330 Hz). Identification of lightninginduced sferics along with lightning location and timing data from the GLD360 network allows for determining the ELF group velocity. In this study lightning source locations were limited to those occurring along the east coast of the Americas as shown in Figure 4.1.
35
Figure 4.1 Location of receivers and propagation paths of the recorded sferics at each receiver
The experimental data which consists of a smoothed average of over 20,000 individual sferics is shown in Figure 4.2.
Figure 4.2 Experimental Data of Propagation Velocity at both Poland and Colorado receivers from lightninginduced sferic propagation from along the east coast of the Americas
The hours 03 in Figure 4.2 correspond to nighttime across the entire path. This figure includes frequencies up to 350 Hz, corresponding to the passband of the ELF receiver. The minimums correspond to daytime across the entire path. When the propagation velocity is at an intermediate between the maximum and minimum, the propagation path overlaps with the daynight terminator. As the day in universal time begins, the path from Colorado to the east coast finishes its overlap with the daynight terminator and transitions to a path propagating through only nighttime ionospheric profiles. The Poland propagation path begins overlapping with the daynight terminator around 3:00 UT. Around 9:00 UT, the propagation path to the Poland
36
receiver no longer overlaps with the daynight terminator. Around 10:00 UT, the propagation path to the Colorado receiver begins to overlap with the daynight terminator one again and transitions to a daytime ionospheric profile from about 13:00 until about 23:00 UT. Around 18:00, the propagation path to the Poland receiver begins to overlap with the daynight terminator until about 23:00 UT when the propagation is through a complete nighttime waveguide configuration.
LWPC Results vs Greifmger Theory
Han and Cummer [Han and Cummer, 2010] developed a method of determining the reference height and the sharpness for any solar zenith angle range while no ionospheric disturbance is occurring. Since the solar zenith angle during anytime can be obtained quite easily, these parameters can be used as inputs for both the LWPC and analytical solutions in order to reproduce Figure 4.2. Table 4.1 is the values for the reference height and the sharpness for all ranges of zenith angles.
37
Zenith Angle (deg) Beta (km'1) h' (km)
030 0.4 71
3050 0.4 72
5060 0.4 73
6070 0.4 75
7080 0.4 76.5
8090 0.4 79
9091.8 0.44 80
91.893.6 0.48 81
93.695.4 0.53 82
95.497.2 0.57 82.9
97.299 0.61 83.9
99180 0.65 84.9
Table 4.1 Values for reference height and sharpness for various ranges of the zenith angle
Greifinger and Greifingerâ€™s [Greifinger and Greifinger, 1978] formulation computes the wave number and attenuation coefficient over the path of propagation for a ELF wave propagating in this Earthionosphere waveguide configuration. Referring back to Chapter 2, this work derives an analytical solution for the wave number, phase velocity and group velocity of ELF waves using the twoparameter ionospheric profile mentioned in Chapter 1 and the collisional frequency profile accepted in most literature as the inputs. Before analyzing LWPC results, analytical solutions for a partial day partial night terminator were developed. Using the power profile mentioned in Chapter 2, the formulation is displayed below.
Â£(/) = ^V2/Wo nfli exp (2fÂ£) exp(i^) (4.1)
Â£(t) = T'lEtfil (4.2)
38
Â£(t) is the computed timedomain electric field. E is the computed frequencydomain electric field. JV is the number of ionospheric profiles segmented into the path. For the LWPC, there is a maximum of ten ionospheric segmentations that can be applied to the propagation path. For the analytical computation, 50 segmentations were used; each of the segmentations were separated by equal distance, p is the propagation distance. The wave number, keff, is explicitly defined in Chapter 2. The effective attenuation is aeff, which is the imaginary component of the phase constant mentioned in Chapter 2 multiplied by the wave number of free space. The quantity Z0 is the effective wave impedance of free space, and P is the power from the power profile corresponding to a certain frequency. Summing these values for the electric field over a frequency range and then applying an inverse Fourier transform, a timedomain solution for the electric field is obtained. The ionospheric parameters used to compute the effective wave number correspond to the zenith angle of the location of the segmentation in the path and time of day in universal time. Figure 4.1 displays the analytical solution over the daynight terminator. The source location for the simulations is 17.7466Â°N, 64.7032Â°W, which is along the east coast within the region for the observational data.
39
0.82
Sferic Group Velocity to Hugo and Poland from Virgin Islands
0 5 10 15 20 25
Hour (UT)
Figure 4.3 Analytical Results for September 28, 2015 for ELF Propagation from Virgin Islands
to Hugo and Poland
All simulations in Figure 4.2 range from 10 Hz to the upper limit of either 350 Hz or 1000 Hz. The dotted lines are simulations from a frequency resolution of 3Hz. The solid lines are simulations from a frequency resolutions of 10 Hz. This figure seems to capture the overall daynight terminator variation reasonably well without capturing the finer more inhomogeneous features corresponding to the path of propagation. The finer features are likely generated by the conductivity of the ground and the orientation of the Earthâ€™s geomagnetic field on the path of propagation. This analytical solution assumes a homogeneous ground conductivity and does not account for variations caused by the geomagnetic dip angle. As mentioned in Chapter 3, the LWPC is a more complete simulation software and formulation than the Greifinger theory for that reason. Using Table 4.1 to compute the sharpness and reference height of the ionospheric profile and segmenting the path equally through ten segments, the LWPC was able to generate a timeseries for the same day.
40
Figure 4.4 LWPCgenerated Results for September 28, 2015 for ELF Propagation from Virgin Islands to Hugo and Poland for the frequency band of 301000hz in lOhz resolution
Figure 4.3 shows the results of LWPC simualtions for the same case of a single source at 17.7466Â°N, 64.7032Â°W. The frequency range and resolution is mentioned in the caption above. The LWPC captures more of the features corresponding to the inhomogeneities of the path than the analytical solution. The upper limit of the frequency range is high at one kilohertz, and this is the reason that the variation between maximum and minimum is not very high. The frequency resolution being so low also means that the variations in the figure are binary rather than being smooth transitions over the entire day. Figure 4.4 displays LWPC computations over the same path for the frequency range of 30330 Hz in 3 Hz resolution.
41
0.9
Time Series (9/28/2015) 30330Hz (3hz resolution)
Figure 4.5 LWPCgenerated Results for September 28, 2015 for ELF Propagation from Virgin Islands to Hugo and Poland for the frequency band of 30330hz in 3hz resolution
Due to a finer resolution, this simulation seems to capture the finer features of the daynight terminator variation in the observed data. The transitions are no longer binary, rather these are more continuous transitions. This is due to a higher frequency resolution. The lowering of the upper limit of the frequency band caused the difference between the maximum and minimum velocity values in the figure to be greater than figure 4.3. This is because lower frequencies yield a lower propagation velocity overall.
Diagnosing A Sudden Ionospheric Disturbance Using ELF Propagation Velocity Variation
A sudden ionospheric disturbance (SID) occurs when the sun emits a sudden Xray flux. This Xray flux can cause secondary ionization in the ionosphere. If the energy is high enough, valence electrons will be released from the ions in the ionosphere increasing the electron density. The increased density then lowers the reflection height and increases the sharpness. According
42
to the theory in Chapter 2, the ELF propagation velocity should increase under the occurrence of a SID.
Here we develop a linear model to estimate the variations in the reference height and sharpness under the occurrence of a SID. Han and Cummer [ Han and Cummer, 2010] show that the change in the reference height is linearly proportional to the logarithm of the magnetic flux density of the Xray.
h'(t) = h'd + dh' (4.3)
dh' = â€” 6.3 log10 Xfluxâ€” Xthreshold (4.4)
P = Pa +  e~dh'lH)
(4.5)
In this model, the Xray flux was inputted from GOES satellite data for September 28, 2015. The sudden ionospheric disturbance occurs around 15:00 UT. h'd and (3d are the undisturbed values for the reference height and sharpness, respectively. The undisturbed values of the reference height and sharpness are the values corresponding to the zenith angle at that time of day for each segment of the path. Each ionospheric segment of the path is then predicted using this linear estimation model above to predict the overall variation of the ionosphere along the path during a SID. Figure 4.5 shows the variation of the LWPCgenerated velocity during a sudden ionospheric disturbance according to this model.
43
Figure 4.6 Top panel displays time series for September 28, 2015 from the Virgin Islands to Hugo including SID calculations. Middle panel displays time series for September 29, 2015, which is a day a SID did not occur. Bottom panel is goes Xray flux data for September 28,
2015.
The ELF propagation velocity variation along the daynight terminator is roughly the same for both time series in Figure 4.5. At 15:00 UT on September 28, 2015, an Mclass sudden ionospheric disturbance occurred. Again, using the model above, this work simulates the ionospheric variation for all segments of the path under the occurrence of a SID. For both paths on September 28, the velocity seems to increase by 2% of the speed of light, roughly. Therefore,
44
it can be assumed that ELF propagation velocity variation can diagnose a sudden ionospheric disturbance.
Referring back to Figure 4.2, during the time of the SID (15:00 UT), the ELF propagation velocity increases by roughly 2% of the speed of light close to the computation from the LWPC. During a sudden ionospheric disturbance, the reference height decreases. Since the correlation between the reference height and the propagation velocity is oscillatory, the sharpness is the dominant factor even when the ionospheric parameters reach abnormal levels like such case. It should be noted that both paths are fully daytime when the SID occurs. In a future study, fully nighttime paths and paths that overlap with the daynight terminator should be studied to develop a better understanding of the variation of the ELF propagation velocity under a SID.
45
y.
Conclusion
The ELF propagation velocity was employed as a diagnostic for the ionosphere in this analysis. From the data and the analytical derivations, the sharpness contributed to the variation of the ELF propagation velocity more than any other quantity. This finding is very novel as in past works it was said that the reference height alters the ELF propagation velocity more than any other factor. This relationship was shown in the variation of the velocity during a SID. During a SID, the sharpness increases and the reference height decreases. Despite the decrease in the reference height, the propagation velocity increased. This finding shows that the ELF propagation velocity may be a good diagnostic for the variation in the composition of the ionosphere during a sudden ionospheric disturbance.
46
BIBLIOGRAPHY
Blaes, P. R. (2015). Remote Sensing and Statistical Analysis of the LightingIonosphere Interaction. (Unpublished doctoral dissertation). Stanford.
Han, F., & Cummer, S. A. (2010). Midlatitude daytime D region variations measured from radio atmospherics. Journal of Geophysical Research: Space Physics, II5(A9). doi: 10.1029/20 lOjaO 15437
Han, F., & Cummer, S. A. (2010). Midlatitude nighttime D region ionosphere variability on hourly to monthly time scales. Journal of Geophysical Research: Space Physics, 115( A9). doi: 10.1029/201 Oj aO 15437
Schmitter, E. D. (2013). Modeling solar flare induced lower ionosphere changes using VLF/LF transmitter amplitude and phase observations at a midlatitude site. Annales Geophvsicae,31(4), 765773. doi: 10.5194/angeo317652013
Greifinger, C., & Greifinger, P. (1978). Approximate method for determining ELF eigenvalues in the earthionosphere waveguide. Radio Science,13(5), 831837. doi: 10.1029/rs013i005p00831
Cummer, S. A. (1997). Lightning and ionospheric remote sensing using VLF/ELF radio atmospherics
Cotts, B. R., Golkowski, M., & Moore, R. (2011). Ionospheric effects of whistler waves from rockettriggered lightning. GEOPHYSICAL RESEARCH LETTERS,38.
Wait, J. R. (n.d.). ,On the propagation ofE.LF. pulses in the Earthionosphere waveguide (Vol. 40, Canadian J. Phys.)
Uman, M. A. (2001). The lightning discharge. Mineola: N.Y.
47
