Oceanography

Oceanography: the thermohaline circulation and the wind-driven gyres.

QuasiGeostrophicSiteElement(label, visual=None, *, zeta=0.0, beta=1.0, Rd=1.0, nu=0.0001, curl=0.0)

Bases: Element

Site of a quasi-geostrophic ocean basin, the relative potential vorticity of the water column at one position.

The site holds the relative potential vorticity zeta = laplacian(psi) - psi / Rd**2, the potential vorticity without the planetary part beta y, with the planetary vorticity gradient beta, the deformation radius Rd, the viscosity nu, and the curl of the wind stress curl at the site.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the basin.	required
visual	VisualSpec | None	How the site is drawn, or None for default_visual_.	None
zeta	Annotated[float, StateSpec()]	The initial relative potential vorticity.	0.0
beta	Annotated[float, ParameterSpec(domain=open(0, oo))]	The gradient of the planetary vorticity.	1.0
Rd	Annotated[float, ParameterSpec(domain=open(0, oo))]	The deformation radius.	1.0
nu	Annotated[float, ParameterSpec(domain=open(0, oo))]	The viscosity.	0.0001
curl	Annotated[float, ParameterSpec()]	The curl of the wind stress.	0.0

Examples:

>>> from skcomplex.models.physics.oceanography import QuasiGeostrophicSiteElement
>>> QuasiGeostrophicSiteElement('site').beta
1.0

Source code in src/skcomplex/models/physics/oceanography/_quasi_geostrophic_site.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    zeta: Annotated[float, StateSpec()] = 0.0,
    beta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    Rd: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    nu: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1e-4,
    curl: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.zeta = zeta
    self.beta = beta
    self.Rd = Rd
    self.nu = nu
    self.curl = curl

StommelBoxesElement(label, visual=None, *, T=3.0, S=0.0, eta1=3.0, eta2=1.2, eta3=0.3)

Bases: Element

Ocean of Stommel's two boxes, the differences of temperature and salinity between the equator and the pole.

Both differences relax to their forcing and are mixed by a flow proportional to the density difference |T - S|: dT/dt = eta1 - T (1 + |T - S|) and dS/dt = eta2 - S (eta3 + |T - S|).

Parameters:

Name	Type	Description	Default
label	str	The label of the ocean.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
T	Annotated[float, StateSpec()]	The initial temperature difference.	3.0
S	Annotated[float, StateSpec()]	The initial salinity difference.	0.0
eta1	Annotated[float, ParameterSpec()]	The thermal forcing.	3.0
eta2	Annotated[float, ParameterSpec()]	The freshwater forcing.	1.2
eta3	Annotated[float, ParameterSpec()]	The ratio of the relaxation times of temperature and salinity.	0.3

Examples:

>>> from skcomplex.models.physics.oceanography import StommelBoxesElement
>>> StommelBoxesElement('ocean').eta2
1.2

Source code in src/skcomplex/models/physics/oceanography/_thermohaline.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    T: Annotated[float, StateSpec()] = 3.0,
    S: Annotated[float, StateSpec()] = 0.0,
    eta1: Annotated[float, ParameterSpec()] = 3.0,
    eta2: Annotated[float, ParameterSpec()] = 1.2,
    eta3: Annotated[float, ParameterSpec()] = 0.3,
) -> None:
    super().__init__(label, visual=visual)
    self.T = T
    self.S = S
    self.eta1 = eta1
    self.eta2 = eta2
    self.eta3 = eta3

double_gyre(n_sites=64, *, length=1.0, beta=1.0, Rd=1.0, nu=0.000125, tau0=1e-05)

Return the quasi-geostrophic double gyre, the wind-driven circulation of a closed ocean basin.

The potential vorticity q = laplacian(psi) - psi / Rd**2 + beta y of a square basin of side length is carried by the geostrophic velocity (-psi_y, psi_x) of the streamfunction psi, driven by the curl of the wind stress -tau0 cos(2 pi y / length), and damped by the viscosity nu. The state is the relative potential vorticity zeta = q - beta y, held at zero on the walls for free slip, so that

