Variables of the Equatorial Ionosphere

Variables of the Equatorial Ionosphere#

By Amadi Brians C.#

bamadi@brianspace.org#

#=========== Import Packages ==============

import os
import glob
import shutil
import matplotlib
import numpy as np
import pandas as pd
import datetime as dt
from pathlib import Path
import cartopy.crs as ccrs
from netCDF4 import Dataset 
from datetime import datetime
import matplotlib.pyplot as plt
import matplotlib.dates as dates
import cartopy.feature as cfeature
from netCDF4 import date2num, num2date
import matplotlib.gridspec as gridspec
from IPython.display import Video, display, HTML

import sys
import requests
import datetime
from urllib.parse import quote

#This user-defined package contains instruction for downloading
#download some files from zenodo, necessary to make 
#some plots such as magnetic equator, etc.
sys.path.append("..")  # path to your helper script folder
from utils.zenodo_tools import get_from_zenodo
from utils.get_dependencies import ensure_dependencies

#This code downloads the files and packages
#into a folder named dependencies
ensure_dependencies()

#Import User-defined packages and files downloaded
sys.path.append("..")  # adjust if needed
from dependencies.convert_waccmx_datesec import *   # or import specific functions
import dependencies.sha as sha

#Read the magnetic equator file
BASE_DIR = os.getcwd()  # current working directory
mag2geo = os.path.join(BASE_DIR, "dependencies", "mag2geo_all.csv")
df2 = pd.read_csv(mag2geo, delim_whitespace = False, header = 0)

# ===== File with WACCM-X EDens =====
wacx = Dataset('../data/WACCMX_subset.nc')

Neutral Winds#

  • Neutral winds; in simple term, are the motions of the atmosphere (uncharged components). They are characterized by oscillations from several sources described subsequently.

  • Sources of Neutral Winds

  • Solar Heating: At any given epoch, one side of Earth faces the Sun (Dayside) while the other side backs the Sun (Nightside). These two sides are separated by the Solar Terminator. The atmosphere on the Dayside expands as a result of solar heating. This creates a pressure gradient that drives atmospheric constituents from the Dayside to Nightside, hence Neutral winds. Watch Video 2.1 to see how winds are generated from solar heating.

#Unmark cell to run code

# ========Display Video========
video = Video("../data/Solheating.mp4", embed=True, width=500, height=400)

#======= Display Source ========
display(HTML(f"""
<div style="display:inline-block; text-align:center;">
    {video._repr_html_()}
    <div style="text-align:right; font-size:12px; color:gray; margin-top:5px;">
        Source: NASA
    </div>
</div>
"""))

  • Tides: These are periodic oscillations observed in the atmosphere due to solar heating, or lunar gravitational pull as shown in Figure 2.1. Tides with period of 24 (Diurnal) and 12 (semi-diurnal) hours are due to solar heating while tides with periods of 24.8 hours are lunar-associated.

# Time in hours (2 days)
t = np.linspace(0, 48, 1000)

# Tidal harmonics
# Diurnal (24h), Semidiurnal (12h), Terdiurnal (8h)
diurnal = 50 * np.sin(2 * np.pi * t / 24)        # amplitude ~50 m/s
semidiurnal = 30 * np.sin(2 * np.pi * t / 12)    # amplitude ~30 m/s
terdiurnal = 10 * np.sin(2 * np.pi * t / 8)      # amplitude ~10 m/s

# Total neutral wind variation (zonal component)
neutral_wind = diurnal + semidiurnal + terdiurnal

# Plot
plt.figure(figsize=(10,5))
plt.plot(t, neutral_wind, label="Total neutral wind", color="k")
plt.plot(t, diurnal, '--', label="Diurnal (24h)")
plt.plot(t, semidiurnal, '--', label="Semidiurnal (12h)")
plt.plot(t, terdiurnal, '--', label="Terdiurnal (8h)")
plt.xlabel("Local time (hours)")
plt.ylabel("Neutral wind (m/s)")
plt.title("Figure 2.1: Synthetic Ionospheric Tides in Neutral Wind")
plt.legend()
plt.grid(True)
plt.show()

  • Don’t be confuse about these nomenclature. These tides are simply defined by the time between any two consecutive peaks (or crest). The time between two neighboring peaks of a diurnal oscillation is 24 hours. This is shown in Figure 2.2 and represented by the blue lines.

