Implement a PotentialTemperatureThermodynamics formulation

docs/make.jl

@@ -18,6 +18,7 @@ mkpath(literated_dir)
 
 example_scripts = [
     "thermal_bubble.jl",
+    "moist_thermal_bubble.jl",
     "cloudy_kelvin_helmholtz.jl",
     # "prescribed_sst.jl", # this is a WIP
 ]
@@ -31,6 +32,7 @@ end
 
 example_pages = Any[
     "Thermal bubble" => "literated/thermal_bubble.md",
+    "Moist and dry thermal bubbles" => "literated/moist_thermal_bubble.md",
     "Cloudy Kelvin-Helmholtz instability" => "literated/cloudy_kelvin_helmholtz.md",
     # "Prescribed SST" => "literated/prescribed_sst.md",
 ]

docs/src/appendix/notation.md

@@ -17,7 +17,7 @@ A few notes about the following table:
   a "reference state", which is an adiabatic, hydrostatic solution to the equations of motion. But there is also an
   "energy reference temperature" and "reference latent heat", which are thermodynamic constants required to define
   the internal energy of moist atmospheric constituents.
-* Mapping to AM fields: `ρe` corresponds to `model.energy_density`, `ρqᵗ` to `model.moisture_density`, and `qᵗ` to `model.specific_moisture`.
+* Mapping to AM fields: `ρe` corresponds to `energy_density(model)`, `ρqᵗ` to `model.moisture_density`, and `qᵗ` to `model.specific_moisture`.
 
 The following table also uses a few conventions that suffuse the source code and which are internalized by wise developers:
 

docs/src/microphysics/mixed_phase_saturation_adjustment.md

@@ -47,7 +47,7 @@ We can compute the saturation specific humidity across the entire range from hom
 ice nucleation up to freezing using the equilibrium model above:
 
 ```@example mixed_phase
-using Breeze.Microphysics: adjustment_saturation_specific_humidity
+using Breeze.Microphysics: equilibrium_saturation_specific_humidity
 
 p = 101325.0
 qᵗ = 0.012
@@ -56,7 +56,7 @@ Tᶠ = eq.freezing_temperature
 Tʰ = eq.homogeneous_ice_nucleation_temperature
 T = range(Tʰ - 10, Tᶠ + 10; length=81) # slightly beyond the mixed-phase range
 
-qᵛ⁺ = [adjustment_saturation_specific_humidity(Tʲ, p, qᵗ, thermo, eq) for Tʲ in T]
+qᵛ⁺ = [equilibrium_saturation_specific_humidity(Tʲ, p, qᵗ, thermo, eq) for Tʲ in T]
 ```
 
 Optionally, we can visualize how `qᵛ⁺` varies with temperature:

docs/src/microphysics/warm_phase_saturation_adjustment.md

@@ -101,10 +101,10 @@ Next, we compute the saturation specific humidity for moist air with
 a carefully chosen moist air mass fraction,
 
 ```@example microphysics
-using Breeze.Microphysics: adjustment_saturation_specific_humidity, WarmPhaseEquilibrium
+using Breeze.Microphysics: equilibrium_saturation_specific_humidity, WarmPhaseEquilibrium
 
 qᵗ = 0.012   # [kg kg⁻¹] total specific humidity
-qᵛ⁺ = Breeze.Microphysics.adjustment_saturation_specific_humidity(T, p, qᵗ, thermo, WarmPhaseEquilibrium())
+qᵛ⁺ = Breeze.Microphysics.equilibrium_saturation_specific_humidity(T, p, qᵗ, thermo, WarmPhaseEquilibrium())
 ```
 
 There are two facts of note. First is that we have identified a situation in which ``qᵗ > qᵛ⁺``,
@@ -164,7 +164,7 @@ To generate a second guess for the secant solver, we start by estimating
 the liquid mass fraction using the guess ``T = T₁``,
 
 ```@example microphysics
-qᵛ⁺₂ = adjustment_saturation_specific_humidity(T₁, p, qᵗ, thermo, WarmPhaseEquilibrium())
+qᵛ⁺₂ = equilibrium_saturation_specific_humidity(T₁, p, qᵗ, thermo, WarmPhaseEquilibrium())
 qˡ₁ = qᵗ - qᵛ⁺₂
 ```
 
