Passive tracers are not conserved on tripolar grid if wind is blowing

[NoraLoose]

I noticed that passive tracers are not conserved in my global simulation (with a tripolar grid and the z-star coordinate) when a prescribed atmosphere is added (i.e., with JRA55 and ClimaOcean). I’ve been able to isolate the issue to an Oceananigans-only MWE, adapted from the tripolar grid tracer conservation test in our suite. The only change is the addition of wind forcing.

Below is the simplified test setup used to reproduce the issue:

function run_tripolar_test_simulation(arch; with_wind = false)

    z_stretched = MutableVerticalDiscretization(collect(-20:0))
    grid = TripolarGrid(arch; size = (20, 20, 20), z = z_stretched)

    function mtn₁(λ, φ)
        λ₁ = 70
        φ₁ = 55
        dφ = 5
        return exp(-((λ - λ₁)^2 + (φ - φ₁)^2) / 2dφ^2)
    end

    function mtn₂(λ, φ)
        λ₁ = 70
        λ₂ = λ₁ + 180
        φ₂ = 55
        dφ = 5
        return exp(-((λ - λ₂)^2 + (φ - φ₂)^2) / 2dφ^2)
    end

    zb = - 20
    h  = - zb + 10
    gaussian_islands(λ, φ) = zb + h * (mtn₁(λ, φ) + mtn₂(λ, φ))

    grid = ImmersedBoundaryGrid(grid, GridFittedBottom(gaussian_islands))
    free_surface = SplitExplicitFreeSurface(grid; substeps=10)
    
    kwargs = (;
        grid,
        free_surface,
        tracers = (:b, :c),
        buoyancy = BuoyancyTracer(),
        vertical_coordinate = ZStar()
    )

    if with_wind
        u₁₀, cᴰ, ρₒ, ρₐ = 10.0, 2.5e-3, 1026.0, 1.225
        τx = -ρₐ / ρₒ * cᴰ * u₁₀ * abs(u₁₀)
        u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τx))
        model = HydrostaticFreeSurfaceModel(; kwargs..., boundary_conditions = (u = u_bcs,))    
    else
        model = HydrostaticFreeSurfaceModel(; kwargs...)
    end
    
    bᵢ(x, y, z) = y < 0 ? 0.06 : 0.01

    # Instead of initializing with random velocities, infer them from a random initial streamfunction
    # to ensure the velocity field is divergence-free at initialization.
    ψ = Field{Center, Center, Center}(grid)
    set!(ψ, rand(size(ψ)...))
    uᵢ = ∂y(ψ)
    vᵢ = -∂x(ψ)

    set!(model, c = (x, y, z) -> rand(), u = uᵢ, v = vᵢ, b = bᵢ)

    simulation = Simulation(model, Δt=2minutes, stop_time=1day)
    
    c = simulation.model.tracers.c
    tracer_integral = Integral(c, dims=(1, 2, 3))

    if with_wind
        prefix = "test_with_wind"
    else
        prefix = "test_no_wind"
    end
    
    integral_ow = JLD2Writer(
        simulation.model,
        (; tracer_integral),
        filename = "$(prefix)_integrals",
        schedule = IterationInterval(1),
        overwrite_existing = true,
    )
    simulation.output_writers[:integral] = integral_ow

    run!(simulation)
    
    return nothing
end

run_tripolar_test_simulation(GPU(); with_wind = false)
run_tripolar_test_simulation(GPU(); with_wind = true)

Visualizing the results shows that the passive tracer c is

• conserved for the simulation without wind, which matches the behavior tested in our CI.
• not conserved for the simulation with wind, despite all other settings being identical.

[simone-silvestri]

It looks like it is a tripolar grid specific problem (I have tried with other grids, immersed and not and tracer is conserved).
Could it play a role with the fact that we need to rotate vectors in a tripolar grid? Maybe there is a bug there...

[glwagner]

is it at all posible there is a bug in the metrics? or we test orthogonality there

[NoraLoose]

If it's a bug in the metrics, why doesn't the problem show up without wind forcing?

Another possibility is that there's a bug in how the halo regions are filled for the top flux boundary conditions.

[glwagner]

I'm not sure --- but it's already mysterious that there is an issue depending on whether there is wind forcing or not, since as far as I can tell the only possible distinction is the nature of the resulting velocity field. The wind forcing does not directly impact tracer advection.

[NoraLoose]

I'd like to better understand how exchanges occur across the tripolar fold, since it's easy to introduce bugs in that region.

A few years ago, I created a schematic showing how tracer and velocity exchanges work on a tripolar B-grid.

@simone-silvestri, on our tripolar C-grid, where exactly is the pivot point located (analogous to the green point in my schematic)? Is it positioned at a vorticity point (as on the B-grid), or at a tracer point?

[glwagner]

I believe it's at a u-point, so that when we make the fold, the two halfs of the grid nest within one another (meaning that one half of the cells are redundant). We resolve the redundancy when we fill the halo regions by overwriting half the values.

[NoraLoose]

