Setting up ice base melting boundary condition

[phyo-wai-thaw]

Hi all,

I have been trying to set up the melting boundary condition at the top surface which is also no-slip.
But I don't think it works. Could you please help me with this?
ρw = 1030.0 # kg/m³ (water density)
ρi = 917.0 # kg/m³ (ice density)
Lf = 3.35e5 # J/kg (latent heat of melting)
cp = 4.18e3 # J/kg/K (heat capacity of seawater)
κT = 1.4e-7 # m²/s (thermal diffusivity)
κS = 1.3e-9 # m²/s (salt diffusivity)

meltrate=-(ρw * cp * κT * dTdz) / (ρi * Lf)
meltrate=-(ρw * κS * dSdz) / (ρi * S_b)
a, b = -0.057, -0.311 # Freezing point parameters
T_b = a * S_b + b # Freezing temperature

Thanks a lot.
I appreciate your time and support.

[glwagner]

Hi @phyo-wai-thaw, can you provide more information, such as a script that can be run, and a description of what you expect to happen versus what you observe to happen?

[phyo-wai-thaw]

Hi @glwagner,
Thanks for your reply.
Actually, I am a new use and I am trying to modify this example (https://clima.github.io/OceananigansDocumentation/stable/literated/ocean_wind_mixing_and_convection/), as ice shelf ocean boundary layer simulation as shown in the fig (https://doi.org/10.1073/pnas.2007541118).

The boundary conditions: top=no-slip and melting, bot=free-slip with sponge layer, others are periodic.
Additionally, there is geostrophic balance forcing.
Here the incomplete code.

using Random
using Printf
using CairoMakie
using NetCDF  # For handling snapshot outputs
using Oceananigans
using Oceananigans.Units: minute, minutes, hour

# Constants
ρw = 1030.0  # kg/m³ (water density)
ρi = 917.0   # kg/m³ (ice density)
Lf = 3.35e5  # J/kg (latent heat of melting)
cp = 4.18e3  # J/kg/K (heat capacity of seawater)
κT = 1.4e-7  # m²/s (thermal diffusivity)
κS = 1.3e-9  # m²/s (salt diffusivity)
f = -1.37e-4     # Coriolis parameter
Cd = 2.4e-3  # Drag coefficient
α = 3.8e-5   # Thermal expansion coefficient
β = 7.8e-4   # Haline contraction coefficient
g = 9.81     # Gravitational acceleration
N0 = 1.75e-3  # Reference buoyancy frequency (s⁻¹)
ν =2e-6   #molecular viscosity

# Case A2 Parameters
U0 = 0.014   # m/s (1.4 cm/s)
T0 = -2.1    # Initial far-field temperature (°C)
S0 = 34.5    # Initial far-field salinity (g/kg)
T0star= 0.0025 # Initial thermal driving (T0star=T0-Tf(S0))
#ustar = sqrt(Cd) * U0

Nx = Ny = 32     # number of points in each of horizontal directions
Nz = 24          # number of points in the vertical direction

Lx = Ly = 64     # (m) domain horizontal extents
Lz = 32          # (m) domain depth

refinement = 2 # controls spacing near surface (higher means finer spaced)
stretching = 12  # controls rate of stretching at bottom

# Normalized height ranging from 0 to 1
h(k) = (k - 1) / Nz

# Linear near-surface generator
ζ₀(k) = 1 + (h(k) - 1) / refinement

# Bottom-intensified stretching function
Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))

# Generating function
z_faces(k) = Lz * (ζ₀(k) * Σ(k) - 1)

grid = RectilinearGrid( size = (Nx, Nx, Nz),
                          x = (0, Lx),
                          y = (0, Ly),
                          z = z_faces,
                        topology = (Periodic, Periodic, Bounded))


# Pressure gradient forcing to maintain geostrophic balance
dPdy = -f* U0 * ρw
geostrophic_forcing_u(x, y, z, t) = dPdy/ρw

# sponge layer
damping_rate=1/100
bottom_mask = GaussianMask{:z}(center=-grid.Lz, width=grid.Lz/10)
u_sponge = Relaxation(rate=damping_rate, mask=bottom_mask, target=U0)
T_sponge = Relaxation(rate=damping_rate, mask=bottom_mask, target=T0)
S_sponge = Relaxation(rate=damping_rate, mask=bottom_mask, target=S0)

buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion = α,
                                                                    haline_contraction = β))