@@ -189,7 +189,7 @@ using CairoMakie
 equilibrium = WarmPhaseEquilibrium()
 T = 230:0.5:320
 r = [saturation_adjustment_residual(Tʲ, 𝒰, thermo, equilibrium) for Tʲ in T]
-qᵛ⁺ = [adjustment_saturation_specific_humidity(Tʲ, p, qᵗ, thermo, equilibrium) for Tʲ in T]
+qᵛ⁺ = [equilibrium_saturation_specific_humidity(Tʲ, p, qᵗ, thermo, equilibrium) for Tʲ in T]
 
 fig = Figure()
 axr = Axis(fig[1, 1], xlabel="Temperature (K)", ylabel="Saturation adjustment residual (K)")
@@ -228,7 +228,7 @@ for (i, qᵗⁱ) in enumerate(qᵗ)
     q = MoistureMassFractions(qᵗⁱ)
     𝒰 = StaticEnergyState(e₀, q, z, p)
     T[i] = compute_temperature(𝒰, microphysics, thermo)
-    qᵛ⁺ = Breeze.Microphysics.adjustment_saturation_specific_humidity(T[i], p, qᵗⁱ, thermo, WarmPhaseEquilibrium())
+    qᵛ⁺ = Breeze.Microphysics.equilibrium_saturation_specific_humidity(T[i], p, qᵗⁱ, thermo, WarmPhaseEquilibrium())
     qˡ[i] = max(0, qᵗⁱ - qᵛ⁺)
 end
 
@@ -295,7 +295,7 @@ for k = 1:grid.Nz
     T[k] = compute_temperature(𝒰, microphysics, thermo)
 
     # Saturation specific humidity via adjustment formula using T[k], pᵣ, and qᵗ
-    qᵛ⁺[k] = Breeze.Microphysics.adjustment_saturation_specific_humidity(T[k], pᵣ, qᵗ, thermo, WarmPhaseEquilibrium())
+    qᵛ⁺[k] = Breeze.Microphysics.equilibrium_saturation_specific_humidity(T[k], pᵣ, qᵗ, thermo, WarmPhaseEquilibrium())
     qˡ[k] = max(0, qᵗ - qᵛ⁺[k])
     rh[k] = 100 * min(qᵗ, qᵛ⁺[k]) / qᵛ⁺[k]
 end

examples/anelastic_bomex.jl

@@ -32,7 +32,8 @@ u_bomex = AtmosphericProfilesLibrary.Bomex_u(FT)
 p₀, θ₀ = 101500, 299.1
 constants = ThermodynamicConstants()
 reference_state = ReferenceState(grid, constants, base_pressure=p₀, potential_temperature=θ₀)
-formulation = AnelasticFormulation(reference_state)
+formulation = AnelasticFormulation(reference_state, thermodynamics=:LiquidIcePotentialTemperature)
+# formulation = AnelasticFormulation(reference_state, thermodynamics=:StaticEnergy)
 
 q₀ = Breeze.Thermodynamics.MoistureMassFractions{eltype(grid)} |> zero
 ρ₀ = Breeze.Thermodynamics.density(p₀, θ₀, q₀, constants)
@@ -41,7 +42,9 @@ Lˡ = constants.liquid.reference_latent_heat
 w′T′, w′q′ = 8e-3, 5.2e-5
 Q = ρ₀ * cᵖᵈ * w′T′ 
 F = ρ₀ * w′q′
+
 ρe_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(Q))
+ρθ_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(ρ₀ * w′T′))
 ρqᵗ_bcs = FieldBoundaryConditions(bottom=FluxBoundaryCondition(F))
 
 u★ = 0.28 # m/s
@@ -69,6 +72,11 @@ end
     return @inbounds - p.ρᵣ[i, j, k] * w_dz_E
 end
 
