Schär mountain wave with terrain-following coordinates

This example simulates the classic Schär mountain wave test case using the fully compressible equations on a terrain-following grid. The test case features a bell-shaped mountain with superimposed small-scale corrugations, generating both propagating and evanescent wave modes.

The Schär mountain wave (Schär et al. (2002)) is a stringent benchmark for terrain-following coordinates because the fine-scale terrain corrugations create steep coordinate-surface slopes that challenge the horizontal pressure gradient computation. We deform the computational grid itself using TerrainFollowingVerticalDiscretization, with a TwoLevelDecay coordinate following Schär et al. (2002). The acoustic modes are integrated with the split-explicit substepper.

References

Physical setup

The simulation initializes an atmosphere with constant Brunt–Väisälä frequency $N$ and uniform background wind $U$. A bell-shaped mountain with superimposed oscillations triggers both propagating gravity waves (for wavenumbers below the critical value $k^*$) and evanescent waves (above $k^*$).

Mountain profile

The terrain follows the Schär mountain profile (Equation 46 by Schär et al. (2002)):

\[h(x) = h_0 \exp\left(-\frac{x^2}{a^2}\right) \cos^2\left(\frac{\pi x}{\lambda}\right)\]

where $h_0$ is the peak height, $a$ is the Gaussian half-width, and $\lambda$ is the wavelength of the terrain corrugations.

Constant stratification base state

For a reference temperature $T_0$, the background atmosphere corresponds to a constant Brunt–Väisälä frequency $N$:

\[N^2 = \frac{g^2}{c_p^d T_0}\]

The density-scale height parameter is:

\[\beta = \frac{g}{R^d T_0}\]

The potential temperature profile that maintains constant $N$ is:

\[\theta(z) = \theta_0 \exp\left(\frac{N^2 z}{g}\right)\]

Linear wave theory

For the linearized mountain wave problem, vertical wavenumber $m$ satisfies the dispersion relation (Appendix A of Klemp et al. (2015)):

\[m^2 = \frac{N^2}{U^2} - \frac{\beta^2}{4} - k^2\]

Waves propagate vertically when $m^2 > 0$ (i.e., $k < k^*$), and decay exponentially when $m^2 < 0$ (i.e., $k > k^*$), where the critical wavenumber is:

\[k^* = \sqrt{\frac{N^2}{U^2} - \frac{\beta^2}{4}}\]

using Breezeusing Oceananigansusing Breeze.TerrainFollowingDiscretization: SlopeInsideInterpolation,                                              TerrainFollowingVerticalDiscretization,                                              TwoLevelDecay,                                              materialize_terrain!using Oceananigans.Grids: xnodes, znodeusing Oceananigans.Unitsusing Oceananigans: Face, Centerusing Breeze.Thermodynamics: dry_air_gas_constantusing Printfusing CairoMakieusing CUDAparse_env(::Type{T}, name, default) where T =    haskey(ENV, name) ? parse(T, ENV[name]) : default

parse_env (generic function with 1 method)

Thermodynamic parameters

We define the base state with surface pressure $p_0 = 1000 \, {\rm hPa}$ and reference temperature $T_0 = 300 \, {\rm K}$. We also set the background wind at $U = 20 \, {\rm m/s}$:

constants = ThermodynamicConstants(Float64)g = constants.gravitational_accelerationcᵖᵈ = constants.dry_air.heat_capacityRᵈ = dry_air_gas_constant(constants)p₀ = 100000                 # Pa - surface pressureT₀ = 300                    # K - reference temperatureθ₀ = T₀                     # K - reference potential temperatureU  = 20                     # m s⁻¹ - uniform background wind

20

Derived atmospheric parameters for an isothermal base state at $T_0$. Note: the standard Schär et al. (2002) parameters use $N = 0.01 \, {\rm s}^{-1}$ (set N² = 1e-4 to match that paper exactly).

