6.1 Free Equatorial Waves

Learning Objectives

After completing this lecture, students should be able to:

understand the role of the equatorial beta approximation
interpret the physical meaning of the equatorial shallow-water equations
solve the equatorial shallow-water equations numerically
analyze numerical solutions of equatorial waves
distinguish among Kelvin, inertia-gravity, and Rossby waves
interpret the meridional structure and propagation of equatorially trapped waves

Concept

The shallow-water equations provide one of the simplest dynamical frameworks for studying large-scale tropical motions. These equations describe the evolution of horizontal winds and pressure anomalies in a single fluid layer and capture many important features of tropical atmosphere-ocean dynamics.

Near the equator, the Coriolis parameter

\[ f = 2\Omega \sin\phi \]

becomes very small and changes sign across the equator. Rather than treating the Coriolis parameter as a constant, we approximate it using a first-order Taylor expansion about the equator,

\[ f \approx \beta y, \]

where

\[ \beta = \frac{2\Omega}{a} \]

is the meridional gradient of the Coriolis parameter, \(a\) is Earth's radius, and \(y\) is the north-south distance from the equator.

This approximation is known as the equatorial beta approximation.

The beta approximation allows the Coriolis force to vary continuously across the equator:

\(f = 0\) at the equator,
\(f > 0\) in the Northern Hemisphere,
\(f < 0\) in the Southern Hemisphere.

As a result, fluid motions experience a restoring effect that tends to trap disturbances near the equator. This mechanism gives rise to a family of equatorially trapped waves, including Kelvin waves, Rossby waves, and inertio-gravity waves.

To investigate these dynamics, we use the linear shallow-water equations on an equatorial beta plane.

Model / Equation

Dimensional Equations

Using the beta approximation, the dimensional linear shallow-water equations can be written as

\[ \frac{\partial u}{\partial t} = \beta y v - \frac{\partial \Phi}{\partial x}, \]

\[ \frac{\partial v}{\partial t} = -\beta y u - \frac{\partial \Phi}{\partial y}, \]

\[ \frac{\partial \Phi}{\partial t} = -c^2 \left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right), \]

where

\(u\) is the zonal velocity,
\(v\) is the meridional velocity,
\(\Phi\) is geopotential,
\(\beta\) is the meridional gradient of the Coriolis parameter,
\(c\) is the gravity-wave speed.

The first equation describes zonal momentum. The term \(\beta y v\) represents the Coriolis force, while the pressure-gradient force is represented by \(-\partial\Phi/\partial x\).

The second equation describes meridional momentum. The term \(-\beta y u\) represents the Coriolis force, while \(-\partial\Phi/\partial y\) represents the meridional pressure-gradient force.

The third equation is the continuity equation, which relates pressure changes to horizontal convergence and divergence.

Nondimensional Form

To simplify the analysis, the equations can be nondimensionalized using characteristic length, time, and velocity scales. In nondimensional form, the equatorial shallow-water equations become

\[ \frac{\partial u}{\partial t} = yv - \frac{\partial p}{\partial x}, \]

\[ \frac{\partial v}{\partial t} = -yu - \frac{\partial p}{\partial y}, \]

\[ \frac{\partial p}{\partial t} = -\frac{\partial u}{\partial x} - \frac{\partial v}{\partial y}. \]

These equations describe freely propagating equatorial waves. The beta effect is represented by the meridional coordinate \(y\), which acts as the nondimensional Coriolis parameter.

Analytical solutions of this system lead to equatorially trapped Kelvin, Rossby, and inertio-gravity waves.

Discrete Form Used in the Numerical Model

The Python model solves the nondimensional shallow-water equations using a forward-Euler time integration scheme. The time derivative is approximated as

\[ \frac{\partial A}{\partial t} \approx \frac{A^{n+1}-A^n}{\Delta t}, \]

where \(A^n\) denotes the value at time step \(n\).

Applying this approximation to the shallow-water equations gives

\[ u^{n+1} = u^{n} + \Delta t \left( yv - \frac{\partial p}{\partial x} \right)^n, \]

\[ v^{n+1} = v^{n} + \Delta t \left( -yu - \frac{\partial p}{\partial y} \right)^n, \]