+@inline function Fρθ_subsidence(i, j, k, grid, clock, fields, p)
+    w_dz_T = ℑzᵃᵃᶜ(i, j, k, grid, w_dz_ϕ, p.wˢ, p.θ_avg)
+    return @inbounds - p.ρᵣ[i, j, k] * w_dz_T
+end
+
 @inline function Fρqᵗ_subsidence(i, j, k, grid, clock, fields, p)
     w_dz_Qᵗ = ℑzᵃᵃᶜ(i, j, k, grid, w_dz_ϕ, p.wˢ, p.qᵗ_avg)
     return @inbounds - p.ρᵣ[i, j, k] * w_dz_Qᵗ
@@ -78,6 +86,7 @@ end
 u_avg_f = Field{Nothing, Nothing, Center}(grid)
 v_avg_f = Field{Nothing, Nothing, Center}(grid)
 e_avg_f = Field{Nothing, Nothing, Center}(grid)
+θ_avg_f = Field{Nothing, Nothing, Center}(grid)
 qᵗ_avg_f = Field{Nothing, Nothing, Center}(grid)
 
 wˢ = Field{Nothing, Nothing, Face}(grid)
@@ -88,17 +97,24 @@ set!(wˢ, z -> w_bomex(z))
 ρu_subsidence_forcing  = Forcing(Fρu_subsidence,  discrete_form=true, parameters=(; u_avg=u_avg_f, wˢ, ρᵣ))
 ρv_subsidence_forcing  = Forcing(Fρv_subsidence,  discrete_form=true, parameters=(; v_avg=v_avg_f, wˢ, ρᵣ))
 ρe_subsidence_forcing  = Forcing(Fρe_subsidence,  discrete_form=true, parameters=(; e_avg=e_avg_f, wˢ, ρᵣ))
+ρθ_subsidence_forcing  = Forcing(Fρθ_subsidence,  discrete_form=true, parameters=(; θ_avg=θ_avg_f, wˢ, ρᵣ))
 ρqᵗ_subsidence_forcing = Forcing(Fρqᵗ_subsidence, discrete_form=true, parameters=(; qᵗ_avg=qᵗ_avg_f, wˢ, ρᵣ))
 
 coriolis = FPlane(f=3.76e-5)
-ρuᵍ = Field{Nothing, Nothing, Center}(grid)
-ρvᵍ = Field{Nothing, Nothing, Center}(grid)
+# ρuᵍ = Field{Nothing, Nothing, Center}(grid)
+# ρvᵍ = Field{Nothing, Nothing, Center}(grid)
+uᵍ = Field{Nothing, Nothing, Center}(grid)
+vᵍ = Field{Nothing, Nothing, Center}(grid)
 uᵍ_bomex = AtmosphericProfilesLibrary.Bomex_geostrophic_u(FT)
 vᵍ_bomex = AtmosphericProfilesLibrary.Bomex_geostrophic_v(FT)
-set!(ρuᵍ, z -> uᵍ_bomex(z))
-set!(ρvᵍ, z -> vᵍ_bomex(z))
-set!(ρuᵍ, ρᵣ * ρuᵍ)
-set!(ρvᵍ, ρᵣ * ρvᵍ)
+set!(uᵍ, z -> uᵍ_bomex(z))
+set!(vᵍ, z -> vᵍ_bomex(z))
+
+ρuᵍ = Field(ρᵣ * uᵍ)
+ρvᵍ = Field(ρᵣ * vᵍ)
+
+# set!(ρuᵍ, ρᵣ * ρuᵍ)
+# set!(ρvᵍ, ρᵣ * ρvᵍ)
 
 @inline Fρu_geostrophic(i, j, k, grid, clock, fields, p) = @inbounds - p.f * p.ρvᵍ[i, j, k]
 @inline Fρv_geostrophic(i, j, k, grid, clock, fields, p) = @inbounds p.f * p.ρuᵍ[i, j, k]
@@ -117,13 +133,19 @@ set!(drying, ρᵣ * drying)
 ρqᵗ_forcing = (ρqᵗ_drying_forcing, ρqᵗ_subsidence_forcing)
 
 Fρe_field = Field{Nothing, Nothing, Center}(grid)
