Public API
NumericalEarth.NumericalEarth — Module
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,
Download Julia (version 1.10 or later).
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},
}Atmospheres
NumericalEarth.Atmospheres.PrescribedAtmosphere — Type
PrescribedAtmosphere(grid, times=[zero(grid)];
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),
tracers = default_atmosphere_tracers(grid, times),
pressure = default_atmosphere_pressure(grid, times),
freshwater_flux = default_freshwater_flux(grid, times))Return a prescribed, time-evolving atmospheric state with data on grid and at given times.
!!! compat Radiation component The downwelling shortwave / longwave radiation part of the top-level radiation component (see PrescribedRadiation, JRA55PrescribedRadiation).
NumericalEarth.Atmospheres.PrescribedPrecipitationFlux — Type
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 freshwater_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.
Bathymetry
NumericalEarth.Bathymetry.ORCAGrid — Function
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., ORCA1() or ORCA12()). 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.
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()orGPU()). Default:CPU().FT: Floating point type. Default:Float64.
Keyword Arguments
dataset: The ORCA dataset to use. Default:ORCA1().ORCA12()is also supported (ORCA1 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., anExponentialDiscretization. Default:(-6000, 0).Nz: Number of vertical levels (only used whenzis a 2-tuple). Default:50.radius: Planet radius. Default:Oceananigans.defaults.planet_radius.with_bathymetry: Iftrue, download the bathymetry and return anImmersedBoundaryGridwithGridFittedBottom. Default:true.active_cells_map: Iftrueandwith_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 viaremove_minor_basins!. Basins are removed from smallest to largest;major_basins = 1keeps 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_maskand bathymetry). Defaults to the dataset scratch cache viadefault_download_directory(dataset).
NumericalEarth.Bathymetry.regrid_bathymetry — Method
regrid_bathymetry(target_grid, metadata;
height_above_water = nothing,
minimum_depth = 0,
major_basins = 1,
interpolation_passes = 1,
cache = true)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_depthis considered land. Default: 0.interpolation_passes: regridding/interpolation passes. The bathymetry is interpolated ininterpolation_passes - 1intermediate 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 byremove_minor_basins!. Basins are removed by order of size: the smallest basins are removed first.major_basins = 1retains only the largest basin. IfInfthen no basins are removed. Default: 1.cache: Iftrue(default), caches the regridded bathymetry to disk and reuses it on subsequent calls with the same grid and parameters. Set tofalseto force recomputation.
NumericalEarth.Bathymetry.regrid_bathymetry — Method
regrid_bathymetry(target_grid; dataset=ETOPO2022(), cache=true, kw...)Regrid bathymetry from dataset onto target_grid. Default: dataset = ETOPO2022().
DataWrangling
NumericalEarth.DataWrangling — Module
Incorporate various datasets to be used for bathymetry, initialization, forcing, restoring, or validation.
NumericalEarth.DataWrangling.BoundingBox — Method
BoundingBox(; longitude=nothing, latitude=nothing, z=nothing)Create a bounding box with latitude, longitude, and z bounds on the sphere.
NumericalEarth.DataWrangling.Column — Type
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.
NumericalEarth.DataWrangling.DatasetRestoring — Type
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 anarchitecture is provided, such asarch_or_grid = CPU()orarch_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, seeinpaint_mask!. Default:NearestNeighborInpainting(Inf).cache_inpainted_data: Iftrue, the data is cached to disk after inpainting for later retrieving. Default:true.
NumericalEarth.DataWrangling.LinearlyTaperedPolarMask — Method
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) : 0NumericalEarth.DataWrangling.Metadata — Method
Metadata(variable_name;
dataset,
dates = all_dates(dataset, variable_name),
dir = default_download_directory(dataset),
region = nothing,
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 areETOPO2022(),ECCO2Monthly(),ECCO2Daily(),ECCO4Monthly(),EN4Monthly(),GLORYSDaily(),GLORYSMonthly(),RepeatYearJRA55(), andMultiYearJRA55().dates: The dates of the dataset (Dates.AbstractDateTimeorCFTime.AbstractCFDateTime). Note thatdatescan either be a range or a vector of dates, representing a time-series. For a single date, useMetadatum.start_date: Ifdates = nothing, we can prescribe the first date of metadata as a date (Dates.AbstractDateTimeorCFTime.AbstractCFDateTime). If outside the date range of the dataset, the first allowable date is chosen. Default: nothing.end_date: Ifdates = nothing, we can prescribe the last date of metadata as a date (Dates.AbstractDateTimeorCFTime.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 aBoundingBoxfor a rectangular region, aColumnfor a single horizontal location, ornothingfor the full domain.filename: The filename(s) for the dataset. Ifnothing, the filename is computed from the dataset type. Can be aString(single file for all dates) or aDatewiseFilename(one file per date).dir: The directory where the dataset is stored.
NumericalEarth.DataWrangling.MetadataSet — Method
MetadataSet(variable_names::Symbol...;
dataset,
dates = all_dates(dataset, first(variable_names)),
date = nothing,
dir = default_download_directory(dataset),
region = nothing,
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 moreSymbols 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 (e.g.ECCO4Monthly(),ERA5HourlyPressureLevels()).dates: shared date axis. Either a singleAbstractDateTime/AbstractCFDateTimeor anAbstractVectorof dates. Defaults toall_dates(dataset, first(variable_names)).date: convenience scalar form; cannot be used together withdates.region: shared spatial region —BoundingBox,Column, ornothing.dir: shared download directory.filenames: an optionalNamedTuplekeyed byvariable_namesoverriding the auto-computed per-variable filenames.start_date,end_date: optional date cropping, matchingMetadata.
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/v4NumericalEarth.DataWrangling.Metadatum — Method
Metadatum(variable_name;
dataset,
region = nothing,
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.
NumericalEarth.DataWrangling.SurfaceFluxRestoring — Type
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)))NumericalEarth.DataWrangling.all_dates — Method
all_dates(metadata)Extract all the dates of the given metadata formatted using the DateTime type.
NumericalEarth.DataWrangling.first_date — Method
first_date(dataset, variable_name)Extracts the first date of the given dataset and variable name formatted using the DateTime type.
NumericalEarth.DataWrangling.last_date — Method
last_date(dataset, variable_name)Extracts the last date of the given dataset and variable name formatted using the DateTime type.
NumericalEarth.DataWrangling.metadata_filename — Function
metadata_filename(dataset, name, date, region)Compute the filename for a single date. Extended by each dataset module.
NumericalEarth.DataWrangling.metadata_filename — Method
metadata_filename(metadata)Return the stored filename(s) of metadata.
NumericalEarth.DataWrangling.native_grid — Function
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.
NumericalEarth.DataWrangling.validate_dataset_coverage — Method
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.
DataWrangling.ECCO
NumericalEarth.DataWrangling.ECCO.ECCOMetadatum — Method
ECCOMetadatum(name;
date = first_date(ECCO4Monthly(), name),
dir = download_ECCO_cache)An alias to construct a Metadatum of ECCO4Monthly().
NumericalEarth.DataWrangling.ECCO.ECCOPrescribedAtmosphere — Function
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.
NumericalEarth.DataWrangling.ECCO.ECCOPrescribedRadiation — Function
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, 0.97),
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.
NumericalEarth.DataWrangling.retrieve_data — Method
retrieve_data(metadata::Metadatum{<:ECCO4DarwinMonthly})Read a ECCO4DarwinMonthly data file and regrid using MeshArrays on to regular lat-lon grid
DataWrangling.EN4
NumericalEarth.DataWrangling.EN4.EN4Metadatum — Method
EN4Metadatum(name;
date = first_date(EN4Monthly()),
dir = download_EN4_cache)An alias to construct a Metadatum of EN4Monthly.
DataWrangling.ERA5
NumericalEarth.DataWrangling.ERA5.ERA5HourlyPressureLevels — Method
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 instantaneous geopotential Φ(λ,φ,p)/g with sub-surface levels clipped at the surface — the right thing over terrain and synoptic Φ-displacement (issue #236).
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).
NumericalEarth.DataWrangling.ERA5.pressure_field — Function
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.
DataWrangling.ETOPO
DataWrangling.GEBCO
NumericalEarth.DataWrangling.GEBCO.GEBCO2024 — Type
GEBCO2024General Bathymetric Chart of the Oceans 2024 release. Global bathymetry and topography at 15 arc-second resolution.
The GEBCO_2024 Grid is a global terrain model for ocean and land, providing elevation data on a 15 arc-second interval grid.
Reference: GEBCO Compilation Group (2024) GEBCO 2024 Grid Data source: https://www.gebco.net/dataandproducts/griddedbathymetrydata/
DataWrangling.GLORYS
DataWrangling.IBCAO
NumericalEarth.DataWrangling.IBCAO.IBCAOv5 — Type
IBCAOv5International 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/
DataWrangling.IBCSO
NumericalEarth.DataWrangling.IBCSO.IBCSOv2 — Type
IBCSOv2International 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/
DataWrangling.JRA55
NumericalEarth.DataWrangling.JRA55.JRA55PrescribedAtmosphere — Function
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 parameterised 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.
NumericalEarth.DataWrangling.JRA55.JRA55PrescribedLand — Function
JRA55PrescribedLand([architecture = CPU()];
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,
other_kw...)Return a PrescribedLand representing JRA55 reanalysis land surface data (river runoff and iceberg calving freshwater fluxes).
NumericalEarth.DataWrangling.JRA55.JRA55PrescribedRadiation — Function
JRA55PrescribedRadiation([architecture = CPU()];
dataset = RepeatYearJRA55(),
start_date = first_date(dataset, :downwelling_shortwave_radiation),
end_date = last_date(dataset, :downwelling_shortwave_radiation),
dir = download_JRA55_cache,
time_indices_in_memory = 10,
time_indexing = Cyclical(),
ocean_surface = SurfaceRadiationProperties(0.05, 0.97),
sea_ice_surface = SurfaceRadiationProperties(0.7, 1.0),
stefan_boltzmann_constant = default_stefan_boltzmann_constant,
region = nothing,
other_kw...)Return a PrescribedRadiation backed by JRA55 downwelling shortwave and longwave FieldTimeSeries. Surface radiative properties (albedo, emissivity) for ocean and sea-ice surfaces default to standard values; pass ocean_surface = nothing (or sea_ice_surface = nothing) to omit a surface.
DataWrangling.ORCA
DataWrangling.OSPapa
NumericalEarth.DataWrangling.OSPapa.OSPapaPrescribedAtmosphere — Function
OSPapaPrescribedAtmosphere(architecture = CPU(), FT = Float32;
start_date = first_date(OSPapaHourly(), :air_temperature),
end_date = last_date(OSPapaHourly(), :air_temperature),
dir = download_OSPapa_cache,
surface_layer_height = 2.5,
max_gap_hours = 72)Construct a PrescribedAtmosphere from Ocean Station Papa buoy observations.
Data is automatically downloaded from the NOAA/PMEL AWS S3 bucket if not already cached locally.
The buoy SW and LW variables are downwelling fluxes. When this atmosphere is used with OceanOnlyModel, ClimaOcean applies its own ocean albedo (default α = 0.05) to compute net absorbed shortwave, and computes upwelling longwave from the model SST via Stefan-Boltzmann. This means the resulting net heat flux will differ from the COARE-computed QNET available via os_papa_prescribed_fluxes. If you need the exact observed net fluxes, use os_papa_prescribed_flux_boundary_conditions instead.
Keyword Arguments
start_date: start of the time rangeend_date: end of the time rangedir: directory for cached data filessurface_layer_height: measurement height in meters (default: 2.5, matching the buoy's temperature/humidity instruments)max_gap_hours: maximum gap size (in hours) to fill by linear interpolation (default: 72)
NumericalEarth.DataWrangling.OSPapa.OSPapaPrescribedRadiation — Function
OSPapaPrescribedRadiation(architecture = CPU(), FT = Float32;
start_date = first_date(OSPapaHourly(), :shortwave_radiation),
end_date = last_date(OSPapaHourly(), :shortwave_radiation),
dir = download_OSPapa_cache,
max_gap_hours = 72,
ocean_surface = SurfaceRadiationProperties(0.05, 0.97),
sea_ice_surface = nothing,
stefan_boltzmann_constant = default_stefan_boltzmann_constant)Construct a PrescribedRadiation from Ocean Station Papa buoy SW/LW observations on a single-column grid.
NumericalEarth.DataWrangling.OSPapa.os_papa_prescribed_flux_boundary_conditions — Method
os_papa_prescribed_flux_boundary_conditions(fluxes;
arch=nothing,
ρ₀=1020.0,
cₚ=3991.0,
u_momentum_flux_correction=no_correction,
v_momentum_flux_correction=no_correction,
temperature_flux_correction=no_correction,
salinity_flux_correction=no_correction)Create Oceananigans FluxBoundaryConditions for u, v, T, S from prescribed OS Papa flux data. Returns a NamedTuple of FieldBoundaryConditions that can be passed directly to ocean_simulation or HydrostaticFreeSurfaceModel via the boundary_conditions keyword argument.
If arch is provided, the flux FieldTimeSeries are first moved to that architecture.
Arguments
fluxes: aNamedTuplereturned byos_papa_prescribed_fluxes
Keyword Arguments
arch: target architecture for the fluxFieldTimeSeries. If provided, each flux FTS is moved toarchviaon_architecture. Defaults tonothing, which keeps them on whatever architecture they were built on.ρ₀: reference ocean density (default: 1020 kg/m³)cₚ: ocean heat capacity (default: 3991 J/(kg·K))u_momentum_flux_correction: discrete-form correction function added to the zonal momentum flux boundary condition (default:no_correction)v_momentum_flux_correction: discrete-form correction function added to the meridional momentum flux boundary condition (default:no_correction)temperature_flux_correction: discrete-form correction function added to the temperature flux boundary condition (default:no_correction)salinity_flux_correction: discrete-form correction function added after converting freshwater flux fromEMPinto the salinity flux boundary condition (default:no_correction)
Each correction function must have the signature (i, j, grid, clock, model_fields, p) and return a value in the same units as the corresponding boundary condition.
Examples
# Basic usage on GPU:
fluxes = os_papa_prescribed_fluxes(GPU(); start_date, end_date)
bcs = os_papa_prescribed_flux_boundary_conditions(fluxes)
ocean = ocean_simulation(grid; Δt=10minutes, boundary_conditions=bcs)
# With a uniform heat flux correction of +5 W/m² to close the heat budget:
heat_correction = (i, j, grid, clock, model_fields, p) -> 5.0 / (p.ρ₀ * p.cₚ)
bcs = os_papa_prescribed_flux_boundary_conditions(fluxes; temperature_flux_correction=heat_correction)NumericalEarth.DataWrangling.OSPapa.os_papa_prescribed_fluxes — Function
os_papa_prescribed_fluxes(architecture = CPU(), FT = Float64;
start_date = first_date(OSPapaFluxHourly(), :net_heat_flux),
end_date = last_date(OSPapaFluxHourly(), :net_heat_flux),
dir = download_OSPapa_cache,
max_gap_hours = 72)Download precomputed air-sea fluxes for Ocean Station Papa from the PMEL ERDDAP ocs_papa_flux dataset (computed with COARE 3.0b) and return them as a NamedTuple of FieldTimeSeries on a Flat, Flat, Flat grid at architecture.
The returned tuple contains:
Qnet: net heat flux (W/m², positive into ocean)Qlat: latent heat flux (W/m², positive upward = ocean cooling)Qsen: sensible heat flux (W/m², positive upward = ocean cooling)SWnet: net shortwave radiation (W/m², positive downward)LWnet: net longwave radiation (W/m², positive upward = ocean cooling)τx,τy: zonal and meridional wind stress (N/m²)evap,rain,EMP: evaporation, precipitation, E-P (mm/hr)Tsk: skin temperature (°C)
Keyword Arguments
start_date: start of the time rangeend_date: end of the time rangedir: directory for cached data filesmax_gap_hours: maximum gap size (in hours) to fill by linear interpolation (default: 72)
DataWrangling.WOA
Diagnostics
NumericalEarth.Diagnostics.MixedLayerDepthField — Method
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.
NumericalEarth.Diagnostics.atmosphere_ocean_freshwater_flux — Method
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.
NumericalEarth.Diagnostics.atmosphere_ocean_heat_flux — Method
atmosphere_ocean_heat_flux(esm::EarthSystemModel)Return the atmosphere-ocean heat flux (W m⁻²) at the atmosphere-ocean interface in a coupled esm.
NumericalEarth.Diagnostics.meridional_heat_transport — Method
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.
The meridional_heat_transport diagnostic currently is only supported only on LongitudeLatitudeGrids.
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 haloNumericalEarth.Diagnostics.net_ocean_freshwater_flux — Method
net_ocean_freshwater_flux(esm::EarthSystemModel)Return the net freshwater mass flux (kg m⁻² s⁻¹) at the ocean's surface in a coupled esm.
NumericalEarth.Diagnostics.net_ocean_heat_flux — Method
net_ocean_heat_flux(esm::EarthSystemModel)Return the net heat flux (W m⁻²) at the ocean's surface in a coupled esm.
NumericalEarth.Diagnostics.sea_ice_ocean_freshwater_flux — Method
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.
NumericalEarth.Diagnostics.sea_ice_ocean_heat_flux — Method
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.
EarthSystemModels
NumericalEarth.EarthSystemModels.EarthSystemModel — Method
EarthSystemModel(radiation, atmosphere, land, sea_ice, ocean;
clock = Clock{Float64}(time=0),
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, ornothingfor a radiatively decoupled surface. Pass aPrescribedRadiation(e.g.JRA55PrescribedRadiation(...)) to enable radiative forcing.atmosphere: A representation of a possibly time-dependent atmospheric state, ornothing. For a prognostic atmosphere, useatmosphere_simulation. For prescribed atmospheric forcing, useJRA55PrescribedAtmosphereorPrescribedAtmosphere.land: Land component, ornothing.sea_ice: A representation of a possibly time-dependent sea ice state, ornothing. For example, the minimalistFreezingLimitedOceanTemperaturerepresents oceanic latent heating during freezing only, but does not evolve sea ice variables. For prognostic sea ice use anOceananigans.SimulationofClimaSeaIce.SeaIceModel.ocean: A representation of a possibly time-dependent ocean state. Currently, onlyOceananigans.Simulations ofOceananigans.HydrostaticFreeSurfaceModelare tested.
Keyword Arguments
clock: Keeps track of time.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 tonothingand will be constructed automatically. To customize the sea ice-ocean heat flux formulation, create interfaces manually usingComponentInterfaces.
NumericalEarth.EarthSystemModels.EarthSystemModel — Method
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 coupledAll keyword arguments accepted by the positional constructor are forwarded.
NumericalEarth.EarthSystemModels.AtmosphereOceanModel — Method
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.
NumericalEarth.EarthSystemModels.OceanOnlyModel — Method
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(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{CPU}(time = 0 seconds, iteration = 0)
├── radiation: Nothing
├── atmosphere: Nothing
├── land: Nothing
├── sea_ice: FreezingLimitedOceanTemperature{ClimaSeaIce.SeaIceThermodynamics.LinearLiquidus{Float64}}
├── ocean: HydrostaticFreeSurfaceModel{CPU, LatitudeLongitudeGrid}(time = 0 seconds, iteration = 0)
└── interfaces: ComponentInterfacesNumericalEarth.EarthSystemModels.OceanSeaIceModel — Method
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
├── ocean: HydrostaticFreeSurfaceModel{CPU, LatitudeLongitudeGrid}(time = 0 seconds, iteration = 0)
└── interfaces: ComponentInterfacesEarthSystemModels.InterfaceComputations
NumericalEarth.EarthSystemModels.InterfaceComputations.BulkTemperature — Type
struct BulkTemperatureA 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.
NumericalEarth.EarthSystemModels.InterfaceComputations.CoefficientBasedFluxes — Type
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 aSimilarityScales, aTuple, or aNamedTuplewith constant or callable entries, or anLargeYeagerTransferCoefficients. Defaults to(1e-3, 1e-3, 1e-3).solver_stop_criteria: Criteria for iterative solver convergence. Ifnothing, creates new criteria usingsolver_toleranceandsolver_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
ComponentInterfacesNumericalEarth.EarthSystemModels.InterfaceComputations.CoefficientBasedFluxes — Type
struct CoefficientBasedFluxes{C, S}A structure for computing turbulent fluxes using bulk transfer coefficients.
transfer_coefficients::Anysolver_stop_criteria::Any
NumericalEarth.EarthSystemModels.InterfaceComputations.ComponentInterfaces — Type
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 pointThreeEquationHeatFlux(): coupled heat/salt/freezing point system (default)
radiation: radiation component. Default:Radiation().radiation: radiation component. Default:nothing.freshwater_density: reference density of freshwater. Default:default_freshwater_density.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. Default:BulkTemperature().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.
NumericalEarth.EarthSystemModels.InterfaceComputations.IceBathHeatFlux — Type
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.
NumericalEarth.EarthSystemModels.InterfaceComputations.IceBathHeatFlux — Type
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 (constantNumberorMomentumBasedFrictionVelocity)
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.002References
- 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.
NumericalEarth.EarthSystemModels.InterfaceComputations.LargeYeagerTransferCoefficients — Type
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
NumericalEarth.EarthSystemModels.InterfaceComputations.LinearStableStabilityFunction — Type
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
NumericalEarth.EarthSystemModels.InterfaceComputations.MomentumBasedFrictionVelocity — Type
MomentumBasedFrictionVelocityA 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)NumericalEarth.EarthSystemModels.InterfaceComputations.MomentumRoughnessLength — Type
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.
NumericalEarth.EarthSystemModels.InterfaceComputations.PolynomialNeutralDragCoefficient — Type
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 otherwiseCalled as a functor: C_dN = drag(U).
References:
- Large, W.G. & Yeager, S.G. (2004): NCAR/TN-460+STR
NumericalEarth.EarthSystemModels.InterfaceComputations.ScalarRoughnessLength — Type
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 to1.6e-4.
NumericalEarth.EarthSystemModels.InterfaceComputations.SimilarityTheoryFluxes — Type
SimilarityTheoryFluxes(FT::DataType = Float64;
gravitational_acceleration = 9.81,
von_karman_constant = 0.4,
turbulent_prandtl_number = 1,
gustiness_parameter = 1.2,
minimum_gustiness = 0.01,
stability_functions = default_stability_functions(FT),
roughness_lengths = default_roughness_lengths(FT),
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.gustiness_parameter: Scaling factor for convective gustiness velocity. Default: 1.2.minimum_gustiness: Minimum gustiness velocity [m/s], used as a floor in stable conditions where convective gustiness is zero. Default: 0.01.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. Default:default_roughness_lengths(FT), formulation taken from Edson et al. (2013).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.
NumericalEarth.EarthSystemModels.InterfaceComputations.SkinTemperature — Type
struct SkinTemperatureA 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.
NumericalEarth.EarthSystemModels.InterfaceComputations.ThreeEquationHeatFlux — Type
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.
NumericalEarth.EarthSystemModels.InterfaceComputations.ThreeEquationHeatFlux — Type
ThreeEquationHeatFlux{FT, U}Three-equation formulation for sea ice-ocean heat flux.
This formulation solves a coupled system for the interface temperature and salinity:
- Heat balance: $\rho c_p \gamma_T (T - T_b) = ℰ q$
- Salt balance: $\gamma_S (S - S_b) = q (S_b - S_i)$
- 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 (constantNumberorMomentumBasedFrictionVelocity)
Example
using NumericalEarth.EarthSystemModels: ThreeEquationHeatFlux
flux = ThreeEquationHeatFlux()
# output
ThreeEquationHeatFlux{Nothing}
├── heat_transfer_coefficient: 0.0095
├── salt_transfer_coefficient: 0.00027142857142857144
└── friction_velocity: 0.002References
- 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.
NumericalEarth.EarthSystemModels.InterfaceComputations.compute_sea_ice_ocean_fluxes! — Method
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.
NumericalEarth.EarthSystemModels.InterfaceComputations.large_yeager_stability_functions — Function
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.
Grids
NumericalEarth.Grids.PressureLevelGrid — Type
PressureLevelGridType alias for any underlying grid whose vertical coordinate is a PressureLevelVerticalDiscretization.
NumericalEarth.Grids.PressureLevelVerticalDiscretization — Type
PressureLevelVerticalDiscretization{G, Geo}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-DVector{Float64}of lengthNz— the representative axis that plot recipes,Lz, and length consumers expect when only the grid is in hand.znodes(field)returns the column-meanVector{Float64}when the field has no horizontal extent (Flattopology orNothinghorizontal locations, e.g. aftermean(...; dims=(1, 2))); otherwise it returns a 3-DFieldof 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.
NumericalEarth.Grids.PressureLevelVerticalDiscretization — Method
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.
InitialConditions
Lands
NumericalEarth.Lands.PrescribedLand — Method
PrescribedLand(freshwater_flux; clock=nothing)Return a PrescribedLand component from a NamedTuple of FieldTimeSeries representing freshwater fluxes (e.g. rivers, icebergs).
If clock is not provided, defaults to a Clock whose time type matches the element type of freshwater_flux.
Oceans
NumericalEarth.Oceans.PrescribedOcean — Type
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:FieldTimeSeriesfor SST.sea_surface_salinity:FieldTimeSeriesfor SSS.velocities:NamedTupleofFieldTimeSeriesfor(u, v).
NumericalEarth.Oceans.SlabOcean — Method
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)whereFTis the float type of the grid.
NumericalEarth.Oceans.ocean_simulation — Method
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 aHydrostaticFreeSurfaceModel-based simulation with a free surface, CATKE vertical mixing, quadratic bottom drag, and a TEOS-10 equation of state. Seehydrostatic_ocean_simulationfor the full kwarg list.:nonhydrostatic— builds aNonhydrostaticModel-based simulation suitable for LES (full 3D pressure, no free surface, no barotropic forcing). Seenonhydrostatic_ocean_simulationfor 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 entriesRadiations
NumericalEarth.Radiations.InterfaceRadiationFlux — Type
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)
NumericalEarth.Radiations.LatitudeDependentAlbedo — Type
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 isFloat64.
Keyword Arguments
diffuse: The diffuse albedo value. Default is0.069.direct: The direct albedo value. Default is0.011.
NumericalEarth.Radiations.PrescribedRadiation — Type
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.
NumericalEarth.Radiations.PrescribedRadiation — Type
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).
NumericalEarth.Radiations.PrescribedRadiation — Method
PrescribedRadiation(downwelling_shortwave, downwelling_longwave;
clock = nothing,
ocean_surface = SurfaceRadiationProperties(0.05, 0.97),
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.
NumericalEarth.Radiations.SeaIceAlbedo — Type
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:
- Base cold albedos: bare ice (0.53) and snow-covered (0.82)
- Temperature reduction within 1C of melting: Δαice = 0.075, Δαsnow = 0.10
- Thin-ice transition to ocean albedo below h_amin = 0.5 m
- 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.
NumericalEarth.Radiations.SeaIceAlbedo — Method
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
NumericalEarth.Radiations.SurfaceRadiationProperties — Method
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.
NumericalEarth.Radiations.TabulatedAlbedo — Type
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.
SeaIces
NumericalEarth.SeaIces.FreezingLimitedOceanTemperature — Type
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.
NumericalEarth.SeaIces.sea_ice_simulation — Function
sea_ice_simulation(grid, ocean=nothing;
Δ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 modelocean: optional ocean simulation to provide surface velocities and salinity for the sea ice
Keyword Arguments
Δt: time step for the sea ice simulationice_salinity: salinity of the sea ice (psu)advection: optional advection scheme for the sea ice model; ifnothing(default), no advection is applied and only thermodynamics evolve the sea ice statetracers: optional tracers to include in the sea ice modelice_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 issea_ice_dynamics(grid, ocean))bottom_heat_boundary_condition: heat boundary condition at the ice-ocean interface (default isIceWaterThermalEquilibriumwith 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 aPhaseTransitionswith 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 aConductiveFluxwith specified conductivity)snow_thermodynamics: thermodynamics for the snow layer (default is a slab thermodynamics with specified conductivity and prescribed temperature)