Public API

NumericalEarth.jl

🌎 A framework for simulating the Earth system at all scales with prescribed or prognostic atmosphere, ocean, sea ice and land components, based on Oceananigans.

Overview

NumericalEarth.jl provides infrastructure for running Earth system model components—ocean, atmosphere, sea ice, and others—coupled together or driven by prescribed datasets. The coupling interface is generic: plug in Oceananigans for ocean dynamics, ClimaSeaIce for sea ice, SpeedyWeather or other atmospheric models, or use reanalysis products like JRA55 and ERA5 as prescribed forcing.

The package handles the complexity of component communication: interpolating between grids, computing air-sea fluxes via similarity theory, managing radiative transfer, and synchronizing time-stepping across components with different temporal resolutions.

NumericalEarth.jl also serves as a sandbox for developing and testing interface parameterizations—bulk flux formulations, roughness length models, albedo schemes, and other boundary layer physics—in a modular environment where they can be validated against observations before deployment in production climate models.

Installation instructions

NumericalEarth is a registered Julia package. So to install it,

1. Download Julia (version 1.10 or later).

2. Launch Julia and type

julia> using Pkg

julia> Pkg.add("NumericalEarth")

This installs the latest version that's compatible with your current environment.

Data Wrangling

Running realistic Earth system simulations requires wrangling gigabytes of observational and reanalysis data into formats your model can ingest. NumericalEarth.jl abstracts away this pain. Point the package at a dataset and a target grid, and it handles the downloading, caching, and regridding automatically.

The Metadatum abstraction provides a unified interface: whether you're initializing ocean temperature from reanalysis or prescribing atmospheric boundary conditions, the workflow is the same.

using NumericalEarth
using Dates

# Load temperature from reanalysis on a specific date
T_init = Metadatum(:temperature; date=DateTime(1993, 1, 1), dataset=ECCO2Daily())

# Build a prescribed atmosphere from JRA55
atmosphere = JRA55PrescribedAtmosphere(arch)

A core abstraction: EarthSystemModel

The coupling infrastructure is anchored by EarthSystemModel, which encapsulates the component models—ocean, sea ice, atmosphere—and specifies how they communicate. Each component can be either prognostic (time-stepped by its own dynamics) or prescribed (interpolated from data). The model handles flux computations at interfaces, grid interpolation between components, and synchronized time-stepping.

We conceive of EarthSystemModel as a model in its own right, not just a container for components. This means it works with all the Oceananigans tools you'd use for any other model—run!(simulation), Callback, Checkpointer, output writers, and the rest.

To illustrate, here's a global ocean simulation driven by prescribed atmospheric reanalysis:

using Oceananigans
using Oceananigans.Units
using Dates
using CUDA
import NumericalEarth

arch = GPU()
grid = LatitudeLongitudeGrid(arch,
                             size = (1440, 560, 10),
                             halo = (7, 7, 7),
                             longitude = (0, 360),
                             latitude = (-70, 70),
                             z = (-3000, 0))

bathymetry = NumericalEarth.regrid_bathymetry(grid)
grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bathymetry))

# Build an ocean simulation initialized to the ECCO state estimate
ocean = NumericalEarth.ocean_simulation(grid)
start_date = DateTime(1993, 1, 1)
set!(ocean.model,
     T=NumericalEarth.Metadatum(:temperature; date=start_date, dataset=NumericalEarth.ECCO2Daily()),
     S=NumericalEarth.Metadatum(:salinity;    date=start_date, dataset=NumericalEarth.ECCO2Daily()))

# Couple the ocean to JRA55 atmospheric forcing
atmosphere = NumericalEarth.JRA55PrescribedAtmosphere(arch)
coupled_model = NumericalEarth.OceanOnlyModel(ocean; atmosphere)
simulation = Simulation(coupled_model, Δt=20minutes, stop_time=30days)
run!(simulation)

This simulation achieves approximately 8 simulated years per day of wall time on an Nvidia H100 GPU.

Since ocean.model is an Oceananigans.HydrostaticFreeSurfaceModel, we can leverage Oceananigans features directly. For example, to plot the surface speed at the end of the simulation:

u, v, w = ocean.model.velocities
speed = Field(sqrt(u^2 + v^2))

using GLMakie
heatmap(view(speed, :, :, ocean.model.grid.Nz), colorrange=(0, 0.5), colormap=:magma, nan_color=:lightgray)

which produces

Installation

NumericalEarth is a registered Julia package. To install from a Julia REPL:

julia> using Pkg

julia> Pkg.add("NumericalEarth")

Use Pkg.add(url="https://github.com/NumericalEarth/NumericalEarth.jl.git", rev="main") to install the latest development version.

Citing

If you use NumericalEarth for your research, teaching, or fun 🤩, everyone in our community will be grateful if you give credit by citing the corresponding Zenodo record, e.g.,

Wagner, G. L. et al. (2026). NumericalEarth/NumericalEarth.jl: v0.1.1 (v0.1.1). Zenodo. https://doi.org/10.5281/zenodo.18665405

and also the recent preprint submitted to the Journal of Advances in Modeling Earth Systems that presents an overview of all the things that make Oceananigans unique:

"High-level, high-resolution ocean modeling at all scales with Oceananigans"

by Gregory L. Wagner, Simone Silvestri, Navid C. Constantinou, Ali Ramadhan, Jean-Michel Campin, Chris Hill, Tomas Chor, Jago Strong-Wright, Xin Kai Lee, Francis Poulin, Andre Souza, Keaton J. Burns, Siddhartha Bishnu, John Marshall, and Raffaele Ferrari

arXiv:2502.14148

bibtex

@article{Oceananigans-overview-paper-2025,
  title = {{High-level, high-resolution ocean modeling at all scales with Oceananigans}},
  author = {G. L. Wagner and S. Silvestri and N. C. Constantinou and A. Ramadhan and J.-M. Campin and C. Hill and T. Chor and J. Strong-Wright and X. K. Lee and F. Poulin and A. Souza and K. J. Burns and S. Bishnu and J. Marshall and R. Ferrari},
  journal = {arXiv preprint},
  year = {2025},
  archivePrefix = {arXiv},
  eprint = {2502.14148},
  doi = {10.48550/arXiv.2502.14148},
  notes = {submitted to the Journal of Advances in Modeling Earth Systems},
}

natural_earth_lines(name; scale = 50)

Return (longitudes, latitudes) for the Natural Earth feature layer name (e.g. "admin_1_states_provinces_lines" or "admin_0_boundary_lines_land") as NaN-separated vectors, ready to draw as disjoint segments with a plotting backend, as in lines!(axis, longitudes, latitudes). scale selects the Natural Earth resolution (10, 50, or 110).

Requires NaturalEarth and GeoInterface to be loaded — the method lives in NumericalEarthNaturalEarthExt.

visualize_nested_domain(grid; parent = nothing, kw...)

Return a Makie.Figure mapping the domain of grid — a topography/bathymetry basemap with the domain drawn as a box — for checking a configuration's geography before running. When parent (a BoundingBox or another grid) is given, its region is drawn as a second box, visualizing the nesting; omit it to map a single domain.

Keyword arguments

• parent: a BoundingBox or grid whose region is drawn as the outer (parent) box. Default: nothing.
• padding: margin (degrees) between the outermost region and the map edge. Default: 2.5.
• resolution: basemap sampling resolution (degrees). Default: 1/30.
• dataset: relief dataset for the basemap. Default: ETOPO2022().
• boundaries: draw Natural Earth state/country boundary lines (requires NaturalEarth and GeoInterface to be loaded). Default: true.
• landmarks: label => (λ, φ) pairs marked with stars, e.g. tuple("ARM SGP" => (-97.485, 36.605)). Default: tuple().
• label, parent_label: legend labels for the two domain boxes. Defaults: "grid", "parent".
• title: axis title. Default: "".

Requires a Makie backend (e.g. CairoMakie) to be loaded — the method lives in NumericalEarthMakieExt.

PrescribedAtmosphere(grid, times=[zero(grid)];
                     source = nothing,
                     clock = Clock{Float64}(time = 0),
                     surface_layer_height = 10, # meters
                     boundary_layer_height = 512, # meters
                     thermodynamics_parameters = AtmosphereThermodynamicsParameters(eltype(grid)),
                     velocities              = default_atmosphere_velocities(grid, times),
                     temperature             = default_atmosphere_temperature(grid, times),
                     specific_humidity       = default_atmosphere_specific_humidity(grid, times),
                     microphysical_variables = default_atmosphere_microphysical_variables(grid, times),
                     pressure                = default_atmosphere_pressure(grid, times),
                     tracers                 = NamedTuple(),
                     precipitation_flux      = default_precipitation_flux(grid, times))

Return a prescribed, time-evolving atmospheric state with data on grid at times.

The state holds velocities, temperature, specific_humidity, an optional microphysical_variables NamedTuple (cloud / precipitation species such as (; qᶜˡ, qʳ, qᶜⁱ, qˢ), empty by default), pressure, gas-species tracers (e.g. CO₂; empty by default), and — for a surface atmosphere — a precipitation_flux (a 3D atmosphere defaults to none, carrying precipitation in microphysical_variables instead).

source records what the atmosphere was built from (a dataset object, e.g. ERA5HourlyPressureLevels(); nothing for a hand-built atmosphere). Dataset constructors set it so that aliases like ERA5PrescribedAtmosphere = PrescribedAtmosphere{<:ERA5Dataset} dispatch on provenance.

Surface (2D, (Center, Center, Nothing)) vs 3D ((Center, Center, Center), adding w) default fields are inferred from the grid's vertical size: Nz == 1 builds a surface atmosphere (ocean / sea-ice coupling), Nz > 1 a 3D atmosphere (e.g. a NestedSimulation parent). Pass any field explicitly to override; pass precipitation_flux = nothing to omit it.

PrescribedPrecipitationFlux(; rain=nothing, snow=nothing)
PrescribedPrecipitationFlux(rain, snow)

Container for prescribed precipitation fluxes. Either component may be nothing to indicate that the corresponding precipitation type is not represented by the atmosphere (e.g. rain-only datasets). Used as the precipitation_flux of a PrescribedAtmosphere; downstream callers query the snow component via surface_snowfall_flux so that prognostic atmospheres with or without snow can dispatch on this type as well.

bulk_drag(model; roughness_length = 0.1, von_karman_constant = 0.4)

Return a callback f(simulation) that fills the model's coupling bottom-stress fields (ρτˣ, ρτʸ) with a bulk neutral surface stress ρτ = −ρ Cᵈ |U| U, with a per-column log-law drag coefficient Cᵈ = (κ / ln(z₁/z₀))² evaluated at the first-cell-center height z₁ above the local surface — a stand-in for land/ocean surface coupling. Methods live in atmosphere-model extensions (e.g. NumericalEarthBreezeExt).

density_from_pressure(
    temperature,
    pressure;
    qᵛ,
    qᶜ,
    qⁱ,
    dry_gas_constant,
    vapor_gas_constant
)

Moist-air density ρ = p / (Rᵐ T) on the grid of temperature, with the moist mixture gas constant Rᵐ = qᵈ Rᵈ + qᵛ Rᵛ, qᵈ = 1 − qᵛ − qᶜ − qⁱ — the same EOS as hydrostatic_pressure_from_surface, so a density built from that pressure is mutually consistent. temperature and pressure are Fields on the same grid; the moisture qᵛ, qᶜ, qⁱ are optional Fields (omitting them gives the dry result ρ = p / (Rᵈ T)). Returns a CenterField. Useful to initialize a compressible model's density from an analysis temperature + (hydrostatic) pressure, e.g. set!(model; ρ = density_from_pressure(T, p; …)).

ORCAGrid(arch = CPU(), FT::DataType = Float64;
         dataset,
         halo = (4, 4, 4),
         z = (-6000, 0),
         Nz = 50,
         radius = Oceananigans.defaults.planet_radius,
         with_bathymetry = true,
         active_cells_map = true,
         major_basins = Inf,
         south_rows_to_remove = default_south_rows_to_remove(dataset),
         dir = default_download_directory(dataset))

Construct an OrthogonalSphericalShellGrid with (Periodic, RightFaceFolded, Bounded) topology using coordinate and metric data from a NEMO eORCA mesh_mask file.

The dataset keyword argument specifies which ORCA configuration to use (e.g., ORCAOne(), ORCAQuarter(), or ORCATwelfth()). The mesh mask and bathymetry files are downloaded automatically via the DataWrangling.ORCA metadata interface.

The horizontal grid (including coordinates, scale factors, and areas) is loaded directly from the mesh_mask NetCDF file. If all staggered NEMO fields are present (T, U, V, F points), they are used directly. If only T and F coordinates are available (glamt/gphit/glamf/gphif), staggered coordinates and metrics are reconstructed approximately using Tripolar-style spherical assumptions. The duplicated columns eORCA carries at its east edge for cyclic exchange are dropped, so Nx is the number of distinct columns: 360 for eORCA1, 1440 for eORCA025 and 4320 for eORCA12.

When with_bathymetry = true (the default), the bathymetry is also downloaded and the grid is returned as an ImmersedBoundaryGrid with a GridFittedBottom.

Positional Arguments

• arch: The architecture (e.g., CPU() or GPU()). Default: CPU().
• FT: Floating point type. Default: Float64.

Keyword Arguments