+Fρθ_field = Field{Nothing, Nothing, Center}(grid)
 dTdt_bomex = AtmosphericProfilesLibrary.Bomex_dTdt(FT)
 set!(Fρe_field, z -> dTdt_bomex(1, z))
 set!(Fρe_field, ρᵣ * cᵖᵈ * Fρe_field)
+set!(Fρθ_field, z -> dTdt_bomex(1, z))
+set!(Fρθ_field, ρᵣ * Fρe_field)
 
 ρe_radiation_forcing = Forcing(Fρe_field)
 ρe_forcing = (ρe_radiation_forcing, ρe_subsidence_forcing)
 
+ρθ_radiation_forcing = Forcing(Fρθ_field)
+ρθ_forcing = (ρθ_radiation_forcing, ρθ_subsidence_forcing)
+
 fig = Figure()
 axe = Axis(fig[1, 1], xlabel="z (m)", ylabel="Fρe (K/s)")
 axq = Axis(fig[1, 2], xlabel="z (m)", ylabel="Fρqᵗ (1/s)")
@@ -144,8 +166,10 @@ closure = AnisotropicMinimumDissipation()
 model = AtmosphereModel(grid; formulation, coriolis, microphysics, closure,
                         scalar_advection = Centered(order=2),
                         momentum_advection = WENO(order=9),
-                        forcing = (ρqᵗ=ρqᵗ_forcing, ρu=ρu_forcing, ρv=ρv_forcing, ρe=ρe_forcing),
-                        boundary_conditions = (ρe=ρe_bcs, ρqᵗ=ρqᵗ_bcs, ρu=ρu_bcs, ρv=ρv_bcs))
+                        forcing = (ρqᵗ=ρqᵗ_forcing, ρu=ρu_forcing, ρv=ρv_forcing, ρθ=ρθ_forcing),
+                        boundary_conditions = (ρθ=ρθ_bcs, ρqᵗ=ρqᵗ_bcs, ρu=ρu_bcs, ρv=ρv_bcs))
+                        # forcing = (ρqᵗ=ρqᵗ_forcing, ρu=ρu_forcing, ρv=ρv_forcing, ρe=ρe_forcing),
+                        # boundary_conditions = (ρe=ρe_bcs, ρqᵗ=ρqᵗ_bcs, ρu=ρu_bcs, ρv=ρv_bcs))
 
 # Values for the initial perturbations can be found in Appendix B
 # of Siebesma et al 2003, 3rd paragraph
@@ -159,32 +183,34 @@ simulation = Simulation(model; Δt=10, stop_time)
 conjure_time_step_wizard!(simulation, cfl=0.7)
 
 # Write a callback to compute *_avg_f
-ρe = energy_density(model)
+e = static_energy(model)
+θ = liquid_ice_potential_temperature(model)
 u_avg = Field(Average(model.velocities.u, dims=(1, 2)))
 v_avg = Field(Average(model.velocities.v, dims=(1, 2)))
-e_avg = Field(Average(ρe / ρᵣ, dims=(1, 2)))
+e_avg = Field(Average(e, dims=(1, 2)))
+θ_avg = Field(Average(θ, dims=(1, 2)))
 qᵗ_avg = Field(Average(model.specific_moisture, dims=(1, 2)))
 
 function compute_averages!(sim)
     compute!(u_avg)
     compute!(v_avg)
     compute!(e_avg)
+    compute!(θ_avg)
     compute!(qᵗ_avg)
     parent(u_avg_f) .= parent(u_avg)
     parent(v_avg_f) .= parent(v_avg)
     parent(e_avg_f) .= parent(e_avg)
+    parent(θ_avg_f) .= parent(θ_avg)
     parent(qᵗ_avg_f) .= parent(qᵗ_avg)
     return nothing
 end
 
 add_callback!(simulation, compute_averages!)
 
-θ = Breeze.AtmosphereModels.PotentialTemperatureField(model)
 qˡ = model.microphysical_fields.qˡ
 qᵛ = model.microphysical_fields.qᵛ
 qᵗ = model.specific_moisture
