Simulatin coastal upwelling circulation

[claudioiturra]

Hello, I would appreciate any example or tutorial link on how to simulate coastal up-welling circulation induced by wind over a stratified section of a shallow coastal ocean (depth of ~100 meters). Currently, I'm tweaking the https://clima.github.io/OceananigansDocumentation/stable/literated/ocean_wind_mixing_and_convection/ example, not been able to add up-welling circulation.

Thanks..

[glwagner]

If you can write down a mathematical description of the problem you want to solve, we can help you implement it.

[ajpeliz]

Hi,
I did a simple upwelling case on a meridional channel with topography. Not exactly what you want but you can build on it.
Upwelling.txt

Note: I don't know why but it doesn't let me add files with .jl extention. Please change Upwelling.txt to Upwelling.jl before you use it

[glwagner]

Here's the code:

using Oceananigans
using GLMakie
using Oceananigans.Units

Δt=0.3minutes
stop_time=5days
kzbackground=5e-2  # Background Vertical viscosity

# Upweling in a periodic channel

Nz=40
Lz=1000
Lx=300e3
Ly=10e3
res=3e3 # about a 3 km resolution grid
Nx=Int(round(Lx/res))
Ny=Int(round(Ly/res))

# Set z-gris streatching 
refinement = 4 # 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)

ZZ=[z_faces(i) for i in 1:40]

my_topo=ones(Nx,Ny)*Lz;
[my_topo[i,j]=Lz-8*i for i in 1:Nx, j in 1:Ny]


# we can do a streatchinhg of the vertical grid!
grid1 = RectilinearGrid(size = (Nx, Ny, Nz),
                          topology=(Bounded,Periodic,Bounded),
                          x = (0, Lx),
                          y = (0, Ly),
                          z = z_faces)
grid = ImmersedBoundaryGrid(grid1, PartialCellBottom(-my_topo))

# later we can download a file with a bathymetry
# to read from a Julia file with a var bathymetry in it
#bathymetry = jlopen("nameofthefile.jld2")["bathymetry"]
# there should be a  field with the same dimensions (Need to know if halo points need to be included)
#
# in this case we name the grid with a different name ex. grid_1 and then
# ImerseBoundaryGrid(grid_1,GridFittedBottom(bathymetry))
# 
# here bathymetry works has a mask of all cells bellow that depth
#

# deafult linear eq (based on a thermal expansion and salinity contractions coefs
buoyancy =SeawaterBuoyancy()

#
#  MIXING !!!
#
# Set horizontal and vertical viscosity (ν) and diffusivity (κ)
#Az=HorizontalScalarDiffusivity(ν=1e5,κ=1e2)
Az=HorizontalScalarDiffusivity(ν=1e-1,κ=1e-1)
Kz=VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization(),ν=kzbackground,κ=1e-3)
# simple vertical mixing to force convection (hydrostatic will not do convection on its own)
convective_adjustment=ConvectiveAdjustmentVerticalDiffusivity(convective_κz= 1.0)
closure=(Az,Kz,convective_adjustment)
closure=(Az,Kz)

#
# BCS  !!!!
#
# In case of boundary conditions check Simone Silvester hands on session at 1h:25min
# we can load fields from a .jld2 file with the same size of the grid
# imagine we need wind stress tau:
#  fileOBC=jldopen("nameofthefile.jld2");
#  tau=fileOBC["tau"]
# Alternatively we need to make tau field here!

# for the moment we go:
tau_back= -3e-4
#tau=tau_back*[cos(j*π/Ny) for i=1:Nx,j=1:Ny]
tau=[tau_back for i=1:Nx,j=1:Ny]
tau=tau_back*ones(Nx,Ny);

# Imposing the Wind Stress
# we do it via a upper boundary condition
u_top_bc = FluxBoundaryCondition(0.0)
v_top_bc = FluxBoundaryCondition(-tau)

# Condition for surface heat flux
# we can impose a static value or a restoring
# See silvester 1h:33 hands on workshop
# here we start with a simple value
sst=ones(Nx,Ny)*20.0 
T_top_bc = ValueBoundaryCondition(sst)

u_bcs=FieldBoundaryConditions(top=u_top_bc)
v_bcs=FieldBoundaryConditions(top=v_top_bc)
T_bcs=FieldBoundaryConditions(top=T_top_bc)
#boundary_conditions=(u = u_bcs, v = v_bcs, T = T_bcs)
boundary_conditions=(u = u_bcs, v = v_bcs)


