Boundary conditions and diffusion with background fields

[liuchihl]

I have modified the simple 1D diffusion example to test how to correct the boundary conditions by using the gradient boundary condition, the script can be found here.

In the experiment, I have specified both the initial temperature (perturbation) and a constant background temperature < T > gradient, so that the total temperature T_total = < T > + T.

The movie below shows T_total and presents two simple experiments with an existing initial temperature. In the blue case, the default no-flux boundary condition is applied (i.e., no gradient in the wall-normal direction), but the gradient at the boundaries is nonzero because the boundary condition does not account for the background temperature. In contrast, the red case includes a boundary condition that forces the gradient of T_total to be 0. With this corrected boundary condition (red curve), the flux at the boundaries is 0, which is physical. However, it remains unclear how to incorporate these corrected fluxes for more complicated configurations.

https://github.com/CliMA/Oceananigans.jl/assets/68127124/69e64cf9-6248-4274-9b6d-5f763827e768

The movie below shows two additional cases without an initial temperature, indicating that the background temperature, < T > defines the entire field. The constant blue line throughout the simulation implies that the background scalar does not diffuse either within the domain or at the boundaries. However, in the red case, despite the absence of diffusion affecting the background temperature, the nonzero flux at the boundaries causes the curve to become smoothed. If diffusion does not affect the background field, would it still make sense if the stratification is not a constant, such as in a Kelvin-Helmholtz instability configuration, e.g., < b > = tanh(z)?

https://github.com/CliMA/Oceananigans.jl/assets/68127124/d350c5ad-2e14-4992-8c5d-f947a1ddf7bb

This brings up the question: what should the default behavior for background fields be? Should they really only appear in the advective terms, or should they also be diffused? And what should the boundary conditions on perturbations be when there are background fields?

[glwagner]

The reason that the advection term is the only term account for is because this is typically what one wants.

If you want diffusion to act on the entire temperature profile (not just the perturbation), then perhaps you shouldn't use a background field at all. Instead you should initialize the temperature to T(z) = N^2 * z + T'(z) where T'(z) is the perturbation and evolve the entire temperature field. This is more correct anyways...

Background fields are useful when one wants to consider either some kind of (implicit) forcing that maintains the background field, preventing it from evolving. If we include a diffusion term associated with the background field, then we defeat much of the purpose of having a background field feature in the first place.

The same is true for boundary conditions. One important use case for background fields is when we can set up a problem where the perturbation is periodic even though the background is not. In that scenario the total field is not periodic, and we could not use a periodic domain. Enter BackgroundFields: decomposing the solution into background + perturbation allows us to model the perturbation in a periodic domain with periodic boundary conditions. If the periodic boundary condition were defined as a boundary condition on the total (background + perturbation), then we wouldn't be able to achieve this important case.

[hdrake]

@glwagner (and @tomchor), I get what you're saying but consider the tilted bottom boundary layer example: https://clima.github.io/OceananigansDocumentation/stable/literated/tilted_bottom_boundary_layer/

In this problem, the buoyancy is decomposed ($b=\bar{b} + b^{\prime}$) into a background component $\bar{b} = N^{2} \hat{z}$ (where $\hat{z}$ is positive in the direction opposite the gravitational vector) and the remaining perturbation. The current example script applies no boundary conditions to the perturbation buoyancy, which I believe defaults to a no-normal-flux condition at the top and bottom boundary. Because no boundary conditions are applied to the background field, there is an implied diffusive flux $\hat{\mathbf{n}} \cdot (-\kappa \nabla \bar{b}) = ±\kappa N^{2} \cos{\theta}$ across these boundaries. This means that the boundary condition on the total buoyancy is a non-zero flux across the seafloor, which does not make much sense for this problem.

