Acoustic wave refraction by wind shear

This example simulates an acoustic pulse propagating through a wind shear layer using the fully compressible Euler equations. When wind speed increases with height, sound waves are refracted: waves traveling with the wind bend downward (trapped near the surface), while waves traveling against the wind bend upward.

The sound speed for a wave traveling in direction $\hat{n}$ is

\[𝕌ˢ = 𝕌ˢⁱ + \mathbf{u} \cdot \hat{n}\]

where $𝕌ˢⁱ$ is the intrinsic wave speed and $\mathbf{u}$ is the wind velocity. This causes wavefronts to tilt toward regions of lower effective sound speed.

This phenomenon explains why distant sounds are often heard more clearly downwind of a source, as sound energy is "ducted" along the surface. For more on this topic, see

Ostashev, V. E. and Wilson, D. K. (2015). Acoustics in Moving Inhomogeneous Media (CRC Press).
Pierce, A. D. (2019). Acoustics: An Introduction to Its Physical Principles and Applications (Springer Cham).

We use stable stratification to suppress Kelvin-Helmholtz instability and a logarithmic wind profile consistent with the atmospheric surface layer.

using Breeze
using Breeze.Thermodynamics: adiabatic_hydrostatic_density
using Oceananigans.Units
using Printf
using CairoMakie

Grid and model setup

Nx, Nz = 128, 64
Lx, Lz = 1000, 200  # meters

grid = RectilinearGrid(size = (Nx, Nz), x = (-Lx/2, Lx/2), z = (0, Lz),
                       topology = (Periodic, Flat, Bounded))

model = AtmosphereModel(grid; dynamics = CompressibleDynamics())

AtmosphereModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×1×64 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── dynamics: CompressibleDynamics
├── formulation: LiquidIcePotentialTemperatureFormulation
├── timestepper: SSPRungeKutta3
├── advection scheme: 
│   ├── momentum: Centered(order=2)
│   ├── ρθ: Centered(order=2)
│   └── ρqᵗ: Centered(order=2)
├── tracers: ()
├── coriolis: Nothing
└── microphysics: Nothing

Background state

We build a hydrostatically balanced reference state using ReferenceState. This provides the background density and pressure profiles.

constants = model.thermodynamic_constants

θ₀ = 300      # Reference potential temperature (K)
p₀ = 101325   # Surface pressure (Pa)

reference = ReferenceState(grid, constants; surface_pressure=p₀, potential_temperature=θ₀)

ReferenceState{Float64}(p₀=101325.0, θ₀=300.0, pˢᵗ=100000.0)

The sound speed at the surface determines the acoustic wave propagation speed.

Rᵈ = constants.molar_gas_constant / constants.dry_air.molar_mass
cᵖᵈ = constants.dry_air.heat_capacity
γ = cᵖᵈ / (cᵖᵈ - Rᵈ)
𝕌ˢⁱ = sqrt(γ * Rᵈ * θ₀)

347.15629052210863

The wind profile follows the classic log-law of the atmospheric surface layer.

U₀ = 20 # Surface velocity (m/s, u★ / κ)
ℓ = 1  # Roughness length [m] -- like, shrubs and stuff

Uᵢ(z) = U₀ * log((z + ℓ) / ℓ)

Uᵢ (generic function with 1 method)

Initial conditions

We initialize a localized Gaussian density pulse representing an acoustic disturbance. For a rightward-propagating acoustic wave, the velocity perturbation is in phase with the density perturbation: $u' = (𝕌ˢ / ρ₀) ρ'$.

δρ = 0.01         # Density perturbation amplitude (kg/m³)
σ = 20            # Pulse width (m)

gaussian(x, z) = exp(-(x^2 + z^2) / 2σ^2)
ρ₀ = interior(reference.density, 1, 1, 1)[]

ρᵢ(x, z) = adiabatic_hydrostatic_density(z, p₀, θ₀, constants) + δρ * gaussian(x, z)
uᵢ(x, z) = Uᵢ(z) #+ (𝕌ˢⁱ / ρ₀) * δρ * gaussian(x, z)

set!(model, ρ=ρᵢ, θ=θ₀, u=uᵢ)

Simulation setup

Acoustic waves travel fast ($𝕌ˢⁱ ≈ 347$ m/s), so we need a small time step. The Courant–Friedrichs–Lewy (CFL) condition is based on the effective sound speed $𝕌ˢ = 𝕌ˢⁱ + \mathrm{max}(U)$.

Δx, Δz = Lx / Nx, Lz / Nz
𝕌ˢ = 𝕌ˢⁱ + Uᵢ(Lz)
Δt = 0.5 * min(Δx, Δz) / 𝕌ˢ
stop_time = 1  # seconds

simulation = Simulation(model; Δt, stop_time)

function progress(sim)
    u, v, w = sim.model.velocities
    msg = @sprintf("Iter: %d, t: %s, max|u|: %.2f m/s, max|w|: %.2f m/s",
                   iteration(sim), prettytime(sim),
                   maximum(abs, u), maximum(abs, w))
    @info msg
end

add_callback!(simulation, progress, IterationInterval(100))

Output

We perturbation fields for density and x-velocity for visualization.

ρ = model.dynamics.density
u, v, w = model.velocities

ρᵇᵍ = CenterField(grid)
uᵇᵍ = XFaceField(grid)

