Sponge layer for regional model

[Nicolas-Theunissen]

Hi! We are participating in the NCAR Open hackathon as team SNL, and have been working on simulating a passive tracer in the Salish Sea using ClimaOcean. I have a few questions about the best way to handle temperature and salinity at the boundaries of our model to insure that they are realistic. So far, we have been initializing our model using ECCO data, roughly following the Mediterranean simulation with ECCO restoring example from ClimaOcean.

Question 1 Looking at the ECCO restoring forcing, since it applies to the whole model space, would it make sense to add an additional sponge layer at the boundary? There does not seem to be any sponge layer in the Mediterranean example.

My understanding of the ECCO restoring function is that it nudges the S and T values back to the interpolated ECCO values on our grid. Implementing a sponge layer would work in a similar way, where the S and T values close to the model border are nudged closer to a known value (in our case the ECCO data). I have attempted to implement this following the sponge layer example here: https://clima.github.io/OceananigansDocumentation/latest/model_setup/forcing_functions/#Relaxation
and the boundary condition documentation here:
https://clima.github.io/OceananigansDocumentation/v0.26.0/model_setup/boundary_conditions/

My implementation is copied below (with a sponge layer implemented on the temperature and ECCO restoring implemented on the salinity). I have also added a figure of the temperature after a 10day spin-up (with ECCO restoring and without a sponge layer). This leads me to Question 2 can there be two forcings applied to the same tracer (ECCOrestoring as well as a sponge layer on the boundaries (which would also be restoring back to ECCO...))? If so, how do I add them together?

Question 3 If it is not possible to have 2 forcings applied to the same tracer, is my implementation of a sponge layer correct?

Question 4 Would changing the mask only apply to the boundaries on ECCOrestoring achieve the same result?

Thank you!

using ClimaOcean
using ClimaOcean.ECCO
using ClimaOcean.ECCO: ECCO4Monthly, NearestNeighborInpainting
using CFTime
using Dates
using Oceananigans
using Oceananigans.Units
using Oceananigans.Units: minute, minutes, hour, hours, day, days, meter, meters, kilometer, kilometers
using Printf


# --------------------------------------------------
# computing architecture
# --------------------------------------------------
arch = GPU()

# --------------------------------------------------
# grid
# --------------------------------------------------
latitude_range = (47, 50)
longitude_range = (234, 239)

# z grid
z_faces = stretched_vertical_faces(;
    depth = 2000, 
    surface_layer_Δz = 2.5,
    surface_layer_height = 25,
    stretching = PowerLawStretching(1.070)
)

# Number of latitude (φ) and longitude (λ) points
# The resolution is 1 / points_per_degree
points_per_degree = 25
Nφ = points_per_degree * (latitude_range[2] - latitude_range[1])
Nλ = points_per_degree * (longitude_range[2] - longitude_range[1])
Nz = length(z_faces) - 1

grid = LatitudeLongitudeGrid(
    arch, 
    size = (Nλ, Nφ, Nz),
    latitude = latitude_range,
    longitude = longitude_range,
    z = z_faces,
    halo = (7, 7, 7)
)

# --------------------------------------------------
# create bathymetry 
# apply ImmersedBoundaryGrid
# --------------------------------------------------

bathymetry = ClimaOcean.regrid_bathymetry(grid;
    interpolation_passes = 40,
    major_basins = 1
)

grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bathymetry)) 

# --------------------------------------------------
# Temperature and salinity restoring force
# --------------------------------------------------
#    :temperature
#    :salinity
#    :u_velocity
#    :v_velocity
#    :free_surface
#    :sea_ice_thickness     
#    :sea_ice_area_fraction
#    :net_heat_flux         

restoring_rate  = 1 / 2days

dates = DateTimeProlepticGregorian(2005, 1, 1) : Month(1) : DateTimeProlepticGregorian(2005, 12, 1)

temperature = ECCOMetadata(:temperature, dates, ECCO4Monthly())
salinity    = ECCOMetadata(:salinity,    dates, ECCO4Monthly())

#FT = ECCORestoring(temperature, arch; rate=restoring_rate, inpainting=NearestNeighborInpainting(50))
FS = ECCORestoring(salinity, arch;    rate=restoring_rate, inpainting=NearestNeighborInpainting(50))
# --------------------------------------------------
# temperature and salinity sponge layer (Only temperature is implemented for now)
# --------------------------------------------------