• dataset: The ORCA dataset to use. Default: ORCAOne(). ORCAQuarter() (eORCA025, quarter-degree) and ORCATwelfth() (eORCA12, twelfth-degree) are also supported (eORCA1 data from Zenodo; https://doi.org/10.5281/zenodo.4436658).
• halo: Halo size tuple (Hx, Hy, Hz). Default: (4, 4, 4).
• z: Vertical coordinate specification. Can be a 2-tuple (z_bottom, z_top), an array of z-interfaces, or, e.g., an ExponentialDiscretization. Default: (-6000, 0).
• Nz: Number of vertical levels (only used when z is a 2-tuple). Default: 50.
• radius: Planet radius. Default: Oceananigans.defaults.planet_radius.
• with_bathymetry: If true, download the bathymetry and return an ImmersedBoundaryGrid with GridFittedBottom. Default: true.
• active_cells_map: If true and with_bathymetry = true, build an active cells map for efficient kernel execution over wet cells only. Default: true.
• major_basins: Number of independent connected ocean basins to retain via remove_minor_basins!. Basins are removed from smallest to largest; major_basins = 1 keeps only the largest. Default: Inf (keep all basins).
• south_rows_to_remove: Number of southern rows to remove from the eORCA grid. The "extended" eORCA grid contains degenerate padding rows near Antarctica that are entirely land. Removing them reduces memory usage and computation.
• dir: Directory to store and look up ORCA files (mesh_mask and bathymetry). Defaults to the dataset scratch cache via default_download_directory(dataset).

regrid_bathymetry(target_grid, metadata;
                  height_above_water = nothing,
                  minimum_depth = 0,
                  major_basins = 1,
                  interpolation_passes = 1,
                  cache = true,
                  overwrite_cache = false)

Return bathymetry that corresponds to metadata onto target_grid.

Arguments

• target_grid: grid to interpolate the bathymetry onto.

Keyword Arguments

• height_above_water: limits the maximum height of above-water topography (where $h > 0$) before interpolating. Default: nothing, which implies that the original topography is retained.

• minimum_depth: minimum depth for the shallow regions, defined as a positive value. h > - minimum_depth is considered land. Default: 0.

• interpolation_passes: regridding/interpolation passes. The bathymetry is interpolated in interpolation_passes - 1 intermediate steps. The more the interpolation steps the smoother the final bathymetry becomes.

  Example

  Interpolating from a 400 x 200 grid to a 100 x 100 grid in 4 passes involves:

  • 400 x 200 → 325 x 175
  • 325 x 175 → 250 x 150
  • 250 x 150 → 175 x 125
  • 175 x 125 → 100 x 100

  If coarsening the original grid, linear interpolation in passes is equivalent to applying a smoothing filter, with more passes increasing the strength of the filter. If refining the original grid, additional passes do not help and no intermediate steps are performed.

• major_basins: Number of "independent major basins", or fluid regions fully encompassed by land, that are retained by remove_minor_basins!. Basins are removed by order of size: the smallest basins are removed first. major_basins = 1 retains only the largest basin. If Inf then no basins are removed. Default: 1.

• cache: If true (default), caches the regridded bathymetry to disk and reuses it on subsequent calls with the same grid, parameters, and dataset file; a re-download of the dataset invalidates the entry. If false, the cache is disabled entirely: nothing is read or written.

• overwrite_cache: If true, skip the cache lookup and overwrite the entry with a freshly regridded result. Default: false.

regrid_bathymetry(target_grid; dataset=ETOPO2022(), cache=true, kw...)

Regrid bathymetry from dataset onto target_grid. Default: dataset = ETOPO2022().

regrid_topography(target_grid, metadata; kw...)

Land surface elevation regridded onto target_grid from metadata, the positive counterpart of regrid_bathymetry. Use this form for region-windowed datasets such as GLO30(), whose metadata carries a BoundingBox region.

regrid_topography(target_grid; dataset = ETOPO2022(), kw...)

Land surface elevation (m, ≥ 0) regridded onto target_grid — the topographic counterpart of regrid_bathymetry for land applications. Returns the positive part of the dataset's bottom height (the elevation over land), with ocean clamped to sea level (0). Accepts the same regridding keywords (interpolation_passes, etc.); there is no depth/minimum_depth notion.

smooth_topography!(elevation; passes) -> Any

Smooth the two-dimensional field elevation in place with passes applications of a two-dimensional binomial (1-2-1) filter, and return it.

Each pass multiplies the amplitude of a Fourier mode with wavenumbers (kˣ, kʸ) by (1 + cos(kˣ Δx)) (1 + cos(kʸ Δy)) / 4: the two-grid-length mode is annihilated exactly, a three-grid-length mode is damped to 1/4 per pass, and well-resolved scales pass through nearly unchanged. Use it to remove the grid-scale orographic roughness left by point-sampled regridding (see regrid_topography), which on a terrain-following coordinate otherwise excites standing grid-scale noise in the near-surface flow.

At bounded edges, the nearest interior elevation is replicated into the stencil so that smoothing does not pull topography toward zero and preserves a constant elevation field.

using Oceananigans
using NumericalEarth

grid = RectilinearGrid(size=(8, 8), extent=(1, 1), topology=(Periodic, Periodic, Flat))
elevation = Field{Center, Center, Nothing}(grid)
set!(elevation, (x, y) -> sinpi(8x) * sinpi(8y))
smooth_topography!(elevation, passes=1)

# output
8×8×1 Field{Center, Center, Nothing} reduced over dims = (3,) on RectilinearGrid on CPU
├── grid: 8×8×1 RectilinearGrid{Float64, Periodic, Periodic, Flat} on CPU with 3×3×0 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: Nothing, top: Nothing, immersed: Nothing
└── data: 14×14×1 OffsetArray(::Array{Float64, 3}, -2:11, -2:11, 1:1) with eltype Float64 with indices -2:11×-2:11×1:1
    └── max=0.0, min=0.0, mean=0.0

Incorporate various datasets to be used for bathymetry, initialization, forcing, restoring, or validation.

BoundingBox(grid; padding = 0)

Create a BoundingBox spanning the horizontal extent of grid, widened on every side by padding (degrees). For a distributed grid the box spans the global domain, identically on every rank.

BoundingBox(; longitude=nothing, latitude=nothing, z=nothing)

Create a bounding box with latitude, longitude, and z bounds on the sphere.

Column(longitude, latitude; z=nothing, interpolation=Linear())

Create a column region at a single horizontal point (longitude, latitude). When used as a Metadata region, native_grid returns a single-column RectilinearGrid and location reduces horizontal dimensions to Nothing.

Keyword Arguments

• z: depth range tuple (z_bottom, z_top) for restricting downloads (used by CopernicusMarine/GLORYS). Default: nothing (full depth).
• interpolation: method for extracting data from the surrounding grid cells. Linear() (default) bilinearly interpolates to the exact point; Nearest() selects the closest grid cell.

DatasetRestoring(metadata::Metadata,
                 arch_or_grid = CPU();
                 rate,
                 mask = 1,
                 time_indices_in_memory = 2, # Not more than this if we want to use GPU!
                 time_indexing = Cyclical(),
                 inpainting = NearestNeighborInpainting(Inf),
                 cache_inpainted_data = true)

Return a forcing term that restores to data from a dataset. The restoring is applied as a forcing on the right hand side of the evolution equations, calculated as:

\[F_ψ = r μ (ψ_{dataset} - ψ)\]

where $μ$ is the mask, $r$ is the restoring rate, $ψ$ is the simulation variable, and $ψ_{dataset}$ is the dataset variable that is linearly interpolated in space and time from the dataset of choice to the simulation grid and time.

Arguments

• metadata: The medatada for a dataset variable to restore. Choices for variables include:

  • :temperature,
  • :salinity,
  • :u_velocity,
  • :v_velocity,
  • :sea_ice_thickness,
  • :sea_ice_area_fraction,
  • :dissolved_inorganic_carbon,
  • :alkalinity,
  • :nitrate,
  • :phosphate,
  • :dissolved_organic_phosphorus,
  • :particulate_organic_phosphorus,
  • :dissolved_iron,
  • :dissolved_silicate,
  • :dissolved_oxygen.

• arch_or_grid: Either the architecture of the simulation, or a grid on which the data is pre-interpolated when loaded. If an architecture is provided, such as arch_or_grid = CPU() or arch_or_grid = GPU(), data is interpolated on-the-fly when the forcing tendency is computed. Default: CPU().

Keyword Arguments

• rate: The restoring rate, i.e., the inverse of the restoring timescale (in s⁻¹).

• mask: The mask value. Can be a function of (x, y, z, time), an array, or a number.

• time_indices_in_memory: The number of time indices to keep in memory. The number is chosen based on a trade-off between increased performance (more indices in memory) and reduced memory footprint (fewer indices in memory). Default: 2.

• time_indexing: The time indexing scheme for the field time series. Default: Cyclical().

• inpainting: inpainting algorithm, see inpaint_mask!. Default: NearestNeighborInpainting(Inf).

• cache_inpainted_data: If true, the data is cached to disk after inpainting for later retrieving. Default: true.

LinearlyTaperedPolarMask(; northern = (70,   75),
                           southern = (-75, -70),
                           z = (-20, 0))

Build a mask that is linearly tapered in latitude between the northern and southern edges. The mask is constant in depth between the z and equals zero everywhere else. The mask is limited to lie between (0, 1). The mask has the following functional form:

n = 1 / (northern[2] - northern[1]) * (φ - northern[1])
s = 1 / (southern[1] - southern[2]) * (φ - southern[2])

valid_depth = (z[1] < z < z[2])

mask = valid_depth ? clamp(max(n, s), 0, 1) : 0

Metadata(variable_name;
         dataset,
         dates = all_dates(dataset, variable_name),
         dir = default_download_directory(dataset),
         region = default_region(dataset),
         filename = nothing,
         start_date = nothing,
         end_date = nothing)

Metadata holding a specific dataset information.

Argument

• variable_name: a symbol representing the name of the variable (for example, :temperature, :salinity, :u_velocity, etc)

Keyword Arguments

• dataset: Supported datasets are returned by supported_datasets.

• dates: The dates of the dataset (Dates.AbstractDateTime or CFTime.AbstractCFDateTime). Note that dates can either be a range or a vector of dates, representing a time-series, or a (start_date, end_date) tuple, which expands to the dataset's native dates in that window (the cadence is the dataset's own). For a single date, use Metadatum.

• start_date: If dates = nothing, we can prescribe the first date of metadata as a date (Dates.AbstractDateTime or CFTime.AbstractCFDateTime). If outside the date range of the dataset, the first allowable date is chosen. Default: nothing.

• end_date: If dates = nothing, we can prescribe the last date of metadata as a date (Dates.AbstractDateTime or CFTime.AbstractCFDateTime). If outside the date range of the dataset, the last allowable date is chosen. Default: nothing.

• region: Specifies the spatial region of the dataset. Can be a BoundingBox for a rectangular region, a Column for a single horizontal location, or nothing for the full domain.

• filename: The filename(s) for the dataset. If nothing, the filename is computed from the dataset type. Can be a String (single file for all dates) or a DatewiseFilename (one file per date).

• dir: The directory where the dataset is stored.

MetadataSet(variable_names::Symbol...;
            dataset,
            dates = all_dates(dataset, first(variable_names)),
            date = nothing,
            dir = default_download_directory(dataset),
            region = default_region(dataset),
            filenames = nothing,
            start_date = nothing,
            end_date = nothing)

A bundle of Metadata for many variables that share dataset, dates, region, and dir — differing only in variable name.

Each element of an mset::MetadataSet, e.g., mset[name] (or equivalently mset.name or mset[i]) is itself a Metadata — or a Metadatum when dates is a single date. Iteration walks the variable axis, yielding one Metadata per variable.

Arguments

• variable_names: one or more Symbols naming the dataset variables to bundle (e.g. :temperature, :salinity). Verbose dataset-internal names — no aliases. It can also be a tuple of Symbols, e.g., (:temperature, :salinity), but not a vector of Symbols.

Keyword Arguments

• dataset: The shared dataset. Supported datasets are returned by supported_datasets.
• dates: Shared date axis. Either a single AbstractDateTime/AbstractCFDateTime or an AbstractVector of dates. Defaults to all_dates(dataset, first(variable_names)).
• date: Convenience scalar form; cannot be used together with dates.
• region: Shared spatial region — BoundingBox, Column, or nothing.
• dir: Shared download directory.
• filenames: An optional NamedTuple keyed by variable_names overriding the auto-computed per-variable filenames.
• start_date, end_date: Optional date cropping, matching Metadata.

Example

using NumericalEarth, Dates

mset = MetadataSet(:temperature, :salinity;
                   dataset = ECCO4Monthly(),
                   date = DateTime(1995, 1, 1))

mset[2] # Metadata for :salinity

using NumericalEarth, Dates

mset = MetadataSet(:temperature, :salinity;
                   dataset = ECCO4Monthly(),
                   date = DateTime(1995, 1, 1))

mset[2] # Metadata for :salinity

# output
Metadatum{ECCO4Monthly, DateTime}:
├── name: salinity
├── dataset: ECCO4Monthly
├── dates: 1995-01-01 00:00:00
├── filename: SALT_1995_01.nc
└── dir: /.julia/scratchspaces/904d977b-046a-4731-8b86-9235c0d1ef02/ECCO/v4

See also Metadata, Metadatum.

Metadatum(variable_name;
          dataset,
          region = default_region(dataset),
          date = first_date(dataset, variable_name),
          filename = nothing,
          dir = default_download_directory(dataset))

A specialized constructor for a Metadata object with a single date, representative of a snapshot in time.

SurfaceFluxRestoring(dataset_restoring::DatasetRestoring)

A thin wrapper around a DatasetRestoring that converts a 3D restoring tendency into a 2D surface flux boundary condition.

When used as a boundary condition (via getbc), the wrapped DatasetRestoring is evaluated at the top cell (k = Nz) and the resulting tendency G is converted to a surface flux as -G * Δz, consistent with the Oceananigans top-flux sign convention (tendency contribution = - J / Δz).

This is intended for use with the additional_surface_fluxes keyword argument of ocean_simulation, allowing a DatasetRestoring to contribute an additional flux at the surface without overwriting the coupled exchange fluxes.

Example

using ..NumericalEarth

restoring = DatasetRestoring(metadata, grid; rate = 1 / 30days)
ocean = ocean_simulation(grid;
    additional_surface_fluxes = (; S = SurfaceFluxRestoring(restoring)))

all_dates(metadata)

Extract all the dates of the given metadata formatted using the DateTime type.

first_date(dataset, variable_name)

Extracts the first date of the given dataset and variable name formatted using the DateTime type.

last_date(dataset, variable_name)

Extracts the last date of the given dataset and variable name formatted using the DateTime type.

metadata_filename(dataset, name, date, region)

Compute the filename for a single date. Extended by each dataset module.

metadata_filename(metadata)

Return the stored filename(s) of metadata.

native_grid(metadata::Metadata, arch=CPU(); halo = (3, 3, 3))

Return the native grid corresponding to metadata with halo size. Returns a LatitudeLongitudeGrid for global or BoundingBox regions, and a column RectilinearGrid for Column regions.

supported_datasets()

Return the dataset constructors currently supported by Metadata.

Currently, these are: ASTERGEDv3(), AVISODaily(), AVISOMonthly(), GLO30(), GLO90(), ECCO2Daily(), ECCO2DarwinMonthly(), ECCO2Monthly(), ECCO4DarwinMonthly(), ECCO4Monthly(), EN4Monthly(), ERA5HourlyLand(), ERA5HourlySingleLevel(), ERA5MonthlyLand(), ERA5MonthlySingleLevel(), ERA5YearlySingleLevel(), ETOPO2022(), GEBCO2026(), GHSBuiltH(), GHSBuiltS(), GLORYSDaily(), GLORYSMonthly(), GLORYSStatic(), GlobalBuildingFootprints3D(), GloFASReanalysis(), IBCAOv5(), IBCSOv2(), MultiYearJRA55(), RepeatYearJRA55(), ORCAOne(), ORCAQuarter(), ORCATwelfth(), OpenLandMapSoilDB(), SeaWiFSMonthly(), SoilGrids2(), WOAAnnual(), WOAMonthly(), and ESAWorldCover().

validate_dataset_coverage(grid, metadata)

Check that grid lies within the spatial coverage of metadata's dataset. Throws an error if the grid extends outside the dataset's domain. The default implementation does nothing (all grids are accepted). Dataset-specific methods can override this to enforce coverage constraints.

ASTERGEDResolution

The two published ASTER GED v3 resolutions: ASTERGEDLow1km (1 km, 36 arcsec) and ASTERGEDHigh100m (100 m, 3.6 arcsec).

ASTERGEDv3 <: AbstractStaticDataset

Global Emissivity Dataset (GED) v3 from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) aboard NASA's Terra satellite: a static 2000–2008 clear-sky mean of land-surface emissivity on a WGS84 lat/lon grid, distributed as HDF5 in 1°×1° tiles. Two resolutions are supported (see ASTERGEDResolution):

• ASTERGEDLow1km — 1 km (36 arcsec), the default.
• ASTERGEDHigh100m — 100 m (3.6 arcsec), for sub-km domains.

ASTER GED provides five narrowband emissivities (TIR bands 10–14). A longwave scheme needs a single broadband value, so the download collapses the five bands using the Ogawa & Schmugge (2004) coefficients. The resulting Field is finite everywhere and can be passed to SurfaceRadiationProperties as its emissivity.

The product also retrieves an emissivity over water (open-ocean broadband ε ≈ 0.98). Genuine gaps appear only where a clear-sky retrieval was never obtained — persistent tropical cloud, summertime snow screened as cloud — and are inpainted when the Field is built, land from land and water from water, so a coastal gap never inherits the ocean value.

ASTER GED is read in regional windows only, so build the Metadatum with a longitude/latitude BoundingBox, most simply derived from the model grid,

region = BoundingBox(grid; padding = default_horizontal_padding(ASTERGEDv3()))

Reading the HDF5 tiles requires ArchGDAL (with the HDF5 driver) and NASA Earthdata credentials (EARTHDATA_USERNAME / EARTHDATA_PASSWORD).

Reference: Hulley et al. (2015), https://doi.org/10.1002/2015GL065564 Data source: https://www.earthdata.nasa.gov/data/catalog/lpcloud-ag100-003

ASTERGEDv3(resolution = ASTERGEDLow1km)

Construct an ASTERGEDv3 dataset. resolution is ASTERGEDLow1km (1 km, default) or ASTERGEDHigh100m (100 m); see ASTERGEDResolution.

julia> using NumericalEarth

julia> ASTERGEDv3()
ASTERGEDv3(resolution = ASTERGEDLow1km)

julia> ASTERGEDv3(resolution = ASTERGEDHigh100m)
ASTERGEDv3(resolution = ASTERGEDHigh100m)

GLO30

Copernicus DEM GLO-30: global 30 m (1 arc-second) Digital Surface Model (DSM), representing the Earth's surface including buildings, infrastructure, and vegetation. Heights are referenced to the EGM2008 geoid; ocean is set to 0.

Because GLO-30 is a global 30 m product (≈ 1.3M × 0.65M cells), it is read in regional windows only: construct the Metadatum with a longitude/latitude BoundingBox and use it with regrid_topography.

Data is read from the Earth Data Hub Zarr store, which requires a (free) DestinE personal access token in the DESTINE_ACCESS_TOKEN environment variable. Register at https://platform.destine.eu/ and create a token at https://earthdatahub.destine.eu/account-settings#my-personal-access-tokens.

Reading the Zarr store requires the Zarr package to be loaded (using Zarr).

Data source: https://earthdatahub.destine.eu/collections/copernicus-dem/datasets/GLO-30

GLO90

Copernicus DEM GLO-90: global 90 m (3 arc-second) Digital Surface Model, the coarser sibling of GLO30. Same source, access, and usage; lower resolution and a correspondingly smaller native grid.

Data source: https://earthdatahub.destine.eu/collections/copernicus-dem/datasets/GLO-90

CopernicusAlbedo(; diffuse_fraction = 0.2)

The Copernicus Global Land Service (CGLS) 1 km surface-albedo dataset, a ten-day time series derived from SPOT/VGT and PROBA-V observations. Provides :albedo, the broadband blue-sky albedo blended from the black-sky (AL_DH_BB) and white-sky (AL_BH_BB) broadband albedos with diffuse fraction diffuse_fraction.