# Time in hours (2 days)
t = np.linspace(0, 48, 1000)

# Tidal harmonics
diurnal = 50 * np.sin(2 * np.pi * t / 24)        # amplitude ~50 m/s
semidiurnal = 30 * np.sin(2 * np.pi * t / 12)    # amplitude ~30 m/s
terdiurnal = 10 * np.sin(2 * np.pi * t / 8)      # amplitude ~10 m/s

# Total neutral wind variation (zonal component)
neutral_wind = diurnal + semidiurnal + terdiurnal

# Plot
fig, ax = plt.subplots(figsize=(12,5))
ax.plot(t, neutral_wind, label="Total neutral wind", color="k")
ax.plot(t, diurnal, '--', label="Diurnal (24h)")
ax.plot(t, semidiurnal, '--', label="Semidiurnal (12h)")
ax.plot(t, terdiurnal, '--', label="Terdiurnal (8h)")

ax.set_xlabel("Local time (hours)")
ax.set_ylabel("Neutral wind (m/s)")
ax.set_title("Figure 2.2: Synthetic Ionospheric Tides in Neutral Wind")
ax.grid(True)
ax.legend()

# Vertical lines marking one diurnal period (you can change the pair to any [t0, t0+24])
t0 = 6
t1 = t0 + 24
ax.vlines([t0, t1], ymin=min(neutral_wind)*1.05, ymax=max(neutral_wind)*1.05,
          colors='blue', linestyles='--', linewidth=2)

# Horizontal double-headed arrow showing the diurnal period
# place it slightly above the plot max
y_arrow = max(neutral_wind) * 1.12
ax.annotate(
    '', xy=(t0, y_arrow), xytext=(t1, y_arrow),
    arrowprops=dict(arrowstyle='<->', lw=2, color='blue')
)
# label the arrow centered
ax.text((t0 + t1) / 2, y_arrow * 1.01, 'Diurnal — 24 h',
        ha='center', va='bottom', fontsize=10, color='blue', weight='bold')

# tighten limits so arrow & lines are visible
ax.set_xlim(t.min(), t.max())
ax.set_ylim(min(neutral_wind) * 1.05, y_arrow * 1.15)

plt.show()

  • Similarly, the time (duration) between any two consecutive peaks (or crests) for a semi-diurnal oscillation is 12 hours. This is shown in Figure 2.3 and indicated by the orange line.

# Time in hours (2 days)
t = np.linspace(0, 48, 1000)

# Tidal harmonics
diurnal = 50 * np.sin(2 * np.pi * t / 24)        # amplitude ~50 m/s
semidiurnal = 30 * np.sin(2 * np.pi * t / 12)    # amplitude ~30 m/s
terdiurnal = 10 * np.sin(2 * np.pi * t / 8)      # amplitude ~10 m/s

# Total neutral wind variation (zonal component)
neutral_wind = diurnal + semidiurnal + terdiurnal

# Plot
fig, ax = plt.subplots(figsize=(12,5))
ax.plot(t, neutral_wind, label="Total neutral wind", color="k")
ax.plot(t, diurnal, '--', label="Diurnal (24h)")
ax.plot(t, semidiurnal, '--', label="Semidiurnal (12h)")
ax.plot(t, terdiurnal, '--', label="Terdiurnal (8h)")

ax.set_xlabel("Local time (hours)")
ax.set_ylabel("Neutral wind (m/s)")
ax.set_title("Figure 2.3: Synthetic Ionospheric Tides in Neutral Wind")
ax.grid(True)
ax.legend()

# Vertical lines marking one diurnal period (you can change the pair to any [t0, t0+24])
t0 = 15
t1 = t0 + 12
ax.vlines([t0, t1], ymin=min(neutral_wind)*1.05, ymax=max(neutral_wind)*1.05,
          colors='orange', linestyles='--', linewidth=2)

# Horizontal double-headed arrow showing the diurnal period
# place it slightly above the plot max
y_arrow = max(neutral_wind) * 1.12
ax.annotate(
    '', xy=(t0, y_arrow), xytext=(t1, y_arrow),
    arrowprops=dict(arrowstyle='<->', lw=2, color='orange')
)
# label the arrow centered
ax.text((t0 + t1) / 2, y_arrow * 1.01, 'Semidiurnal — 12 h',
        ha='center', va='bottom', fontsize=10, color='blue', weight='bold')