N² = g^2 / (cᵖᵈ * T₀)       # s⁻² - Brunt–Väisälä frequency squared (≈ 3.2e-4)N  = sqrt(N²)               # s⁻¹ - Brunt–Väisälä frequencyβ  = g / (Rᵈ * T₀)          # m⁻¹ - density scale parameter

0.0001139362871088201

Schär mountain parameters

The mountain profile with the parameters used by Schär et al. (2002) is:

h₀ = 250          # m - peak mountain height (use 25 m for strict linearity)a  = 5000         # m - Gaussian half-width parameterλ  = 4000         # m - terrain corrugation wavelengthK  = 2 * pi / λ   # rad m⁻¹ - terrain wavenumberhill(x) = h₀ * exp(-(x / a)^2) * cos(pi * x / λ)^2

hill (generic function with 1 method)

Grid setup

The domain spans 200 km horizontally (±100 km) and 20 km vertically. We use a TerrainFollowingVerticalDiscretization with uniform spacing in the computational coordinate; the terrain-following transformation deforms these surfaces to follow the mountain, so no vertical grid stretching near the surface is needed. The defaults below are a quick preview; for the standard corrugated Schär validation, use BREEZE_SCHAR_NX=400 BREEZE_SCHAR_NZ=200 so the 4 km terrain wavelength is resolved by roughly eight cells.

arch_name = lowercase(get(ENV, "BREEZE_SCHAR_ARCH", CUDA.functional() ? "gpu" : "cpu"))arch = arch_name == "gpu" ? GPU() : CPU()Nx = parse_env(Int, "BREEZE_SCHAR_NX", 201)Nz = parse_env(Int, "BREEZE_SCHAR_NZ", 50)domain_width = parse_env(Float64, "BREEZE_SCHAR_L", 100kilometers)const domain_height = parse_env(Float64, "BREEZE_SCHAR_H", 20kilometers)z_faces = TerrainFollowingVerticalDiscretization(collect(range(0, domain_height, length=Nz+1));                                                 formulation = TwoLevelDecay(large_scale_height = domain_height / 2,                                                                     small_scale_height = domain_height / 8))grid = RectilinearGrid(arch, size = (Nx, Nz),                       halo = (5, 5),                       x = (-domain_width/2, domain_width/2), z = z_faces,                       topology = (Periodic, Flat, Bounded))

201×1×50 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CUDAGPU with 5×0×5 halo
├── Periodic x ∈ [-50000.0, 50000.0) regularly spaced with Δx=497.512
├── Flat y                           
└── Bounded  z ∈ [0.0, 20000.0]      variably spaced with terrain-following min(Δr)=400.0, max(Δr)=400.0

Terrain

Apply the Schär mountain profile to the grid. This materializes the terrain components inside the grid's TerrainFollowingVerticalDiscretization. The pressure-gradient stencil (TerrainMetrics) is built automatically by CompressibleDynamics from the grid; the default is SlopeOutsideInterpolation and is overridden below with slope_stencil = SlopeInsideInterpolation().

materialize_terrain!(grid, hill)

201×1×50 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CUDAGPU with 5×0×5 halo
├── Periodic x ∈ [-50000.0, 50000.0) regularly spaced with Δx=497.512
├── Flat y                           
└── Bounded  z ∈ [0.0, 20000.0]      variably spaced with terrain-following min(Δr)=400.0, max(Δr)=400.0

Plot: Terrain-following coordinate surfaces

Visualize how the computational grid deforms to follow the corrugated terrain. The terrain-following coordinate exactly captures the mountain profile regardless of vertical resolution.