#
# Set model
#
model =HydrostaticFreeSurfaceModel(; grid,
		        		closure,
                                        coriolis = FPlane(f=1e-4),
       		        		boundary_conditions,
		        		buoyancy,
		        		tracers = (:T, :S))

                 	              # coriolis,
                                      # coriolis = FPlane(f=1e-4),

# initial fields;

Tini=ones(Nx,Ny,Nz);
for k=1:Nz 
 Tini[:,:,k].=k;
end

#set!(model, T = 15.0, S=35.0)
set!(model, T = Tini, S=35.0)

#
# define the simulation
#
#

simulation = Simulation(model; Δt, stop_time)
progress(sim) = @info "Time : $(prettytime(sim.model.clock.time)), 
                       max(|u|): $(maximum(abs, sim.model.velocities.u)), 
                       max(|v|): $(maximum(abs, sim.model.velocities.v)), 
                       max(T) : $(maximum(model.tracers.T))"

simulation.callbacks[:progress] =Callback(progress, IterationInterval(100))
run!(simulation)

println(" ")
println(" ")
println(" Sim end ")
println(" ")
println(" ")
println(" Plotting !!! ")
println(" ")
println(" ")

# Extract surface fields
uu=Array(interior(model.velocities.u, :,:,grid.Nz));
vv=Array(interior(model.velocities.v, :,:,grid.Nz));
tt=Array(interior(model.tracers.T   , :,:,grid.Nz));

# make Makie fig
fig = Figure(size = (500, 900));
#def axis
axis_kwargs = (xlabel="x km",
               ylabel="y km",
               aspect = AxisAspect(grid.Lx/grid.Ly),
               limits = ((0,grid.Lx*1e-3), (0, grid.Ly*1e-3)));
# colormaping ranges (same as caxis)
taulims=(-0.001,0.001);
ulims  =(-0.5,0.5);
tlims  =(15,20);
# logica da divisão dos eixos ax=Axis(fig[2,1]) duas linhs uma coluna
ax_tau  = Axis(fig[1, 1]; title = "τₓ", axis_kwargs...);
ax_u    = Axis(fig[2, 1]; title = "Surface u", axis_kwargs...);
ax_v    = Axis(fig[3, 1]; title = "Surface v", axis_kwargs...);
ax_t    = Axis(fig[4, 1]; title = "Surface T", axis_kwargs...);
XX=Array(1:Nx)*res*1e-3; #in km
YY=Array(1:Ny)*res*1e-3; #in km
hm_tau = heatmap!(ax_tau, XX,YY,tau; colormap = :balance, colorrange = taulims);
Colorbar(fig[1, 2], hm_tau; label = "N m⁻2")
hm_uu = heatmap!(ax_u, XX,YY,uu; colormap = :balance, colorrange = ulims);
Colorbar(fig[2, 2], hm_uu; label = "m s⁻1")
hm_vv = heatmap!(ax_v, XX,YY,vv; colormap = :balance, colorrange = ulims);
Colorbar(fig[3, 2], hm_vv; label = "m s⁻1")
hm_tt = heatmap!(ax_t, XX,YY,tt; colormap = :balance, colorrange = tlims);
Colorbar(fig[4, 2], hm_tt; label = "C")
# display it!
fig



bottom = interior(grid.immersed_boundary.bottom_height, :, 1, 1)

# section
figsec = Figure(size = (500, 900));
#def axis
axis_kwargs = (xlabel="x km",
               ylabel="z m",
               limits = ((0,grid.Lx*1e-3), (ZZ[1],0)));

# logica da divisão dos eixos ax=Axis(fig[2,1]) duas linhs uma coluna
axsect    = Axis(figsec[1, 1]; title = "Section t", axis_kwargs...);
tt=Array(interior(model.tracers.T      , :,3,:));
vv=Array(interior(model.velocities.v   , :,3,:));
for i in 1:Nx 
    for j in 1:Nz
       if(tt[i,j]<1.0) 
         tt[i,j] = NaN 
         vv[i,j] = NaN 
      end
    end
end
contourf!(axsect,XX,ZZ,vv)
contour!(axsect,XX,ZZ,tt, color=:red)


# quiver
figquiv = Figure(size = (500, 500));
#def axis
axis_kwargs = (xlabel="x km",
               ylabel="y km",
      #        aspect = AxisAspect(grid.Lx/grid.Ly),
               limits = ((0,grid.Lx*1e-3), (0, grid.Ly*1e-3))  );
