Langmuir Example: Understanding Differences in BuoyancyTracer() and SeawaterBuoyancy()

[anpa5279]

Hi, I am hoping to get feedback as to what I am doing wrong when comparing these two cases below. For my first case, I took the Langmuir example and ran it. For my second case, I took the Langmuir example and converted buoyancy to temperature via T = b/(g*beta) where beta is the thermal expansion coefficient and changed BuoyancyTracer() to SeawaterBuoyancy(). For both runs, I used the same randomization in the ICs to compare. I thought that these runs would produce the same results (at least close to floating point precision) as one another, but they diverge quickly in time. Below is time t = 0 for both cases and their differences:
BuoyancyTracer():

SeawaterBuoyancy():

Differences:

Below is the last time step for both cases and their differences:
BuoyancyTracer():

SeawaterBuoyancy():

Differences:

Am I wrong to assume that they should match? If not, how can I change my cases so they do match? Thanks!

BuoyancyTracer() code:

grid = RectilinearGrid(; size=(128, 128, 64), extent=(128, 128, 64))

 amplitude = 0.8 # m
wavelength = 60  # m
wavenumber = 2π / wavelength # m⁻¹
 frequency = sqrt(g_Earth * wavenumber) # s⁻¹

# The vertical scale over which the Stokes drift of a monochromatic surface wave
# decays away from the surface is `1/2wavenumber`, or
const vertical_scale = wavelength / 4π

# Stokes drift velocity at the surface
const Uˢ = amplitude^2 * wavenumber * frequency # m s⁻¹
uˢ(z) = Uˢ * exp(z / vertical_scale)
∂z_uˢ(z, t) = 1 / vertical_scale * Uˢ * exp(z / vertical_scale)

τx = -3.72e-5 # m² s⁻², surface kinematic momentum flux
u_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(τx))

Jᵇ = 2.307e-8 # m² s⁻³, surface buoyancy flux
N² = 1.936e-5 # s⁻², initial and bottom buoyancy gradient

b_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Jᵇ),
                                                bottom = GradientBoundaryCondition(N²))

coriolis = FPlane(f=1e-4) # s⁻¹

model = NonhydrostaticModel(; grid, coriolis,
                            advection = WENO(order=9),
                            tracers = :b,
                            buoyancy = BuoyancyTracer(),
                            stokes_drift = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
                            boundary_conditions = (u=u_boundary_conditions, b=b_boundary_conditions))

Ξ(x, y, z) = randn(Xoshiro(1234), (grid.Nx + grid.Ny + grid.Nz+3))[Int(1 + round(grid.Nx*x/grid.Lx+grid.Ny*y/grid.Ly-grid.Nz*z/grid.Lz))] * exp(z / 4)
initial_mixed_layer_depth = 33 # m
stratification(z) = z < - initial_mixed_layer_depth ? N² * z : N² * (-initial_mixed_layer_depth)

bᵢ(x, y, z) = stratification(z) + 1e-1 * Ξ(x, y, z) * N² * model.grid.Lz
u★ = sqrt(abs(τx))
uᵢ(x, y, z) = u★ * 1e-1 * Ξ(x, y, z)
wᵢ(x, y, z) = u★ * 1e-1 * Ξ(x, y, z)

set!(model, u=uᵢ, w=wᵢ, b=bᵢ)

simulation = Simulation(model, Δt=45.0,stop_time=2hours)

conjure_time_step_wizard!(simulation, cfl=1.0, max_Δt=1minute)

function progress(simulation)
    u, v, w = simulation.model.velocities

    # Print a progress message
    msg = @sprintf("i: %04d, t: %s, Δt: %s, umax = (%.1e, %.1e, %.1e) ms⁻¹, wall time: %s\n",
                   iteration(simulation),
                   prettytime(time(simulation)),
                   prettytime(simulation.Δt),
                   maximum(abs, u), maximum(abs, v), maximum(abs, w),
                   prettytime(simulation.run_wall_time))

    @info msg

    return nothing
end

simulation.callbacks[:progress] = Callback(progress, IterationInterval(20))

output_interval = 5minutes

fields_to_output = merge(model.velocities, model.tracers)

simulation.output_writers[:fields] =
    JLD2Writer(model, fields_to_output,
               schedule = TimeInterval(output_interval),
               filename = "local_outputs/original example/langmuir_turbulence_fields.jld2",
               overwrite_existing = true)

u, v, w = model.velocities
b = model.tracers.b

 U = Average(u, dims=(1, 2))
 V = Average(v, dims=(1, 2))
 B = Average(b, dims=(1, 2))
wu = Average(w * u, dims=(1, 2))
wv = Average(w * v, dims=(1, 2))

simulation.output_writers[:averages] =
    JLD2Writer(model, (; U, V, B, wu, wv),
               schedule = AveragedTimeInterval(output_interval, window=2minutes),
               filename = "local_outputs/original example/langmuir_turbulence_averages.jld2",
               overwrite_existing = true)
run!(simulation)

SeawaterBuoyancy() code:

grid = RectilinearGrid(; size=(128, 128, 64), extent=(128, 128, 64))

 amplitude = 0.8 # m
wavelength = 60  # m
wavenumber = 2π / wavelength # m⁻¹
 frequency = sqrt(g_Earth * wavenumber) # s⁻¹

# The vertical scale over which the Stokes drift of a monochromatic surface wave
# decays away from the surface is `1/2wavenumber`, or
const vertical_scale = wavelength / 4π

# Stokes drift velocity at the surface
const Uˢ = amplitude^2 * wavenumber * frequency # m s⁻¹
uˢ(z) = Uˢ * exp(z / vertical_scale)
∂z_uˢ(z, t) = 1 / vertical_scale * Uˢ * exp(z / vertical_scale)