\[ p^{n+1} = p^{n} + \Delta t \left( -\frac{\partial u}{\partial x} - \frac{\partial v}{\partial y} \right)^n. \]

These equations provide a simple numerical method for integrating the equatorial shallow-water equations forward in time.

The Python code uses centered finite differences for the spatial derivatives and a forward-Euler scheme for the time integration. Although this approach is straightforward to implement, numerical stability depends on the choice of time step and grid spacing.

Free Equatorial Waves

When

\[ Q = 0, \qquad r = 0, \qquad \epsilon = 0, \]

the equatorial shallow-water equations support freely propagating wave solutions. Following Matsuno (1966), solutions can be expressed in the form

\[ [u,v,p] = [\hat{u}(y),\hat{v}(y),\hat{p}(y)] e^{i(kx-\omega t)}, \]

where

\(k\) is the zonal wavenumber,
\(\omega\) is the wave frequency,
\(\hat{u}(y)\) is the meridional structure of the zonal velocity,
\(\hat{v}(y)\) is the meridional structure of the meridional velocity,
\(\hat{p}(y)\) is the meridional structure of the pressure anomaly.

The exponential factor

\[ e^{i(kx-\omega t)} \]

describes the zonal propagation and temporal oscillation of the wave. The sign of \(\omega\) determines the propagation direction, while the magnitude of \(\omega\) determines the oscillation frequency.

Using Euler's formula,

\[ e^{i(kx-\omega t)} = \cos(kx-\omega t) + i\sin(kx-\omega t), \]

the physical solution is obtained by taking the real part of the complex expression.

In the analytical solutions presented below, the meridional velocity \(\hat{v}\) contains the imaginary unit \(i\). This indicates that the meridional velocity is 90 degrees out of phase with the zonal velocity and pressure fields.

For example, if the zonal velocity and pressure are represented by

\[ u,p \propto \sin(kx-\omega t), \]

then the meridional velocity is proportional to

\[ v \propto \cos(kx-\omega t). \]

Therefore, in the Python simulations used in this course, the zonal velocity and pressure are initialized using a sine function, while the meridional velocity is initialized using a cosine function. This approach produces the same phase relationship as the analytical complex solutions without explicitly using complex numbers in the numerical model.

The functions

\[ \hat{u}(y), \qquad \hat{v}(y), \qquad \hat{p}(y) \]

describe how the wave amplitude varies with latitude. These meridional structures determine whether the solution corresponds to a Kelvin wave, Rossby wave, or inertio-gravity wave.

Substituting the wave solution into the equatorial shallow-water equations transforms the problem into an eigenvalue problem. The resulting equations have solutions that are trapped near the equator and can be expressed using Hermite functions.

Meridional Structures

The meridional structure of equatorial waves is described by Hermite functions

\[ \psi_n(y) = H_n(y)e^{-y^2/2}, \]

where

\(H_n(y)\) is the Hermite polynomial,
\(n\) is the meridional mode number.

The Gaussian factor

\[ e^{-y^2/2} \]

acts as an equatorial trapping function. It has its maximum amplitude at the equator (\(y=0\)) and decreases exponentially with distance from the equator. As a result, wave amplitudes become progressively weaker at higher latitudes.

The Hermite polynomial \(H_n(y)\) determines the meridional structure of the wave by specifying the number and location of meridional nodes. Different values of \(n\) therefore produce different wave structures while remaining confined to the equatorial region by the Gaussian envelope.

The lowest-order Hermite functions are

\[ \psi_0=e^{-y^2/2}, \]

\[ \psi_1=2ye^{-y^2/2}, \]

\[ \psi_2=(4y^2-2)e^{-y^2/2}. \]

These functions correspond to increasingly complex meridional structures. The \(n=0\) mode has a single maximum at the equator, while higher-order modes develop additional meridional nodes and extend farther away from the equator.

For each wave family, the meridional structures

\[ \hat{u}(y), \qquad \hat{v}(y), \qquad \hat{p}(y) \]

can be expressed as combinations of Hermite functions. The meridional mode number \(n\) determines the wave structure, while the dispersion relationship determines the relationship between zonal wavenumber \(k\) and frequency \(\omega\).

Different equatorial wave families correspond to different values of \(n\) and satisfy different dispersion relationships.

Kelvin Wave (\(n=-1\))

Dispersion Relationship

