geophysix: physical constants, Earth geometry and atmospheric state for JAX

geophysix is the stack’s library of physical facts about the Earth: its constants, its size and shape, its rotation, the thermodynamics of its air and the equation of state of its sea water, and the units a gas column is reported in. It sits at GeoNumerics tier 0, imports nothing from the stack, and gives every library above it one value for each constant where today there are several. Status: proposed, at design-sketch depth.

At a glance¶

geophysix sits at the bottom of GeoNumerics because everything above it needs the same Earth. Core cannot import it, so Core functions that need the Earth’s size or shape take it as an argument, and geophysix supplies the value. This page is the design; core types, algorithms and the code survey are in geophysix: design details.

The problem¶

Turning a satellite methane plume into an emission rate needs a handful of physical facts about the Earth and its air, and each library that needs one writes its own copy:

The copies do not agree. The same constant appears with different values across the stack, so two libraries given the same inputs can return different answers:

The Boltzmann constant and the standard atmosphere are written out four times (radtrax, plumax twice, xrops), and three repositories hold the same constants.py, copied from one to the next (finitevolX, somax and sesmax-core). Each difference is small, but each makes a cross-library test fail at a tolerance it should pass, and none of them is a choice anyone made.

One Earth, passed as data¶

Everything geophysix holds is a closed-form function of a few constants, and the design rests on that. The constants live in one record. Every function takes the record, or the one constant it needs, as an argument whose default is that record. Nothing reads a module-level global inside a traced function.

The ellipsoid. The Earth is an ellipsoid of revolution with semi-major axis a and flattening f; WGS84 has a = 6,378,137 m and f = 1 / 298.257223563. With e² = f (2 − f) the squared eccentricity, a point at geodetic longitude λ, latitude φ and height h sits at

The area of a cell between latitudes φ₁ and φ₂, spanning Δλ of longitude, has a closed form, with b = a (1 − f):

On a sphere (e → 0) it reduces to A = R² Δλ (sin φ₂ − sin φ₁). A sphere is the same record with f = 0 and a chosen mean radius, so every function works on both.

Rotation. With Ω the rotation rate and R the radius, the Coriolis parameter and its meridional gradient are f = 2Ω sin φ and β = 2Ω cos φ / R. A β-plane linearises them about a reference latitude: f(y) = f₀ + β₀ (y − y₀).

The air. Dry air is an ideal gas, p = ρ R_d T_v with T_v the virtual temperature, and its number density is n = p / (k_B T). Potential temperature is θ = T (p₀ / p)^κ with κ = R_d / c_p, and the thickness of a layer between two pressure levels is

Sea water. Density from temperature T, salinity S and pressure: the linear equation of state ρ = ρ₀ (1 − α (T − T₀) + β_S (S − S₀)) first, TEOS-10 later.

Columns. A gas column is reported in several units. With X the enhancement in ppm·m, V_m = R T / p the molar volume of air and M the molar mass of the gas:

None of these needs a grid, a carrier or data. They are plain functions of arrays and constants, which is why they sit at the bottom of GeoNumerics.

Modules and boundaries¶

Six modules, one per kind of fact. geophysix owns closed-form physics at a point; anything on a grid, a field or a carrier belongs to the libraries above it.

What geophysix owns, and what it hands off¶

The test for every name: is it a closed-form fact about the Earth, its air or its sea water, evaluated at a point? Then geophysix owns it.

What moves in¶

In use¶

The methane example once more: pixel areas on the ellipsoid, the enhancement converted to a mass, the dry-air column for XCH₄, the heights of the wind levels, and the Coriolis field a transport model needs. Every geophysix name is proposed; each is a function of arrays with the Earth passed in.

# Shapes: Ny × Nx plume pixels on a lon/lat grid; Z pressure levels

import einx
import geophysix as gpx

# proposed (geophysix)
# WGS84: a = 6,378,137 m, f = 1 / 298.257223563
earth = gpx.earth.WGS84

