Bulk formula implementation

[yairchn]

Description

In 'momentum_bc.jl' the bulk formula is operating on 'state_int⁻' which is the first DG point above z=0. The computed flux from the formula is applied however at the lower level 'state⁻'. This flux is too large for this point and this creating problems that are mostly evident if turbulence is shut off and +1m/s horizontal wind at the surface can turn into -200m/s winds within 30min simulation.
A consistent implementation as suggested by this branch

original bulk formula code

u1⁻ = state_int⁻.ρu / state_int⁻.ρ
u_int⁻_tan = u1⁻ - dot(u1⁻, n) .* SVector(n)
normu_int⁻_tan = norm(u_int⁻_tan)
-- NOTE: difference from design docs since normal points outwards
C = bc_momentum.drag.fn(state⁻, aux⁻, t, normu_int⁻_tan)
τn = C * normu_int⁻_tan * u_int⁻_tan

Consistent implementation (in consistent with the assumptions of bulk formula though)

u1⁻ = state⁻.ρu / state⁻.ρ
u⁻_tan = u1⁻ - dot(u1⁻, n) .* SVector(n)
normu⁻_tan = norm(u⁻_tan)
# NOTE: difference from design docs since normal points outwards
C = bc_momentum.drag.fn(state⁻, aux⁻, t, normu⁻_tan)
τn = C * normu⁻_tan * u⁻_tan
fluxᵀn.ρu += state⁻.ρ * τn
fluxᵀn.energy.ρe += state⁻.ρu' * τn

See also

ClimateMachine.jl/src/Atmos/Model/bc_momentum.jl

Line 160 in 6c11bec

# original version

does ensure the vanishing of the horizontal velocities at the surface but is somewhat inconsistent with the approach of a bulk formula that should operate at some finite height above the surface.

[yairchn]

The issue can be reproduced by calling
include("test/Atmos/EDMF/stable_bl_anelastic1d.jl") in a Julia RPEL
and after its is done calling in the same RPEL
""
include("test/Atmos/EDMF/compute_mse.jl")
dons_arr = get_dons_arr(diag_arr)
using Plots
z = dons_arr[end]["z"]
ρu_1 = dons_arr[end]["prog_ρu_1"]
plot(ρu_1,z,xlabel = "ρu_1",ylabel = "z (m)")
""
The setup is now a non rotating, single column of dry air with some stable temperature profile without any heat fluxes at the surface a sponge at the top and with initial u=1m/s and v=0. In this setup I have no turbulence (i.e. turbulence = ConstantKinematicViscosity(FT(0)),). I would expect in this setup the profiles of u, v, e to stay unchanged in time.

Note - Since it is a 1D Anelastic SingleStack there are no soundwaves and time steps can be large (1sec). Thus I can easily run 0.5h simulation which I am showing here. After 0.5h I see a value of u=-200m/s.

[tapios]

@akshaysridhar suggested raising the lower boundary to a nominal height (e.g., \Delta z = 10 m). Then we can use the corresponding \Delta z for the calculation of fluxes in the similarity theory and apply them as numerical fluxes on the face of the mesh element facing the boundary.

[yairchn]

If I understand correctly this is a change in the manner in which the transfer coefficient 'C' is computed:
C = bc_momentum.drag.fn(state⁻, aux⁻, t, normu_int⁻_tan)

but does it entail a change in of how 'τn' is computed as well:
τn = C * normu_int⁻_tan * u_int⁻_tan

namely in this expression, where are 'normu_int⁻_tan' and 'u_int⁻_tan' computed? at the z=0 interior point or at the z>0 interior point?

[tapios]

If I understand correctly this is a change in the manner in which the transfer coefficient 'C' is computed:
C = bc_momentum.drag.fn(state⁻, aux⁻, t, normu_int⁻_tan)

but does it entail a change in of how 'τn' is computed as well:
τn = C * normu_int⁻_tan * u_int⁻_tan

namely in this expression, where are 'normu_int⁻_tan' and 'u_int⁻_tan' computed? at the z=0 interior point or at the z>0 interior point?

Both would be computed at the interface point (z=0, but we'd nominally interpret as, e.g., z=10 m).

[yairchn]

Both would be computed at the interface point (z=0, but we'd nominally interpret as, e.g., z=10 m).

sounds like a good plan! I am happy to help if needed (@akshaysridhar).

[yairchn]

@akshaysridhar if you are changing the bulk formula notice that the same issues arises in the energy_bc.jl

ClimateMachine.jl/src/Atmos/Model/bc_energy.jl

Line 204 in 6c11bec

fluxᵀn.energy.ρe -= C_h * ρ_avg * normu_int⁻_tan * (MSE - MSE_int)

and moisture_bc.jl

ClimateMachine.jl/src/Atmos/Model/bc_moisture.jl

Line 115 in 6c11bec

fluxᵀn.ρu -= C_q * normu_int⁻_tan * (q_tot - q_tot_int) .* state_int⁻.ρu