Bubble plume modeling

[samanthasrkim]

Hey all,

I would like to model the transportation of methane gas from the sea bottom to the surface (for example, originating from a leak). Ideally, the model result will give the diluted gas concentration as a function of the depth.
I see that there is a GasExchange model, which I could eventually adapt to methane. Is this a good approach, or would anyone recommend another function/model?
I am also wondering if creating a gas inflow at the sea floor is possible instead of an exchange with the air.

For now, it is fine to start simple, and later add more complexity (everything can always be more complex...).

Thanks :)

[jagoosw]

Hi,

I'm glad to hear you're interest in OceanBioME!

If the methane is dissolved or the bubbles are sufficiently small then you could use a forcing or flux boundary condition (in the Oceananigans documentation) to release them as a tracer from the bottom, and you can use the GasExchange bc at the top.

You would only need to copy the parameters for methane (I will find the paper link later) and it should work otherwise.

[samanthasrkim]

Hey Jago,

Thanks for your answer! Yes, it looks very promising for our application because I started testing with COMSOL Multiphysics, but we want to combine some machine learning or data assimilation later, which seems inadequate with most commercial software (I am from data assimilation, so I was also very happy to see a module for it).

Let me know if you have found the paper. I will give more time for this project later this month, so no hurry.
I will also try to contribute by sharing it, so if you have any additional inputs I need to think of, let me know.

[jagoosw]

Yeah, making data assimilation possible is one of the CliMA people's priorities, so hopefully, we've got the pieces in place to make it relatively easy. I have used it a little bit so if you'd like any more examples of how that should work in OceanBioME.

The paper is Wanninkhof, 2014, and I think the only change you'd need to make is to fill in the parameter values when you make the boundary condition, and tell it how to get the CH4 concentration from the model (where I define MethaneSurfaceConcentration):

import OceanBioME.Models.GasExchangeModel: surface_value

struct MethaneSurfaceConcentration end

@inline surface_value(::MethaneSurfaceConcentration, i, j, grid, clock, model_fields) = @inbounds model_fields.CH₄[i, j, grid.Nz]

transfer_velocity = SchmidtScaledTransferVelocity(schmidt_number = PolynomialParameterisation{4}(FT; coefficients = (a, b, c, d, e))) # (a, b, c, d, e) from Table 1

water_concentration = MethaneSurfaceConcentration()

air_concentration = PartiallySolubleGas(air_concentration = 9352.7, # change to whatever CH4 is 
                                        solubility = Wanninkhof92Solubility(A1, A2, A3, B1, B2, B3)) #(A1, A2, A3, B1, B2, B3) from Table 2

CH₄_surface_bc = GasExchangeBoundaryCondition(water_concentration, air_concentration, transfer_velocity)

[samanthasrkim]

Super, thank you very much. I'll try!

[samanthasrkim]

Hey all!
I started working on this model. Thanks again for the guidance; it was very helpful to get acquainted with Oceanigans!
So, I have another question (which may seem trivial...): how to make the bubbles rise?

I used a forcing to create a bubble inlet at the bottom. I also use the example of convecting plankton to understand the different model settings in Oceanigans. This example seems similar to my case, but with the bubble properties and dynamics instead of the plankton dynamics.
Long story short, now I am at the step where I would like to simulate a tracer that rises simply because of the density difference.

• I tried using the buoyancy as a tracer in a Nonhydrostaticmodel: buoyancy=BuoyancyTracer():

model = NonhydrostaticModel(
grid = grid, -> simple grid
closure = ScalarDiffusivity(ν=1e-4, κ=1e-4), -> some diffusivity
tracers = (:b, :G), -> a tracer b for the buoyancy and a tracer G for the gas
buoyancy = BuoyancyTracer(),
forcing = (; G = gas_forcing) -> gas_forcing is the forcing that creates the gas inlet at the bottom
)

Is it the correct approach, or is there a straightforward way to do it that I didn't see? Also, is this approach suitable if, later, I want to specify the gas parameters, including initial velocity and bubble/particle size, specific density, for example?

Looking forward to hearing from you!

[glwagner]