# --- pixel areas, on the ellipsoid ---
# proposed (geophysix)
# Aᵢⱼ = (b² Δλⱼ / 2) [q(φ)] from φᵢ to φᵢ₊₁
area = gpx.earth.cell_area(
    lat_edges, lon_edges, earth=earth
)  # (Ny+1,), (Nx+1,) → (Ny, Nx)

# --- the enhancement as a mass ---
# proposed (geophysix)
# Ω_kg = M_CH₄ · 10⁻⁶ X / V_m, with V_m = R T / p
column = gpx.columns.convert(  # (Ny, Nx) → (Ny, Nx)
    enhancement,
    units_in="ppm_m",
    units_out="kg_m2",
    gas="CH4",
    pressure=surface_pressure,
    temperature=surface_temperature,
)
# IME = Σᵢⱼ Ω_kg Aᵢⱼ
mass = einx.dot("y x, y x ->", column, area)  # → ()

# --- dry-air column, for XCH₄ ---
# proposed (geophysix)
# N_dry = N_A (p_s / g − W) / M_dry
dry_air = gpx.columns.dry_air_column(  # (Ny, Nx) → (Ny, Nx)
    surface_pressure, water_vapour=water_vapour_column
)

# --- heights of the wind levels ---
# proposed (geophysix)
# zₖ₊₁ − zₖ = (R_d T̄_v / g) ln(pₖ / pₖ₊₁)
thickness = gpx.thermo.hypsometric_thickness(  # (Z,), (Z,) → (Z−1,)
    pressure_levels, virtual_temperature
)

# --- the Coriolis field of a transport model ---
# proposed (geophysix)
# f = 2 Ω sin φ
coriolis = gpx.rotation.coriolis(lat_centres)  # (Ny,) → (Ny,)

Each call takes its constants from the defaults, so a caller that passes nothing gets WGS84, CODATA 2018 and standard gravity, and a caller that wants a sphere of radius 6,371,000 m passes earth=gpx.earth.sphere(6_371_000.0). The sum for the mass is shown with einx; inside a pipeline it is integrix.rule with these areas as weights.

Integration across the stack¶

Twelve libraries would import geophysix, across three layers. For most the change is an import and the end of a local constants.py; for xrops it is a thin DataArray layer over geophysix functions.

The change at each site is an import. The constants module of finitevolX, somax and sesmax-core is the clearest case:

# Today:
# today (finitevolX)
from finitevolx._src.utils.constants import GRAVITY, OMEGA, R_EARTH

# After:
# proposed (geophysix)
from geophysix.constants import STANDARD

# proposed (geophysix)
from geophysix.earth import SPHERE

The values change with the import: R_EARTH = 6,371,200 m becomes the record’s mean radius. That is the one place the switch moves a number, so each consumer’s switch lands with its own tolerance check (decision 4).

Decisions and open questions¶

Two decided on the layer map, four leaning, three open. Each entry gives the position first, then the trade-off.

1. Layer: GeoNumerics tier 0¶

Decided (5 Nov 2026, on the layer map). geophysix imports nothing from the stack, and every GeoNumerics library, three GeoStack packages and three models import it. Core sits below and cannot import it: a Core function that needs the Earth takes the radius, or a and f, as an argument.

2. Scope: the Earth, its air and its sea water¶

Decided (5 Nov 2026, on the layer map). Constants, the ellipsoid and sphere metrics, Coriolis and β-planes, geodesics, atmospheric thermodynamics, seawater equations of state and column units. Model scales and non-dimensionalisation stay in the models; diagnostics of fields go to diagnostix.

3. Constants as records passed in¶

Leaning: frozen records, passed as arguments with defaults. gpx.constants.STANDARD and gpx.earth.WGS84 are immutable records, and every function takes the one it needs as a keyword argument. A model that wants another planet, a sphere or a different convention passes its own record; nothing patches a global. Plain module-level floats would be shorter, but they are the reason three libraries hold three radii today.

4. Which value for each constant¶