\[ \omega = k. \]

The phase speed is always positive, so Kelvin waves propagate only eastward.

Meridional Structure

\[ \hat{u}= \psi_0, \]

\[ \hat{v} = 0, \]

\[ \hat{p} = \psi_0. \]

Characteristics

eastward propagation only
no meridional velocity
pressure and zonal wind are in phase
strongest amplitude at the equator

Inertio-Gravity Wave (\(n=0\))

Dispersion Relationship

\[ \omega^2 = k^2+1. \]

Therefore,

\[ \omega = \pm\sqrt{k^2+1}, \]

giving eastward and westward propagating solutions.

Meridional Structure

\[ \hat{u} = \psi_1, \]

\[ \hat{v} = 2 i\, (\omega+k)\psi_0, \]

\[ \hat{p} = \psi_1. \]

Characteristics

eastward and westward propagation
strongest gravity-wave character
symmetric \(v\) structure
antisymmetric \(u\) and \(p\) structures

Rossby and Inertio-Gravity Waves (\(n \ge 1\))

For higher meridional modes, the solutions have the general form

\[ \hat{v} = i\, (\omega^2-k^2)\psi_n, \]

\[ \hat{u} = \frac{\omega-k}{2}\psi_{n+1} + n(\omega+k)\psi_{n-1}, \]

\[ \hat{p} = \frac{\omega-k}{2}\psi_{n+1} - n(\omega+k)\psi_{n-1}. \]

The dispersion relationship determines whether the solution corresponds to a Rossby wave or an inertio-gravity wave.

Rossby Wave (\(n \ge 1\))

Dispersion Relationship

\[ \omega = -\frac{k} {k^2+2n+1}. \]

The negative sign indicates westward propagation.

Example: \(n=1\)

\[ \hat{v} = i\, \psi_1, \]

\[ \hat{u} = \frac{\omega-k}{2}\psi_2 + (\omega+k)\psi_0, \]

\[ \hat{p} = \frac{\omega-k}{2}\psi_2 - (\omega+k)\psi_0. \]

Example: \(n=2\)

\[ \hat{v} = i\, \psi_2, \]

\[ \hat{u} = \frac{\omega-k}{2}\psi_3 + 2(\omega+k)\psi_1, \]

\[ \hat{p} = \frac{\omega-k}{2}\psi_3 - 2(\omega+k)\psi_1. \]

Characteristics

westward propagation
slower than gravity waves
strongly influenced by the beta effect
increasingly complex meridional structure as \(n\) increases

Inertio-Gravity Wave (\(n \ge 1\))

Dispersion Relationship

\[ \omega^2 = k^2+2n+1. \]

Therefore,

\[ \omega = \pm\sqrt{k^2+2n+1}. \]

Example: \(n=1\)

\[ \hat{v} = i\, \psi_1, \]

\[ \hat{u} = \frac{\omega-k}{2}\psi_2 + (\omega+k)\psi_0, \]

\[ \hat{p} = \frac{\omega-k}{2}\psi_2 - (\omega+k)\psi_0. \]

Example: \(n=2\)

\[ \hat{v} = i\, \psi_2, \]

\[ \hat{u} = \frac{\omega-k}{2}\psi_3 + 2(\omega+k)\psi_1, \]

\[ \hat{p} = \frac{\omega-k}{2}\psi_3 - 2(\omega+k)\psi_1. \]

Characteristics

eastward and westward propagation
higher frequencies than Rossby waves
increasingly complex meridional structures for larger \(n\)

Summary

Wave Type	Mode Number	Dispersion Relationship
Kelvin	\(n=-1\)	\(\omega=k\)
Inertio-Gravity	\(n=0\)	\(\omega=\pm\sqrt{k^2+1}\)
Rossby	\(n\ge1\)	\(\omega=-k/(k^2+2n+1)\)
Inertio-Gravity	\(n\ge1\)	\(\omega=\pm\sqrt{k^2+2n+1}\)

The Python simulations in this module use these analytical solutions as initial conditions to demonstrate the propagation and structure of equatorial waves.

Connection to the Python Model

The Python model solves the nondimensional shallow-water equations using a forward-Euler time integration scheme.

For example, the zonal momentum equation is implemented as

dudt = Y * v - ddx(p) - r * u
u_new = u + DT * dudt