@simone-silvestri recently put in a feature that allows you to augment particle velocities with additional components, eg a rising velocity:

CliMA/Oceananigans.jl#4448

[samanthasrkim]

Hey, thank you, it is exactly what I mean.
For the velocity, it was just for the initial condition that v=0 at the sea bottom, but I think it is not necessary like this.

1. I went through the code, and it is very clear! Also nice that you added G_escape. Based on your previous answers, I tried to add the GasExchangeModel as a boundary condition for the top, using:

CH₄_surface_bc = GasExchangeBoundaryCondition(; water_concentration, air_concentration, transfer_velocity, wind_speed)

I am not sure how to understand water_concentration (I understand what it means in the model, but not how it is built). I've checked the glossary in OceanBioME:

"water_concentration should usually be a [Tracer]Concentration, where is the name of the tracer (you will have to build your own if this is not OxygenConcentration), or a CarbonDioxideConcentration, which diagnoses the partial pressure of CO₂ in the water."

So I also checked in GasExchange.jl to try to understand, based on OxygenConcentration (

OceanBioME.jl/src/Models/GasExchange/GasExchange.jl

Lines 122 to 135 in a84ab98

	"""
	OxygenGasExchangeBoundaryCondition(FT = Float64;
	transfer_velocity = SchmidtScaledTransferVelocity(schmidt_number = OxygenPolynomialSchmidtNumber(FT)),
	water_concentration = OxygenConcentration(),
	air_concentration = 9352.7, # mmolO₂/m³
	wind_speed = 2,
	kwagrs...)
	Returns a `FluxBoundaryCondition` for the gas exchange between oxygen dissolved in the water
	specified by the the `OxygenConcentration` in the base model, and `air_concentration` with `transfer_velocity`
	(see `GasExchangeBoundaryCondition` for details).
	`kwargs` are passed on to `GasExchangeBoundaryCondition`.
	"""

).

If you have more documentation about it, let me know.

I end up having a "working" script, but the simulation stops at the first time step because of negative values in a log(). I am unsure where this happens, but it might come from the parameter values. I'll copy the script, so if you have time, feel free to see if you see a mistake.

2. As an alternative, I considered using the water_concentration (which is, to be sure that I understand correctly, the gas concentration at the interface water/air) as a flux boundary condition in:

G_escape = FluxBoundaryCondition((x, t, G, flux_rate) -> flux_rate * G, parameters = 1/hour, field_dependencies = :G)

using water_concentration to set "all gas reaching the surface can escape". But I still don't understand the water_concentration format.

Thanks again for your help, and it is very helpful to understand Oceanigans and OceanBioME better!

Here is the script:

using Oceananigans, Oceananigans.Units, Printf
using Oceananigans.BuoyancyFormulations: AbstractBuoyancyFormulation
import Oceananigans.BuoyancyFormulations: required_tracers, buoyancy_perturbationᶜᶜᶜ, ∂x_b, ∂y_b, ∂z_b, top_buoyancy_flux, bottom_buoyancy_flux
import OceanBioME.Models.GasExchangeModel: surface_value, PolynomialParameterisation, SchmidtScaledTransferVelocity, GasExchangeBoundaryCondition, PartiallySolubleGas, Wanninkhof92Solubility

struct GasDrivenBuoyancy{FT} <: AbstractBuoyancyFormulation{Nothing}
    density_anomaly :: FT # gΔρ/ρₒ
end
struct MethaneSurfaceConcentration end

required_tracers(::GasDrivenBuoyancy) = (:G,)

## Setting the buoyancy of the water to be b=-rho beta G, where G is the gas fraction of the water <<1
buoyancy_perturbationᶜᶜᶜ(i, j, k, grid, buoyancy::GasDrivenBuoyancy, C) = buoyancy.density_anomaly * C.G[i, j, k]