The normal way of implementing this problem is to apply a gradient boundary condition, $\frac{\partial b^{\prime}}{\partial z} = -N^{2}\cos{\theta}$, at the bottom, to impose the no-flux condition. (In my MITgcm configurations, I verified that implementing the boundary condition this way recovered the classical analytical solutions of Garrett et al. 1993. Here is another example using Dedalus to replicate some of Jacob Wenegrat's results, where I also include a non-zero gradient boundary condition on the perturbation buoyancy to impose no flux on the total buoyancy.)

I understand why one might want background fields to be constructed a certain way for idealized problems, but it should be communicated clearly (and consistently) to the user exactly how this works. Is the tilted bottom boundary layer example incorrectly implemented or do I still not understand how the background fields work?

[hdrake]

The same is true for boundary conditions. One important use case for background fields is when we can set up a problem where the perturbation is periodic even though the background is not. In that scenario the total field is not periodic, and we could not use a periodic domain. Enter BackgroundFields: decomposing the solution into background + perturbation allows us to model the perturbation in a periodic domain with periodic boundary conditions. If the periodic boundary condition were defined as a boundary condition on the total (background + perturbation), then we wouldn't be able to achieve this important case.

I understand the reasoning for making boundary conditions on perturbations periodic, but what about for bounded dimensions, where one might want to model one (or both) of the boundaries as a solid wall that is impermeable to all fluxes?

[hdrake]

@glwagner (and @tomchor), I get what you're saying but consider the tilted bottom boundary layer example: https://clima.github.io/OceananigansDocumentation/stable/literated/tilted_bottom_boundary_layer/

In this problem, the buoyancy is decomposed (b=b¯+b′) into a background component b¯=N2z^ (where z^ is positive in the direction opposite the gravitational vector) and the remaining perturbation. The current example script applies no boundary conditions to the perturbation buoyancy, which I believe defaults to a no-normal-flux condition at the top and bottom boundary. Because no boundary conditions are applied to the background field, there is an implied diffusive flux n^⋅(−κ∇b¯)=±κN2cos⁡θ across these boundaries. This means that the boundary condition on the total buoyancy is a non-zero flux across the seafloor, which does not make much sense for this problem.

@francispoulin, I am now realizing that this may be related to some of the differences you are seeing in https://github.com/CliMA/Oceananigans.jl/discussions/3529

[francispoulin]

The paper that I cited used periodic boundary conditions for the deviations in the horizontal. Not physically meaningful but it is one way to avoid the walls.

[hdrake]

The paper that I cited used periodic boundary conditions for the deviations in the horizontal. Not physically meaningful but it is one way to avoid the walls.

Agreed. But I am arguing that there is a problem with the boundary conditions in the vertical in your script (and the tilted bottom boundary layer script in the docs). Rearranging the boundary conditions in WT2020,

we get $\frac{\partial b}{\partial z} = - N^{2}_{\infty} \cos {\theta}$ at $z=0$.

[glwagner]

Because no boundary conditions are applied to the background field, there is an implied diffusive flux $\hat{\mathbf{n}} \cdot (-\kappa \nabla \bar{b}) = ±\kappa N^{2} \cos{\theta}$ across these boundaries.

How can this be? There is no $\nabla \cdot \kappa \nabla \bar{b}$ term in the equations --- all terms that solely involve the mean background field are neglected.

Note, I think the total buoyancy is conserved if we use a no-flux condition on the perturbation. The total perturbation is conserved and the background, if prescribed constant, remains constant.

However, the pointwise values of the perturbation field would change if different boundary conditions are used, for sure.

[glwagner]

I understand why one might want background fields to be constructed a certain way for idealized problems, but it should be communicated clearly (and consistently) to the user exactly how this works.

This doesn't really need to be said! Of course, we want to be both clear, and consistent. Otherwise what are we doing here?

I'm not sure where the communication is inconsistent but please fix anything that's wrong.

As for clarity, we can add a line to the page on background fields:

https://clima.github.io/OceananigansDocumentation/stable/model_setup/background_fields/

that says something like "note that boundary conditions apply to the perturbation fields".

As for a convenience feature that applies boundary conditions to the total rather than the perturbation fields: this would be a bit of work. It has seemed to me unnecessary, since if you know what the background fields are, you can apply any boundary condition you like to the perturbations. But convenience is laudable and friendly of course so we can consider it. I think we have to pass the background fields deep into fill_halo_regions! so that, eventually when halos are computed for Gradient or Value boundary conditions, the contribution of the background field can be included.

[hdrake]

Because no boundary conditions are applied to the background field, there is an implied diffusive flux n^⋅(−κ∇b¯)=±κN2cos⁡θ across these boundaries.

How can this be? There is no ∇⋅κ∇b¯ term in the equations --- all terms that solely involve the mean background field are neglected.

The slope-aligned equations in the tilted bottom boundary equations are, I think, derived and formulated differentially than equation sets in other idealized problems with a "background field". The problem is first formulated in terms of conservation equations for total buoyancy (and momentum). Because the background field has no time tendency, the buoyancy evolution equation becomes an evolution equation for just the perturbations but include advection and diffusion of total buoyancy still (including both perturbations and background). There actually should be a $\nabla \cdot (\kappa \nabla \bar{b})$ term in the equations, but with a constant diffusivity and background stratification this term drops out in the interior. This is how all of the dozens of papers on the tilted bottom boundary layer formulate the problem.

[hdrake]

As for clarity, we can add a line to the page on background fields:

https://clima.github.io/OceananigansDocumentation/stable/model_setup/background_fields/

that says something like "note that boundary conditions apply to the perturbation fields".

Yes, I think that would suffice (even though it might be obvious to some).

[glwagner]

There actually should be a term in the equations, but with a constant diffusivity and background stratification this term drops out in the interior. This is how all of the dozens of papers on the tilted bottom boundary layer formulate the problem.

Should we conflate the mean buoyancy with the background field in Oceananigans? All of those papers are interested in computing the evolution of the mean buoyancy field, right? This would have to be treated separately if we try to conflate the mean buoyancy and the Oceananigans background field?

[glwagner]

As for clarity, we can add a line to the page on background fields: https://clima.github.io/OceananigansDocumentation/stable/model_setup/background_fields/ that says something like "note that boundary conditions apply to the perturbation fields".

Yes, I think that would suffice (even though it might be obvious to some).

Very little is obvious I think.

This is a great addition. Just to give a little context --- background fields were easy to add, so I put them in! But I don't use them for my work. I wrote up some basic documentation, but it is not meant to be the Final World by any means. Please feel free to contribute docs. I think an example could be pretty helpful too, showing the decomposition of a field into background and perturbation, then deriving boundary conditions on the perturbation and implementing that. That would be an awesome pedagogical example for the docs. The tilted bottom boundary layer example should also be made clearer if that's possible.

We can also change the functionality, ie how boundary conditions work. But some feature additions might be really invasive to the source code, so there's always a balance to be weighed. It seems at the end of the day it's always pretty important to be hyperintentional about setting up a particular problem, studying the docs and studying the equations you are trying to solve in order to make sure that the simulated problem is what's intended.

[hdrake]

I propose two changes: 1) An example that re-implements the existing 1D diffusion example but with a background field (modifying the boundary conditions on the perturbations to make the two exactly equivalent). @liuchihl, this would be a great first contribution for you 2) Correction to the boundary condition in the present tilted geometry example. I am working on this now.