For free equatorial waves, the damping and forcing terms are set to zero:

r = 0.0
eps = 0.0
Q = 0.0

The model then evolves the specified initial condition forward in time.

Sample Code

Python scripts for this lecture are provided separately.

Template file:

Template: 6_1-shallow_water_template.py

Example command:

python 6_1-shallow_water_template.py

Exercise

Follow the steps below to build the numerical solutions gradually.

Step 1. Run the Template Model and Understand the Code Structure

Start from the original shallow-water template script. The purpose of this step is to confirm that the model runs correctly and to understand how the numerical model is organized.

First, create a copy of the original script:

cp 6_1-shallow_water_template.py 6_1-shallow_water_step1.py

Then run the script:

python 6_1-shallow_water_step1.py

After the script finishes, examine the generated output files:

6_1-shallow_water_step1.mp4
6_1-shallow_water_step1.png

In this default configuration, the model uses zero initial conditions and zero forcing. Therefore, no wave propagation or circulation response should develop.

Step 1.1 Explore the Grid System

Find the grid section in the script:

DX = 0.2
DY = 0.2

x = np.arange(-10 + DX/2, 10 - DX/2, DX)
y = np.arange(-5 - DY/2, 5 + DY/2, DY)

X = x[None, :]
Y = y[:, None]

Experiment with different grid spacings:

DX = 0.5
DY = 0.5

and

DX = 0.1
DY = 0.1

Run the model after each modification.

Focus on the following questions:

How does changing DX and DY affect the number of grid points?
How does changing the grid spacing affect the computational cost?
How do the dimensions of X, Y, u, v, and p change?

After exploring these sensitivities, restore the original values:

DX = 0.2
DY = 0.2

Step 1.2 Explore the Time Parameters

Find the time parameters:

DT = 0.02
NT = 600

Experiment with the following values:

DT = 0.05
NT = 200

and

DT = 0.01
NT = 1200

Run the model after each modification.

Focus on the following questions:

How does changing DT affect the time resolution?
How does changing NT affect the total integration length?
What is the total simulation time given by

\[ T = DT \times NT ? \]

How does changing these parameters affect the animation length and computational cost?

After exploring these sensitivities, restore the original values:

DT = 0.02
NT = 600

Step 1.3 Examine the Initial Conditions

Find the function:

def initial_condition(X, Y):

This function defines the initial values of \(u\), \(v\), and \(p\). In the default configuration, the amplitude is set to zero:

A = 0.0
p = A * env_x * env_y * np.sin(phase)
u = np.zeros_like(p)
v = np.zeros_like(p)

Focus on the following questions:

What variables are returned by initial_condition()?
Why does setting A = 0.0 produce zero pressure anomaly?
Why are u and v also set to zero?
What would happen if the Kelvin-wave block were activated?

Step 1.4 Examine the Forcing Term

Find the function:

def forcing_term(X, Y, t):

In the default configuration, the forcing is set to zero:

Q = np.zeros_like(X + Y)

Focus on the following questions:

What variable is returned by forcing_term()?
Why does np.zeros_like(X + Y) create a two-dimensional forcing field?
What would happen if the Gaussian forcing block were activated?

Step 1.5 Examine the Model Equation

Find the function:

def model_step(u, v, p, t):

This function computes one forward-Euler time step:

dudt = Y * v - ddx(p) - r * u
dvdt = -Y * u - ddy(p) - r * v
dpdt = -(ddx(u) + ddy(v)) - eps * p + Q

u_new = u + DT * dudt
v_new = v + DT * dvdt
p_new = p + DT * dpdt

Focus on the following questions:

Which terms represent the Coriolis force?
Which terms represent the pressure-gradient force?
Which terms represent convergence or divergence?
What is the role of r and eps?
Why does the solution remain zero when both the initial condition and forcing are zero?

This step provides the baseline model configuration. In the next step, you will modify the initial condition to generate free equatorial waves.

Step 2. Generate Equatorial Waves

In this step, we introduce nonzero initial conditions to generate freely propagating equatorial waves.

First, create a working copy of the template script:

cp 6_1-shallow_water_template.py 6_1-shallow_water_step2.py

Use the following configuration throughout Step 2:

DX = 0.2
DY = 0.2