-qᵛ⁺ = Breeze.AtmosphereModels.SaturationSpecificHumidityField(model)
-θ_avg = Average(θ, dims=(1, 2)) |> Field
+qᵛ⁺ = SaturationSpecificHumidityField(model)
 qˡ_avg = Average(qˡ, dims=(1, 2)) |> Field
 
 fig = Figure()
@@ -212,7 +238,8 @@ ax_ρe = Axis(fig_lower[1, 1], xlabel="Energy density", ylabel="z (m)")
 ax_e  = Axis(fig_lower[1, 2], xlabel="Specific energy", ylabel="z (m)")
 ax_θ  = Axis(fig_lower[1, 3], xlabel="Potential temperature (K)", ylabel="z (m)")
 
-e = specific_energy(model)
+e = static_energy(model)
+ρe = ρᵣ * e
 ρe_avg = Average(ρe, dims=(1, 2)) |> Field
 e_avg = Average(e, dims=(1, 2)) |> Field
 
@@ -247,7 +274,7 @@ function progress(sim)
     qᵗ = sim.model.specific_moisture
     qᵗmax = maximum(qᵗ)
 
-    ρe = energy_density(sim.model)
+    ρe = static_energy_density(sim.model)
     ρemin = minimum(ρe)
     ρemax = maximum(ρe)
 

examples/atmos_convection.jl

@@ -21,7 +21,7 @@ set!(model, θ=θᵢ, u=Ξᵢ, v=Ξᵢ)
 
 simulation = Simulation(model, Δt=10, stop_time=20minutes)
 conjure_time_step_wizard!(simulation, cfl=0.7)
-ρe = energy_density(model)
+ρe = static_energy_density(model)
 ∫ρe = Field(Integral(ρe))
 
 function progress(sim)

examples/cloudy_kelvin_helmholtz.jl

@@ -28,11 +28,8 @@ using Printf
 # Grid resolution is modest but enough to clearly resolve the Kelvin-Helmholtz billows and
 # rolled-up moisture filament.
 
-Nx = 384   # horizontal resolution
-Nz = 128   # vertical resolution
-
-Lx = 10e3  # domain length
-Lz =  3e3  # domain height
+Nx, Nz = 384, 128   # resolution
+Lx, Lz = 10e3, 3e3  # domain extent
 
 grid = RectilinearGrid(; size = (Nx, Nz), x = (0, Lx), z = (0, Lz),
                          topology = (Periodic, Flat, Bounded))
@@ -53,14 +50,14 @@ model = AtmosphereModel(grid; advection=WENO(order=5), microphysics)
 # N² = \frac{g}{θ₀} \frac{∂θ}{∂z} ,
 # ```
 #
-# We initialize the potential temperature that gives constant Brunt–Väisälä frequency,
+# We initialize with a potential temperature that gives constant Brunt–Väisälä frequency,
 # representative of mid-tropospheric stability. The (dry) Brunt–Väisälä frequency is
 #
 # ```math
 # N² = \frac{g}{θ} \frac{∂θ}{∂z}
 # ```
 #
-# and thus, for constant ``N²`` the above implies ``θ = θ₀ \exp{(N² z / g)}``.
+# and thus, for constant ``N²``, the above implies that ``θ = θ₀ \exp{(N² z / g)}``.
 
 thermo = ThermodynamicConstants()
 g = thermo.gravitational_acceleration
@@ -81,17 +78,29 @@ N = 0.01                  # target dry Brunt–Väisälä frequency (s⁻¹)
 
 # First, we set up the shear layer using a ``\tanh`` profile:
 
-z₀    = 1e3     # center of shear & moist layer (m)
-Δzᶸ   = 150     # shear layer half-thickness (m)
-U_top = 25      # upper-layer wind (m/s)
-U_bot =  5      # lower-layer wind (m/s)
-uᵇ(z) = U_bot + (U_top - U_bot) * (1 + tanh((z - z₀) / Δzᶸ)) / 2
+z₀  = 1e3  # center of shear & moist layer (m)
+Δzᵘ = 150  # shear layer half-thickness (m)
+U₀  =  5   # base wind speed (m/s)
+ΔU  = 20   # upper-layer wind (m/s)
+uᵇ(z) = U₀ + ΔU * (1 + tanh((z - z₀) / Δzᵘ)) / 2
 
 # For the moisture layer, we use a Gaussian in ``z`` centered at ``z₀``:
 