Leaning: the published conventions, one each. The IUGG mean radius R₁ = 6,371,008.8 m for the default sphere; standard gravity g₀ = 9.80665 m s⁻²; Ω = 7.292115 × 10⁻⁵ s⁻¹ from WGS84; CODATA 2018 for k_B, N_A and R. finitevolX, somax and sesmax-core use 6,371,200 m today, IGRF’s geomagnetic reference radius rather than a mean Earth radius, so their results move by about 3 × 10⁻⁵ relative when they switch; each switch carries a test of that size, and a model that must reproduce an old run passes the old radius explicitly.

5. Geodesics: the sphere in Core, the ellipsoid here¶

Leaning. Great-circle distance on a sphere of given radius stays where Core needs it (gaussx’s localisation, kernellib’s kernels), with its default radius equal to geophysix’s and a test that says so. Geodesics on the ellipsoid (distance and azimuth, after Karney) are geophysix’s. The private haversines in geopatcher and xrproducts go.

6. Units: SI floats, no unit library¶

Leaning. Every argument and return is in SI, with the unit in the docstring and, where two are common, in the name (pressure_pa). The column module is the one place units are converted, and it does so explicitly by name. A unit library such as pint would not trace under JAX.

7. TEOS-10: port or wrap?¶

Open. xrops calls gsw, which is NumPy and C, so it neither traces nor differentiates. Leaning: port the subset the stack uses (density, α, β_S, conservative temperature, absolute salinity) to JAX, with gsw as the test reference, as integrix does with quadax.

8. Saturation vapour pressure and moist thermodynamics¶

Open. No library in the stack computes it today, and the formulas disagree with one another (Magnus, Goff–Gratch, Murphy–Koop). Leaning: one formula, Murphy–Koop, chosen by an argument, and added only when sesmax-atmos or plumax needs it.

9. What to take from pyproj and Shapely¶

Open, to be investigated. pyproj (PROJ and GeographicLib) and Shapely (GEOS) are C libraries, so nothing they compute traces under jit or differentiates. The candidates to port as maths, not code, are:

• from pyproj: the ellipsoid table; the auxiliary latitudes (geocentric, parametric, authalic, conformal, rectifying); Karney’s direct and inverse geodesic and the geodesic polygon area.

• from Shapely: point-in-polygon masks, planar polygon area and distance to a polygon. These are generic geometry, so they would go to gridax or functax, not here.

The working position: map projections stay with pyproj on the host, through geoarray’s pyproj.CRS, until something has to differentiate through one; pyproj’s Geod is the test reference for whatever is ported. To settle: which consumers need ellipsoidal geodesics and polygon areas under jit (swath footprints and basin-scale budgets are the likely ones; plume-scale work is served by a local tangent plane); how Karney’s iterative inverse behaves under jit and grad, including at coincident and antipodal points; and the licences, since PROJ, GeographicLib and pyproj are MIT and Shapely is BSD-3, but GEOS is LGPL, so nothing is copied from its source.

Roadmap¶

Five waves, each a GitHub issue epic in the geophysix repository once it exists. The constants and the Earth come first: every other module reads them, and they end the most copies.

Waves 2 and 3 can run side by side. Wave 4 is one issue per consumer, filed in that consumer’s repository and linked to the geophysix epic; somax needs none, since sesmax-core replaces it.

Sources¶

• The GeoStack layer map: geophysix’s layer, tier, imports and moves.

• National Geospatial-Intelligence Agency (2014): WGS84’s a, f and Ω.

• CODATA recommended values: k_B, N_A, R.

• Karney (2013): geodesics on the ellipsoid.

• TEOS-10: the thermodynamic equation of seawater, and the gsw reference code.

• Code read on 5 Nov 2026 at geostack’s submodule pins and sesmax main: finitevolX _src/utils/constants.py, operators/diagnostics.py, operators/eos.py; somax _src/constants/; sesmax-core _src/constants.py, _src/scales.py; xrops grad/_src/constants.py, atm/_src/, atm/gas/ch4/_src/column.py, ocn/_src/kinematics.py; geonnax _frames.py; kernellib functional/_geo.py; gaussx _inference/_localization.py; geoops plume/_src/array.py, geom/_src/operators.py; geoproducts himawari/constants.py; radtrax atmosphere/gas.py; plumax coupled/rtm.py, hapi_lut/beers.py, met/wrf.py.