We should test whether wind forcing away from the tripolar fold leads to tracer non-conservation. That would help confirm whether the issue is actually related to the fold region.

[NoraLoose]

The tripolar fold seams to be the issue, because wind forcing applied only in the Southern Hemisphere does not lead to tracer non-conservation.

Here is the script. Same as above, just with slightly modified wind forcing.

function run_tripolar_test_simulation(arch; wind_away_from_tripolar_fold = false)
    z_stretched = MutableVerticalDiscretization(collect(-20:0))
    grid = TripolarGrid(arch; size = (20, 20, 20), z = z_stretched)

    function mtn₁(λ, φ)
        λ₁ = 70
        φ₁ = 55
        dφ = 5
        return exp(-((λ - λ₁)^2 + (φ - φ₁)^2) / 2dφ^2)
    end

    function mtn₂(λ, φ)
        λ₁ = 70
        λ₂ = λ₁ + 180
        φ₂ = 55
        dφ = 5
        return exp(-((λ - λ₂)^2 + (φ - φ₂)^2) / 2dφ^2)
    end

    zb = - 20
    h  = - zb + 10
    gaussian_islands(λ, φ) = zb + h * (mtn₁(λ, φ) + mtn₂(λ, φ))

    grid = ImmersedBoundaryGrid(grid, GridFittedBottom(gaussian_islands))
    free_surface = SplitExplicitFreeSurface(grid; substeps=10)
    
    kwargs = (;
        grid,
        free_surface,
        tracers = (:b, :c),
        buoyancy = BuoyancyTracer(),
        vertical_coordinate = ZStar()
    )

    τx = wind_away_from_tripolar_fold ?
        ((x, y, z) -> y < 0 ? -1e-4 : 0.0) :
        ((x, y, z) -> -1e-4)

    prefix = wind_away_from_tripolar_fold ?
        "wind_southern_hemisphere_only" :
        "wind_global"
    
    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τx))
    model = HydrostaticFreeSurfaceModel(; kwargs..., boundary_conditions = (u = u_bcs,))    
    
    bᵢ(x, y, z) = y < 0 ? 0.06 : 0.01

    ψ = Field{Face, Face, Center}(grid)
    set!(ψ, rand(size(ψ)...))
    uᵢ = ∂y(ψ)
    vᵢ = -∂x(ψ)

    set!(model, c = (x, y, z) -> rand(), u = uᵢ, v = vᵢ, b = bᵢ)

    simulation = Simulation(model, Δt=2minutes, stop_time=1day)
    
    c = simulation.model.tracers.c
    tracer_integral = Integral(c, dims=(1, 2, 3))

    integral_ow = JLD2Writer(
        simulation.model,
        (; tracer_integral),
        filename = "$(prefix)_integrals",
        dir = target_dir,
        schedule = IterationInterval(1),
        overwrite_existing = true,
    )
    simulation.output_writers[:integral] = integral_ow

    run!(simulation)
    
    return nothing
end

run_tripolar_test_simulation(GPU(); wind_away_from_tripolar_fold = false)
run_tripolar_test_simulation(GPU(); wind_away_from_tripolar_fold = true)

[glwagner]

Ok I have a thought. I don't think we are interpreting vector components correctly on the tripolar grid.

@NoraLoose if you form the zonal wind stress from a streamfunction instead, and then compute both τx and τy from the wind stress streamfunction --- does that solve the issue?

Basically I think that we are forcing the u component but this is not consistently a zonal forcing. As a result, there is an inconsistent forcing on either side of the fold.

Need to think about the best way to solve this. We have to decide how to interpret user input (we may need to add more information somewhere, either globally to the model constructor or locally to BoundaryCondition)

[glwagner]

I guess this is something like

ψτ = Field{Face, Face, Center}(grid)
set!(ψτ, (x, y, z) -> -1e-4 * y)
τx = + ∂y(ψτ)
τy = - ∂x(ψτ)

[NoraLoose]

I just tried what you suggested, and swapped in

    ψτ = Field{Face, Face, Center}(grid)
    set!(ψτ, (x, y, z) -> -1e-4 * y)
    τx = + ∂y(ψτ)
    τy = - ∂x(ψτ)
    
    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τx))
    v_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τy))
    model = HydrostaticFreeSurfaceModel(; kwargs..., boundary_conditions = (u = u_bcs, v = v_bcs, ))

and tracers are conserved!

[NoraLoose]

Need to think about the best way to solve this. We have to decide how to interpret user input (we may need to add more information somewhere, either globally to the model constructor or locally to BoundaryCondition)

And we will have to handle JRA55 forcing in ClimaOcean correctly, and regrid in a way that the vector components are correctly interpreted. (Applying JRA55 forcing is how the issue arose for me originally.)

[glwagner]

ok so at least we have a way to correctly set momentum boundary conditions...

the next question is how we should write the user interface. There's a few options. There's also the slightly annoying situation that u, v are not zonal, meridional on TripolarGrid. But it seems like we do want to write things like set!(model, u=something) and something is zonal.

I believe JRA55 is rotated, but @simone-silvestri can answer more conclusively.