From 5b1ed6f43b8ee93cf9b75913c6694767e7d74da8 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 5 Mar 2026 11:34:57 -0900 Subject: [PATCH 1/3] Replace GM with area-scaled biharmonic viscosity in coarse examples MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Remove IsopycnalSkewSymmetricDiffusivity (Gent-McWilliams) from the one-degree and global climate simulation examples. Add a new `area_scaled_biharmonic_viscosity()` function that computes ν = Az² / λ, where Az is the horizontal cell area and λ is a damping timescale (default 15 days). This resolution-dependent viscosity adapts to local grid cell size, providing appropriate grid-scale dissipation at all latitudes. Co-Authored-By: Claude Opus 4.6 --- examples/global_climate_simulation.jl | 7 +++---- examples/one_degree_simulation.jl | 13 +++++-------- src/NumericalEarth.jl | 1 + src/Oceans/Oceans.jl | 2 +- src/Oceans/ocean_simulation.jl | 21 ++++++++++++++++++++- 5 files changed, 30 insertions(+), 14 deletions(-) diff --git a/examples/global_climate_simulation.jl b/examples/global_climate_simulation.jl index 5a40eeac2..7f8f91916 100644 --- a/examples/global_climate_simulation.jl +++ b/examples/global_climate_simulation.jl @@ -1,7 +1,7 @@ # # Global climate simulation # # This example configures a global ocean--sea ice simulation at 1.5ᵒ horizontal resolution with -# realistic bathymetry and a few closures including the "Gent-McWilliams" `IsopycnalSkewSymmetricDiffusivity`. +# realistic bathymetry and an area-scaled biharmonic viscosity for grid-scale dissipation. # The atmosphere is represented by a 4-layer [SpeedyWeather](https://github.com/SpeedyWeather/SpeedyWeather.jl) # simulation on the T63 spectral grid (this grid has approximately 1.875ᵒ resolution). # @@ -48,9 +48,8 @@ momentum_advection = VectorInvariant() tracer_advection = WENO(order=5) free_surface = SplitExplicitFreeSurface(grid; substeps=40) catke_closure = NumericalEarth.Oceans.default_ocean_closure() -eddy_closure = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1e3, κ_symmetric=1e3) -viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarBiharmonicDiffusivity(ν=1e12) -closures = (catke_closure, eddy_closure, viscous_closure) +viscous_closure = area_scaled_biharmonic_viscosity() +closures = (catke_closure, viscous_closure) nothing #hide # The ocean simulation, complete with initial conditions for temperature and salinity from ECCO. diff --git a/examples/one_degree_simulation.jl b/examples/one_degree_simulation.jl index 9a2dbc005..cc3805b10 100644 --- a/examples/one_degree_simulation.jl +++ b/examples/one_degree_simulation.jl @@ -1,7 +1,7 @@ # # [One-degree global ocean--sea ice simulation](@id one-degree-ocean-seaice) # # This example configures a global ocean--sea ice simulation at 1ᵒ horizontal resolution with -# realistic bathymetry and a few closures including the "Gent-McWilliams" `IsopycnalSkewSymmetricDiffusivity`. +# realistic bathymetry and an area-scaled biharmonic viscosity for grid-scale dissipation. # The simulation is forced by repeat-year JRA55 atmospheric reanalysis # and initialized by temperature, salinity, sea ice concentration, and sea ice thickness # from the ECCO state estimate. @@ -47,13 +47,10 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); # ### Closures # -# We include a Gent-McWilliams isopycnal diffusivity as a parameterization for the mesoscale -# eddy fluxes. For vertical mixing at the upper-ocean boundary layer we include the CATKE -# parameterization. +# We use an area-scaled biharmonic viscosity for grid-scale dissipation and +# the CATKE parameterization for vertical mixing at the upper-ocean boundary layer. -using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity, AdvectiveFormulation - -eddy_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew=1e3, κ_symmetric=1e3, skew_flux_formulation=AdvectiveFormulation()) +horizontal_viscosity = area_scaled_biharmonic_viscosity() vertical_mixing = NumericalEarth.Oceans.default_ocean_closure() # ### Ocean simulation @@ -65,7 +62,7 @@ momentum_advection = WENOVectorInvariant(order=5) tracer_advection = WENO(order=5) ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surface, - closure=(eddy_closure, vertical_mixing)) + closure=(horizontal_viscosity, vertical_mixing)) @info "We've built an ocean simulation with model:" @show ocean.model diff --git a/src/NumericalEarth.jl b/src/NumericalEarth.jl index 8b23f354a..ff98195be 100644 --- a/src/NumericalEarth.jl +++ b/src/NumericalEarth.jl @@ -47,6 +47,7 @@ export LinearlyTaperedPolarMask, DatasetRestoring, ocean_simulation, + area_scaled_biharmonic_viscosity, ORCAGrid, sea_ice_simulation, atmosphere_simulation, diff --git a/src/Oceans/Oceans.jl b/src/Oceans/Oceans.jl index 93b1d2727..e26b26a2c 100644 --- a/src/Oceans/Oceans.jl +++ b/src/Oceans/Oceans.jl @@ -1,6 +1,6 @@ module Oceans -export ocean_simulation, SlabOcean +export ocean_simulation, SlabOcean, area_scaled_biharmonic_viscosity using Oceananigans using Oceananigans.Units diff --git a/src/Oceans/ocean_simulation.jl b/src/Oceans/ocean_simulation.jl index 88627bb1b..15c551511 100644 --- a/src/Oceans/ocean_simulation.jl +++ b/src/Oceans/ocean_simulation.jl @@ -4,7 +4,8 @@ using Oceananigans.BoundaryConditions: DefaultBoundaryCondition using Oceananigans.ImmersedBoundaries: immersed_peripheral_node, inactive_node, MutableGridOfSomeKind using Oceananigans.OrthogonalSphericalShellGrids -using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization +using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization, + HorizontalScalarBiharmonicDiffusivity using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity, @@ -84,6 +85,24 @@ function default_ocean_closure(FT=Oceananigans.defaults.FloatType) return CATKEVerticalDiffusivity(VerticallyImplicitTimeDiscretization(), FT; mixing_length, turbulent_kinetic_energy_equation) end +@inline _area_scaled_biharmonic_viscosity(i, j, k, grid, ℓx, ℓy, ℓz, clock, fields, λ) = + Az(i, j, k, grid, ℓx, ℓy, ℓz)^2 / λ + +""" + area_scaled_biharmonic_viscosity(FT=Oceananigans.defaults.FloatType; timescale=15days) + +Return a `HorizontalScalarBiharmonicDiffusivity` with a grid-area-scaled viscosity +coefficient `ν = Az² / timescale`, where `Az` is the horizontal area of each grid cell. +This produces a viscosity that adapts to local grid cell size, providing appropriate +damping at all latitudes. +""" +function area_scaled_biharmonic_viscosity(FT=Oceananigans.defaults.FloatType; timescale=15days) + return HorizontalScalarBiharmonicDiffusivity(FT; + ν = _area_scaled_biharmonic_viscosity, + discrete_form = true, + parameters = timescale) +end + function default_radiative_forcing(grid) ϵʳ = 0.6 # red fraction λʳ = 1 # red decay scale From 047fcdeeaf6a0425b0529120a7511369e2d2a21f Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 5 Mar 2026 11:45:56 -0900 Subject: [PATCH 2/3] Pass explicit timescale to area_scaled_biharmonic_viscosity in examples Co-Authored-By: Claude Opus 4.6 --- examples/global_climate_simulation.jl | 2 +- examples/one_degree_simulation.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/global_climate_simulation.jl b/examples/global_climate_simulation.jl index 7f8f91916..b0eb14604 100644 --- a/examples/global_climate_simulation.jl +++ b/examples/global_climate_simulation.jl @@ -48,7 +48,7 @@ momentum_advection = VectorInvariant() tracer_advection = WENO(order=5) free_surface = SplitExplicitFreeSurface(grid; substeps=40) catke_closure = NumericalEarth.Oceans.default_ocean_closure() -viscous_closure = area_scaled_biharmonic_viscosity() +viscous_closure = area_scaled_biharmonic_viscosity(timescale=15days) closures = (catke_closure, viscous_closure) nothing #hide diff --git a/examples/one_degree_simulation.jl b/examples/one_degree_simulation.jl index cc3805b10..b1039c791 100644 --- a/examples/one_degree_simulation.jl +++ b/examples/one_degree_simulation.jl @@ -50,7 +50,7 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); # We use an area-scaled biharmonic viscosity for grid-scale dissipation and # the CATKE parameterization for vertical mixing at the upper-ocean boundary layer. -horizontal_viscosity = area_scaled_biharmonic_viscosity() +horizontal_viscosity = area_scaled_biharmonic_viscosity(timescale=15days) vertical_mixing = NumericalEarth.Oceans.default_ocean_closure() # ### Ocean simulation From 336b890099607b7ffe1537f2240d448177b150a6 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Thu, 5 Mar 2026 13:35:58 -0900 Subject: [PATCH 3/3] Add IsopycnalSkewSymmetricDiffusivity with symmetric-only diffusivity to examples MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Use κ_symmetric=1000, κ_skew=0 for isopycnal lateral tracer mixing without the GM skew flux. Co-Authored-By: Claude Opus 4.6 --- examples/global_climate_simulation.jl | 5 ++++- examples/one_degree_simulation.jl | 8 ++++++-- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/examples/global_climate_simulation.jl b/examples/global_climate_simulation.jl index b0eb14604..3a1c4d002 100644 --- a/examples/global_climate_simulation.jl +++ b/examples/global_climate_simulation.jl @@ -47,9 +47,12 @@ nothing #hide momentum_advection = VectorInvariant() tracer_advection = WENO(order=5) free_surface = SplitExplicitFreeSurface(grid; substeps=40) +using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity + catke_closure = NumericalEarth.Oceans.default_ocean_closure() viscous_closure = area_scaled_biharmonic_viscosity(timescale=15days) -closures = (catke_closure, viscous_closure) +isopycnal_closure = IsopycnalSkewSymmetricDiffusivity(κ_skew=0, κ_symmetric=1e3) +closures = (catke_closure, viscous_closure, isopycnal_closure) nothing #hide # The ocean simulation, complete with initial conditions for temperature and salinity from ECCO. diff --git a/examples/one_degree_simulation.jl b/examples/one_degree_simulation.jl index b1039c791..a77ed8887 100644 --- a/examples/one_degree_simulation.jl +++ b/examples/one_degree_simulation.jl @@ -47,10 +47,14 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(bottom_height); # ### Closures # -# We use an area-scaled biharmonic viscosity for grid-scale dissipation and +# We use an area-scaled biharmonic viscosity for grid-scale dissipation, +# an isopycnal symmetric diffusivity for lateral tracer mixing, and # the CATKE parameterization for vertical mixing at the upper-ocean boundary layer. +using Oceananigans.TurbulenceClosures: IsopycnalSkewSymmetricDiffusivity + horizontal_viscosity = area_scaled_biharmonic_viscosity(timescale=15days) +isopycnal_diffusivity = IsopycnalSkewSymmetricDiffusivity(κ_skew=0, κ_symmetric=1e3) vertical_mixing = NumericalEarth.Oceans.default_ocean_closure() # ### Ocean simulation @@ -62,7 +66,7 @@ momentum_advection = WENOVectorInvariant(order=5) tracer_advection = WENO(order=5) ocean = ocean_simulation(grid; momentum_advection, tracer_advection, free_surface, - closure=(horizontal_viscosity, vertical_mixing)) + closure=(horizontal_viscosity, isopycnal_diffusivity, vertical_mixing)) @info "We've built an ocean simulation with model:" @show ocean.model