DT = 0.02
NT = 600

r = 0.00
eps = 0.00

Q = np.zeros_like(X + Y)

These settings remove damping and forcing so that the model simulates freely propagating equatorial waves.

Step 2.1 Kelvin Wave

A Kelvin wave corresponds to the special case \(n=-1\), for which the meridional velocity is zero:

\[ v = 0. \]

Create a new script:

cp 6_1-shallow_water_step2.py 6_1-shallow_water_Kelvin.py

Open the script and locate the initial_condition() function. Activate the Kelvin-wave initial condition:

p = A * env_x * env_y * np.sin(phase)
u = p.copy()
v = np.zeros_like(p)

Then disable the zero-initial-condition block.

Run the model:

python 6_1-shallow_water_Kelvin.py

Examine the generated files:

6_1-shallow_water_Kelvin.mp4
6_1-shallow_water_Kelvin.png

Questions:

Does the disturbance propagate eastward or westward?
Why is \(v=0\) for a Kelvin wave?
Why are \(u\) and \(p\) in phase?
Why is the wave amplitude largest near the equator?

Step 2.2 Inertio-Gravity Wave (\(n=0\))

Next, generate an \(n=0\) inertio-gravity wave. Unlike Kelvin waves, inertio-gravity waves can propagate both eastward and westward.

Create a new script:

cp 6_1-shallow_water_step2.py 6_1-shallow_water_IG_n0.py

Replace the initial_condition() function with:

def initial_condition(X, Y):
    A = 1.5
    x0 = 0.0
    sigma_x = 4.0
    k = 0.5

    omega = np.sqrt(k**2 + 1.0)          # Eastward
    #omega = -1.0 * np.sqrt(k**2 + 1.0)  # Westward

    env_x = np.exp(-0.5 * ((X - x0) / sigma_x)**2)
    env_y = np.exp(-0.5 * Y**2)

    phase = k * (X - x0)

    psi0 = hermiteH(0, Y) * env_y
    psi1 = hermiteH(1, Y) * env_y

    u = A * env_x * psi1 * np.sin(phase)
    v = A * env_x * 2.0 * (omega + k) * psi0 * np.cos(phase)
    p = A * env_x * psi1 * np.sin(phase)

    return u, v, p

First, run the eastward-propagating case:

python 6_1-shallow_water_IG_n0.py

Then modify the sign of omega:

#omega = np.sqrt(k**2 + 1.0)          # Eastward
omega = -1.0 * np.sqrt(k**2 + 1.0)   # Westward

Run the script again and compare the two simulations.

Questions:

How does changing the sign of \(\omega\) affect the propagation direction?
Why do \(u\) and \(p\) use sin(phase) while \(v\) uses cos(phase)?
How does the inertio-gravity wave differ from the Kelvin wave?
How does the westward-propagating branch differ from the eastward-propagating branch?

Step 2.3 General Equatorial Wave Solution (\(n \ge 1\))

For \(n \ge 1\), both Rossby waves and inertio-gravity waves share the same meridional structures. The difference between the two wave families is determined by the dispersion relationship used to calculate the frequency \(\omega\).

Create a new script:

cp 6_1-shallow_water_step2.py 6_1-shallow_water_general_wave.py

Replace the initial_condition() function with:

def initial_condition(X, Y):
    A = 2.0
    x0 = 0.0
    sigma_x = 4.0

    n = 1
    k = 0.5

    # Rossby wave
    omega = k / (k**2 + 2.0*n + 1.0)

    # Eastward inertio-gravity wave
    #omega = np.sqrt(k**2 + 2.0*n + 1.0)

    # Westward inertio-gravity wave
    #omega = -np.sqrt(k**2 + 2.0*n + 1.0)

    env_x = np.exp(-0.5 * ((X - x0) / sigma_x)**2)
    env_y = np.exp(-0.5 * Y**2)

    phase = k * (X - x0)

    psi_nm1 = hermiteH(n - 1, Y) * env_y
    psi_n   = hermiteH(n,     Y) * env_y
    psi_np1 = hermiteH(n + 1, Y) * env_y

    u = A * env_x * (
        0.5 * (omega - k) * psi_np1
        + n * (omega + k) * psi_nm1
    ) * np.sin(phase)

    v = A * env_x * (
        (omega**2 - k**2) * psi_n
    ) * np.cos(phase)

    p = A * env_x * (
        0.5 * (omega - k) * psi_np1
        - n * (omega + k) * psi_nm1
    ) * np.sin(phase)

    return u, v, p