τx = -3.72e-5 # m² s⁻², surface kinematic momentum flux
u_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(τx))

Jᵇ = 2.307e-8 # m² s⁻³, surface buoyancy flux
N² = 1.936e-5 # s⁻², initial and bottom buoyancy gradient
β = 2.0e-4     # 1/K, thermal expansion coefficient

T_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Jᵇ/(g_Earth*β)),
                                                bottom = GradientBoundaryCondition(N²/(g_Earth*β)))

coriolis = FPlane(f=1e-4) # s⁻¹
buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion = β), constant_salinity = 35.0)
model = NonhydrostaticModel(; grid, coriolis,
                            advection = WENO(order=9),
                            tracers = :T,
                            buoyancy = buoyancy,
                            stokes_drift = UniformStokesDrift(∂z_uˢ=∂z_uˢ),
                            boundary_conditions = (u=u_boundary_conditions, T=T_boundary_conditions))

Ξ(x, y, z) = randn(Xoshiro(1234), (grid.Nx + grid.Ny + grid.Nz+3))[Int(1 + round(grid.Nx*x/grid.Lx+grid.Ny*y/grid.Ly-grid.Nz*z/grid.Lz))] * exp(z / 4)
initial_mixed_layer_depth = 33 # m
stratification(z) = z < - initial_mixed_layer_depth ? N² * z/(g_Earth*β) : N² * (-initial_mixed_layer_depth)/(g_Earth*β)

Tᵢ(x, y, z) = stratification(z) + 1e-1 * Ξ(x, y, z) * N² * model.grid.Lz/(g_Earth*β)
u★ = sqrt(abs(τx))
uᵢ(x, y, z) = u★ * 1e-1 * Ξ(x, y, z)
wᵢ(x, y, z) = u★ * 1e-1 * Ξ(x, y, z)

set!(model, u=uᵢ, w=wᵢ, T=Tᵢ)

simulation = Simulation(model, Δt=45.0,stop_time=2hours)

conjure_time_step_wizard!(simulation, cfl=1.0, max_Δt=1minute)

function progress(simulation)
    u, v, w = simulation.model.velocities

    # Print a progress message
    msg = @sprintf("i: %04d, t: %s, Δt: %s, umax = (%.1e, %.1e, %.1e) ms⁻¹, wall time: %s\n",
                   iteration(simulation),
                   prettytime(time(simulation)),
                   prettytime(simulation.Δt),
                   maximum(abs, u), maximum(abs, v), maximum(abs, w),
                   prettytime(simulation.run_wall_time))

    @info msg

    return nothing
end

simulation.callbacks[:progress] = Callback(progress, IterationInterval(20))

output_interval = 5minutes

u, v, w = model.velocities
T = model.tracers.T
b = T * g_Earth * β

simulation.output_writers[:fields] =
    JLD2Writer(model, (; u, v, w, b, T),
               schedule = TimeInterval(output_interval),
               filename = "langmuir_turbulence_fields.jld2",#local_outputs/original example with temperature/
               overwrite_existing = true)

u, v, w = model.velocities

 U = Average(u, dims=(1, 2))
 V = Average(v, dims=(1, 2))
 T = Average(T, dims=(1, 2))
 B = T * g_Earth * β
wu = Average(w * u, dims=(1, 2))
wv = Average(w * v, dims=(1, 2))

simulation.output_writers[:averages] =
    JLD2Writer(model, (; U, V, B, wu, wv, T),
               schedule = AveragedTimeInterval(output_interval, window=2minutes),
               filename = "langmuir_turbulence_averages.jld2",#local_outputs/original example with temperature/
               overwrite_existing = true)
run!(simulation)

[glwagner]

The two results look very close to me --- they are a fraction of a percent in the horizontal average, right? What happens if you run with the same buoyancy model, but different random seeds? Are the differences also comparable then?

[anpa5279]

Sorry, I was not clear in my post. I did not show the percent difference, but the absolute difference w_d = w -w1. In the velocity fields, the differences are the same order of magnitude as the velocity fields. Below I attached the percent difference in the horizontal average at the last time step.
Percent difference relative to the horizontal average:

I took also the horizontal average of the absolute values since some of the averages were very close to zero on the planes in the velocity fields, making the error look worse than it is.
Percent difference relative to the horizontal absolute average:

I just ran these cases with Ξ(x, y, z) = randn(Xoshiro(4321), (grid.Nx+grid.Hx, grid.Ny+grid.Hy, grid.Nz+grid.Hz))[Int(1 + round(grid.Nx*x/grid.Lx)), Int(1 + round(grid.Ny*y/grid.Ly)), Int(1 + round(-grid.Nz*z/grid.Lz))] * exp(z / 4) to change my random seed. Here are the results at the initial time step:
BuoyancyTracer():

SeawaterBuoyancy():

Percent difference relative to the horizontal absolute average:

The last time step:
BuoyancyTracer():

SeawaterBuoyancy():

Percent difference relative to the horizontal absolute average:

[glwagner]

Sorry for being dense, but the buoyancy profiles look identical at the last time step. What exactly is the difference --- can you describe in words? Is the mixed layer deeper, shallower, more energetic, less energetic?

I'm confused what's going on with the horizontally-averaged temperature in the SeawaterBuoyancy case. You are finding a net drift in the temperature below the mixed layer? How is that possible? For example the horizontally-averaged temperature at z=-60 m is about T_xy = -0.6 in the last time step, but -0.9 in the first time-step? So there is net warming below the mixed layer?

It may help to plot the horizontally-averaged results on top of each other.