# tighten limits so arrow & lines are visible
ax.set_xlim(t.min(), t.max())
ax.set_ylim(min(neutral_wind) * 1.05, y_arrow * 1.15)

plt.show()

  • There is also a special class of tides called terdiurnal. They have a time range of 8 hours between any two consecutive peaks as shown in Figure 2.4.

# Time in hours (2 days)
t = np.linspace(0, 48, 1000)

# Tidal harmonics
diurnal = 50 * np.sin(2 * np.pi * t / 24)        # amplitude ~50 m/s
semidiurnal = 30 * np.sin(2 * np.pi * t / 12)    # amplitude ~30 m/s
terdiurnal = 10 * np.sin(2 * np.pi * t / 8)      # amplitude ~10 m/s

# Total neutral wind variation (zonal component)
neutral_wind = diurnal + semidiurnal + terdiurnal

# Plot
fig, ax = plt.subplots(figsize=(12,5))
ax.plot(t, neutral_wind, label="Total neutral wind", color="k")
ax.plot(t, diurnal, '--', label="Diurnal (24h)")
ax.plot(t, semidiurnal, '--', label="Semidiurnal (12h)")
ax.plot(t, terdiurnal, '--', label="Terdiurnal (8h)")

ax.set_xlabel("Local time (hours)")
ax.set_ylabel("Neutral wind (m/s)")
ax.set_title("Figure 2.4: Synthetic Ionospheric Tides in Neutral Wind")
ax.grid(True)
ax.legend()

# Vertical lines marking one diurnal period (you can change the pair to any [t0, t0+24])
t0 = 18
t1 = t0 + 8
ax.vlines([t0, t1], ymin=min(neutral_wind)*1.05, ymax=max(neutral_wind)*1.05,
          colors='green', linestyles='--', linewidth=2)

# Horizontal double-headed arrow showing the diurnal period
# place it slightly above the plot max
y_arrow = max(neutral_wind) * 1.12
ax.annotate(
    '', xy=(t0, y_arrow), xytext=(t1, y_arrow),
    arrowprops=dict(arrowstyle='<->', lw=2, color='blue')
)
# label the arrow centered
ax.text((t0 + t1) / 2, y_arrow * 1.01, 'Terdiurnal — 8 h',
        ha='center', va='bottom', fontsize=10, color='blue', weight='bold')

# tighten limits so arrow & lines are visible
ax.set_xlim(t.min(), t.max())
ax.set_ylim(min(neutral_wind) * 1.05, y_arrow * 1.15)

plt.show()

  • Gravity Waves: Consider a bulk of atmospheric gases hitting an elevated feature such as a mountain. They will be displaced upwards, and then downwards due to gravity. This process gives rise to an upward and downward oscillation thus producing a unique class of waves known as Atmospheric Gravity Waves (AGWs). These AGWs produced by mountain or other elevations are specifically termed Mountain (or Orographic) Waves. A snapshot of these waves is shown in Figure 2.5.

  • Rossby (Planetary) Waves: As Earth rotates, does it cause the neutral atmosphere to move? Yes, it does. This motion is slow but can effectively modify the global circulation of the atmosphere and hence, neutral winds.

# Grid
x = np.linspace(0, 2000, 200)   # horizontal distance [km]
z = np.linspace(100, 300, 100)     # altitude [km]
X, Z = np.meshgrid(x, z)

# Wave parameters
A = 20             # amplitude (m/s or K)
kx = 2*np.pi/500   # horizontal wavelength 500 km
mz = 2*np.pi/20    # vertical wavelength 20 km
omega = 2*np.pi/8  # period = 8 hours
c_phase = omega/kx # phase speed

# Time axis
t = np.linspace(0, 24, 200)   # 24 hours

# Function to generate wave field (e.g. wind perturbation)
def wave_field(tt):
    return A * np.sin(kx*X + mz*Z - omega*tt)

# --- Static snapshot ---
plt.figure(figsize=(8,5))
snapshot = wave_field(6)  # snapshot at t=6 hr
plt.contourf(X, Z, snapshot, levels=20, cmap="RdBu_r")
plt.colorbar(label="Perturbation")
plt.xlabel("Horizontal distance (km)")
plt.ylabel("Altitude (km)")
plt.title("Figure 2.5: Atmospheric wave snapshot (t=6 hr)")
plt.show()