-q_max = 0.012  # peak specific humidity (kg/kg)
-Δz_q = 200     # moist layer half-width (m)
-qᵇ(z) = q_max * exp(-(z - z₀)^2 / 2Δz_q^2)
+qᵗ₀ = 0.012  # peak specific humidity (kg/kg)
+Δzᵗ = 200     # moist layer half-width (m)
+qᵇ(z) = qᵗ₀ * exp(-(z - z₀)^2 / 2Δzᵗ^2)
+
+# We initialize the model via Oceananigans `set!`, adding also a bit of random noise.
+
+δθ = 0.01
+δu = 1e-3
+δq = 0.05 * qᵗ₀
+
+θᵢ(x, z) = θᵇ(z) + δθ * rand()
+qᵗᵢ(x, z) = qᵇ(z) + δq * rand()
+uᵢ(x, z) = uᵇ(z) + δu * rand()
+
+set!(model; u=uᵢ, qᵗ=qᵗᵢ, θ=θᵢ)
 
 # ## The Kelvin-Helmholtz instability
 #
@@ -107,12 +116,11 @@ qᵇ(z) = q_max * exp(-(z - z₀)^2 / 2Δz_q^2)
 #
 # Let's plot the initial state as well as the Richardson number.
 
-z = znodes(grid, Center())
-
-dudz = @. (U_top - U_bot) * sech((z - z₀) / Δzᶸ)^2 / 2Δzᶸ
-Ri = N^2 ./ dudz.^2
+U = Field(Average(model.velocities.u, dims=(1, 2)))
+Ri = N^2 / ∂z(U)^2
 
-using CairoMakie
+Qᵗ = Field(Average(model.specific_moisture, dims=1))
+θ = Field(Average(liquid_ice_potential_temperature(model), dims=1))
 
 fig = Figure(size=(800, 500))
 
@@ -121,11 +129,12 @@ axq = Axis(fig[1, 2], xlabel = "qᵇ (kg/kg)", title="Total liquid")
 axθ = Axis(fig[1, 3], xlabel = "θᵇ (K)", title="Potential temperature")
 axR = Axis(fig[1, 4], xlabel = "Ri", ylabel="z (m)", title="Richardson number")
 
-lines!(axu, uᵇ.(z), z)
-lines!(axq, qᵇ.(z), z)
-lines!(axθ, θᵇ.(z), z)
-lines!(axR, Ri, z)
+lines!(axu, U)
+lines!(axq, Qᵗ)
+lines!(axθ, θ)
+lines!(axR, Ri)
 lines!(axR, [1/4, 1/4], [0, Lz], linestyle = :dash, color = :black)
+
 xlims!(axR, 0, 0.8)
 axR.xticks = 0:0.25:1
 
@@ -137,20 +146,6 @@ end
 
 fig
 
-# ## Define initial conditions
-#
-# We initialize the model via Oceananigans `set!`, adding also a bit of random noise.
-
-δθ = 0.01
-δu = 1e-3
-δq = 0.05 * q_max
-
-θᵢ(x, z) = θᵇ(z) + δθ * rand()
-qᵗᵢ(x, z) = qᵇ(z) + δq * rand()
-uᵢ(x, z) = uᵇ(z) + δu * rand()
-
-set!(model; u=uᵢ, qᵗ=qᵗᵢ, θ=θᵢ)
-
 # ## Set up and run the simulation
 #
 # We construct a simulation and use the time-step wizard to keep the CFL number under control.
@@ -176,7 +171,7 @@ add_callback!(simulation, progress, TimeInterval(1minute))
 # the potential temperatures and the specific humidities (vapour, liquid, ice).
 u, v, w = model.velocities
 ξ = ∂z(u) - ∂x(w)
-θ = PotentialTemperatureField(model)
+θ = liquid_ice_potential_temperature(model)
 outputs = merge(model.velocities, model.microphysical_fields, (; ξ, θ))
 
 filename = "wave_clouds.jld2"