@inline ∂x_b(i, j, k, grid, buoyancy::GasDrivenBuoyancy, C) = buoyancy.density_anomaly * ∂xᶠᶜᶜ(i, j, k, grid, C.G)
@inline ∂y_b(i, j, k, grid, buoyancy::GasDrivenBuoyancy, C) = buoyancy.density_anomaly * ∂yᶜᶠᶜ(i, j, k, grid, C.G)
@inline ∂z_b(i, j, k, grid, buoyancy::GasDrivenBuoyancy, C) = buoyancy.density_anomaly * ∂zᶜᶜᶠ(i, j, k, grid, C.G)

@inline    top_buoyancy_flux(i, j, grid, ::GasDrivenBuoyancy, top_tracer_bcs, clock, fields) = getbc(top_tracer_bcs.G, i, j, grid, clock, fields)
@inline bottom_buoyancy_flux(i, j, grid, ::GasDrivenBuoyancy, bottom_tracer_bcs, clock, fields) = getbc(bottom_tracer_bcs.G, i, j, grid, clock, fields)

## Setting the MethaneSurfaceConcentration
@inline surface_value(::MethaneSurfaceConcentration, i, j, grid, clock, C) = @inbounds C.G[i, j, grid.Nz]

grid = RectilinearGrid(topology = (Periodic, Flat, Bounded), 
                                   size = (32, 32), 
                                   extent = (20, 10))

buoyancy = BuoyancyForce(GasDrivenBuoyancy(0.902*9.81))

G_injection = Forcing((x, z, t, injection_rate) -> ifelse((9.5<x<10.5) & (z < -9.5), injection_rate, 0), parameters = 1e-3/hour) # G is a volume fraction that must be << 1

# For the bc at interface air/water:
transfer_velocity = SchmidtScaledTransferVelocity(schmidt_number = PolynomialParameterisation{4}(coefficients = (2101.2, -131.54, 4.4931, -0.08676, 0.00070663)))
water_concentration = MethaneSurfaceConcentration()
air_concentration = PartiallySolubleGas(air_concentration = 2000, # average mole fraction in ppb
                                        solubility = Wanninkhof92Solubility(-68.8862, 101.4956, 28.7314, -0.076146, 0.043970, -0.006872))
wind_speed = 1

CH₄_surface_bc = GasExchangeBoundaryCondition(; water_concentration,
                                                air_concentration,
                                                transfer_velocity,
                                                wind_speed)

# Changing G_bcs for the GasExchangeBoundaryCondition
G_bcs = FieldBoundaryConditions(top = CH₄_surface_bc)
#G_bcs = FieldBoundaryConditions(top = G_escape)

advection = WENO()

model = NonhydrostaticModel(; grid, 
                            buoyancy, 
                            advection, 
                            tracers = (:G, :T, :S), 
                            forcing = (; G = G_injection), 
                            boundary_conditions = (; G = G_bcs), 
                            hydrostatic_pressure_anomaly=nothing)

buoyancy_timescale = √(minimum_zspacing(grid) / (1e-3 * 0.902 * 9.81))

set!(model, T=10.0, S=35.0)

simulation = Simulation(model, Δt = buoyancy_timescale, stop_time = 12hours)

[samanthasrkim]

@glwagner cool, I'll have a look! Thanks

[jagoosw]

I think there are two problems here: the gas exchange boundary conditions are for dissolved gases, not bubbles, and how I wrote the buoyancy model assumes that the gas is a volume fraction $\ll1$, whereas the boundary condition is expecting the concentration in some other units (I guess mmol/m^3 or something similar). This is causing the model to fail because there is a huge downward flux of G, making a very strong buoyancy.

It might be most straightforward to work out the factor converting between volume fraction and methane concentration, and then put that into the buoyancy model, but for the boundary condition, I think you will need to let G escape in some way, like the boundary condition I wrote.

G should also dissolve into the water and the dissolved form of the gas can have the GasExchangeBoundaryCondition but I don't know how the dissolution is modelled etc.

As for your question about water_concentration, for "simple" gases like oxygen, the concentration involved in the surface exchange is just the concentration of dissolved gas, but for e.g. carbon dioxide we have to compute it from the DIC and alkalinity. I guess you can just use the dissolved concentration for methane but I don't know if it is involved in some seawater chemistry that means you would need todo something else.

[samanthasrkim]

Okay, I see. This is a good start, I will continue with it! Thanks