plt.show()

  • Video 2.2 shows gravity waves tracked by NASA scientists. Observe that these waves exhibit circular wavefronts. Hence, they are called concentric gravity waves (CGWs).

# ========Display Video========
video = Video("../data/GravityWave.mp4", embed=True, width=500, height=400)

#======= Display Source ========
display(HTML(f"""
<div style="display:inline-block; text-align:center;">
    {video._repr_html_()}
    <div style="text-align:right; font-size:12px; color:gray; margin-top:5px;">
        Source: NASA
    </div>
</div>
"""))

Magnetic Field#

  • An important variable in equatorial ionosphere and many other areas in space is the Magnetic field. In simple terms, it is the vector field used in describing the Earth’s magnetic influence or any other magnetized or magnetic material.

    • Sources of Earth’s Magnetic Field

      • Geodynamo: The motion of molten iron in Earth’s outer core is the main source of Earth’s magnetic Field. It produces this field by driving electric current.

      • Space Current: A small but significant source of magnetic field is the variation in magnetospheric and ionospheric currents, commonly referred to as External currents since they are formed beyond Earth’s surface unlike the geodynamo current generated within Earth.

Electric Field#

  • Since Electric Field is the field around a charge, hence, charge distribution or density creates electric fields as represented in Poisson equation

\begin{equation} \nabla\cdot \textbf{E} = \rho_{c}/\epsilon_{0} \end{equation}

Considering electrostatics, \(\textbf{E}\) can be written as: \begin{equation} \textbf{E} = -\nabla{\phi} \end{equation}

where \(\phi\) is the potential and a function of charge distribution. Hence, the sources of electric field in the equatorial ionosphere are the same as agents responsible for charge distribution.

  • Sources of Electric Fields

    • Neutral Winds: Under the influence of atmospheric motion, ions and electrons drift differentially, thus creating charge separation. This leads to potential gradient, known as Electric Field.

      • An implication of this polarization electric field is the redistribution of plasma to ensure current continuity.

    • Magnetic and gravitational Fields: Similar to neutral winds, pressure gradient of Magnetic and Gravitational fields cause ions and electrons to exhibit differential motion that subsequently create electric fields.

Orientation/Description of Variables#

  • As clearly seen, these variables have direction and magnitude. How do we describe these variables with respect to the globe. For ease of illustration, let’s use neutral winds for this description.

  • Winds moving towards the North or South are said to be moving poleward, while winds crossing the equator are specifically referred to as transequatorial (or meridional) as shown in Figure 2.6.

#create a variable for datesec
variable = wacx['datesec'][:]
wac_time = time_conv(variable)

fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(projection=ccrs.Orthographic(central_longitude = -40))
ax.gridlines(lw = 2, color = 'gray', ls = '--')

im = ax.pcolormesh(wacx['lon'][:],  wacx['lat'][:], wacx['EDens'][6, 12, :, :], 
                   transform=ccrs.PlateCarree(), cmap='plasma')

# scatter points (use ax.scatter so transform is honored)
size = 50
i = 0
lons = df2['lon__' + str(i)]
lats = df2['lat__' + str(i)]
ax.scatter(lons, lats, transform=ccrs.PlateCarree(), 
            marker='o', color = 'r', s = size,)

# Draw coastlines on top
ax.coastlines(color='black', linewidth=1)

ax.coastlines()
ax.set_global()
ax.set_ylabel('Latitude (deg)')
ax.set_xlabel('Longitude (deg)')

# ---------- Meridional ----------
# xy = arrow tip location, xytext = label/arrow start location
ax.annotate("", xy=(-40, -40), xytext=(-40, 40),
            xycoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            textcoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            arrowprops=dict(arrowstyle='<->', facecolor='black', lw=2))

# Add the label at the middle of the arrow
ax.text(-36, 0, "Meridional", transform=ccrs.PlateCarree(),
        ha="center", va="center", fontsize=16, color="black",
         rotation=90, rotation_mode="anchor",
        bbox=dict(facecolor="white", alpha=0.6, edgecolor="none"))

# ---------- North arrow ----------
# xy = arrow tip location, xytext = label/arrow start location
ax.annotate("Northward", xy=(-50, 40), xytext=(-60, 10),
            xycoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            textcoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            arrowprops=dict(facecolor='black', width=2, headwidth=8),
           fontsize=15)

ax.annotate("Southward", xy=(-50, -40), xytext=(-60, -10),
            xycoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            textcoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            arrowprops=dict(facecolor='black', width=2, headwidth=8),
           fontsize=15)