x_grid = xnodes(grid, Center())fig = Figure(size=(900, 400))ax = Axis(fig[1, 1],          xlabel = "x (m)",          ylabel = "Height (m)",          title = "Terrain-following grid near the Schär mountain")# Plot coordinate surfaces at selected vertical levelscoordinate_surface_indices = unique(round.(Int, range(1, Nz + 1, length = 7)))CUDA.@allowscalar for k in coordinate_surface_indices    z_surface = [znode(i, 1, k, grid, Center(), Center(), Face()) for i in 1:Nx]    lines!(ax, x_grid, z_surface, color = :gray, linewidth = 0.5)end# The bottom coordinate surface (terrain)z_bottom = CUDA.@allowscalar [znode(i, 1, 1, grid, Center(), Center(), Face()) for i in 1:Nx]lines!(ax, x_grid, z_bottom, color = :brown, linewidth = 2, label = "Terrain surface")band!(ax, x_grid, zero(z_bottom), z_bottom, color = (:brown, 0.2))# Analytical terrain for comparisonx_fine = range(-domain_width/6, domain_width/6, length=2000)lines!(ax, collect(x_fine), [hill(x) for x in x_fine],       color = :black, linestyle = :dash, linewidth = 1, label = "Analytical h(x)")axislegend(ax, position = :rt)xlims!(ax, -domain_width/6, domain_width/6)ylims!(ax, -100, 3500)fig

Potential temperature profile

The potential temperature profile that maintains constant Brunt–Väisälä frequency is:

\[\theta(z) = \theta_0 \exp\left(\frac{N^2 z}{g}\right)\]

The profile depends only on height. We give it a two-argument method as well so the same function works both as reference_potential_temperature (called as θ(z) by the reference-state builder) and in set! (called as θ(x, z) on this Flat-in-y grid).

potential_temperature_profile(z) = θ₀ * exp(N² * z / g)potential_temperature_profile(x, z) = potential_temperature_profile(z)

potential_temperature_profile (generic function with 2 methods)

Rayleigh damping layer

A sponge layer near the domain top absorbs upward-propagating waves and prevents spurious reflections from the rigid lid. The split-explicit acoustic substepper applies this sponge as part of the acoustic update.

const sponge_depth = domain_height / 4const sponge_damping_rate = 0.1

0.1

Model construction

Build a compressible model with split-explicit acoustic substepping, 9th-order WENO advection, and terrain corrections. On a TerrainFollowingVerticalDiscretization grid, CompressibleDynamics automatically activates the terrain-following physics — contravariant vertical velocity, corrected pressure gradient, terrain-aware divergence — and builds the pressure-gradient stencil with the slope_stencil kwarg. The reference_potential_temperature enables a perturbation pressure approach for the horizontal pressure gradient that reduces the truncation error inherent in terrain-following coordinates (Klemp (2011)).

time_discretization = SplitExplicitTimeDiscretization(acoustic_cfl = 0.5,                                                      sponge = UpperSponge(damping_rate = sponge_damping_rate,                                                                           depth = sponge_depth))dynamics = CompressibleDynamics(time_discretization;                                slope_stencil = SlopeInsideInterpolation(),                                surface_pressure = p₀,                                reference_potential_temperature = potential_temperature_profile)model = AtmosphereModel(grid; dynamics, advection = WENO(order=9))

