Coriolis

The Coriolis option determines whether the fluid experiences the effect of the Coriolis force, or rotation. Several Coriolis approximations are available, from the simple -plane to full spherical Coriolis.

Coriolis vs. rotation

If you are wondering why this option is called "Coriolis" it is because rotational effects could include the Coriolis and centripetal forces, both of which arise in non-inertial reference frames. But here the model only considers the Coriolis force.

No rotation

By default there is no rotation. This can be made explicit by passing coriolis = nothing to a model constructor.

Traditional -plane

To set up an -plane with, for example, Coriolis parameter  

julia> coriolis = FPlane(f=1e-4)
FPlane{Float64}(f=0.0001)

An -plane can also be specified at some latitude on a spherical planet with a planetary rotation rate. For example, to specify an -plane at a latitude of    on Earth that has a rotation rate of   

julia> coriolis = FPlane(rotation_rate=7.292115e-5, latitude=45)
FPlane{Float64}(f=0.000103126)

in which case the value of is given by .

Coriolis term for constant rotation in a Cartesian coordinate system

One can use ConstantCartesianCoriolis to set up a Coriolis acceleration term where the Coriolis parameter is constant and the rotation axis is arbitrary. For example, with    ,

julia> coriolis = ConstantCartesianCoriolis(fx=0, fy=2e-4, fz=1e-4)
ConstantCartesianCoriolis{Float64}: fx = 0.00e+00, fy = 2.00e-04, fz = 1.00e-04

Or alternatively, the same result can be achieved by specifying the magnitude of the Coriolis frequency f and the rotation_axis. So another way to get a Coriolis acceleration with the same values is:

julia> using LinearAlgebra: norm

julia> (fx, fy, fz) = (0, 2e-4, 1e-4)
(0, 0.0002, 0.0001)

julia> f = norm((fx, fy, fz))
0.00022360679774997898

julia> rotation_axis = (fx, fy, fz) ./ f # rotation_axis has to be a unit vector
(0.0, 0.8944271909999159, 0.4472135954999579)

julia> coriolis = ConstantCartesianCoriolis(; f, rotation_axis)
ConstantCartesianCoriolis{Float64}: fx = 0.00e+00, fy = 2.00e-04, fz = 1.00e-04

An -plane with non-traditional Coriolis terms can also be specified at some latitude on a spherical planet with a planetary rotation rate. For example, to specify an -plane at a latitude of    on Earth which has a rotation rate of   

julia> coriolis = ConstantCartesianCoriolis(rotation_rate=7.292115e-5, latitude=45)
ConstantCartesianCoriolis{Float64}: fx = 0.00e+00, fy = 1.03e-04, fz = 1.03e-04

in which case   and  .

Traditional -plane

To set up a -plane the background rotation rate and the parameter must be specified. For example, a -plane with   and    can be set up with

julia> coriolis = BetaPlane(f₀=1e-4, β=1.5e-11)
BetaPlane{Float64}(f₀=0.0001, β=1.5e-11)

Alternatively, a -plane can also be set up at some latitude on a spherical planet with a planetary rotation rate and planetary radius. For example, to specify a -plane at a latitude of   on Earth which has a rotation rate of    and a radius of  

julia> coriolis = BetaPlane(rotation_rate=7.292115e-5, latitude=-10, radius=6371e3)
BetaPlane{Float64}(f₀=-2.53252e-5, β=2.25438e-11)

in which case   and  .

Non-traditional -plane