Files are global on a regular 1/112° latitude-longitude grid; build the Metadata with a lon/lat BoundingBox to window a region at read time. Downloads come from the C3S satellite-albedo catalog entry and require the CDSAPI backend: using CDSAPI with ~/.cdsapirc credentials, the same setup as ERA5 (see this module's README).

julia> using NumericalEarth

julia> CopernicusAlbedo()
CopernicusAlbedo{Float64}(0.2)

CopernicusAlbedoClimatology(; diffuse_fraction = 0.2, years = 2019:2019)

A 12-month climatology of the CopernicusAlbedo blue-sky broadband albedo, built by averaging all ten-day files of years month by month (NaN-aware, per pixel). The monthly-mean files are generated on demand by build_monthly_climatology! (triggered automatically on download) and yield a 12-slot FieldTimeSeries with Cyclical() time indexing.

julia> using NumericalEarth

julia> CopernicusAlbedoClimatology()
CopernicusAlbedoClimatology{Float64, UnitRange{Int64}}(0.2, 2019:2019)

build_monthly_climatology!(dataset::CopernicusAlbedoClimatology;
                           name = :albedo,
                           months = 1:12,
                           dir = default_download_directory(dataset),
                           latitude_chunk = 1120)

For each calendar month in months, download every ten-day albedo file of dataset.years falling in that month, average the black-sky and white-sky broadband albedos pixel by pixel (NaN-aware: pixels missing in a ten-day file are excluded from its mean), and write one monthly-mean NetCDF to dir under the name computed by metadata_filename. Months whose file already exists are skipped. Returns the paths of the 12 (or length(months)) monthly files.

ECCOMetadatum(name;
              date = first_date(ECCO4Monthly(), name),
              dir = download_ECCO_cache)

An alias to construct a Metadatum of ECCO4Monthly().

ECCOPrescribedAtmosphere([architecture = CPU()];
                          dataset = ECCO4Monthly(),
                          start_date = first_date(dataset, :air_temperature),
                          end_date = last_date(dataset, :air_temperature),
                          dir = default_download_directory(dataset),
                          time_indices_in_memory = 10,
                          time_indexing = Cyclical(),
                          surface_layer_height = 2,  # meters
                          other_kw...)

Return a PrescribedAtmosphere representing ECCO state estimate data. The atmospheric data will be held in FieldTimeSeries objects containing

• velocities: u, v
• air temperature and humidity: T, q
• surface pressure: p
• freshwater flux: rain

Note: downwelling shortwave / longwave radiation is now part of the top-level radiation component. Use ECCOPrescribedRadiation to load ECCO SW/LW into a PrescribedRadiation.

ECCOPrescribedRadiation([architecture = CPU()];
                        dataset = ECCO4Monthly(),
                        start_date = first_date(dataset, :downwelling_shortwave),
                        end_date = last_date(dataset, :downwelling_shortwave),
                        dir = default_download_directory(dataset),
                        time_indices_in_memory = 10,
                        time_indexing = Cyclical(),
                        ocean_surface = SurfaceRadiationProperties(0.05, default_water_emissivity),
                        sea_ice_surface = SurfaceRadiationProperties(0.7, 1.0),
                        stefan_boltzmann_constant = default_stefan_boltzmann_constant,
                        other_kw...)

Return a PrescribedRadiation backed by ECCO downwelling shortwave and longwave fields.

retrieve_data(metadata::Metadatum{<:ECCO4DarwinMonthly})

Read a ECCO4DarwinMonthly data file and regrid using MeshArrays on to regular lat-lon grid

EN4Metadatum(name;
             date = first_date(EN4Monthly()),
             dir = download_EN4_cache)

An alias to construct a Metadatum of EN4Monthly.

Hourly ERA5-Land reanalysis (reanalysis-era5-land): a global 0.1° land-surface product covering skin/soil temperature, soil moisture, 2 m temperature/dewpoint, and snow depth/water-equivalent. Ocean cells are masked and never inpainted.

ERA5HourlyPressureLevels(; pressure_levels = ERA5_all_pressure_levels, z = nothing)

Hourly ERA5 pressure-levels dataset metadata. By default (z = nothing) the native grid's vertical coordinate is a 3-D PressureLevelVerticalDiscretization built from the time-evolving geopotential Φ(λ,φ,p,t)/g with sub-surface levels clipped at the surface — the right thing over terrain and synoptic Φ-displacement (issue #236). Each pressure level's height follows the reanalysis in time as the atmosphere's clock advances.

To pin the z-coordinate to something static — e.g. for a quick test without an extra Φ download — pass a precomputed vector: z = standard_atmosphere_z_interfaces(pressure_levels).

Monthly-mean counterpart of ERA5HourlyLand.

ERA5PrescribedAtmosphere([architecture = CPU()];
                         dataset = ERA5HourlySingleLevel(),
                         start_date = first_date(dataset, :temperature),
                         end_date = last_date(dataset, :temperature),
                         dir = download_ERA5_cache,
                         time_indices_in_memory = 24,
                         time_indexing = Cyclical(),
                         surface_layer_height = 10,    # meters
                         boundary_layer_height = 512,  # meters
                         thermodynamics_parameters = nothing,
                         region = nothing,
                         other_kw...)

Return a PrescribedAtmosphere representing ERA5 single-level reanalysis, suitable for regional hindcast forcing. Eastward/northward 10 m winds, 2 m temperature, and surface pressure are loaded directly; specific humidity is derived from the 2 m dewpoint and surface pressure (qᵛ = qᵛ⁺(Tᵈ, pˢ)); total precipitation is converted from hourly-accumulated depth (m) to a mass flux (kg m⁻² s⁻¹) at load time and wrapped in a PrescribedPrecipitationFlux.

region (a BoundingBox) restricts the download and the native grid to a sub-domain; the coupled model interpolates the native-resolution atmosphere onto the exchange grid. Pass thermodynamics_parameters to share a specific thermodynamics with the rest of the model (defaults to AtmosphereThermodynamicsParameters at the data's float type).

ERA5PrescribedAtmosphere(bounding_box::BoundingBox, dates;
                         architecture = CPU(),
                         dataset = ERA5HourlyPressureLevels(),
                         dir = download_ERA5_cache,
                         time_indices_in_memory = nothing,
                         thermodynamics_parameters = nothing,
                         other_kw...)

Return a 3-D PrescribedAtmosphere built from ERA5 pressure-level reanalysis over bounding_box at the requested dates — a range or vector of dates, or a (start_date, end_date) tuple that expands to the dataset's native (hourly or monthly) cadence — on ERA5's native grid: a PressureLevelGrid at the reanalysis' native horizontal resolution with a time-varying geopotential-height vertical (each pressure level's height follows the reanalysis as the atmosphere's clock advances). Each variable loads natively and raw (per pressure level, no vertical remap); a downstream model (e.g. a NestedSimulation child) interpolates the parent onto its own grid on the fly, at the current heights. time_indices_in_memory defaults to all dates.

The atmosphere holds eastward/northward velocities, temperature, specific_humidity, microphysical_variables = (; qᶜˡ, qʳ, qᶜⁱ, qˢ) (cloud liquid/ice + rain/snow water content), and pressure (the level coordinate, via pressure_level_field). Use as the parent of a NestedSimulation.

ERA5PrescribedRadiation([architecture = CPU()]; dataset = ERA5HourlySingleLevel(),
                        time_indices_in_memory = 24, other_kw...)

ERA5 downwelling radiation, suitable for regional hindcast forcing. ERA5 stores these as energy accumulated over the previous hour (J m⁻²); the load-time conversion_units divides by the accumulation interval to recover the mean flux (W m⁻²). Remaining keyword arguments go to prescribed_radiation.

pressure_field(metadata::ERA5PressureMetadatum, arch=CPU(); halo=(3,3,3))

Return a Field{Nothing, Nothing, Center} on the native grid of metadata holding the pressure value (Pa) at each vertical level. Levels are ordered bottom-to-top (k=1 is the highest pressure level). The Nothing horizontal locations make this field broadcast against full 3-D fields without copying.

GEBCO2026

General Bathymetric Chart of the Oceans 2026 release. Global bathymetry and topography at 15 arc-second resolution.

Data source: https://www.gebco.net/dataandproducts/griddedbathymetrydata/

GHSBuiltH

GHS-BUILT-H R2023A average net building height (ANBH, epoch 2018, 100 m, World Mollweide ESRI:54009). Provides :building_height — the mean height in meters of the built-up pixels within each cell. A height of 0 over non-built land is a valid value, not missing, so the field is not inpainted. Only the product no-data is masked to NaN.

Because it is a global 100 m product on a 1000 km Mollweide tile grid, it is read in regional windows only: construct the Metadatum with a longitude/latitude BoundingBox. The windowed tile download and Mollweide→EPSG:4326 warp require using ArchGDAL.

Reference: Pesaresi, M. & Politis, P. (2023), GHS-BUILT-H R2023A, European Commission JRC, doi:10.2905/85005901-3A49-48DD-9D19-6261354F56FE.

julia> using NumericalEarth.DataWrangling.GHSL

julia> GHSBuiltH()
GHSBuiltH()

GHSBuiltS(; resolution = GHSBuiltS100m, epoch = 2020)

GHS-BUILT-S R2023A built-up surface (World Mollweide ESRI:54009). Provides :built_up_fraction — the plan-area index λᵖ ∈ [0, 1] of each cell covered by buildings, obtained from the built-up surface (m² per cell) by dividing by the native cell area and clamping to [0, 1]. A fraction of 0 over non-built land is a valid value, not missing.

resolution selects the native pixel size, GHSBuiltS10m or GHSBuiltS100m (see GHSBuiltSResolution); epoch the reference year. The 10 m product is only published for epoch = 2018; the 100 m product covers 1975–2030 in 5-year steps. The default GHSBuiltS100m/2020 pairs naturally with a ~100 m model grid and the epoch-2018 GHSBuiltH height; pass resolution = GHSBuiltS10m for the finest built fraction (much larger tiles).

Read in regional windows only (see GHSBuiltH); the tile download and Mollweide→EPSG:4326 warp require using ArchGDAL.

Reference: Pesaresi, M. & Politis, P. (2023), GHS-BUILT-S R2023A, European Commission JRC, doi:10.2905/9F06F36F-4B11-47EC-ABB0-4F8B7B1D72EA.

julia> using NumericalEarth.DataWrangling.GHSL

julia> GHSBuiltS()
GHSBuiltS(resolution = GHSBuiltS100m, epoch = 2020)

julia> GHSBuiltS(resolution = GHSBuiltS10m)
GHSBuiltS(resolution = GHSBuiltS10m, epoch = 2018)

GHSBuiltSResolution

The two published GHS-BUILT-S native pixel sizes: GHSBuiltS10m (10 m) and GHSBuiltS100m (100 m).

GlobalBuildingFootprints3D(; resolution = 3)

3D-GloBFP building footprints: ~1.3 billion building polygons (LoD1), each carrying an estimated height, distributed globally as per-tile shapefiles in EPSG:4326.

The heights are rasterized onto a fine grid of resolution meters, giving a single variable :building_height (m, 0 where unbuilt). building_morphometry reduces that raster onto a coarser target grid.

Being a global vector product, it is read in regional windows only: build the Metadatum with a longitude/latitude BoundingBox. The tile download and rasterization require using ArchGDAL.

Heights are machine-learning estimates (RMSE 1.9–14.6 m) and biased low. The product is licensed CC BY 4.0 over footprint geometry from Microsoft Building Footprints, OpenStreetMap (ODbL), and Google–Microsoft Open Buildings.

Reference: Che et al. (2024), https://doi.org/10.5194/essd-16-5357-2024

julia> using NumericalEarth

julia> GlobalBuildingFootprints3D()
GlobalBuildingFootprints3D(resolution = 3 m)

building_morphometry(target_grid; dataset = GlobalBuildingFootprints3D(), region,
                     maximum_raster_cells = 400_000_000)

Per-cell building morphometry on target_grid (a LatitudeLongitudeGrid, coarser than the dataset rasterization resolution), aggregated from the fine 3D-GloBFP building-height raster over region. Returns a NamedTuple of Fields:

• plan_area_index λᵖ — fraction of fine cells that are built.
• mean_building_height h — mean height over the built fine cells.
• building_height_deviation σʰ — standard deviation of height over the built fine cells.
• maximum_building_height hᵐᵃˣ — maximum height.
• gross_building_height — mean height over all fine cells (= λᵖ·h), the digital surface lift.
• frontal_area_index λᶠ — windward wall area from height steps, direction-averaged: (Σₓ|δh|·dy + Σᵧ|δh|·dx) / (4·A), with A the cell area.

A region whose native raster exceeds maximum_raster_cells (default 400_000_000 cells, 3.2 GB of Float64) is reduced in latitude bands sized to the limit, reproducing the single pass exactly while memory stays bounded regardless of the region size; band raster files are deleted once reduced (the downloaded footprint tiles stay cached).

Downloading and rasterizing the footprints requires using ArchGDAL.

GloFASReanalysis(; system_version="version_4_0")

The GloFAS-ERA5 historical river-discharge reanalysis (cems-glofas-historical, product_type="consolidated", hydrological_model="lisflood"). Daily river discharge in m³ s⁻¹ on a global 0.05° grid (v4), 1979–present.

GloFASPrescribedLand(grid;
                     dataset = GloFASReanalysis(),
                     start_date = first_date(dataset, :river_discharge),
                     end_date = last_date(dataset, :river_discharge),
                     dir = download_GloFAS_cache,
                     time_indices_in_memory = 10,
                     time_indexing = Cyclical(),
                     region = nothing,
                     freshwater_density = 1000,
                     maximum_search_radius = 5,
                     spread_radius = 1.2,
                     maximum_spread_cells = 8,
                     other_kw...)

Return a PrescribedLand that forces the ocean with GloFAS river discharge routed onto the coastline of grid (the target ocean grid).

GloFAS provides river discharge (m³ s⁻¹) already accumulated downstream to river mouths by the LISFLOOD channel-routing model, driven by ERA5 runoff. The mouths are located from the dataset's land/ocean boundary and mapped to the active ocean cells of grid surrounding each mouth; the discharge is then deposited as a freshwater mass flux that conserves volume (see build_river_routing).

Keyword Arguments

• freshwater_density: density converting the discharge (m³ s⁻¹) to a mass flux (kg m⁻² s⁻¹). Default: 1000.
• maximum_search_radius: search distance in grid cells for the ocean cell receiving a mouth. Default: 5.
• spread_radius: radius in degrees over which each mouth's discharge is divided equally. Default: 1.2.
• maximum_spread_cells: cap on that footprint, nearest first. Default: 8.

See also JRA55PrescribedLand.

IBCAOv5

International Bathymetric Chart of the Arctic Ocean Version 5.1 (2025). 100m resolution bathymetry for the Arctic Ocean (north of 64°N), including Greenland ice sheet surface elevation. Data is provided in Polar Stereographic projection (EPSG:3996) and reprojected to WGS84 geographic coordinates at 0.01° resolution on first use.

Reference: Jakobsson et al. (2024), https://doi.org/10.1038/s41597-024-04278-w Data source: https://www.gebco.net/data-products/gridded-bathymetry-data/arctic-ocean/

IBCSOv2

International Bathymetric Chart of the Southern Ocean Version 2 (2024 Annual Release). High-resolution (500m) bathymetry for the Southern Ocean (south of 50°S).

Reference: Dorschel et al. (2022), https://doi.org/10.1594/PANGAEA.937574 Data source: https://ibcso.org/ibcso-2024-annual-release/

JRA55PrescribedAtmosphere([architecture = CPU()];
                          dataset = RepeatYearJRA55(),
                          start_date = first_date(dataset, :temperature),
                          end_date = last_date(dataset, :temperature),
                          dir = download_JRA55_cache,
                          time_indices_in_memory = 10,
                          time_indexing = Cyclical(),
                          surface_layer_height = 10,  # meters
                          region = nothing,
                          other_kw...)

Return a PrescribedAtmosphere representing JRA55 reanalysis data. Each atmospheric field is constructed via FieldTimeSeries(::JRA55Metadata), which uses a DatasetBackend parameterized by JRA55 metadata so that the JRA55-specific set! (chunked-yearly NetCDF) is dispatched. The region keyword restricts the atmosphere to a sub-domain of the global JRA55 grid.

Note: downwelling shortwave / longwave radiation is now part of the top-level radiation component. Use JRA55PrescribedRadiation to load JRA55 SW/LW into a PrescribedRadiation.

JRA55PrescribedLand(grid;
                    dataset = RepeatYearJRA55(),
                    start_date = first_date(dataset, :river_freshwater_flux),
                    end_date = last_date(dataset, :river_freshwater_flux),
                    dir = download_JRA55_cache,
                    time_indices_in_memory = 10,
                    time_indexing = Cyclical(),
                    region = nothing,
                    maximum_search_radius = 5,
                    spread_radius = 1.2,
                    maximum_spread_cells = nothing,
                    other_kw...)

Return a PrescribedLand holding the JRA55-do river runoff and iceberg calving fluxes, routed onto the coastline of the target ocean grid by build_flux_routing.

Keyword Arguments

• maximum_search_radius: search distance in grid cells for the ocean cell receiving a mouth. Default: 5.
• spread_radius: radius in degrees over which each mouth's discharge is divided equally. Default: 1.2.
• maximum_spread_cells: cap on that footprint, nearest first. Default: nothing (uncapped).

See also GloFASPrescribedLand.

JRA55PrescribedRadiation([architecture = CPU()]; dataset = RepeatYearJRA55(),
                         time_indices_in_memory = 10, other_kw...)

JRA55 downwelling radiation. Remaining keyword arguments go to prescribed_radiation.

OpenLandMapSoilDB <: AbstractStaticDataset

OpenLandMap-soilDB global soil properties at 30 m (Hengl et al., 2026), predicted from spatiotemporal machine learning over Landsat/MODIS/Sentinel covariates.

Delivers static soil texture mass fractions (:sand_fraction, :silt_fraction, :clay_fraction, in kg/kg) and fine-earth bulk density (:bulk_density, in kg/m³) over the three native depth intervals 0–30, 30–60, and 60–100 cm — stored as a three-dimensional field whose vertical axis carries the depths (deepest first). Data are plain geographic EPSG:4326 so no reprojection is needed.

Because the global grid is ~1.44M × 528k cells, this dataset is read in regional windows only: construct the Metadatum with a longitude/latitude BoundingBox. Reading the cloud-optimized GeoTIFFs requires ArchGDAL to be loaded (using ArchGDAL); access is anonymous, no credentials needed.

Coverage spans latitudes −56° to 76° (Antarctica excluded); permanent ice and sand deserts are masked to NaN (pass inpainting to Field to fill them).

Data source: https://stac.openlandmap.org/ (CC-BY 4.0).

using NumericalEarth

metadatum = Metadatum(:clay_fraction;
                      dataset = OpenLandMapSoilDB(),
                      region = BoundingBox(longitude = (-112.3, -111.9),
                                           latitude  = (36.0, 36.4)))

metadatum.filename

# output
"OpenLandMap_clay_fraction_lon_-112.3_-111.9_lat_36.0_36.4.nc"

SeaWiFSMonthly(; resolution = 1)

Monthly composites of near-surface chlorophyll-a from SeaWiFS, reprocessing R2018.0, covering the mission record September 1997 to December 2010. Provides :chlorophyll in mg m⁻³.

The native product is mapped at 1/12°; resolution is the degree spacing actually requested. Files come from the NOAA CoastWatch ERDDAP server, which needs no credentials.

A FieldTimeSeries built over a twelve-month window of Metadata(:chlorophyll; dataset=SeaWiFSMonthly(), dates) drives a run of any length with a repeating seasonal cycle, under the default Cyclical() time indexing.

ESAWorldCover(; version = WorldCoverV200, aggregation_factor = 12)

ESA WorldCover global 10 m land-cover classification.

version selects the release, either WorldCoverV200 (2021, the default) or WorldCoverV100 (2020); see WorldCoverVersion.

aggregation_factor is the integer number of 10 m pixels reduced per side into one aggregated cell (12 → ~110 m, the default, so each aggregated cell still samples ~144 sub-pixels at ~100 m model resolution; 120 → ~1 km, cheaper for a large region). Class codes are never averaged: the fine raster is counted block-wise into a majority class field, a vegetation-fraction field, and one per-class area-fraction field. The integer factor keeps the aggregated cell boundaries on fine-pixel boundaries, so the categorical field is never reprojected.

The Map band is a UInt8 class code and 0 is no-data. The legend (ESA_WORLDCOVER_CLASS_NAMES) is

Because the source is a 10 m categorical raster, it is read in regional windows only: build the Metadatum with a longitude/latitude BoundingBox. Available variables are :landcover_class (the majority class), a per-class :<class>_fraction for each name above (e.g. :cropland_fraction), and :vegetation_fraction, the mosaic weight f_veg summing tree cover, shrubland, grassland, cropland, herbaceous wetland, and mangroves — see ESA_WORLDCOVER_VEGETATED_CLASSES for why that set is a modeling choice.

Reading the anonymous Cloud-Optimized GeoTIFF tiles from the public s3://esa-worldcover/ bucket requires the ArchGDAL package to be loaded (using ArchGDAL).

using NumericalEarth

ESAWorldCover()

# output
ESAWorldCover(version = WorldCoverV200, aggregation_factor = 12)

Data source: https://esa-worldcover.org/en/data-access ; DOI v200 10.5281/zenodo.7254221 (Zanaga et al. 2022), license CC-BY 4.0.

WorldCoverVersion

The published ESA WorldCover releases: WorldCoverV100 (2020) and WorldCoverV200 (2021).

The two use different classification algorithms and ESA warns they are not comparable for change detection.

MixedLayerDepthField(bm, grid, tracers;
                     difference_criterion = 2.87e-4,
                     reference_depth = 10.0)

Mixed-layer-depth diagnostic using the de Boyer Montégut et al. (2004) DR003 buoyancy-difference criterion. The default difference_criterion = 2.87e-4 m s⁻² corresponds to Δσ = 0.03 kg m⁻³ at ρ₀ = 1025 kg m⁻³. The default reference_depth = 10 m is the standard dBM reference; buoyancy at that depth is obtained by linear interpolation between adjacent cell centers, so the diagnostic is insensitive to the vertical resolution near the surface.

References

• de Boyer Montégut, C., G.Madec, A. S.Fischer, A.Lazar, and D.Iudicone (2004), Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology, J. Geophys. Res., 109, C12003, doi:10.1029/2004JC002378.

atmosphere_ocean_freshwater_flux(esm::EarthSystemModel)

Return the atmosphere-ocean freshwater mass flux (kg m⁻² s⁻¹) at the atmosphere-ocean interface in a coupled esm.

atmosphere_ocean_heat_flux(esm::EarthSystemModel)

Return the atmosphere-ocean heat flux (W m⁻²) at the atmosphere-ocean interface in a coupled esm.

frazil_heat_flux(esm::EarthSystemModel)

Return the two-dimensional frazil heat flux (W m⁻²) in a coupled esm.

meridional_heat_transport(esm::EarthSystemModel;
                          reference_temperature = 0)

Return the meridional heat transport for the coupled esm::EarthSystemModel by computing the meridional heat flux.

The meridional heat transport is computed via:

\[\mathrm{MHT} ≡ ρᵒᶜ cᵒᶜ ∫ v (T - T_{\rm ref}) \, \mathrm{d}x \, \mathrm{d}z\]

Above, $T_{\rm ref}$ is a reference temperature and $ρᵒᶜ$ and $cᵒᶜ$ are the ocean reference density and specific heat capacity respectively.

Arguments

• esm: An EarthSystemModel.

Keyword Arguments

• reference_temperature: The reference temperature (in ᵒC) used for the calculation; default: 0 ᵒC.

  Reference temperature

  The reference temperature is only relevant when we compute the meridional heat transport over a section where there is a net volume transport. If we are computing the diagnostic globally, i.e., around a whole latitude circle, then by necessity there is no net volume transport and thus the reference temperature value is irrelevant. Section-averaged transport could also be considered as a reference temperature to remove residual barotropic volume fluxes in basin-scale/regional analyses where a net volume transport is present.

Example

using ..NumericalEarth
using Oceananigans

grid = RectilinearGrid(size = (4, 5, 2), extent = (1, 1, 1),
                       topology = (Periodic, Bounded, Bounded))

ocean = ocean_simulation(grid;
                         momentum_advection = nothing,
                         tracer_advection = nothing,
                         closure = nothing,
                         coriolis = nothing)

sea_ice = sea_ice_simulation(grid, ocean)

atmosphere = PrescribedAtmosphere(grid, [0.0])
radiation = PrescribedRadiation(grid)

esm = OceanSeaIceModel(ocean, sea_ice; atmosphere, radiation)

mht = meridional_heat_transport(esm)

# output

Integral of BinaryOperation at (Center, Face, Center) over dims (1, 3)
└── operand: BinaryOperation at (Center, Face, Center)
    └── grid: 4×5×2 RectilinearGrid{Float64, Periodic, Bounded, Bounded} on CPU with 3×3×2 halo

net_ocean_freshwater_flux(esm::EarthSystemModel)

Return the net freshwater mass flux (kg m⁻² s⁻¹) at the ocean's surface in a coupled esm.

net_ocean_heat_flux(esm::EarthSystemModel)

Return the net heat flux (W m⁻²) at the ocean's surface in a coupled esm.

sea_ice_ocean_freshwater_flux(esm::EarthSystemModel)

Return the sea ice-ocean freshwater mass flux (kg m⁻² s⁻¹) at the sea ice-ocean interface in a coupled esm.

sea_ice_ocean_heat_flux(esm::EarthSystemModel)

Return the sea ice-ocean heat flux (W m⁻²) at the sea ice-ocean interface in a coupled esm.

EarthSystemModel(radiation, atmosphere, land, sea_ice, ocean;
                 clock = default_earth_system_clock(atmosphere),
                 ocean_reference_density = reference_density(ocean),
                 ocean_heat_capacity = heat_capacity(ocean),
                 sea_ice_reference_density = reference_density(sea_ice),
                 sea_ice_heat_capacity = heat_capacity(sea_ice),
                 interfaces = nothing)

Construct a coupled earth system model. Components are passed in struct order (top to bottom): radiation, atmosphere, land, sea_ice, ocean. Pass nothing for components that are absent. For simpler configurations, see OceanOnlyModel, OceanSeaIceModel, and AtmosphereOceanModel.

Arguments

• radiation: Radiation component, or nothing for a radiatively decoupled surface. Pass a PrescribedRadiation (e.g. JRA55PrescribedRadiation(...)) to enable radiative forcing.
• atmosphere: A representation of a possibly time-dependent atmospheric state, or nothing. For a prognostic atmosphere, use atmosphere_simulation. For prescribed atmospheric forcing, use JRA55PrescribedAtmosphere or PrescribedAtmosphere.
• land: Land component, or nothing.
• sea_ice: A representation of a possibly time-dependent sea ice state, or nothing. For example, the minimalist FreezingLimitedOceanTemperature represents oceanic latent heating during freezing only, but does not evolve sea ice variables. For prognostic sea ice use an Oceananigans.Simulation of ClimaSeaIce.SeaIceModel.
• ocean: A representation of a possibly time-dependent ocean state. Currently, only Oceananigans.Simulations of Oceananigans.HydrostaticFreeSurfaceModel are tested.

Keyword Arguments

• clock: Keeps track of time. Defaults to a clock whose time type matches the atmosphere's (Float64 without an atmosphere).
• ocean_reference_density: Reference density for the ocean. Defaults to value from ocean model.
• ocean_heat_capacity: Heat capacity for the ocean. Defaults to value from ocean model.
• sea_ice_reference_density: Reference density for sea ice. Defaults to value from sea ice model.
• sea_ice_heat_capacity: Heat capacity for sea ice. Defaults to value from sea ice model.
• interfaces: Component interfaces for coupling. Defaults to nothing and will be constructed automatically. To customize the sea ice-ocean heat flux formulation, create interfaces manually using ComponentInterfaces.

EarthSystemModel(; radiation = nothing,
                   atmosphere = nothing,
                   land = nothing,
                   sea_ice = nothing,
                   ocean = nothing,
                   kw...)

Keyword-only constructor for EarthSystemModel. Equivalent to the positional form, but lets you pass only the components you actually have:

EarthSystemModel(; atmosphere, ocean)                       # ocean + atmosphere
EarthSystemModel(; atmosphere, sea_ice, ocean, radiation)   # full coupled

All keyword arguments accepted by the positional constructor are forwarded.

AtmosphereLandModel(atmosphere, land; radiation=nothing, kw...)

Construct a coupled atmosphere–land model. Convenience constructor for EarthSystemModel with an atmosphere and land but no ocean or sea ice. All keyword arguments are forwarded to EarthSystemModel.

The atmosphere–land turbulent fluxes are computed via SimilarityTheoryFluxes using land-side roughness and a β-reduced surface specific humidity (qᵃᵗ + β·[qₛ - qᵃᵗ]).

AtmosphereOceanModel(atmosphere, ocean; kw...)

Construct a coupled atmosphere–ocean model. Convenience constructor for EarthSystemModel with an atmosphere and ocean but no sea ice. All keyword arguments are forwarded to EarthSystemModel.

OceanOnlyModel(ocean; atmosphere=nothing, radiation=nothing, land=nothing, kw...)

Construct an ocean-only model without a sea ice component. This is a convenience constructor for EarthSystemModel that sets sea_ice to FreezingLimitedOceanTemperature (a simple freezing limiter that does not evolve sea ice variables).

The atmosphere, radiation, and land keywords can be used to specify prescribed components (e.g., JRA55PrescribedAtmosphere). All other keyword arguments are forwarded to EarthSystemModel.

using NumericalEarth, Oceananigans

grid = LatitudeLongitudeGrid(GPU();
                             size = (20, 20, 4),
                             z = (-100, 0),
                             latitude = (-80, 80),
                             longitude = (0, 360),
                             halo = (6, 6, 3))

ocean = ocean_simulation(grid, closure=nothing)
set!(ocean.model, T=20, S=35, u=0.01, v=-0.005)

ocean = OceanOnlyModel(ocean)
# output

EarthSystemModel{GPU}(time = 0 seconds, iteration = 0)
├── radiation: Nothing
├── atmosphere: Nothing
├── land: Nothing
├── sea_ice: FreezingLimitedOceanTemperature
├── ocean: HydrostaticFreeSurfaceModel{CUDAGPU, LatitudeLongitudeGrid}(time = 0 seconds, iteration = 0)
└── interfaces: ComponentInterfaces

OceanSeaIceModel(ocean, sea_ice; atmosphere=nothing, radiation=nothing, land=nothing, kw...)

Construct a coupled ocean–sea ice model.

This is a convenience constructor for EarthSystemModel with explicit ocean and sea ice components and optional prescribed atmosphere, prescribed radiation, and prescribed land.

Positional arguments are ocean then sea_ice.

Example

using NumericalEarth, Oceananigans

grid = LatitudeLongitudeGrid(size = (20, 20, 4),
                             z = (-100, 0),
                             latitude = (-80, 80),
                             longitude = (0, 360),
                             halo = (6, 6, 3))

ocean = ocean_simulation(grid, closure=nothing)
set!(ocean.model, T=20, S=35, u=0.01, v=-0.005)

sea_ice = sea_ice_simulation(grid, ocean)

hi(λ, φ) = φ > 70 || φ < -70
set!(sea_ice.model, h=hi, ℵ=hi)

coupled_model = OceanSeaIceModel(ocean, sea_ice)

# output

EarthSystemModel{CPU}(time = 0 seconds, iteration = 0)
├── radiation: Nothing
├── atmosphere: Nothing
├── land: Nothing
├── sea_ice: SeaIceModel{CPU, LatitudeLongitudeGrid}(time = 0 seconds, iteration = 0)
├── ocean: HydrostaticFreeSurfaceModel{CPU, LatitudeLongitudeGrid}(time = 0 seconds, iteration = 0)
└── interfaces: ComponentInterfaces

AltitudeCorrection(surface_elevation, atmosphere_elevation; lapse_rate = 6.5e-3)

A moist-environmental lapse-rate correction of the near-surface atmosphere state for the mismatch between the desired surface elevation (surface_elevation, the first argument) and the elevation the atmosphere data assumes (atmosphere_elevation). The elevation difference Δz = surface_elevation − atmosphere_elevation is materialized on the exchange grid under the hood; the user just supplies the two elevations (anything set! accepts — a Field, function, number, or exchange-grid array).

The gravitational acceleration and dry-air gas constant used in the hydrostatic pressure adjustment are not supplied here — they are pulled from the atmosphere's thermodynamics via thermodynamic_constants when the atmosphere ComponentExchanger materializes the correction, so they aren't duplicated.

Applied in place to the regridded exchange state each step,

\[T ← T − Γ\,Δz, \qquad p ← p\,\exp\!\left(\frac{-g\,Δz}{Rᵈ\,\bar{T}}\right), \qquad \bar{T} = T − Γ\,Δz/2,\]

with specific humidity q conserved (adiabatic lifting conserves q). Δz is naturally ≈ 0 over ocean / sea-ice (sea level), so a single correction is correct across surfaces, and it applies to any atmosphere (prescribed or online).

struct BulkTemperature

A type to represent the interface temperature used in fixed-point iteration for interface fluxes following similarity theory. The interface temperature is not calculated but instead provided by either the ocean or the sea ice model.

CoefficientBasedFluxes(FT = Oceananigans.defaults.FloatType;
                       transfer_coefficients = (1e-3, 1e-3, 1e-3),
                       solver_stop_criteria = nothing,
                       solver_tolerance = 1e-8,
                       solver_maxiter = 20)

Return the structure for computing turbulent fluxes using bulk transfer coefficients. Used in bulk flux calculations to determine the exchange of momentum, heat, and moisture between the ocean/ice surface and the atmosphere.

Arguments

• FT: (optional) Float type for the coefficients, defaults to Oceananigans.defaults.FloatType

Keyword Arguments

• transfer_coefficients: Transfer coefficients for momentum, heat, and moisture. Can be a SimilarityScales, a Tuple, or a NamedTuple with constant or callable entries, or an LargeYeagerTransferCoefficients. Defaults to (1e-3, 1e-3, 1e-3).
• solver_stop_criteria: Criteria for iterative solver convergence. If nothing, creates new criteria using solver_tolerance and solver_maxiter.
• solver_tolerance: Tolerance for solver convergence when creating new stop criteria, defaults to 1e-8.
• solver_maxiter: Maximum iterations for solver when creating new stop criteria, defaults to 20

Example

using Oceananigans
using NumericalEarth

grid = RectilinearGrid(size=3, z=(-1, 0), topology=(Flat, Flat, Bounded))
ocean = ocean_simulation(grid; timestepper = :QuasiAdamsBashforth2)

ao_fluxes = CoefficientBasedFluxes(transfer_coefficients = (1e-2, 1e-3, 1e-3))

interfaces = ComponentInterfaces(nothing, ocean; atmosphere_ocean_fluxes=ao_fluxes)

# output
ComponentInterfaces

struct CoefficientBasedFluxes{C, S}

A structure for computing turbulent fluxes using bulk transfer coefficients.

• transfer_coefficients::Any

• solver_stop_criteria::Any

ComponentInterfaces(atmosphere, ocean, sea_ice=nothing; kwargs...)

Construct component interfaces for atmosphere-ocean-sea ice coupling.

Keyword Arguments

• sea_ice_ocean_heat_flux: formulation for sea ice-ocean heat flux. Options are:

  • IceBathHeatFlux(): bulk heat flux with interface at freezing point
  • ThreeEquationHeatFlux(): coupled heat/salt/freezing point system (default)

• radiation: radiation component. Default: nothing.

• freshwater_density: reference density of freshwater. Default: default_freshwater_density.

• latent_heat_of_fusion: latent heat [J kg⁻¹] the ocean supplies to melt snowfall and icebergs. Default: default_latent_heat_of_fusion.

• atmosphere_ocean_fluxes: flux formulation for atmosphere-ocean interface. Default: SimilarityTheoryFluxes().

• atmosphere_sea_ice_fluxes: flux formulation for atmosphere-sea ice interface. Default: SimilarityTheoryFluxes().

• atmosphere_ocean_interface_temperature: temperature formulation for atmosphere-ocean interface. Options are BulkTemperature() (default) and SkinTemperature(internal_flux), where for the ocean internal_flux is a DiffusiveFlux(κ, δ) with either a prescribed diffusivity κ or an InteriorDiffusivity() assessed from the ocean turbulence closure.

• atmosphere_ocean_interface_specific_humidity: specific humidity formulation. Default: default_ao_specific_humidity(ocean).

• atmosphere_sea_ice_interface_temperature: temperature formulation for atmosphere-sea ice interface. Default: default_ai_temperature(sea_ice).

• ocean_reference_density: reference density for the ocean. Default: reference_density(ocean).

• ocean_heat_capacity: heat capacity for the ocean. Default: heat_capacity(ocean).

• ocean_temperature_units: temperature units for the ocean. Default: DegreesCelsius().

• sea_ice_temperature_units: temperature units for sea ice. Default: DegreesCelsius().

• sea_ice_reference_density: reference density for sea ice. Default: reference_density(sea_ice).

• sea_ice_heat_capacity: heat capacity for sea ice. Default: heat_capacity(sea_ice).

• gravitational_acceleration: gravitational acceleration. Default: default_gravitational_acceleration.

ConstantTortuosity()

Trivial tortuosity model: Dᵛ_eff = Dᵛ₀. Used by DryLayerVaporPistonVelocity when the soil air space is not modeled explicitly.

ConvectiveGustiness{FT}(; gustiness_parameter = 1.2, minimum_gustiness = 0.01)

Beljaars (1995)-type convective gustiness: in unstable conditions $(Jᵇ > 0)$ the bulk velocity is enhanced by $Uᴳ = β (Jᵇ h_{bℓ})^{1/3}$, where $Jᵇ = -u★b★$ is the surface buoyancy flux, $h_{bℓ}$ the boundary layer height, and $β$ the gustiness_parameter. In stable conditions the enhancement falls back to the minimum_gustiness floor [m/s].

struct CriticalSaturation

Evaporation efficiency after Manabe (1969): the surface is saturated (β = 1) above a critical saturation 𝒮ᶜ, and the efficiency falls off linearly below it,

\[β(𝒮) = \min(𝒮 / 𝒮ᶜ, 1), 𝒮 = Mˡᵃ / Mˡᵃ⁺.\]

Used as the efficiency of FractionalHumidity. The type declares its land-state dependency (the saturation 𝒮); the interface materializes exactly that into the land interface state.

DiffusiveFlux(κ, δ)

Internal flux $J = - κ (Tₛ - Tᵢ) / δ$ between the interior temperature $Tᵢ$, located a distance $δ$ below the interface (typically half the spacing of the topmost interior cell), and the interface temperature $Tₛ$. The diffusivity κ (m² s⁻¹) is either a prescribed constant or an InteriorDiffusivity assessed from the interior model.

DryLayerHumidity(FT = Oceananigans.defaults.FloatType,
                 phase = AtmosphericThermodynamics.Liquid();
                 dry_layer_depth,
                 vapor_exchange,
                 thermal_exchange_depth,
                 porosity)

Surface specific-humidity formulation for the dry-layer model: qⁱⁿ is solved from a vapor-flux balance between a Fick flux through an unresolved dry surface layer and the atmospheric vapor flux, following Ye and Pielke (1993) with the dry-layer (DSL) resistance of Yamanaka et al. (1997) and Swenson and Lawrence (2014). The formulation plugs into the compute_interface_state solver exactly where SkinHumidity does, and reduces to a saturated-skin boundary condition when the slab is wet enough (𝒮 ≥ 𝒮ᶜ).

• dry_layer_depth — depth diagnostic, e.g. StorageBasedDryLayerDepth.
• vapor_exchange — δᵛ_min, Dᵛ₀, tortuosity (a DryLayerVaporPistonVelocity).
• thermal_exchange_depth — ℓᵀ (m), the same depth used by SkinTemperature(DiffusiveFlux) on the temperature side. Controls the interpolation Tᵉ = Tⁱⁿ + χ(Tˡᵃ − Tⁱⁿ) with χ = clip(δᵛ/ℓᵀ, 0, 1).
• porosity — ν, soil porosity (matches the hydrology closure; needed for the Millington–Quirk tortuosity).

The dry-layer source humidity is the saturation value at the front temperature, qᵉ = qᵛ⁺(Tᵉ), with no matric-suction (Kelvin-equation) reduction of the pore relative humidity (the hₛ of Ye and Pielke (1993), after Philip 1957) — that factor departs appreciably from 1 only at extreme dryness.

DryLayerVaporPistonVelocity(FT = Oceananigans.defaults.FloatType;
                            minimum_dry_layer_depth,
                            molecular_diffusivity,
                            wet_transition_width = 5 * minimum_dry_layer_depth,
                            tortuosity = ConstantTortuosity())

Parameters of the dry-layer vapor piston velocity wᵈ = Dᵛ_eff / max(δᵛ, δᵛ_min), the reciprocal of the dry-surface-layer soil resistance r_soil = δᵛ/Dᵛ_eff of Yamanaka et al. (1997) and Swenson and Lawrence (2014). The tortuosity model is a singleton type — ConstantTortuosity or PowerLawTortuosity, after Millington and Quirk (1961) — dispatched on by effective_vapor_diffusivity. The piston velocity feeds the DryLayerHumidity flux balance.

wet_transition_width (m) is the width over which the saturated-skin (wet) humidity transitions to the dry-layer series solution: the two are combined with a logistic weight in δᵛ centered at δᵛ_min + wet_transition_width/2, so the transition is infinitely differentiable (see Kavetski and Kuczera (2007)) and ≈99 % complete across δᵛ ∈ [δᵛ_min, δᵛ_min + wet_transition_width]. Pass 0 to recover a sharp switch at δᵛ = δᵛ_min.

struct FractionalHumidity

Surface specific humidity as a fraction of saturation at the surface temperature,

\[qˢ = β · qᵛ⁺(Tₛ),\]

where the evaporation efficiency β is set by efficiency — a CriticalSaturation (Manabe, 1969) or a constant Number. Unlike SkinHumidity, the saturation is taken at the skin temperature: β is a surface evaporation efficiency, not a deep reservoir. BulkHumidity is the 𝒮ᶜ → 0 corner (saturated wherever 𝒮 > 0).

IceBathHeatFlux(FT::DataType = Oceananigans.defaults.FloatType;
                heat_transfer_coefficient = 0.006,
                friction_velocity = 0.02)

Construct an IceBathHeatFlux with the specified parameters.

Keyword Arguments

• heat_transfer_coefficient: turbulent heat exchange coefficient. Default: 0.006.
• friction_velocity: friction velocity value or formulation. Default: 0.02.

IceBathHeatFlux{FT, U}

Bulk formulation for sea ice-ocean heat flux.

The interface temperature is fixed at the freezing point of the surface salinity, and the heat flux is computed using bulk transfer:

\[Q = \rho_o c_o \alpha_h u_* (T - T_m)\]

where $\alpha_h$ is the heat transfer coefficient and $u_*$ is the friction velocity.

Fields

• heat_transfer_coefficient::FT: turbulent heat exchange coefficient $\alpha_h$ (dimensionless)
• friction_velocity::U: friction velocity value or formulation (constant Number or MomentumBasedFrictionVelocity)

Example

using NumericalEarth.EarthSystemModels: IceBathHeatFlux

flux = IceBathHeatFlux(heat_transfer_coefficient = 0.006, friction_velocity = 0.002)

# output
IceBathHeatFlux{Float64}
├── heat_transfer_coefficient: 0.006
└── friction_velocity: 0.002

References

• Holland and Jenkins (1999): Holland, D. M., & Jenkins, A. (1999). Modeling thermodynamic ice–ocean interactions at the base of an ice shelf. Journal of Physical Oceanography, 29(8), 1787-1800.

InteriorDiffusivity(FT = Oceananigans.defaults.FloatType; minimum_diffusivity = 1.4e-7)

Diffusivity for a DiffusiveFlux that is assessed from the interior model (for example the near-surface vertical diffusivity predicted by the ocean turbulence closure) instead of being prescribed. The assessed value is floored by minimum_diffusivity, which defaults to the molecular thermal diffusivity of seawater, guarding stably-stratified conditions in which modeled diffusivities vanish.

LargeYeagerTransferCoefficients{FT, D, SF}

Wind-dependent bulk transfer coefficients following the NCAR/Large & Yeager (2004) algorithm as used in OMIP-2 simulations.

Computes all three transfer coefficients (drag, heat, moisture) from a neutral drag polynomial with stability corrections via Monin-Obukhov similarity functions. The empirical heat and moisture coefficients correspond to the Stanton and Dalton numbers of L&Y eq. 6c-6d.

Pass as transfer_coefficients to CoefficientBasedFluxes:

ly = LargeYeagerTransferCoefficients()
fluxes = CoefficientBasedFluxes(transfer_coefficients = ly,
                                solver_stop_criteria = FixedIterations(5))

References:

• Large, W.G. & Yeager, S.G. (2004): NCAR/TN-460+STR

LinearStableStabilityFunction{FT}

A simple linear stability function for stable conditions: $ψ = -c ζ$, bounded at $|ζ| ≤ ζ_{max}$.

Used by the NCAR/Large-Yeager (2004) bulk formulae with $c = 5$ and $ζ_{max} = 10$.

References:

• Large, W.G. & Yeager, S.G. (2004): NCAR/TN-460+STR

MomentumBasedFrictionVelocity

A friction velocity formulation that computes the friction velocity from momentum stresses.

The friction velocity is computed as:

\[u_* = \sqrt{\frac{|\boldsymbol{\tau}|}{\rho_o}}\]

where τ is the magnitude of the momentum stress vector and ρᵒᶜ is the ocean reference density.

Example

using NumericalEarth.EarthSystemModels: MomentumBasedFrictionVelocity

fv = MomentumBasedFrictionVelocity()

# output
MomentumBasedFrictionVelocity (computed from momentum stresses)

MomentumRoughnessLength(FT = Float64;
                        gravitational_acceleration = default_gravitational_acceleration,
                        maximum_roughness_length = 1.0,
                        air_kinematic_viscosity = 1.5e-5,
                        wave_formulation = 0.011,
                        smooth_wall_parameter = 0.11)

Construct a MomentumRoughnessLength object that represents the momentum roughness length that regulates the exchange of momentum, heat, and water vapor between the ocean and the atmosphere.

Keyword Arguments

• gravitational_acceleration: The gravitational acceleration. Default: default_gravitational_acceleration.
• maximum_roughness_length: The maximum roughness length. Default: 1e-1.
• air_kinematic_viscosity: The air kinematic viscosity. Default: 1.5e-5.
• wave_formulation: The gravity wave parameter. Default: 0.011.
• smooth_wall_parameter: The smoothwallparameter parameter. Default: 0.11.

PolynomialNeutralDragCoefficient{FT}

Neutral 10-m drag coefficient from the Large & Yeager (2004) polynomial:

C_dN = (a/U + b + c*U - d*U^6) × 10⁻³   for U < high_wind_speed_threshold
C_dN = high_wind_drag_coefficient         otherwise

Called as a functor: C_dN = drag(U).

References:

• Large, W.G. & Yeager, S.G. (2004): NCAR/TN-460+STR

PowerLawTortuosity()

Millington–Quirk tortuosity: Dᵛ_eff = Dᵛ₀ · θᵍ^(10/3) / ν² where θᵍ = ν − θˡ is the gas-filled pore fraction. Reduces vapor diffusivity in near-saturated soils.

ScalarRoughnessLength(FT = Float64;
                      air_kinematic_viscosity = temperature_dependent_viscosity,
                      reynolds_number_scaling_function = empirical_scaling_function,
                      maximum_roughness_length = 1.6e-4)

Construct a ScalarRoughnessLength object that represents the scalar roughness length that regulates the exchange of heat and water vapor between the ocean and the atmosphere.

Keyword Arguments

• air_kinematic_viscosity::Function: The function to compute the air kinematic viscosity.
• reynolds_number_scaling_function::Function: The function to compute the Reynolds number scaling factor.
• maximum_roughness_length::Float: The maximum roughness length value. Defaults to 1.6e-4.

SimilarityTheoryFluxes(FT::DataType = Float64;
                       gravitational_acceleration = 9.81,
                       von_karman_constant = 0.4,
                       turbulent_prandtl_number = 1,
                       subgrid_velocities = ConvectiveGustiness{FT}(),
                       stability_functions = default_stability_functions(FT),
                       roughness_lengths = default_roughness_lengths(FT),
                       zero_plane_displacement = 0,
                       similarity_form = LogarithmicSimilarityProfile(),
                       solver_stop_criteria = nothing,
                       solver_tolerance = 1e-8,
                       solver_maxiter = 100)

SimilarityTheoryFluxes contains parameters and settings to calculate air-interface turbulent fluxes using Monin–Obukhov similarity theory.

Keyword Arguments

• von_karman_constant: The von Karman constant. Default: 0.4.
• turbulent_prandtl_number: The turbulent Prandtl number. Default: 1.
• subgrid_velocities: Empirical subgrid velocity enhancement of the bulk velocity: nothing, a Number [m/s], a formulation implementing vsgs², or a SubgridVelocityCorrection composing several. Default: ConvectiveGustiness{FT}().
• stability_functions: The stability functions. Default: default_stability_functions(FT) that follow the formulation of Edson et al. (2013).
• roughness_lengths: The roughness lengths used to calculate the characteristic scales for momentum, temperature and water vapor. Each may be a formulation, a Number, or — at the atmosphere–land interface only — a Field{Center, Center, Nothing} of per-cell values. Default: default_roughness_lengths(FT), formulation taken from Edson et al. (2013).
• zero_plane_displacement: The zero-plane displacement d [m] of surfaces with tall roughness elements (buildings, plant canopy): the similarity profiles are evaluated at the height Δh - d above the interface. A Number, or — at the atmosphere–land interface only — a Field{Center, Center, Nothing} of per-cell displacements. Default: 0 (undisplaced).
• similarity_form: The type of similarity profile used to relate the atmospheric state to the interface fluxes / characteristic scales.
• solver_tolerance: The tolerance for convergence. Default: 1e-8.
• solver_maxiter: The maximum number of iterations. Default: 100.

struct SkinHumidity

Surface specific humidity qˢ solved from a vapor-flux balance at the land surface, the humidity analog of SkinTemperature.

Vapor reaches the surface by diffusing up from saturated soil at the saturation depth d (the surface_thickness), where the soil air is saturated at qᵛ⁺(Tᵢ) — evaluated at the interior (bulk land) temperature, since the reservoir sits at depth below the surface. Fick's law across d gives the internal (soil) vapor flux

\[J^q = - κ^q/d \, (qˢ - qᵛ⁺)\]

with soil vapor diffusivity κ^q (vapor_diffusivity). The surface is massless, so qˢ is the value for which this soil flux balances the atmospheric vapor flux carried away by turbulence — qˢ is solved inside the interface fixed-point iteration (see compute_interface_humidity), exactly as SkinTemperature solves Tₛ from a surface energy balance.

surface_thickness is a Number (fixed d). A future WetnessDependentSurfaceThickness will let d grow as the soil dries, making evaporation self-limiting.

struct SkinTemperature

A type to represent the interface temperature used in the flux calculation. The interface temperature is calculated from the flux balance at the interface. In particular, the interface temperature $Tₛ$ is the root of:

\[F(Tₛ) - Jᵀ = 0\]

where $Jᵀ$ are the fluxes at the top of the interface (turbulent + radiative), and $F$ is the internal diffusive flux dependent on the interface temperature itself.

Note that all fluxes positive upwards.

StorageBasedDryLayerDepth(FT = Oceananigans.defaults.FloatType;
                          maximum_dry_layer_depth,
                          dry_layer_onset_saturation,
                          dry_layer_exponent = 2)

Diagnostic dry-layer depth δᵛ as a function of land saturation 𝒮:

\[δᵛ(𝒮) = δᵛ_{max} \left[1 - \min(𝒮 / 𝒮ᶜ, 1)\right]^η\]

δᵛ = 0 when 𝒮 ≥ 𝒮ᶜ (wet branch), growing toward δᵛ_max as the slab dries.

Here dry_layer_onset_saturation (𝒮ᶜ) is the slab saturation at which the dry layer first appears.

Diagnosing the dry-surface-layer thickness from near-surface moisture follows the CLM5 scheme of Swenson and Lawrence (2014) (also used by ClimaLand, with maximum depth ≈ 15 mm after Shokri and Or 2011); this closure differs only in using the slab saturation 𝒮 as the moisture variable and a power-law shape.

SubgridVelocityCorrection(FT = Float64;
                          convective = ConvectiveGustiness{FT}(),
                          mesoscale = nothing)

Composition of the empirical subgrid velocity enhancements applied to the bulk velocity, $U² = Δu² + Δv² + Uᶜ² + Uᵐ²$:

• convective: a state-dependent convective gustiness, by default ConvectiveGustiness.
• mesoscale: a static velocity scale [m/s] representing mesoscale variability unresolved on coarse grids, e.g. mahrt_sun_subgrid_velocity(Δx). Default nothing (no contribution).

Either slot may be nothing, a Number [m/s], or a formulation implementing vsgs²(correction, u★, b★, h_bℓ).

ThreeEquationHeatFlux(FT::DataType = Oceananigans.defaults.FloatType;
                      heat_transfer_coefficient = 0.0095,
                      salt_transfer_coefficient = heat_transfer_coefficient / 35,
                      friction_velocity = 0.002)

Construct a ThreeEquationHeatFlux with the specified parameters.

Default values follow Shi et al. (2021) with $R = \alpha_h / \alpha_s = 35$.

Keyword Arguments

• heat_transfer_coefficient: turbulent heat exchange coefficient $\alpha_h$. Default: 0.0095.
• salt_transfer_coefficient: turbulent salt exchange coefficient $\alpha_s$. Default: $\alpha_h / 35 \approx 0.000271$.
• friction_velocity: friction velocity value or formulation. Default: 0.002.

ThreeEquationHeatFlux{FT, U}

Three-equation formulation for sea ice-ocean heat flux.

This formulation solves a coupled system for the interface temperature and salinity:

1. Heat balance: $\rho c_p \gamma_T (T - T_b) = ℰ q$
2. Salt balance: $\gamma_S (S - S_b) = q (S_b - S_i)$
3. Freezing point: $T_b = T_m(S_b)$

where $T_b$ and $S_b$ are the interface temperature and salinity, $\gamma_T = \alpha_h u_*$ and $\gamma_S = \alpha_s u_*$ are turbulent exchange velocities, $L$ is the latent heat of fusion, and $q$ is the melt rate (computed, not input).

Fields

• heat_transfer_coefficient::FT: turbulent heat exchange coefficient $\alpha_h$ (dimensionless)
• salt_transfer_coefficient::FT: turbulent salt exchange coefficient $\alpha_s$ (dimensionless)
• internal_heat_flux::FT: diffusive flux inside the sea ice (ConductiveFlux)
• friction_velocity::U: friction velocity value or formulation (constant Number or MomentumBasedFrictionVelocity)

Example

using NumericalEarth.EarthSystemModels: ThreeEquationHeatFlux

flux = ThreeEquationHeatFlux()

# output
ThreeEquationHeatFlux{Nothing}
├── heat_transfer_coefficient: 0.0095
├── salt_transfer_coefficient: 0.00027142857142857144
└── friction_velocity: 0.002

References

• Holland and Jenkins (1999): Holland, D. M., & Jenkins, A. (1999). Modeling thermodynamic ice–ocean interactions at the base of an ice shelf. Journal of Physical Oceanography, 29(8), 1787-1800.
• Shi et al. (2021): Shi, X., Notz, D., Liu, J., Yang, H., & Lohmann, G. (2021). Sensitivity of Northern Hemisphere climate to ice-ocean interface heat flux parameterizations. Geosci. Model Dev., 14, 4891-4908.

atmosphere_land_interface(grid, atmosphere, land;
                          fluxes               = default_atmosphere_land_fluxes(land, eltype(grid)),
                          temperature          = BulkTemperature(),
                          velocity_difference  = RelativeVelocity(),
                          specific_humidity    = default_al_specific_humidity(land))

Build the atmosphere–land interface on grid from atmosphere and land with the given turbulent-flux closure, interface-temperature model, atmosphere-relative velocity model, and specific-humidity formulation. Pass the result as atmosphere_land_interface = ... to ComponentInterfaces / AtmosphereLandModel to override the default.

The flux closure's roughness lengths and zero-plane displacement may be per-cell Fields at (Center, Center, Nothing) on grid — for example from urban_roughness or a canopy roughness closure — localized to each cell before the Monin–Obukhov solve.

compute_sea_ice_ocean_fluxes!(coupled_model)

Compute heat, salt, and momentum fluxes at the sea ice-ocean interface.

This function computes:

• Frazil heat flux: heat released when ocean temperature drops below freezing (all formulations)
• Interface heat flux: heat flux from ocean to ice, computed using the specified formulation
• Salt flux: salt exchange due to ice growth/melt
• Momentum stresses: ice-ocean momentum transfer

The interface heat flux formulation is determined by coupled_model.interfaces.sea_ice_ocean_interface.flux_formulation.

large_yeager_stability_functions(FT = Float64)

NCAR/Large-Yeager (2004) stability functions combining:

• Unstable: Paulson (1970) with γ = 16
• Stable: linear ψ = -5ζ, bounded at |ζ| ≤ 10

Used for OMIP-2 protocol compliance.

mahrt_sun_subgrid_velocity(Δx; threshold = 5e3)

Return the mesoscale subgrid velocity [m/s] for grid spacing Δx [m] following Mahrt and Sun (1995), as implemented in the revised MM5 surface layer scheme of Jiménez et al. (2012):

\[V_{sg} = 0.32 \, [\max(Δx / 5000 - 1, 0)]^{0.33}\]

The enhancement is zero for Δx ≤ threshold (default 5 km). Pass the result as the mesoscale slot of SubgridVelocityCorrection.

surface_layer_diagnostics(model; ...)
surface_layer_diagnostics(
    model,
    interface;
    velocity_height,
    temperature_height,
    specific_humidity_height
)

Return the named tuple (; u, v, T, q) of lazy KernelFunctionOperations that evaluate the Monin–Obukhov profile of interface at the requested heights above the interface: the wind components at velocity_height, the in-situ temperature at temperature_height, and the specific humidity at specific_humidity_height, all in meters. The heights default to the WMO observation heights, velocity_height = 10 and temperature_height = specific_humidity_height = 2, and interface defaults to model.interfaces.atmosphere_land_interface.

Each variable is the similarity profile evaluated with the scales that interface stored at the last flux solve,

\[φ(z) = φₛ + \frac{φ_★}{ϰ} Π(z - d, ℓ_φ, L_★) ,\]

with the wind direction taken from the stored stress and, over the ocean, the surface current added back. The profiles are NaN inside the roughness sublayer z - d ≤ ℓ_φ, and evaluating them at or below the zero-plane displacement d is an error.

The temperature is in situ: θ_★ scales a potential temperature referenced to z = 0, so the adiabatic increment g z / cᵖᵐ is subtracted at the geometric height z.

Where the flux solve enhances the bulk velocity with subgrid velocities (convective gustiness, a mesoscale scale), that enhancement is carried by u_★ and is not undone here, so u and v exceed the resolved wind. T and q are unaffected: their scales are built from the resolved increments.

PressureLevelGrid

Type alias for any underlying grid whose vertical coordinate is a PressureLevelVerticalDiscretization.

PressureLevelVerticalDiscretization

A vertical discretization for pressure-level reanalysis data on a LatitudeLongitudeGrid. Per-cell heights are geopotential[i, j, k] / gravitational_acceleration; znode(i, j, k, grid) returns that directly.

Two znodes paths:

• znodes(grid, ℓ...) returns the column-mean z profile as a 1-D Vector{Float64} of length Nz — the representative axis that plot recipes, Lz, and length consumers expect when only the grid is in hand.
• znodes(field) returns the column-mean Vector{Float64} when the field has no horizontal extent (Flat topology or Nothing horizontal locations, e.g. after mean(...; dims=(1, 2))); otherwise it returns a 3-D Field of per-cell heights at the field's own location.

geopotential (units m²/s²) is a 3-D Field or a TimeSeriesInterpolation over a FieldTimeSeries. The former gives a static z-coordinate; the latter gives a time-evolving one driven by an attached Clock.

The LatitudeLongitudeGrid constructor needs a value for Lz; we compute it as extrema(geopotential) / g inside generate_coordinate.

PressureLevelVerticalDiscretization(geopotential;
                                    gravitational_acceleration,
                                    surface_geopotential = nothing)

Build a discretization backed by per-column geopotential (m²/s²). znode divides by gravitational_acceleration at read time.

If surface_geopotential is provided (a 2-D Field, m²/s²), columns are clipped so that geopotential[i,j,k] ≥ surface_geopotential[i,j]. Required when the source is ERA5 pressure-level data, because sub-surface levels are filled with non-physical extrapolations that would break the column-monotonicity assumed by _fractional_indices. Clipping fixes only the coordinate: the raw sub-surface data stays on those levels, so column_fractional_z_index also clamps interpolation to the first above-ground level (never sampling it). The clip source is retained on the discretization and exposed through surface_elevation.

surface_elevation(grid)

Return the surface elevation (m) of the orography underlying grid as a two-dimensional (Center, Center, Nothing) field on the source's native horizontal grid — for a PressureLevelGrid, the clip-source surface geopotential divided by the gravitational acceleration. Return nothing when the surface elevation is unknown (non-pressure-level grids, or a discretization built without surface_geopotential).

abstract type AbstractLand end

Marker abstract type for prognostic NumericalEarth land components (SlabLand and any future multi-layer variants). PrescribedLand is a prescribed-forcing component and subtypes AbstractPrescribedComponent instead.

abstract type AbstractUrbanRoughness end

Supertype of the urban morphometric roughness closures. A closure maps the plan-area index λᵖ and mean building height h to the momentum roughness length ℓᵐ and zero-plane displacement d via aerodynamic_parameters; it is also callable as closure(λᵖ, h).

BucketHydrology(FT = Oceananigans.defaults.FloatType;
                maximum_water_storage = 150)

Single-bucket land hydrology. The land water mass per area Mˡᵃ (water_storage) evolves under P − E, clamped to [0, Mˡᵃ⁺] where Mˡᵃ⁺ is the bucket capacity (maximum_water_storage, the soil-science "field capacity"; 150 kg m⁻² ≈ 15 cm of equivalent liquid water).

The atmosphere-facing surface_saturation is the continuous saturation 𝒮 = clamp(Mˡᵃ/Mˡᵃ⁺, 0, 1). The interface surface-humidity models derive their own availability from it: BulkHumidity is the binary limit (saturated where 𝒮 > 0), FractionalHumidity scales saturation by β(𝒮), and SkinHumidity solves a vapor-flux balance. The availability parameterization belongs on the interface, not here.

maximum_water_storage may be a scalar or a per-cell Field.

CosbyConductivity(FT = Oceananigans.defaults.FloatType;
                  intercept = -0.884,
                  sand_coefficient = 0.0153)

Saturated hydraulic conductivity from sand, after Cosby et al. (1984) Table 5:

\[\log_{10} K^{+} = a + b\,S,\]

with sand S in % and K⁺ in inch hour⁻¹. Read with saturated_conductivity, which takes a mass fraction (kg/kg) and returns m s⁻¹. The defaults are the published fit, regressed on 1,448 US samples across 11 texture classes.

DarcyDeepLiquidFlux(exchange_length, liquid_density)

Darcy exchange between the slab bottom and a deep reservoir held at pressure head Πᵈ, set via VariablySaturatedHydrology(deep_pressure_head = …). The exchange length ℓ_l (m) separates the slab bottom from the deep reservoir; the deep height is z_D = z_b − ℓ_l. With hydraulic heads h_b = z_b + Π_b and h_D = z_D + Πᵈ,

\[J^{lb} = \rho^l K_b\,\frac{h_D - h_b}{\ell_l}.\]

struct EmpiricalFrontalArea{FT} <: NumericalEarth.Lands.AbstractFrontalAreaEstimator

Take λᶠ = q·λᵖ² + l·λᵖ, regressed by Kanda et al. (2013) on the building databases of Tokyo, Nagoya, Toulouse, Berlin, Salt Lake City and Los Angeles (their eq. 1, fitted over 0.05 < λᵖ < 0.45). Unlike IsotropicFrontalArea this assumes no element geometry, and it respects the λᶠ < 2λᵖ envelope they observe across real cities. The default, since a mean-height dataset never measures λᶠ.

ForceRestoreEnergy(FT = Oceananigans.defaults.FloatType;
                   dry_heat_capacity = 1480 * 1500 * 0.10,
                   liquid_heat_capacity = 4186,
                   deep_temperature = 280.0,
                   deep_time_scale = 12 * 3600)

Slab energy model that relaxes the bulk land temperature toward a prescribed deep climatological temperature deep_temperature (math Tᵈᵉᵉᵖ) on the time scale deep_time_scale (math τ), in addition to the surface energy flux.

deep_temperature may be a Number, per-cell AbstractField, or any state-indexable climatology (e.g. FieldTimeSeries).

FreeDrainageFlux(liquid_density)

Free-drainage bottom boundary. ∂z Π = 0 so ∂z h = 1 and Jˡ = -ρˡ K. ρˡ is the intrinsic liquid-water density (kg m⁻³).

HYPRESPedotransfer(FT = Oceananigans.defaults.FloatType;
                   organic_matter = 1,
                   topsoil_depth = 0.3,
                   residual_liquid_fraction = 0.01,
                   pore_connectivity_exponent = 0.5,
                   regression_coefficients = HYPRES_REGRESSION)

The Wösten et al. (1999) continuous pedotransfer function (their Table 5), fitted to HYPRES, the database of HYdraulic PRoperties of European Soils: closed-form regressions in clay %, silt %, organic matter %, bulk density (g/cm³) and a topsoil/subsoil flag. Clay is the < 2 μm fraction and silt the 2–50 μm fraction, matching the USDA split the soil datasets report. Wösten's class regressions, which bin soil into eleven texture classes, are not implemented.

Soil shallower than topsoil_depth (m) is topsoil. organic_matter (%) is a uniform value. residual_liquid_fraction (θʳ) and pore_connectivity_exponent (ηᴷ) are constants, since HYPRES supplies no usable regression for either: the defaults are Wösten's Table 4 class fit for θʳ and the Mualem exponent an unreduced matching point pairs with. regression_coefficients is a HYPRESRegression.

Predicted θs is returned as the porosity ν. Units are converted to model units: αᵃᵉ (cm⁻¹) → m⁻¹, K₀ (cm day⁻¹) → m s⁻¹.

Texture and organic matter are floored at 1 % and bulk density held within 0.5–2.0 g/cm³, because they enter as 1/x and ln x. Wösten warns the regressions "should not be used for the assignment of hydraulic properties to soils outside Europe".

InfiltrationCapacityRunoff(infiltration_capacity)

Cap the downward infiltration rate at infiltration_capacity (kg m⁻² s⁻¹, positive magnitude). Any precipitation exceeding the cap becomes surface runoff:

\[J^{ls} = \max(-P^l, -J^l_{cap}), \qquad R^{\mathrm{sfc}} = J^{ls} - (-P^l) \ge 0.\]

No subsurface runoff.

IsotropicFrontalArea()

Take λᶠ = λᵖ, exact for cubical elements, whose frontal and plan area densities are equal — the shape the array constant A of MorphometricRoughness was fitted to.

LinearReservoirDrainage(drainage_time_scale, equilibrium_storage)

Cheap linear-reservoir drainage:

\[J^{lb} = -\frac{\max(M^{la} - M_{eq}, 0)}{\tau_l}.\]

Storage above equilibrium_storage drains out the bottom on time scale τ_l.

MorphometricRoughness(FT = Oceananigans.defaults.FloatType; kw...)

Construct a MorphometricRoughness closure.

• array_constant (A) — default 4.43, the staggered-array fit. Under UniformHeight the square-array fit is 3.59, but it pairs with correction_factor = 0.55: sheltering is not captured by A alone, and 3.59 on its own overpredicts ℓᵐ by about a factor of two. Leave both at their defaults under VariableHeight, whose constants were regressed with A = 4.43 and β = 1 held fixed for either array type.
• drag_coefficient (Cᴰ) — building drag coefficient. Default 1.2, for a sharp-edged cube; rounded or cylindrical elements are lower by up to a factor of two.
• correction_factor (β) — lumped drag correction, folding in the velocity-profile shape, turbulence intensity and length scale, incidence angle and corner rounding. Default 1 reproduces the staggered array; 0.55 is the square-array value.
• von_karman_constant (ϰ) — default 0.4.
• frontal_area_estimator — frontal-area estimator: EmpiricalFrontalArea (default), IsotropicFrontalArea for the idealized cube array.
• height_distribution — height-spread treatment: VariableHeight (default) or UniformHeight for the idealized equal-height array.
• bare_soil_roughness (ℓˢᵒⁱˡ, m) — momentum roughness length where the built fraction vanishes. Default 0.03.
• minimum_built_fraction — built-fraction floor below which the cell reduces to bare soil. Default 0.01.
• maximum_displacement_ratio — ceiling on the packing displacement ratio d/h that feeds the drag partition. Default 1, the value at which d/h → 1 and ℓᵐ → 0 as the array closes into a smooth surface. Capping below 1 freezes the (1 − d/h) suppression and makes ℓᵐ turn back up in the skimming limit. Under VariableHeight the returned displacement is bounded by hᵐᵃˣ rather than by this ratio.

julia> using NumericalEarth.Lands

julia> MorphometricRoughness()
MorphometricRoughness{Float64} with EmpiricalFrontalArea() and VariableHeight()

julia> MorphometricRoughness(height_distribution = UniformHeight())
MorphometricRoughness{Float64} with EmpiricalFrontalArea() and UniformHeight()

struct MorphometricRoughness{FT, E, H} <: NumericalEarth.Lands.AbstractUrbanRoughness

Morphometric roughness closure. Maps (λᵖ, h) to (ℓᵐ, d) from the plan-area index and mean building height alone: d from the plan-area packing, ℓᵐ from a drag partition over the frontal area λᶠ estimated by frontal_area_estimator, both then corrected for the spread of element heights by height_distribution. ℓᵐ is non-monotonic in λᵖ (isolated → wake-interference → skimming).

Formulated by Macdonald et al. (1998) for idealized obstacle arrays and extended to height-heterogeneous real cities by Kanda et al. (2013).

NoDeepLiquidFlux()

Zero deep liquid flux. The slab is closed at the bottom.

NoRunoff()

No runoff. All precipitation infiltrates (Jˡˢ = −Pˡ), no subsurface export.

PrescribedLand(freshwater_flux; clock=nothing, river_routing=nothing)

Return a PrescribedLand from a NamedTuple of FieldTimeSeries freshwater fluxes (e.g. rivers, icebergs).

Keyword Arguments

• clock: the component's Clock. Default: nothing, a Clock matching the element type of freshwater_flux.
• river_routing: NamedTuple of RiverRouting, one per component, scattering each onto coastal ocean cells conserving volume (see GloFASPrescribedLand, JRA55PrescribedLand). Default: nothing, which interpolates the flux pointwise onto the ocean grid.

RiverRouting

A static map from river-mouth cells on a forcing dataset's native grid to the active (wet) cells of a target ocean grid, used to deposit river discharge as a conservative freshwater mass flux (kg m⁻² s⁻¹).

Contributions are grouped by destination ocean cell so the scatter writes each ocean cell exactly once (no atomics). For destination cell c, the contributing river mouths are contribution_outlet_{i,j}[offsets[c]:offsets[c+1]-1] with contribution_weight = outlet_weight / Aᵒᶜᵉᵃⁿ (see build_river_routing), chosen so the area integral of the deposited flux equals the total mass delivered.

SlabLand{FT, G, Clk, T, W, B, F, D, E, H}

A composable slab land-surface component. The default configuration — SlabEnergy + BucketHydrology — is the classic bucket slab introduced by Manabe (1969). Replace either axis to swap in a different energy or hydrology closure. Aerodynamic roughness lengths are a property of the atmosphere-land flux closure (atmosphere_land_fluxes), not of the land model.

Fields

• grid : Oceananigans grid (typically (Nx, Ny, Flat)).
• clock : Oceananigans.TimeSteppers.Clock.
• temperature : prognostic bulk land temperature T (K).
• water_storage : prognostic land water mass per area Mˡᵃ (kg m⁻²).
• saturation : diagnostic surface saturation 𝒮 = Mˡᵃ/Mˡᵃ⁺ ∈ [0, 1] (–).
• fluxes : NamedTuple of flux/forcing Fields the coupler writes.
• diagnostics : NamedTuple of closure-owned diagnostic Fields.
• energy : an AbstractEnergyBalance (parameters).
• hydrology : an AbstractHydrology (parameters).

SlabLand(grid;
         energy    = SlabEnergy(eltype(grid)),
         hydrology = BucketHydrology(eltype(grid)),
         clock     = Clock{eltype(grid)}(time = 0))

Construct a SlabLand with the chosen closures. The prognostic temperature and water_storage fields and the diagnostic saturation field are allocated on grid; the flux accumulators the coupler writes are sized from each closure's flux_variables declaration.

UniformHeight()

Treat every roughness element as being of the mean height h, so (ℓᵐ, d) follow from the plan-area and frontal-area indices alone. The idealized obstacle-array limit of Macdonald et al. (1998).

VanGenuchtenConductivity(FT = Oceananigans.defaults.FloatType;
                         matching_point_conductivity, pore_size_uniformity,
                         pore_connectivity_exponent = 1//2,
                         water_viscosity = WaterViscosity(FT))

Unsaturated hydraulic conductivity as a function of saturation 𝒮, combining the Mualem (1976) pore-bundle model with the van Genuchten (1980) retention shape:

\[K(\mathcal S) = K_0\,\mathcal S^{\eta^K}\left[1 - (1 - \mathcal S^{1/\mathscr{m}})^{\mathscr{m}}\right]^2, \qquad \mathscr{m} = 1 - 1/\mathscr{n}.\]

matching_point_conductivity (K₀, m s⁻¹) is the conductivity this curve reaches at 𝒮 = 1, pore_size_uniformity (𝓃) must match the retention curve's, and pore_connectivity_exponent (ηᴷ, –) is the exponent on saturation: a larger value throttles conductivity more steeply as the soil drains. Each may be a scalar or a Field (see soil_hydraulic_properties).

water_viscosity is a WaterViscosity scaling K with the slab temperature; pass nothing for an isothermal conductivity.

VanGenuchtenRetention(FT = Oceananigans.defaults.FloatType;
                      inverse_air_entry_head, pore_size_uniformity)

Empirical soil-water retention curve mapping liquid pore fraction θˡ (or saturation 𝒮) to soil matric pressure head Π_m (m, negative in unsaturated soil), following van Genuchten (1980):

\[\Pi_m(\mathcal S) = -\frac{1}{\alpha^{\mathrm{ae}}}\left[\mathcal S^{-1/\mathscr{m}} - 1\right]^{1/\mathscr{n}}, \qquad \mathscr{m} = 1 - 1/\mathscr{n}.\]

inverse_air_entry_head (αᵃᵉ, m⁻¹) is the reciprocal of the air-entry pressure head, the suction at which the largest pores begin to empty. pore_size_uniformity (𝓃, –) sets how narrow the band of suctions is over which the soil drains. Each may be a scalar or a Field (see soil_hydraulic_properties).

VariableHeight(FT = Oceananigans.defaults.FloatType; kw...)

Construct a VariableHeight distribution. All four defaults are the Method 1 regressions of Kanda et al. (2013).

• displacement_constants (a₀, b₀, c₀) — displacement fit, their eq. 10. Default (1.29, 0.36, -0.17); their Method 2 gives (0.86, 0.28, -0.18).
• roughness_constants (a₁, b₁, c₁) — roughness rescaling, their eq. 12. Default (0.71, 20.21, -0.77); their Method 2 gives (0.93, 8.93, 4.68).
• height_spread_constants (s₁, s₀) — σʰ = s₁·h + s₀ (m), their eq. 2. Default (1.05, -3.7), whose zero crossing near h = 3.5 m is about one story.
• maximum_height_constants (m₁, m₂) — hᵐᵃˣ = m₁·(σʰ)^m₂ (m), their eq. 3. Default (12.51, 0.77). The most scattered of the four, since a city's tallest building is often an outlying landmark.

julia> using NumericalEarth.Lands

julia> VariableHeight()
VariableHeight()

struct VariableHeight{FT} <: NumericalEarth.Lands.AbstractHeightDistribution

Correct the uniform-height (ℓᵐ, d) for the spread of roughness-element heights, after Kanda et al. (2013). The tallest elements absorb the most momentum, so d is referenced to hᵐᵃˣ rather than to h, and ℓᵐ is rescaled by a polynomial in λᵖ σʰ / h.

A mean-height dataset measures neither σʰ nor hᵐᵃˣ, so both are derived from h by the regressions the same study fitted to real cities (their eqs. 2 and 3). Pass measured per-cell statistics instead wherever a footprint-level dataset supplies them.

VariablySaturatedHydrology(FT = Oceananigans.defaults.FloatType;
                           slab_depth,
                           porosity,
                           residual_liquid_fraction = 0,
                           storage_height,
                           liquid_density = 1000,
                           retention_curve,
                           hydraulic_conductivity,
                           deep_liquid_flux = NoDeepLiquidFlux(),
                           deep_pressure_head = 0,
                           runoff = NoRunoff())

Conservative slab hydrology with an augmented liquid fraction ϑˡ. The storage variable Mˡᵃ evolves under signed fluxes (positive upward at the surface; positive upward at the bottom).

slab_depth, porosity, residual_liquid_fraction, and storage_height each accept either a scalar (uniform) or a Field that varies grid point by grid point (e.g. from a soil-texture map).

• slab_depth (hˡᵃ, m) — slab thickness.
• porosity (ν) — total pore fraction.
• residual_liquid_fraction (θʳ) — minimum pore liquid (default 0).
• storage_height (hˢˢ, m) — saturated storage height; the reciprocal of the specific storage (1/Sₛ), i.e. the pressure head built per unit fractional over-saturation.
• retention_curve — e.g. VanGenuchtenRetention.
• hydraulic_conductivity — e.g. VanGenuchtenConductivity.
• deep_liquid_flux — bottom-boundary closure: NoDeepLiquidFlux, FreeDrainageFlux, DarcyDeepLiquidFlux, or LinearReservoirDrainage.
• deep_pressure_head — the deep-reservoir pressure head (m), passed to the deep-flux closure as Πᵈ. Default 0.
• runoff — runoff closure: NoRunoff or InfiltrationCapacityRunoff.

WaterCoupledEnergy(FT = Oceananigans.defaults.FloatType;
                   dry_heat_capacity = 1480 * 1500 * 0.10,
                   liquid_heat_capacity = 4186,
                   reference_temperature = 273.15,
                   deep_temperature = 280,
                   deep_conductance = nothing,
                   deep_time_scale = nothing,
                   advect_deep_liquid_energy = true,
                   advect_surface_liquid_energy = false)

Force-restore energy closure with Mˡᵃ-dependent heat capacity and conservative treatment of water-energy advection across the slab boundaries.

Use either deep_conductance (Λᵈᵉᵉᵖ, W m⁻² K⁻¹) or deep_time_scale (τᵈᵉᵉᵖ, s) to set the deep restoring strength; exactly one must be supplied. With deep_time_scale, the effective conductance is cˡᵃ(Mˡᵃ)/τᵈᵉᵉᵖ (matches ForceRestoreEnergy).

dry_heat_capacity accepts a Number (uniform) or a per-cell AbstractField — for example pavement carrying building thermal mass next to soil that does not.

reference_temperature is the reference for internal energy eˡ(T) = cˡ(T − Tᵣ). The temperature update is invariant to Tᵣ: fluxes whose advective switch is disabled enter and leave at the slab temperature, so Tᵣ cancels exactly (up to the positivity floor on Mˡᵃ). A standard choice is the triple point 273.15 K.

WaterViscosity(FT = Oceananigans.defaults.FloatType;
               activation_temperature = 507.88,
               pole_temperature = 149.3,
               reference_temperature = 293)

Temperature dependence of the dynamic viscosity of water, as the factor by which it scales the hydraulic conductivity of soil:

\[\Theta(T) = \exp\!\left[\frac{T_1}{T_0 - T_2} - \frac{T_1}{T - T_2}\right],\]

with T₁ the activation_temperature and T₂ the pole_temperature, the law of Deck et al. (2026) (their Equations A19–A20) after Hansson et al. (2004). Read with viscosity_correction.

Θ is unity at reference_temperature (T₀), the temperature at which matching_point_conductivity holds; laboratory K₀ is measured near 20 °C.

WeynantsPedotransfer(FT = Oceananigans.defaults.FloatType;
                     organic_carbon = 0.58,
                     regression_coefficients = WEYNANTS_REGRESSION)

The Weynants et al. (2009) pedotransfer function: closed-form regressions mapping clay %, sand %, bulk density and organic carbon % to all six Mualem–van Genuchten parameters, fitted in a single inversion against measured retention and conductivity curves.

organic_carbon (% by weight, per the Weihermüller et al. (2017) erratum) is a uniform value, since the 30 m texture datasets do not carry it; it enters αᵃᵉ and K₀ only. 0.58 % is the same soil as HYPRESPedotransfer's 1 % organic matter under the van Bemmelen convention that organic matter is 58 % carbon. regression_coefficients is a WeynantsRegression.

Predicted θs is returned as the porosity ν and θʳ is zero. Units are converted to model units: αᵃᵉ (cm⁻¹) → m⁻¹, K₀ (cm day⁻¹) → m s⁻¹. K₀ is a matrix value: the fit excluded measurements wetter than 6 cm suction, so macropore flow is not in it.

Predictors are clamped to the fit's own ranges — clay ≤ 54.5 %, sand 5.6–97.8 %, ρᵇ 0.89–1.77 g/cm³, OC 0.01–6.6 % — so a soil past those bounds receives the parameters of the boundary.

surface_temperature(energy::AbstractEnergyBalance, land) -> AbstractField

Return the field of surface temperatures the atmosphere reads as a BC.

SlabEnergy(FT = Oceananigans.defaults.FloatType;
           dry_heat_capacity = 1480 * 1500 * 0.10,
           liquid_heat_capacity = 4186)

Pure slab energy balance — the τ → ∞ limit of ForceRestoreEnergy, with no deep restoring term. Returns a ForceRestoreEnergy configured so the restoring contribution vanishes: deep_time_scale = Inf (so (Tᵈ − T)/τ = 0) and deep_temperature is set to a finite sentinel that is never used.

aerodynamic_parameters(
    c::NumericalEarth.Lands.AbstractUrbanRoughness,
    cell
) -> Tuple{Any, Any}

Evaluate closure from a per-cell cell NamedTuple, reading its plan_area_index and mean_building_height keys.

aerodynamic_parameters(
    c::MorphometricRoughness{FT},
    λᵖ,
    h
) -> Tuple{Any, Any}

Momentum roughness length ℓᵐ and zero-plane displacement d (meters) of a built-up surface with plan-area index λᵖ and mean building height h, under closure. Returns (ℓᵐ, d). Endpoints reduce cleanly: λᵖ → 0 returns the bare-soil (ℓˢᵒⁱˡ, 0); λᵖ → 1 is the skimming limit. An invalid (NaN/negative) input returns (NaN, NaN). Under VariableHeight the displacement is bounded by the tallest element rather than by h, and routinely exceeds the mean building height.

aerodynamic_parameters(
    c::MorphometricRoughness{FT},
    λᵖ,
    h,
    σʰ,
    hᵐᵃˣ,
    λᶠ
) -> Tuple{Any, Any}

Momentum roughness length ℓᵐ and zero-plane displacement d (meters) from measured per-cell morphometry: the plan-area index λᵖ, mean building height h, height standard deviation σʰ, tallest element hᵐᵃˣ and frontal-area index λᶠ, as a footprint-level dataset supplies them. The frontal_area_estimator and the σʰ/hᵐᵃˣ regressions of VariableHeight drop out; the drag partition and the height-spread corrections themselves remain. hᵐᵃˣ is floored at h, and the displacement parameter X = (σʰ + h)/hᵐᵃˣ keeps its [0, 1] clamp, so a low-biased maximum height cannot invert the correction. Under UniformHeight σʰ and hᵐᵃˣ are ignored and the result is the obstacle-array closure with measured λᶠ.

build_river_routing(target_grid, outlet_i, outlet_j, outlet_λ, outlet_φ, outlet_weight;
                    maximum_search_radius = 5, spread_radius = 1.2, maximum_spread_cells = nothing)

Map each mouth at (outlet_λ, outlet_φ) onto the active ocean cells of target_grid, returning a RiverRouting that deposits outlet_weight[n] * value[outletₙ] / Aᵒᶜᵉᵃⁿ. The weight is the freshwater density for a volumetric discharge (m³ s⁻¹), the source-cell area for a per-area mass flux (kg m⁻² s⁻¹).

coastal_outlet_indices(discharge)

Return (outlet_i, outlet_j, outlet_λ, outlet_φ) for the river-mouth cells of a discharge Field whose ocean cells are NaN (the GloFAS convention). A river mouth is a finite (land/river) cell with at least one NaN (ocean) horizontal neighbor — i.e. the point where the routed river network meets the coast.

compute_aerodynamic_roughness!(
    ℓᵐ,
    d,
    closure::NumericalEarth.Lands.AbstractUrbanRoughness,
    properties,
    grid
) -> Tuple{Any, Any}

Fill the momentum roughness length ℓᵐ and zero-plane displacement d (meters) in place by applying an urban closure (AbstractUrbanRoughness) over every cell of grid. properties is a NamedTuple of the closure's per-cell inputs — scalars, arrays or Fields; the urban closures expect plan_area_index and mean_building_height.

routable_grid(grid) -> Bool

Whether grid has a coastline to route river mouths onto.

saturated_conductivity(c::CosbyConductivity, sand) -> Any

Macropore-inclusive saturated hydraulic conductivity K⁺ (m s⁻¹) from sand mass fraction (kg/kg), for an infiltration cap such as InfiltrationCapacityRunoff.

using NumericalEarth

## sand, loam and clay, in mm hour⁻¹
K⁺ = [saturated_conductivity(CosbyConductivity(), sand) for sand in (0.92, 0.43, 0.20)]
round.(K⁺ .* 3.6e6, digits = 1)

# output
3-element Vector{Float64}:
 84.8
 15.1
  6.7

soil_hydraulic_properties(
    sand,
    silt,
    clay,
    bulk_density;
    slab_depth,
    ptf,
    matching_heads
)

Reduce the 3-D texture (sand, silt, clay, kg/kg) and bulk_density (kg/m³) Fields to a NamedTuple of 2-D Field{Center, Center, Nothing} van Genuchten properties

(; porosity, residual_liquid_fraction, inverse_air_entry_head,
   pore_size_uniformity, matching_point_conductivity, pore_connectivity_exponent)

whose keys match the keyword arguments of VariablySaturatedHydrology, VanGenuchtenRetention and VanGenuchtenConductivity.

The pedotransfer function ptf is applied to each depth layer of the inputs' grid, at that layer's depth, and the layers inside [-slab_depth, 0] are combined so that the slab reproduces the thickness-weighted mean retention curve of the column: ν, θʳ and ηᴷ are thickness-weighted arithmetic means, K₀ is the harmonic mean of layers in series, and αᵃᵉ and 𝓃 are solved so the effective curve passes through the mean curve at the two suction heads matching_heads (m). The defaults are field capacity and the permanent wilting point, with field capacity at 1 m after Balsamo et al. (2009).

A cell that is NaN in any input is NaN in every predicted output; fill the texture first, e.g. with Field(metadatum; inpainting = NearestNeighborInpainting(n)).

urban_roughness(
    h,
    λᵖ,
    σʰ,
    hᵐᵃˣ,
    λᶠ;
    closure
) -> Tuple{Oceananigans.Fields.Field{Center, Center, Nothing, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}, Oceananigans.Fields.Field{Center, Center, Nothing, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}}

Momentum roughness length ℓᵐ and zero-plane displacement d (as Fields on the grid of h) from measured per-cell morphometry: the mean building height h, plan-area index λᵖ, height standard deviation σʰ, maximum height hᵐᵃˣ and frontal-area index λᶠ — the fields a footprint-level dataset such as GlobalBuildingFootprints3D aggregates. The closure (default MorphometricRoughness) is fed the measured height heterogeneity in place of its frontal-area estimator and σʰ/hᵐᵃˣ regressions. Where λᵖ → 0 the result reduces to a bare-soil roughness.

urban_roughness(
    h,
    λᵖ;
    closure
) -> Tuple{Oceananigans.Fields.Field{Center, Center, Nothing, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}, Oceananigans.Fields.Field{Center, Center, Nothing, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B}}

Momentum roughness length ℓᵐ and zero-plane displacement d (as Fields on the grid of h) from a mean building-height field h and a plan-area index field λᵖ, under closure (default MorphometricRoughness). Where λᵖ → 0 the result reduces to a bare-soil roughness.

NestedModel(parent, child::AbstractModel)

Pair a child Oceananigans model with a parent state (e.g. a 3D PrescribedAtmosphere). The result is itself an AbstractModel whose property access and integrator-protocol methods forward to the child, except for time_step!, which advances the child and then ticks the parent to match the child's new clock time.

Construct via child_model plus this wrapper, or use the NestedSimulation convenience that bundles model construction and Simulation wrapping.

NestedSimulation

Type alias for Simulation{<:NestedModel} — an Oceananigans Simulation whose model is a NestedModel. Useful for ::NestedSimulation dispatch.

NestedSimulation(parent, child_model::AbstractModel; simulation_kwargs...)

Convenience for Simulation(NestedModel(parent, child_model); simulation_kwargs...).

simulation_kwargs are forwarded to the Simulation constructor (Δt, stop_time, stop_iteration, verbose, etc.). Parent-sync is handled by NestedModel's time_step!, so no callback is installed — the result is a plain Oceananigans Simulation whose model happens to be a NestedModel.

blend_parent_terrain!(elevation, parent_elevation; width)

Blend the child elevation (a two-dimensional field) toward parent_elevation (the parent's surface elevation on the same grid) over the outermost width cells: the blend weight ramps linearly from the parent's elevation at the boundary to the child's width cells inward, so the terrain at the open boundaries matches the orography the parent state was produced with.

parent_boundary_conditions(grid;
                           variables,
                           sides    = (:west, :east, :south, :north),
                           schemes  = NamedTuple(),
                           bc_types = NamedTuple())

Build a NamedTuple of FieldBoundaryConditions driving the child's lateral boundaries directly from a parent's FieldTimeSeries.

Arguments

• grid: the child grid. (Passed for symmetry with the rest of the API; not used directly — Oceananigans picks up the boundary node coordinates during regularization.)

• variables: a NamedTuple mapping child field names to the parent FieldTimeSeries that should drive them, e.g. (u = parent.velocities.u, v = parent.velocities.v).

• sides: a tuple of Symbols naming which boundaries to drive. Choices are :west, :east, :south, :north, :bottom, :top.

• schemes: optional NamedTuple keyed by child field name giving the NormalFlowBoundaryCondition scheme (e.g. a PerturbationAdvection). Fields not listed default to scheme = nothing (stiff Dirichlet). Only consulted when the BC type for that field is NormalFlowBoundaryCondition.

• bc_types: optional NamedTuple keyed by child field name giving the BC constructor to use for that field. Each entry is either a single BC constructor applied to every side, or a per-side NamedTuple (keyed by :west/:east/:south/:north) selecting the BC type per side. Defaults to NormalFlowBoundaryCondition. For cell-centered scalars (e.g. ρ, ρθ in a compressible LAM), pass ValueBoundaryCondition to avoid the NormalFlowBoundaryCondition-on-Center asymmetric-halo behavior. For a momentum component, use per-side types — NormalFlowBoundaryCondition on the wall-normal side and ValueBoundaryCondition on the tangential side — so the tangential momentum is set (Dirichlet) rather than left to the normal-flow halo fill.

parent_forcings(; variables, rate, mask = 1)

Build a NamedTuple of Relaxation forcings driving the child interior toward the parent's FieldTimeSeries values.

Arguments

• variables: a NamedTuple mapping child field names to the parent FieldTimeSeries that should drive them, e.g. (u = parent.velocities.u, v = parent.velocities.v) — same convention as parent_boundary_conditions.

• rate: the relaxation rate (s⁻¹). May be a scalar (applied to all variables) or a NamedTuple keyed by child variable.

• mask: a callable (x, y, z) -> [0, 1], a Field, a scalar, or an Oceananigans-supplied mask object (e.g. GaussianMask). May be a NamedTuple for per-variable masks. Default 1 relaxes uniformly across the domain.

Returns

A NamedTuple of Oceananigans.Forcings.Relaxation ready to splat into the forcing kwarg of an Oceananigans / Breeze model constructor.

ChlorophyllOptics(FT = Oceananigans.defaults.FloatType;
                  clear_water_attenuation = 0.0232,
                  chlorophyll_scaling = 0.074,
                  chlorophyll_exponent = 0.674)

Parameters relating the blue-green absorption coefficient of sea water to its chlorophyll concentration, following Manizza et al. (2005). Defaults are the values given there. See absorption_coefficient.

PrescribedOcean(grid, times=[zero(grid)];
                density = 1025,
                heat_capacity = 4000,
                clock = Clock{FT}(time=0),
                sea_surface_temperature = default_prescribed_sst(grid, times),
                sea_surface_salinity = default_prescribed_sss(grid, times),
                velocities = default_prescribed_velocities(grid, times))

A prescribed ocean component for EarthSystemModel with sea surface temperature, salinity, and velocities prescribed as FieldTimeSeries.

The ocean state does not evolve in response to surface fluxes — it follows the prescribed data. Surface fluxes are still computed so the atmosphere feels the ocean.

Arguments

• grid: An Oceananigans grid for the ocean domain.
• times: Time instances for the prescribed data. Default: [zero(grid)] (constant).

Keyword Arguments

• density: Seawater density in kg/m³. Default: 1025.
• heat_capacity: Seawater specific heat capacity in J/(kg·K). Default: 4000.
• clock: Clock for tracking ocean time. Default: Clock{FT}(time=0).
• sea_surface_temperature: FieldTimeSeries for SST.
• sea_surface_salinity: FieldTimeSeries for SSS.
• velocities: NamedTuple of FieldTimeSeries for (u, v).

SlabOcean(grid;
          FT = eltype(grid),
          depth = 50,
          density = 1025,
          heat_capacity = 4000,
          clock = Clock{FT}(time=0))

A slab ocean component for EarthSystemModel. Represents the ocean as a single well-mixed layer of fixed depth, with temperature evolving in response to net surface heat fluxes:

∂T/∂t = Q / (ρ cₚ H)

where Q is the net downward surface heat flux (W/m²), ρ is the seawater density, cₚ is the heat capacity, and H is the slab depth.

Internally, the EarthSystemModel coupling assembles a temperature flux Jᵀ (in K m/s), so the tendency is computed as ∂T/∂t = -Jᵀ / H.

Use set!(ocean, T=value) to initialize the slab temperature after construction.

Arguments

• grid: An Oceananigans grid for the slab ocean domain.

Keyword Arguments

• FT: Floating point type for the slab ocean fields. Default: eltype(grid).
• depth: Depth of the slab in meters. Default: 50.
• density: Seawater density in kg/m³. Default: 1025.
• heat_capacity: Seawater specific heat capacity in J/(kg·K). Default: 4000.
• clock: Clock for tracking slab ocean time. Default: Clock{FT}(time=0) where FT is the float type of the grid.

TwoColorRadiation(grid; first_color_fraction = 0.58,
                        first_decay_scale = 0.35,
                        chlorophyll_optics = ChlorophyllOptics(eltype(grid)),
                        chlorophyll = equivalent_chlorophyll(chlorophyll_optics, 1 / 23))

Return TwoColorRadiation that computes the radiative flux divergence associated with a two-color radiation flux that decays according to Beer's law,

\[I(z) = ϵ₁ I₀ \exp(κ₁ z) + (1 - ϵ₁) I₀ \exp(κ₂ z)\]

where $I₀$ is the surface flux, $ϵ₁$ is the first_color_fraction, and $κ₁$, $κ₂$ are the absorption coefficients of the two colors. The first color is red and is absorbed within first_decay_scale of the surface; the second is blue-green and its absorption coefficient $κ₂$ follows the water's chlorophyll through absorption_coefficient.

chlorophyll is anything stateindex resolves: a number for globally uniform optics, a surface Field or array for a fixed spatial pattern, a function of (λ, φ, z, t), or a FieldTimeSeries, which is interpolated in time and so carries a seasonal cycle. Cyclical time indexing turns a twelve-month climatology, such as SeaWiFSMonthly, into optics that repeat every year.

The defaults are the red band of Manizza et al. (2005) and the chlorophyll whose decay scale is 23 m, which is the Jerlov Type I value of Paulson and Simpson (1977) for the clearest open-ocean water.

using NumericalEarth
using Oceananigans

grid = LatitudeLongitudeGrid(size=(4, 4, 4), longitude=(0, 360), latitude=(-60, 60), z=(-100, 0))

TwoColorRadiation(grid)

# output

TwoColorRadiation with red fraction 0.58 decaying over 0.35 m
└── blue-green chlorophyll: 0.147 mg m⁻³, decaying over 23.0 m

absorption_coefficient(optics::ChlorophyllOptics, chlorophyll)

Return the blue-green absorption coefficient in m⁻¹ implied by a chlorophyll concentration in mg m⁻³,

\[κ = a + b \, C^{n}\]

with $a$ the clear-water attenuation and $b$, $n$ the chlorophyll scaling and exponent of optics. Clear subtropical water, $C ≈ 0.05$, decays over about 30 m; a subpolar bloom, $C ≈ 1$, over about 10 m. equivalent_chlorophyll inverts the relation.

using NumericalEarth

round(1 / absorption_coefficient(ChlorophyllOptics(), 1.0), digits=1)

# output

10.3

equivalent_chlorophyll(optics::ChlorophyllOptics, absorption_coefficient)

Return the chlorophyll concentration in mg m⁻³ whose blue-green absorption coefficient is absorption_coefficient, the inverse of absorption_coefficient. Use it to express a Jerlov water type, whose optics are quoted as a decay scale, as a chlorophyll.

using NumericalEarth

round(equivalent_chlorophyll(ChlorophyllOptics(), 1 / 23), digits=3)

# output

0.147

ocean_simulation(grid; model = :hydrostatic, kwargs...)

Construct and return an ocean simulation tailored to grid. The model keyword selects the underlying Oceananigans model formulation:

• :hydrostatic (default) — builds a HydrostaticFreeSurfaceModel-based simulation with a free surface, CATKE vertical mixing, quadratic bottom drag (semi-implicit unless implicit_bottom_drag = false), and a TEOS-10 equation of state. See hydrostatic_ocean_simulation for the full kwarg list.

• :nonhydrostatic — builds a NonhydrostaticModel-based simulation suitable for LES (full 3D pressure, no free surface, no barotropic forcing). See nonhydrostatic_ocean_simulation for the full kwarg list.

Remaining kwargs are forwarded to the per-model builder; an unknown kwarg for the selected model raises the usual MethodError.

Examples

julia> using NumericalEarth, Oceananigans

julia> grid = RectilinearGrid(size = (10, 10, 6), extent = (1, 1, 1), halo=(6, 6, 5));

julia> ocean = ocean_simulation(grid)
Simulation of HydrostaticFreeSurfaceModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 30 minutes
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: Inf days
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│   └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
└── output_writers: OrderedDict with no entries

julia> les = ocean_simulation(grid; model = :nonhydrostatic, Δt = 2)
Simulation of NonhydrostaticModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── Next time step: 2 seconds
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: NaN days
├── stop_time: Inf days
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 4 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   ├── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
│   └── nan_checker => Callback of NaNChecker for u on IterationInterval(100)
└── output_writers: OrderedDict with no entries

river_mouth_vertical_diffusivity(
    grid,
    river_routing;
    river_mouth_diffusivity,
    river_mouth_mixing_depth
) -> Union{Oceananigans.TurbulenceClosures.VerticalScalarDiffusivity{Oceananigans.Utils.VerticallyImplicitTimeDiscretization, _A, _B, typeof(NumericalEarth.Oceans.masked_river_mouth_diffusivity)} where {_A, _B}, Oceananigans.TurbulenceClosures.VerticalScalarDiffusivity{Oceananigans.Utils.VerticallyImplicitTimeDiscretization, _A, _B, K} where {_A, _B, K<:(Oceananigans.TurbulenceClosures.DiscreteDiffusionFunction{_A, _B, _C, P, typeof(NumericalEarth.Oceans.masked_river_mouth_diffusivity)} where {_A, _B, _C, P<:(Oceananigans.Fields.Field{Center, Center, Center, Nothing, G, I, D, T, B, Nothing} where {G, I, D, T, B})})}}

A VerticalScalarDiffusivity adding river_mouth_diffusivity over the top river_mouth_mixing_depth metres of the cells river_routing discharges into, so a fresh plume trapped in the surface layer cannot drive the salinity to zero.

InterfaceRadiationFlux{F}

Container for the diagnostic radiative fluxes at an air–surface interface. The same struct type is instantiated per surface (ocean, sea ice, snow, ...).

Fields

• upwelling_longwave :: F ϵσT⁴
• downwelling_longwave :: F ϵℐꜜˡʷ (absorbed by the surface)
• downwelling_shortwave :: F (1−α)ℐꜜˢʷ (transmitted into the surface)

LatitudeDependentAlbedo([FT::DataType=Float64]; diffuse = 0.069, direct = 0.011)

Constructs a LatitudeDependentAlbedo object. The albedo of the ocean surface is assumed to be a function of the latitude, obeying the following formula (Large and Yeager, 2009):

α(φ) = α.diffuse - α.direct * cos(2φ)

where φ is the latitude, α.diffuse is the diffuse albedo, and α.direct is the direct albedo.

Arguments

• FT::DataType: The data type of the albedo values. Default is Float64.

Keyword Arguments

• diffuse: The diffuse albedo value. Default is 0.069.
• direct: The direct albedo value. Default is 0.011.

PrescribedRadiation(grid, times = [zero(grid)]; kwargs...)

Construct a PrescribedRadiation with zero downwelling shortwave and longwave fields on grid. Useful for emission-only mode (the surface radiates ϵσT⁴ but no incoming radiation is absorbed). All other keyword arguments are forwarded to the FTS-form constructor.

PrescribedRadiation{G, T, FT, SW, LW, S, IF, TI}

Top-level radiation component holding prescribed downwelling shortwave and longwave radiation as FieldTimeSeries, plus per-surface radiative properties (albedo, emissivity) and the Stefan–Boltzmann constant. Diagnostic radiative fluxes (one InterfaceRadiationFlux per surface) are populated by the apply_air_sea_*_radiative_fluxes! kernels at every step; interface_fluxes is nothing until the radiation is paired with an EarthSystemModel (which allocates the per-surface buffers on the exchange grid).

PrescribedRadiation(downwelling_shortwave, downwelling_longwave;
                    clock = nothing,
                    ocean_surface = SurfaceRadiationProperties(0.05, default_water_emissivity),
                    sea_ice_surface = SurfaceRadiationProperties(0.7, 1.0),
                    snow_surface = nothing,
                    land_surface = nothing,
                    stefan_boltzmann_constant = default_stefan_boltzmann_constant)

Construct a PrescribedRadiation component from FieldTimeSeries of downwelling shortwave and longwave radiation. Grid + times are inferred from the shortwave FTS.

Pass *_surface = nothing to omit that surface from surface_properties.

SeaIceAlbedo{FT, HI, HS, TS}

Sea ice albedo parameterization following the CCSM3 scheme (Briegleb et al. 2004).

Computes broadband albedo as a function of ice thickness, snow depth, and surface temperature. The scheme blends between bare ice and snow-covered albedos, with a temperature-dependent reduction near the melting point to implicitly represent melt pond formation.

Algorithm:

1. Base cold albedos: bare ice (0.53) and snow-covered (0.82)
2. Temperature reduction within 1C of melting: Δαice = 0.075, Δαsnow = 0.10
3. Thin-ice transition to ocean albedo below h_amin = 0.5 m
4. Snow cover interpolation: full snow albedo at hsnow > hsmin = 0.02 m

References:

• Briegleb, B.P., C.M. Bitz, E.C. Hunke, W.H. Lipscomb, and M.M. Schramm (2004): Scientific description of the sea ice component in CCSM3. NCAR Tech Note.
• Briegleb, B.P. and B. Light (2007): NCAR/TN-472+STR.

SeaIceAlbedo(ice_thickness, snow_thickness, surface_temperature;
                  ice_albedo = 0.54,
                  snow_albedo = 0.83,
                  ice_melt_reduction = 0.075,
                  snow_melt_reduction = 0.10,
                  melting_temperature = 0.0,
                  temperature_range = 1.0,
                  ocean_albedo = 0.06,
                  minimum_ice_thickness = 0.5,
                  minimum_snow_depth = 0.02)

Construct a CCSM3 sea ice albedo parameterization. Requires references to the sea ice model's ice thickness, snow thickness, and surface temperature fields.

Broadband albedos are approximate averages of the visible and near-IR bands weighted by solar spectrum (52% visible, 48% near-IR):

• ice_albedo ≈ 0.52 x 0.73 + 0.48 x 0.33 ≈ 0.54
• snow_albedo ≈ 0.52 x 0.96 + 0.48 x 0.68 ≈ 0.83

SurfaceRadiationProperties(albedo, emissivity)

Bundle the radiative properties of a single surface (ocean, sea ice, snow, land) that participate in radiative flux computation: shortwave reflectivity (albedo) and longwave emissivity (emissivity).

albedo may be a Number, a LatitudeDependentAlbedo, a TabulatedAlbedo, or any other object for which stateindex is defined. emissivity may be a Number or any other stateindex-able object.

TabulatedAlbedo(arch = CPU(), FT = Float64;
                S₀ = convert(FT, 1365),
                α_table  = α_payne,
                φ_values = (0:2:90) ./ 180 * π,
                𝓉_values = 0:0.05:1)

Constructs a TabulatedAlbedo object that interpolated the albedo from a value table α_table that is function of latitude φ and atmospheric transimissivity 𝓉.

Note: TabulatedAlbedo assumes that the latitude and the transissivity in the table are uniformly spaced.

The transmissivity of the atmosphere is calculated as the ratio of the downwelling solar radiation to the maximum possible downwelling solar radiation for a transparent atmosphere, function of hour of the day, latitude, and day in the year.

Arguments

• arch: The architecture to use. Default: CPU().
• FT: The floating-point type to use. Default: Float64.

Keyword Arguments

• S₀: The solar constant. Default: convert(FT, 1365).
• α_table: The table of albedo values. Default: α_payne.
• φ_values: The latitude values for the table. Default: (0:2:90) ./ 180 * π.
• 𝓉_values: The transmissivity values for the table. Default: 0:0.05:1.

FreezingLimitedOceanTemperature(FT=Float64; liquidus=LinearLiquidus(FT))

The minimal possible sea ice representation, clipping the temperature below to the freezing point. Not really a "model" per se, however, it is the most simple way to make sure that temperature does not dip below freezing.

The melting temperature is a function of salinity and is controlled by the liquidus.

sea_ice_simulation(grid, ocean=nothing;
                   clock = Clock(grid),
                   stop_time = default_stop_time(grid, clock),
                   Δt = 5minutes,
                   ice_salinity = 4, # psu
                   advection = nothing,
                   tracers = (),
                   ice_heat_capacity = 2100, # J kg⁻¹ K⁻¹
                   ice_consolidation_thickness = 0.05, # m
                   sea_ice_density = 900, # kg m⁻³
                   snow_density = 330, # kg m⁻³
                   dynamics = sea_ice_dynamics(grid, ocean),
                   bottom_heat_boundary_condition = nothing,
                   top_heat_boundary_condition = nothing,
                   timestepper = :SplitRungeKutta3,
                   phase_transitions = PhaseTransitions(eltype(grid);
                                                        heat_capacity=ice_heat_capacity,
                                                        density=sea_ice_density),
                   conductivity = 2, # W m⁻¹ K⁻¹
                   internal_heat_flux = ConductiveFlux(; conductivity),
                   snow_thermodynamics = default_snow_thermodynamics(grid))

Construct a sea ice simulation with the given grid and optional ocean simulation. The sea ice model is configured with a slab thermodynamics, Elasto-Visco-Plastic rheology, and a SplitExplicit Runge-Kutta 3rd order time stepper by default. The thermodynamics include conductive internal heat flux, and the option to specify top and bottom heat boundary conditions. The dynamics include a semi-implicit ocean stress formulation, with the option to specify a free drift velocity.

Arguments

• grid: the grid on which to build the sea ice model
• ocean: optional ocean simulation to provide surface velocities and salinity for the sea ice

Keyword Arguments

• clock: Clock for the underlying model. Defaults to Clock(grid), a numeric clock starting at time = 0. Pass a DateTime-based clock to step the simulation in calendar time (e.g. when coupling).
• stop_time: Stop time for the simulation. Defaults to Inf for numeric clocks, or DateTime(9999, 12, 31, 23, 59, 59) for DateTime clocks. On Reactant architectures it defaults to nothing, since Reactant does not support stop_time.
• Δt: time step for the sea ice simulation
• ice_salinity: salinity of the sea ice (psu)
• advection: optional advection scheme for the sea ice model; if nothing (default), no advection is applied and only thermodynamics evolve the sea ice state
• tracers: optional tracers to include in the sea ice model
• ice_heat_capacity: heat capacity of the sea ice (J kg⁻¹ K⁻¹)
• ice_consolidation_thickness: thickness threshold for sea ice consolidation (m)
• sea_ice_density: density of the sea ice (kg m⁻³)
• snow_density: density of the snow (kg m⁻³)
• dynamics: sea ice dynamics model to use (default is sea_ice_dynamics(grid, ocean))
• bottom_heat_boundary_condition: heat boundary condition at the ice-ocean interface (default is IceWaterThermalEquilibrium with ocean surface salinity)
• top_heat_boundary_condition: heat boundary condition at the ice-atmosphere interface (default is a prescribed temperature calculated in the flux computation)
• timestepper: time stepper to use for the sea ice model (default is :SplitRungeKutta3)
• phase_transitions: phase transition properties for the sea ice (default is a PhaseTransitions with specified heat capacity and density)
• conductivity: thermal conductivity for the internal heat flux (W m⁻¹ K⁻¹)
• internal_heat_flux: internal heat flux formulation for the sea ice (default is a ConductiveFlux with specified conductivity)
• snow_thermodynamics: thermodynamics for the snow layer (default is a slab thermodynamics with specified conductivity and prescribed temperature)