AtmosphereModel{GPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 201×1×50 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CUDAGPU with 5×0×5 halo
├── dynamics: CompressibleDynamics{SplitExplicitTimeDiscretization}
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float64}
├── timestepper: AcousticRungeKutta3
├── advection scheme: 
│   ├── momentum: WENO{5, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│   ├── ρθ: WENO{5, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
│   └── ρqᵛ: WENO{5, Float64, Oceananigans.Utils.BackendOptimizedDivision}(order=9)
├── forcing: @NamedTuple{ρᵈ::Returns{Float64}, ρu::Returns{Float64}, ρv::Returns{Float64}, ρw::Returns{Float64}, ρθ::Returns{Float64}, ρqᵛ::Returns{Float64}, ρe::Returns{Float64}}
├── tracers: ()
├── coriolis: Nothing
└── microphysics: Nothing

Initial conditions

We initialize the atmosphere in discrete hydrostatic balance using Exner function integration. This is essential for compressible models on terrain-following grids: the equation of state alone does not produce a pressure field in discrete hydrostatic balance. Because we passed reference_potential_temperature to CompressibleDynamics, the model has already computed this discrete reference state for us in reference_state.density!

set!(model,     ρ = model.dynamics.reference_state.density,     θ = potential_temperature_profile,     u = U,     v = 0,     w = 0,     enforce_mass_conservation = false)

Note: we set w = 0, yet the flow must follow the terrain at the surface (w = u ∂h/∂x). On a terrain-following grid ρw carries a kinematic bottom boundary condition that enforces this automatically — update_state! fills it in, so the contravariant w̃ vanishes at the ground from the first step. No manual bottom-face initialization is required.

Oceananigans.TimeSteppers.update_state!(model)

Simulation

The split-explicit scheme uses an advective outer time step; the acoustic waves are handled by inner substeps.

horizontal_spacing = domain_width / Nxtime_step = parse_env(Float64, "BREEZE_SCHAR_DT", 0.5 * horizontal_spacing / U)stop_time = parse_env(Float64, "BREEZE_SCHAR_STOP_TIME", 2hours)simulation = Simulation(model; Δt = time_step, stop_time = stop_time)Oceananigans.Diagnostics.erroring_NaNChecker!(simulation)

Progress callback to monitor simulation health:

wall_clock = Ref(time_ns())function progress(sim)    elapsed = 1e-9 * (time_ns() - wall_clock[])    u, v, w = sim.model.velocities    msg = @sprintf("Iter: %d, time: %s, wall time: %s, max|u|: %.2f, max|w|: %.4f m/s",                   iteration(sim), prettytime(sim), prettytime(elapsed),                   maximum(abs, u), maximum(abs, w))    wall_clock[] = time_ns()    @info msg    return nothingendadd_callback!(simulation, progress, name=:progress, IterationInterval(500))

Output

Save vertical velocity and contravariant vertical velocity for post-processing:

w = model.velocities.wcontravariant_w = model.dynamics.contravariant_vertical_velocityfilename = get(ENV, "BREEZE_SCHAR_OUTPUT", "mountain_waves")output_directory = dirname(filename)output_directory == "." || mkpath(output_directory)outputs = (; w, contravariant_w)simulation.output_writers[:jld2] = JLD2Writer(model, outputs;                                              filename,                                              schedule = TimeInterval(2minutes),                                              overwrite_existing = true)run!(simulation)

[ Info: Initializing simulation...
[ Info: Iter: 0, time: 0 seconds, wall time: 11.844 seconds, max|u|: 20.00, max|w|: 3.5611 m/s
[ Info:     ... simulation initialization complete (15.108 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (1.268 minutes).
[ Info: Iter: 500, time: 1.667 hours, wall time: 1.768 minutes, max|u|: 26.99, max|w|: 3.8023 m/s
[ Info: Simulation is stopping after running for 1.991 minutes.
[ Info: Simulation time 2 hours equals or exceeds stop time 2 hours.

Analytical solution

The linear analytical solution for mountain waves provides a validation benchmark. Following Appendix A of Klemp et al. (2015), the vertical velocity field is computed via Fourier integration over wavenumber space.

Fourier transform of terrain

The Fourier transform of the Schär mountain profile (Equation A8):

\[\hat{h}(k) = \frac{\sqrt{\pi} h_0 a}{4} \left[ e^{-a^2(K+k)^2/4} + e^{-a^2(K-k)^2/4} + 2e^{-a^2 k^2/4} \right]\]

ĥ(k) = sqrt(pi) * h₀ * a / 4 * (exp(-a^2 * (K + k)^2 / 4) +                                exp(-a^2 * (K - k)^2 / 4) +                                2 * exp(-a^2 * k^2 / 4))

ĥ (generic function with 1 method)

Dispersion relation

Vertical wavenumber squared (Equation A5) and critical wavenumber (Equation A11) by Klemp et al. (2015):

m²(k) = (N² / U^2 - β^2 / 4) - k^2k★ = sqrt(N² / U^2 - β^2 / 4)

0.00089147756144519

Linear vertical velocity

Compute the analytical vertical velocity $w(x, z)$ via Equation A10 by Klemp et al. (2015):

\[w(x, z) = -\frac{U}{\pi} e^{\beta z/2} \left[ \int_0^{k^*} k \hat{h}(k) \sin(m z + k x) \, \mathrm{d}k + \int_{k^*}^{\infty} k \hat{h}(k) e^{-|m| z} \sin(k x) \, \mathrm{d}k \right]\]

Above, the first integral represents propagating waves and the second represents evanescent waves.

"""    w_linear(x, z; nk=100)Compute the 2-D linear vertical velocity `w(x,z)` from the analytical solution(Appendix A, Equation A10 of Klemp et al., 2015)."""function w_linear(x, z; nk=100)    k = range(0, 10 * k★; length=nk)    m_abs = @. sqrt(abs(m²(k)))    integrand = @. k * ĥ(k) * ifelse(m²(k) >= 0,                                   sin(m_abs * z + k * x),                                   exp(-m_abs * z) * sin(k * x))    # Numerical integration using trapezoidal rule:    wavenumber_spacing = step(k)    integral = wavenumber_spacing * (sum(integrand) - (first(integrand) + last(integrand)) / 2)    return -(U / pi) * exp(β * z / 2) * integralend

Results: Comparison with analytical solution

We compare the simulated vertical velocity field at the final time with the linear analytical solution. The terrain-following grid exactly represents the corrugated mountain profile, avoiding the staircase artifacts of the immersed boundary method. The heatmaps display fields in physical $(x, z)$ coordinates, with the deformed grid automatically mapped to the true terrain geometry.

First, we compute analytical vertical velocity $w$ on the same grid as the simulation.

w_analytical = Field{Center, Nothing, Face}(grid)set!(w_analytical, w_linear)fig = Figure(size=(900, 800), fontsize=14)ax1 = Axis(fig[1, 1], ylabel = "z (m)", title = "Simulated w at final time")ax2 = Axis(fig[2, 1], xlabel = "x (m)", ylabel = "z (m)", title = "Linear Analytical w")hidexdecorations!(ax1, grid = false)hm1 = heatmap!(ax1, w, colormap = :balance, colorrange = (-1, 1))hm2 = heatmap!(ax2, w_analytical, colormap = :balance, colorrange = (-1, 1))Colorbar(fig[1:2, 2], hm1, label = "w (m s⁻¹)")for ax in (ax1, ax2)    ax.limits = ((-30e3, 30e3), (0, 10e3))endfig

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfoversioninfo()

Julia Version 1.12.6
Commit 15346901f00 (2026-04-09 19:20 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_PKG_SERVER_REGISTRY_PREFERENCE = eager
  JULIA_VERSION = 1.12.6
  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 PkgPkg.status()

Status `/__w/Breeze.jl/Breeze.jl/docs/Project.toml`
⌃ [86bc3604] AtmosphericProfilesLibrary v0.1.8
  [660aa2fb] Breeze v0.8.0 `.`
  [052768ef] CUDA v6.2.1
  [13f3f980] CairoMakie v0.15.13
⌅ [6a9e3e04] CloudMicrophysics v0.37.1
  [e30172f5] Documenter v1.17.0
  [daee34ce] DocumenterCitations v1.4.1
  [b6400b83] DocumenterCodeBlocks v1.1.0
  [7da242da] Enzyme v0.13.198
  [46192b85] GPUArraysCore v0.2.0
  [63c18a36] KernelAbstractions v0.9.42
  [98b081ad] Literate v2.21.0
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.110.14
⌅ [a01a1ee8] RRTMGP v0.21.9
  [3c362404] Reactant v0.2.278
  [b77e0a4c] InteractiveUtils v1.11.0
  [44cfe95a] Pkg v1.12.1
  [9a3f8284] Random v1.11.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.