zeta_t + divergence(u zeta) = -beta psi_x - tau0 (2 pi / length) sin(2 pi y / length) + nu laplacian(zeta),

with psi = inverse_helmholtz(zeta, 1 / Rd, ('fixed', 0)). In the viscous steady regime the interior follows the Sverdrup balance beta psi_x = curl(tau), and the return flow closes in Munk western boundary layers of width (nu / beta) ** (1 / 3), so the streamfunction forms a subtropical and a subpolar gyre of opposite sign whose extrema lie close to the western wall. The defaults keep the wind weak, so the boundary layers stay viscous rather than inertial, and the circulation is steady after a few hundred units of time. The catalog entry streams the streamfunction as the observable psi.

Parameters:

Name	Type	Description	Default
n_sites	int	The number of sites along each side of the basin.	64
length	float	The side of the basin.	1.0
beta	float	The gradient of the planetary vorticity.	1.0
Rd	float	The deformation radius.	1.0
nu	float	The viscosity.	0.000125
tau0	float	The amplitude of the wind stress.	1e-05

Returns:

Name	Type	Description
system	System	The model, with the relative potential vorticity basin.zeta.

Examples:

>>> from skcomplex.models.physics.oceanography import double_gyre
>>> double_gyre(16).elements[0].values['beta']
1.0

Source code in src/skcomplex/models/physics/oceanography/_gyres.py

@register_model(
    field='physics/oceanography',
    tags=('quasi-geostrophic', 'gyre', 'western boundary current', 'field', 'inversion', 'pde'),
    references=('Munk 1950', 'arXiv:1910.00583'),
    defaults={'until': 400.0, 'dt': 0.4, 'solver': RK4Solver(dt=0.4)},
    observables=lambda system: [ObservableSpec('psi', _streamfunction, group='basin', reduce=_sites)],
)
def double_gyre(
    n_sites: int = 64,
    *,
    length: float = 1.0,
    beta: float = 1.0,
    Rd: float = 1.0,
    nu: float = 1.25e-4,
    tau0: float = 1e-5,
) -> System:
    """Return the quasi-geostrophic double gyre, the wind-driven circulation of a closed ocean basin.

    The potential vorticity `q = laplacian(psi) - psi / Rd**2 + beta y` of a square basin of side `length` is carried by
    the geostrophic velocity `(-psi_y, psi_x)` of the streamfunction `psi`, driven by the curl of the wind stress
    `-tau0 cos(2 pi y / length)`, and damped by the viscosity `nu`. The state is the relative potential vorticity
    `zeta = q - beta y`, held at zero on the walls for free slip, so that

        zeta_t + divergence(u zeta) = -beta psi_x - tau0 (2 pi / length) sin(2 pi y / length) + nu laplacian(zeta),

    with `psi = inverse_helmholtz(zeta, 1 / Rd, ('fixed', 0))`. In the viscous steady regime the interior follows the
    Sverdrup balance `beta psi_x = curl(tau)`, and the return flow closes in Munk western boundary layers of width
    `(nu / beta) ** (1 / 3)`, so the streamfunction forms a subtropical and a subpolar gyre of opposite sign whose
    extrema lie close to the western wall. The defaults keep the wind weak, so the boundary layers stay viscous rather
    than inertial, and the circulation is steady after a few hundred units of time. The catalog entry streams the
    streamfunction as the observable `psi`.

    Args:
        n_sites:
            The number of sites along each side of the basin.

        length:
            The side of the basin.

        beta:
            The gradient of the planetary vorticity.

        Rd:
            The deformation radius.

        nu:
            The viscosity.

        tau0:
            The amplitude of the wind stress.

    Returns:
        system:
            The model, with the relative potential vorticity `basin.zeta`.

    Examples:
        >>> from skcomplex.models.physics.oceanography import double_gyre
        >>> double_gyre(16).elements[0].values['beta']
        1.0
    """
    spacing = length / n_sites
    latitudes = (np.arange(n_sites) + 0.5) * spacing
    curl = -tau0 * (2 * np.pi / length) * np.sin(2 * np.pi * latitudes / length)
    values = {
        'zeta': 0.0,
        'beta': beta,
        'Rd': Rd,
        'nu': nu,
        'curl': np.broadcast_to(curl, (n_sites, n_sites)).copy(),
    }
    basin = FieldElement(
        'basin',
        element=QuasiGeostrophicSiteElement,
        space=LatticeSpace((n_sites, n_sites), spacing=spacing, boundary='fixed'),
        boundary={'zeta': ('fixed', 0.0)},
        values=values,
    )
    transport = AdvectionInteraction(
        'transport',
        target='basin',
        state='zeta',
        drift=_velocity,
        source=_beta_drift,
        limiter=None,
        riemann=None,
    )
    wind = RateInteraction('wind', target='basin', state='zeta', rate=_wind)
    viscosity = DiffusionInteraction('viscosity', target='basin', state='zeta', coefficient=_viscosity)
    return System(elements=[basin], interactions=[transport, wind, viscosity])