no_slip_bc = ValueBoundaryCondition(0.0)
velocity_bcs=FieldBoundaryConditions(top=no_slip_bc)

**# To set up T_bc and and S_bc according to the melting boundary condition
# here I have no idea how to set up, I only know three equations
meltrate=-(ρw * cp * κT * dTdz) / (ρi * Lf)
meltrate=-(ρw * κS * dSdz) / (ρi * S_b)
T_b = -0.057 * S_b - 0.311      # Freezing temperature at ice base**


# Model Setup
model = NonhydrostaticModel(
    grid = grid,
    advection = UpwindBiased(order=5),
    coriolis = FPlane(f=f),
    tracers = (:T, :S),
    buoyancy = buoyancy,
    forcing=(u = (geostrophic_forcing_u, u_sponge), v= (geostrophic_forcing_u, u_sponge), T=T_sponge, S=S_sponge),
    boundary_conditions = (; u=velocity_bcs, v=velocity_bcs),
    closure = AnisotropicMinimumDissipation()
)


# Random noise damped at top and bottom
Ξ(z) = randn() * z / grid.Lz * (1 + z / grid.Lz)  # Noise function

# Temperature initial condition: a stable density gradient with random noise superposed.
Tᵢ(x, y, z) = T0 + 9.5e-4 * z + 9.5e-4 * grid.Lz * 1e-6 * Ξ(z)  # Temperature profile with noise

# Salinity initial condition: linear gradient with noise
Sᵢ(x, y, z) = S0 + 4.0e-4 * z + 4.0e-4 * grid.Lz * 1e-6 * Ξ(z)  # Salinity profile with noise

# Set initial conditions using `set!` function
set!(model, T=Tᵢ, S=Sᵢ)

simulation = Simulation(model, Δt=10.0, stop_time=120minutes)

wizard = TimeStepWizard(cfl=0.9, max_change=1.1, max_Δt=1minute)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

# Print a progress message
progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, max(|w|) = %.1e ms⁻¹, wall time: %s\n",
                                iteration(sim), prettytime(sim), prettytime(sim.Δt),
                                maximum(abs, sim.model.velocities.w), prettytime(sim.run_wall_time))

add_callback!(simulation, progress_message, IterationInterval(20))


# Create a NamedTuple with eddy viscosity
eddy_viscosity = (; νₑ = model.diffusivity_fields.νₑ)

# Define Diagnostics
stratification = (; N²=g / ρw * (α * ∂z(model.tracers.T) - β * ∂z(model.tracers.S)))  # Stratification (Brunt-Väisälä frequency squared)


# Define the filename for the output
filename = "melting"

# Set up the output writer for slices
simulation.output_writers[:slices] =
    JLD2OutputWriter(model, merge(model.velocities, model.tracers, eddy_viscosity, stratification, #meltrate),
                     filename = filename * ".jld2",
                     
                     schedule = TimeInterval(1minute),
                     overwrite_existing = true)

run!(simulation)

Thanks again.

[bethanlily]

Hi @phyo-wai-thaw.

It seems like you're asking how to use the three equation model for the melting to define your boundary conditions for temperature and salinity.
To do this you'll need to first solve the system of equations for the melt rate, the temperature at the ice base and the salinity at the ice base. You'll probably have to write the temperature/salt gradient as a function of the interface value and the value within your domain e.g. $(T_b - T(z_1)) / (z_b - z_1)$. You can use a package such as NonLinearSolve.jl to do the solve.

You can then assign the boundary conditions using a field dependent boundary condition. An example can be found in the documentation for a flux boundary condition. In this case, I think you want to assign value boundary conditions for T and S based on the solutions to your system of equation. You will want to write field and parameter dependent functions that solve the equations and extract the relevant solution for both T and S.

[bethanlily]

Unfortunately it seems like RootSolvers is just for solving single variable equations, rather than the system of (4) equations I need to solve. If you have an other leads for doing this on the GPU that would be great!

[glwagner]

Ah, very interesting. I think we should develop something.

Can you use fixed point iteration? This involves starting with some initial guess, and then updating the guess using a partial evaluation of your nonlinear system. It's often used for the multivariable problem of estimating surface fluxes from Monin-Obukhov theory.

[phyo-wai-thaw]

Thank you all.

I tried to set up, but still not working.
Could you please check it (test_melt_bcs.ipynb)?

I appreciate your time.

[glwagner]

Your boundary condition functions don't have the right signature. Read these docs carefully:

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

Note the difference between discrete_form=true and not.