set!(ρᵇᵍ, (x, z) -> adiabatic_hydrostatic_density(z, p₀, θ₀, constants))
set!(uᵇᵍ, (x, z) -> Uᵢ(z))

ρ′ = Field(ρ - ρᵇᵍ)
u′ = Field(u - uᵇᵍ)

U = Average(u, dims = 1)
R = Average(ρ, dims = 1)
W² = Average(w^2, dims = 1)

filename = "acoustic_wave.jld2"
outputs = (; ρ′, u′, w, U, R, W²)

simulation.output_writers[:jld2] = JLD2Writer(model, outputs; filename,
                                              schedule = TimeInterval(0.01),
                                              overwrite_existing = true)

run!(simulation)

[ Info: Initializing simulation...
[ Info: Iter: 0, t: 0 seconds, max|u|: 105.91 m/s, max|w|: 0.00 m/s
[ Info:     ... simulation initialization complete (22.536 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (3.794 seconds).
[ Info: Iter: 100, t: 333.448 ms, max|u|: 105.91 m/s, max|w|: 0.49 m/s
[ Info: Iter: 200, t: 666.895 ms, max|u|: 106.25 m/s, max|w|: 0.49 m/s
[ Info: Simulation is stopping after running for 29.720 seconds.
[ Info: Simulation time 1 second equals or exceeds stop time 1 second.
[ Info: Iter: 300, t: 1 second, max|u|: 106.26 m/s, max|w|: 0.32 m/s

Visualization

Load the saved perturbation fields and create a snapshot.

ρ′ts = FieldTimeSeries(filename, "ρ′")
u′ts = FieldTimeSeries(filename, "u′")
wts = FieldTimeSeries(filename, "w")
Uts = FieldTimeSeries(filename, "U")
Rts = FieldTimeSeries(filename, "R")
W²ts = FieldTimeSeries(filename, "W²")

times = ρ′ts.times
Nt = length(times)

fig = Figure(size = (900, 600), fontsize = 12)

axρ = Axis(fig[1, 2]; aspect = 5, ylabel = "z (m)")
axw = Axis(fig[2, 2]; aspect = 5, ylabel = "z (m)")
axu = Axis(fig[3, 2]; aspect = 5, xlabel = "x (m)", ylabel = "z (m)")
axR = Axis(fig[1, 1]; xlabel = "⟨ρ⟩ (kg/m³)")
axW = Axis(fig[2, 1]; xlabel = "⟨w²⟩ (m²/s²)")
axU = Axis(fig[3, 1]; xlabel = "⟨u⟩ (m/s)")

hidexdecorations!(axρ)
hidexdecorations!(axw)
colsize!(fig.layout, 1, Relative(0.1))

n = Observable(Nt)
ρ′n = @lift ρ′ts[$n]
u′n = @lift u′ts[$n]
wn = @lift wts[$n]
Un = @lift Uts[$n]
Rn = @lift Rts[$n]
W²n = @lift W²ts[$n]

ρlim = δρ / 4
ulim = 1

hmρ = heatmap!(axρ, ρ′n; colormap = :balance, colorrange = (-ρlim, ρlim))
hmw = heatmap!(axw, wn; colormap = :balance, colorrange = (-ulim, ulim))
hmu = heatmap!(axu, u′n; colormap = :balance, colorrange = (-ulim, ulim))

lines!(axU, Un)
lines!(axR, Rn)
lines!(axW, W²n)

Colorbar(fig[1, 3], hmρ; label = "ρ′ (kg/m³)")
Colorbar(fig[2, 3], hmw; label = "w (m/s)")
Colorbar(fig[3, 3], hmu; label = "u′ (m/s)")

title = @lift "Acoustic wave in log-layer shear — t = $(prettytime(times[$n]))"
fig[0, :] = Label(fig, title, fontsize = 16, tellwidth = false)

CairoMakie.record(fig, "acoustic_wave.mp4", 1:Nt, framerate = 18) do nn
    n[] = nn
end

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.12.2
Commit ca9b6662be4 (2025-11-20 16:25 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD EPYC 7R13 Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver3)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 8 virtual cores)
Environment:
  JULIA_GPG = 3673DF529D9049477F76B37566E3C7DC03D6E495
  JULIA_LOAD_PATH = :@breeze
  JULIA_VERSION = 1.12.2
  JULIA_DEPOT_PATH = /usr/local/share/julia:
  JULIA_PATH = /usr/local/julia
  JULIA_PROJECT = @breeze

These were the top-level packages installed in the environment:

import Pkg
Pkg.status()

Status `/__w/Breeze.jl/Breeze.jl/docs/Project.toml`
  [86bc3604] AtmosphericProfilesLibrary v0.1.7
  [660aa2fb] Breeze v0.2.1 `.`
  [052768ef] CUDA v5.9.6
  [13f3f980] CairoMakie v0.15.8
  [6a9e3e04] CloudMicrophysics v0.29.1
  [e30172f5] Documenter v1.16.1
  [daee34ce] DocumenterCitations v1.4.1
  [98b081ad] Literate v2.21.0
  [85f8d34a] NCDatasets v0.14.10
  [9e8cae18] Oceananigans v0.103.1
  [a01a1ee8] RRTMGP v0.21.6
  [b77e0a4c] InteractiveUtils v1.11.0
  [44cfe95a] Pkg v1.12.0
  [9a3f8284] Random v1.11.0

This page was generated using Literate.jl.