function create_boundary_mask(nx, ny, nz) #Creates a mask with the boundaries set to 1
    # Create a 3D mask initialized to 0
    mask = zeros(Int, nx, ny, nz)
    
    # Apply the mask on the x-boundaries
    mask[1, :, :] .= 1
    mask[end, :, :] .= 1
    
    # Apply the mask on the y-boundaries
    mask[:, 1, :] .= 1
    mask[:, end, :] .= 1
    
    return mask
end

damping_rate = 1/100 # relax fields on a 100 second time-scale
temperature_gradient = 0.001 # ⁰C m⁻¹
target_temperature = temperature


edge_mask = create_boundary_mask(Nφ, Nλ, Nz)

T_sponge = Relaxation(rate=damping_rate, mask=edge_mask, target=target_temperature)




# --------------------------------------------------
# boundary conditions
# --------------------------------------------------
c_bcs = FieldBoundaryConditions()

# --------------------------------------------------
# ocean model
# --------------------------------------------------
#momentum_advection = VectorInvariant()
#tracer_advection   = Centered(order=2)

#free_surface = SplitExplicitFreeSurface(grid; substeps=30)

ocean = ocean_simulation(grid;
#	momentum_advection, 
#	tracer_advection, 
#	free_surface,
    forcing = (S = FS, T = T_sponge), #Can I have two forcings for the same tracer?
    tracers = (:T, :S, :c),
    boundary_conditions = (c = c_bcs,)
) 

[glwagner]

On Question 1, the restoring can be applied to a portion of the domain using the mask argument to ECCORestoring. Check the docstring to see more, here is an example:

ClimaOcean.jl/experiments/one_degree_simulation/one_degree_simulation.jl

Lines 60 to 69 in 745dbb0

	mask = LinearlyTaperedPolarMask(southern=(-80, -70), northern=(70, 90), z=(z_below_surface, 0))
	dates = DateTimeProlepticGregorian(1993, 1, 1) : Month(1) : DateTimeProlepticGregorian(1993, 11, 1)
	temperature = ECCOMetadata(:temperature; dates, version=ECCO4Monthly(), dir="./")
	salinity = ECCOMetadata(:salinity; dates, version=ECCO4Monthly(), dir="./")
	# inpainting = NearestNeighborInpainting(30) should be enough to fill the gaps near bathymetry
	FT = ECCORestoring(arch, temperature; grid, mask, rate=restoring_rate, inpainting=NearestNeighborInpainting(50))
	FS = ECCORestoring(arch, salinity; grid, mask, rate=restoring_rate, inpainting=NearestNeighborInpainting(50))
	forcing = (T=FT, S=FS)

Since you may not want just to mask at the north and south (so a "polar" mask won't work), you may need to use a field.

[glwagner]

Question 2: you can implement multiple forcings on a tracer using a tuple. There is some recent discussion about that here:

CliMA/Oceananigans.jl#4127

Although, I do not think you will need to do this.

[glwagner]

Looking at questions 3 and 4 I think I have answered them but let me know if there is more to say.

[Nicolas-Theunissen]

The default boundary conditions are not necessarily periodic for a LatitudeLongitudeGrid. It is inferred from the longitude. If the longitude spans 360 degrees. Can you help me understand where you got this impression from? I'm interested in making this clear from the docs.

I was reading about boundary conditions here: https://clima.github.io/OceananigansDocumentation/v0.26.0/model_setup/boundary_conditions/

That makes a lot more sense! Thank you!

[glwagner]

Oof! I am not sure where you found that link, but it is very old, version 0.26! Current Oceananigans is version 0.95.14 (also it is good to ensure that you are running the latest version in your code...)

[glwagner]

Look to the bottom left side of the screen on the page too:

changing it to "stable" produces this screen:

among other things that have happened over the many years between 0.26 and 0.95, Documenter started supporting dark mode!

[glwagner]

Here's a link to the latest

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

[navidcy]

Oof! I am not sure where you found that link, but it is very old, version 0.26!

Google often gives this as a result! Others have bumped onto these archaeological versions of the docs and were asking me about it....