# ---------- Transequatorial ----------
# xy = arrow tip location, xytext = label/arrow start location
ax.annotate("", xy=(-15, -20), xytext=(-15, 20),
            xycoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            textcoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            arrowprops=dict(arrowstyle='<->', facecolor='black', lw=2))

ax.annotate("", xy=(-10, -20), xytext=(-10, 20),
            xycoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            textcoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            arrowprops=dict(arrowstyle='<->', facecolor='black', lw=2))

# Add the label at the middle of the arrow
ax.text(-8, 0, "Transequatorial", transform=ccrs.PlateCarree(),
        ha="center", va="center", fontsize=16, color="black",
         rotation=90, rotation_mode="anchor",
        bbox=dict(facecolor="white", alpha=0.6, edgecolor="none"))

plt.show()

  • On the other hand, winds in the west-east or east-west direction are said to be zonal.

  • The terms Northward, Southward, Eastward or Westward are also valid for these winds if headed North, South, East or West respectively.

  • If winds are moving from the poles (North or South) to the equator, they can be called equatorward.

fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(projection=ccrs.Orthographic(central_longitude = -40))
ax.gridlines(lw = 2, color = 'gray', ls = '--')

im = ax.pcolormesh(wacx['lon'][:],  wacx['lat'][:], wacx['EDens'][6, 12, :, :], 
                   transform=ccrs.PlateCarree(), cmap='plasma')

# scatter points (use ax.scatter so transform is honored)
size = 50
i = 0
lons = df2['lon__' + str(i)]
lats = df2['lat__' + str(i)]
ax.scatter(lons, lats, transform=ccrs.PlateCarree(), 
            marker='o', color = 'r', s = size,)

# Draw coastlines on top
ax.coastlines(color='black', linewidth=1)

ax.coastlines()
ax.set_global()
ax.set_ylabel('Latitude (deg)')
ax.set_xlabel('Longitude (deg)')

# ---------- Zonal ----------
# ---------- East–West arrow (double-headed with rotated centered text) ----------
# Draw the arrow only (from lon=-20 to lon=20 at lat=0)
ax.annotate("", xy=(45, 0), xytext=(-120, 0),
            xycoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            textcoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            arrowprops=dict(arrowstyle='<->', color='black', lw=2))

# Add the label at the middle of the arrow
ax.text(-36, -5, "Zonal", transform=ccrs.PlateCarree(),
        ha="center", va="center", fontsize=16, color="black",
         rotation=0, rotation_mode="anchor",
        bbox=dict(facecolor="white", alpha=0.6, edgecolor="none"))

ax.annotate("", xy=(45, -15), xytext=(-20, -15),
            xycoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            textcoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            arrowprops=dict(facecolor='black', width=2, headwidth=8),
           fontsize=15)

# Add the label at the middle of the arrow
ax.text(-8, -20, "Eastward", transform=ccrs.PlateCarree(),
        ha="center", va="center", fontsize=16, color="black",
         rotation=0, rotation_mode="anchor",
        bbox=dict(facecolor="white", alpha=0.6, edgecolor="none"))

ax.annotate("", xy=(-120, -15), xytext=(-60, -15),
            xycoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            textcoords=ccrs.PlateCarree()._as_mpl_transform(ax),
            arrowprops=dict(facecolor='black', width=2, headwidth=8),
           fontsize=15)

# Add the label at the middle of the arrow
ax.text(-75, -20, "Westward", transform=ccrs.PlateCarree(),
        ha="center", va="center", fontsize=16, color="black",
         rotation=0, rotation_mode="anchor",
        bbox=dict(facecolor="white", alpha=0.6, edgecolor="none"))

plt.show()

These descriptions are used for neutral winds and electric fields. A similar description is used for the magnetic field but there are some technicalities to be aware of.

  • The magnetic field, at a given location and epoch, has the vertical (positive downwards), and horizontal (Northward/Eastward)components. There is a displacement between the geographic and geomagnetic components known as the declination. This angle varies from one longitudinal sector to another.

  • Additionally, the imaginary lines representing this field emerges in the magnetic South and re-enters earth in the magnetic North pole. As earlier stated, it is horizontal at the equator but vertical at the poles. The angle formed by magnetic field lines with respect to the horizontal is known as Inclination or angle of Dip.

Contents