A non-traditional -plane requires either 5 parameters (by default Earth's radius and rotation rate are used):

julia> NonTraditionalBetaPlane(fz=1e-4, fy=2e-4, β=4e-11, γ=-8e-11)
NonTraditionalBetaPlane{Float64}(fz = 1.00e-04, fy = 2.00e-04, β = 4.00e-11, γ = -8.00e-11, R = 6.37e+06)

or the rotation rate, radius, and latitude:

julia> NonTraditionalBetaPlane(rotation_rate=5.31e-5, radius=252.1e3, latitude=10)
NonTraditionalBetaPlane{Float64}(fz = 1.84e-05, fy = 1.05e-04, β = 4.15e-10, γ = -1.46e-10, R = 2.52e+05)

Spherical Coriolis

For simulations on latitude-longitude grids, use HydrostaticSphericalCoriolis (for hydrostatic models) or SphericalCoriolis (for nonhydrostatic models). These evaluate the Coriolis parameter from the grid's latitude:

julia> coriolis = HydrostaticSphericalCoriolis()
SphericalCoriolis
├─ rotation rate: 7.29e-05 s⁻¹ = 1.00 Ω_Earth
├─ formulation: HydrostaticFormulation
└─ scheme: EnstrophyConserving

A custom rotation rate can be specified:

julia> coriolis = HydrostaticSphericalCoriolis(rotation_rate=1e-4)
SphericalCoriolis
├─ rotation rate: 1.00e-04 s⁻¹ = 1.37 Ω_Earth
├─ formulation: HydrostaticFormulation
└─ scheme: EnstrophyConserving

Discretization schemes

The Coriolis term   requires interpolating velocities on the C-grid, since and are not collocated. The scheme keyword argument controls how this interpolation is performed. Different schemes have different conservation properties and different behavior near immersed boundaries.

Five schemes are available:

Scheme	Conserves	Immersed boundary correction
EnstrophyConserving()	Potential enstrophy	No
EnergyConserving()	Kinetic energy	No
TriadScheme()	Both energy and enstrophy	No
ActiveWeightedEnstrophyConserving()	Potential enstrophy	Yes
ActiveWeightedEnergyConserving()	Kinetic energy	Yes

The default scheme is EnstrophyConserving() for HydrostaticSphericalCoriolis, FPlane, and BetaPlane.

Selecting a scheme

Pass the scheme keyword to any Coriolis constructor:

julia> using Oceananigans.Advection: EnergyConserving, EnstrophyConserving

julia> coriolis = FPlane(f=1e-4, scheme=EnergyConserving())
FPlane{Float64}(f=0.0001)

julia> coriolis = FPlane(f=1e-4, scheme=EnstrophyConserving())
FPlane{Float64}(f=0.0001)

julia> coriolis = HydrostaticSphericalCoriolis(scheme=EnergyConserving())
SphericalCoriolis
├─ rotation rate: 7.29e-05 s⁻¹ = 1.00 Ω_Earth
├─ formulation: HydrostaticFormulation
└─ scheme: EnergyConserving

Active-weighted schemes for immersed boundaries

When using immersed boundaries, the standard 4-point averaging of velocities in the Coriolis term includes masked (land) points where velocity is zero. This dilutes the Coriolis force near boundaries, creating spurious numerical boundary layers (Jamart and Ozer, 1986).

The ActiveWeightedEnstrophyConserving and ActiveWeightedEnergyConserving schemes correct this by dividing the interpolated result by the number of active (non-masked) nodes in the stencil, rather than the full stencil size:

julia> coriolis = FPlane(f=1e-4, scheme=Oceananigans.Coriolis.ActiveWeightedEnstrophyConserving())
FPlane{Float64}(f=0.0001)

julia> coriolis = HydrostaticSphericalCoriolis(scheme=Oceananigans.Coriolis.ActiveWeightedEnergyConserving())
SphericalCoriolis
├─ rotation rate: 7.29e-05 s⁻¹ = 1.00 Ω_Earth
├─ formulation: HydrostaticFormulation
└─ scheme: ActiveWeightedEnergyConserving

When to use active-weighted schemes

The active-weighted correction can reduce spurious numerical boundary layers along simple, flat immersed boundaries. However, for complex topography (narrow passages, sharp capes, jagged coastlines) or large values, the amplification factor can inject energy and produce grid-scale checkerboard artifacts along coastlines. For this reason, the standard (non-active-weighted) schemes are the default. Users should test the active-weighted schemes carefully before adopting them in production simulations.

Triad (Energy- and Enstrophy-Conserving) scheme

The TriadScheme scheme is based on the triad formulation of Arakawa and Lamb (1981). It uses a 12-point stencil that conserves both kinetic energy and potential enstrophy in the limit of horizontally non-divergent flow:

julia> coriolis = FPlane(f=1e-4, scheme=Oceananigans.Coriolis.TriadScheme())
FPlane{Float64}(f=0.0001)

The default scheme for FPlane, BetaPlane, and HydrostaticSphericalCoriolis is EnstrophyConserving.