First, run the Rossby-wave case:

python 6_1-shallow_water_general_wave.py

Note

The analytical Rossby-wave dispersion relationship is often written as

\[ \omega = -\frac{k} {k^2+2n+1}, \]

which indicates westward propagation. In this exercise, however, the sign of \(\omega\) depends on the phase convention used to initialize the numerical solution. To obtain the correct westward-propagating Rossby wave in the provided code, use

\[ \omega = \frac{k} {k^2+2n+1}. \]

The propagation direction remains westward because the sign convention is absorbed into the sinusoidal phase used to construct the initial condition.

Next, replace the Rossby-wave frequency with the eastward inertio-gravity frequency:

#omega = k / (k**2 + 2.0*n + 1.0)
omega = np.sqrt(k**2 + 2.0*n + 1.0)
#omega = -np.sqrt(k**2 + 2.0*n + 1.0)

Run the model again and compare the result.

Finally, replace the frequency with the westward inertio-gravity frequency:

#omega = k / (k**2 + 2.0*n + 1.0)
#omega = np.sqrt(k**2 + 2.0*n + 1.0)
omega = -np.sqrt(k**2 + 2.0*n + 1.0)

Run the model and compare all three simulations.

After comparing the three wave families for \(n=1\), repeat the experiments using

n = 2

Questions:

How does the propagation direction differ between Rossby and inertio-gravity waves?
How do the eastward and westward inertio-gravity waves differ?
How does the propagation speed compare among the three wave types?
How does the meridional structure change when \(n\) increases from 1 to 2?
Do the disturbances remain trapped near the equator?
Why do Rossby waves propagate much more slowly than inertio-gravity waves?

Step 2 Summary

After completing Step 2, compare the simulations:

Script	Wave Type	Mode Number	Frequency
`2_3-shallow_water_Kelvin.py`	Kelvin wave	\(n=-1\)	\(\omega=k\)
`2_3-shallow_water_IG_n0.py`	Inertio-gravity wave	\(n=0\)	\(\omega=\pm\sqrt{k^2+1}\)
`2_3-shallow_water_general_wave.py`	Rossby or inertio-gravity wave	\(n \ge 1\)	depends on selected \(\omega\)

Focus on the following comparisons:

propagation direction
propagation speed
meridional structure
phase relationship between \(u\), \(v\), and \(p\)
degree of equatorial trapping

Key Points

The equatorial beta approximation describes how the Coriolis parameter varies near the equator.
The equatorial shallow-water equations provide a simple framework for studying tropical climate dynamics.
Free equatorial waves include Kelvin waves, inertio-gravity waves, and Rossby waves.
The meridional structures of equatorial waves are described by Hermite functions and are trapped near the equator by a Gaussian envelope.
Different wave types are distinguished by their dispersion relationships and propagation characteristics.
Numerical integration of the shallow-water equations provides a useful tool for visualizing tropical wave dynamics.

Homework

Complete the following tasks:

Generate one example of a free equatorial wave simulation.
Choose one wave type: Kelvin, inertio-gravity, or Rossby wave.
Submit the corresponding MP4 animation file.
Provide a brief discussion, approximately one paragraph, describing: - the propagation direction, - the meridional structure, - the phase relationship among \(u\), \(v\), and \(p\).

Optional Further Experiments

Students interested in further exploration may investigate one or more of the following:

Test different zonal wavenumbers (\(k\)) and examine changes in phase speed and wavelength.
Investigate higher-order meridional modes (\(n \ge 2\)).
Compare eastward and westward inertio-gravity waves.
Compare Rossby and inertio-gravity waves with the same meridional mode number.
Examine how grid spacing and time step affect the numerical solution.

Discussion Questions

Why are equatorial waves trapped near the equator?
What physical processes distinguish Kelvin, inertio-gravity, and Rossby waves?
Why does the Kelvin wave have no meridional velocity?
Why are Rossby waves slower than inertio-gravity waves?
What aspects of the real tropical atmosphere are not represented by this idealized shallow-water model?