stommel(T=3.0, S=0.0, *, eta1=3.0, eta2=1.2, eta3=0.3)

Return Stommel's two-box thermohaline circulation, bistable between a thermal and a haline mode.

The dimensionless differences of temperature T and salinity S between an equatorial and a polar box relax to their forcing and are mixed by a flow proportional to the density difference |T - S|: dT/dt = eta1 - T (1 + |T - S|) and dS/dt = eta2 - S (eta3 + |T - S|). For the default parameters a strong circulation driven by temperature and a weak one driven by salinity are both stable.

Parameters:

Name	Type	Description	Default
T	float	The initial temperature difference.	3.0
S	float	The initial salinity difference.	0.0
eta1	float	The thermal forcing.	3.0
eta2	float	The freshwater forcing.	1.2
eta3	float	The ratio of the relaxation times of temperature and salinity.	0.3

Returns:

Name	Type	Description
system	System	The model, with the differences ocean.T and ocean.S.

Examples:

>>> from skcomplex.models.physics.oceanography import stommel
>>> stommel(eta2=1.0).elements[0].eta2
1.0

Source code in src/skcomplex/models/physics/oceanography/_thermohaline.py

@register_model(
    field='physics/oceanography',
    tags=('ode', 'bistability', 'thermohaline circulation', 'tipping'),
    references=('Stommel 1961', 'Cessi 1994', 'Ghil and Lucarini 2020'),
    defaults={'until': 50.0, 'dt': 0.1},
)
def stommel(
    T: float = 3.0,
    S: float = 0.0,
    *,
    eta1: float = 3.0,
    eta2: float = 1.2,
    eta3: float = 0.3,
) -> System:
    """Return Stommel's two-box thermohaline circulation, bistable between a thermal and a haline mode.

    The dimensionless differences of temperature `T` and salinity `S` between an equatorial and a polar box relax to
    their forcing and are mixed by a flow proportional to the density difference `|T - S|`:
    `dT/dt = eta1 - T (1 + |T - S|)` and `dS/dt = eta2 - S (eta3 + |T - S|)`. For the default parameters a strong
    circulation driven by temperature and a weak one driven by salinity are both stable.

    Args:
        T:
            The initial temperature difference.

        S:
            The initial salinity difference.

        eta1:
            The thermal forcing.

        eta2:
            The freshwater forcing.

        eta3:
            The ratio of the relaxation times of temperature and salinity.

    Returns:
        system:
            The model, with the differences `ocean.T` and `ocean.S`.

    Examples:
        >>> from skcomplex.models.physics.oceanography import stommel
        >>> stommel(eta2=1.0).elements[0].eta2
        1.0
    """
    return System(elements=[StommelBoxesElement('ocean', T=T, S=S, eta1=eta1, eta2=eta2, eta3=eta3)])