[hdrake]

There actually should be a term in the equations, but with a constant diffusivity and background stratification this term drops out in the interior. This is how all of the dozens of papers on the tilted bottom boundary layer formulate the problem.

Should we conflate the mean buoyancy with the background field in Oceananigans? All of those papers are interested in computing the evolution of the mean buoyancy field, right? This would have to be treated separately if we try to conflate the mean buoyancy and the Oceananigans background field?

Maybe not, but my point is that these are being conflated currently in https://clima.github.io/OceananigansDocumentation/v0.90.13/literated/tilted_bottom_boundary_layer/

[glwagner]

Hmm, I think perhaps the mean buoyancy actually has to include an additional component that evolves (in addition to the static background component) which is independent of $x$ (upslope coordinate) but maybe varies in $z$ (slope-normal coordinate), right?

It sounds like what that example needs is to include a derivation of the equations we want to solve + boundary conditions, and then to impose intended boundary conditions (whatever those need to be) consistent with the boundary conditions that we'd like to impose on the total (non-decomposed) problem?

Edit: I noticed that you suggested modifying the boundary conditions in the tilted bottom boundary layer problem. My only suggestion is to also include a write up of the equations and a reference to the docs equations, that will make it clear the precise connection between the "original" equations and the perturbation equations that have to be solved by Oceananigans to achieve the desired goal.

[glwagner]

I propose two changes:

1. An example that re-implements the existing 1D diffusion example but with a background field (modifying the boundary conditions on the perturbations to make the two exactly equivalent). @liuchihl, this would be a great first contribution for you
2. Correction to the boundary condition in the present tilted geometry example. I am working on this now.

Do we need both of those examples? I only ask because examples are a little expensive, both computationally (they have to run in CI) and also because they have to be maintained. So we want every "example" to present significant value...

An alternative is to put the 1D example somewhere else, perhaps just a very simple illustration in the docs for background field (that could just be a short write up, then, rather than a full blown example). The 1D and 2D cases are basically the same in terms of boundary conditions right?

Open to various things but just wanted to raise that point since we have tried to be very selective about examples in the past (and we are still very short on examples for they hydrostatic model + immersed boundaries). I also want to build out a much larger selection of simulation examples here: https://github.com/glwagner/AwesomeOceananigans but progress is slow...

[hdrake]

Good point. A simple illustration in the docs for background field will suffice.