axquiv    = Axis(figquiv[1, 1]; title = "Section t", axis_kwargs...);
vv=Array(interior(model.velocities.v  ,:,:,grid.Nz));
uu=Array(interior(model.velocities.u  ,:,:,grid.Nz));
quiver!(axquiv,XX,YY,uu[2:end,:],vv)

# plot bottom mask
bottom_boundary = interior(grid.immersed_boundary.bottom_height, :, 1, 1)
bot = 0 * XX
bot = bot.-grid.Lz
fig2 = Figure(size = (700, 200))
ax = Axis(fig2[1, 1],
                 xlabel="x [km]",
                 ylabel="z [m]",
                 limits=((0, grid.Lx/1e3), (-grid.Lz, 0)))
band!(ax, XX, bot, bottom_boundary, color = :mediumblue)

# DISPLAY PLOTS !!!!!!!!!!!!!
window = display(GLMakie.Screen(), fig2)
window = display(GLMakie.Screen(), figquiv)
window = display(GLMakie.Screen(), figsec)

[claudioiturra]

Hello, first, thank you so much for taking the time to create the code and explain it. The video by Simone Silvestre is excellent (https://www.youtube.com/watch?v=3bhDFI6RaJM&t=3601s) and helped me a lot to understand the basics of Oceananigans. The example is similar to the upwelling example in ROMS, thanks for that! Just one question: even when I completely set the wind stress to 0, I still see a cold temperature near the coast and not a strong vertical thermal gradient in the ocean under no wind and strong thermal stratification. I would appreciate any help understanding why, under no wind, we don't see a strong stratification throughout the entire ocean. I have to tell you that the code contains a lot of information, and I'm still trying to figure out some parts. Again, thanks.

[glwagner]

@claudioiturra I have to say, parsing all of that code is a bit beyond me. I strongly suggest starting from scratch rather than working with existing code, however. Basically it is quite important that you understand every piece of code and select the parameters of your simulation intentionally. As you build up the code, you can test the effect of various pieces. This will help ensure there is not some hidden component of an experiment that affects your results. This exercise will also help cement your understand of Oceananigans. My opinion is that new simulations should be constructed from scratch as much as possible, just using existing simulations as a guide but not copy and pasting code that we do not understand.

Good luck and don't hesitate to continue asking questions, the more specific, the better!

[ajpeliz]

Hello Claudio,
This code is part of a set of numerical labs I just started developing for students. I find Oceananigans very useful as students only have to install Julia. For research I use ROMS and I played with the upwelling example many times during my PhD in Rutgers. Anyway, as for Oceananigans, I am also a newcomer and I cannot advice you. I'm still not aware of the details of the model physics, in particular the bottom, the partial cells functioning and the handling of the mixing. I need to research it myself. I suggest you try with a flat topography and check if you have the same problem. Next, probably, play with the mixing ...
Please report, I'll check it myself as soon as I've time.

[glwagner]

I can help you understand your own script if you like! This is the mixing you prescribed:

Az=HorizontalScalarDiffusivity(ν=1e-1,κ=1e-1)
Kz=VerticalScalarDiffusivity(VerticallyImplicitTimeDiscretization(),ν=kzbackground,κ=1e-3)
# simple vertical mixing to force convection (hydrostatic will not do convection on its own)
convective_adjustment=ConvectiveAdjustmentVerticalDiffusivity(convective_κz= 1.0)
closure=(Az,Kz,convective_adjustment)
closure=(Az,Kz)

This combines horizontal and vertical Laplacian diffusion; the vertical diffusion is treated implicitly in time. The convective adjustment is also treated implicitly in time, which I verified at the REPL:

julia> convective_adjustment=ConvectiveAdjustmentVerticalDiffusivity(convective_κz= 1)
ConvectiveAdjustmentVerticalDiffusivity{VerticallyImplicitTimeDiscretization}(background_κz=0.0 convective_κz=1 background_νz=0.0 convective_νz=0.0)

Partial cells adjust the vertical size of the bottommost cell. This leaves the cells rectangular but adjusts their volume, height, and lateral areas . This can be inspected with

dz_op = zspacings(grid, Center(), Center(), Center())
dz = Field(dz_op)
compute!(dz)

I used your code along with

Oceananigans.ImmersedBoundaries.mask_immersed_field!(dz, NaN)
heatmap(interior(dz, :, 2, :))

to get

it doesn't look like the partial cells are having much of an effect, not sure why. But maybe looking at the difference between a GridFittedBottom and PartialCellBottom would reveal more.

[ajpeliz]

Also, forgot to tell that the filled contours in the plot is for along channel velocity.not temperature.