Public Documentation

<!– Title –> <h1 align="center"> ClimaOcean.jl </h1>

<!– description –> <p align="center"> <strong>🌎 A framework for realistic ocean-only and coupled ocean + sea-ice simulations driven by prescribed atmospheres and based on <a href=https://github.com/CliMA/Oceananigans.jl>Oceananigans</a> and <a href=https://github.com/CliMA/ClimaSeaIce.jl>ClimaSeaIce</a></strong>. </p>

Installation

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

julia> using Pkg

julia> Pkg.add("ClimaOcean")

julia> Pkg.instantiate()

Use Pkg.add(url="https://github.com/CliMA/ClimaOcean.jl.git", rev="main") to install the latest version of ClimaOcean. For more information, see the documentation for Pkg.jl.

Why? What's the difference between ClimaOcean and Oceananigans?

Oceananigans is a general-purpose library for ocean-flavored fluid dynamics. ClimaOcean implements a framework for driving realistic Oceananigans simulations with prescribed atmospheres, and coupling them to prognostic sea ice simulations.

A core abstraction: ClimaOcean.OceanSeaIceModel

Our system for realistic modeling is anchored by ClimaOcean.OceanSeaIceModel, which encapsulates the ocean simulation, sea ice simulation, prescribed atmospheric state, and specifies how the three communicate. To illustrate how OceanSeaIceModel works we set up a simulation on a grid with 10 vertical levels and 1/4-degree horizontal resolution:

using Oceananigans
using Oceananigans.Units
using Dates, CFTime
import ClimaOcean

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

bathymetry = ClimaOcean.regrid_bathymetry(grid) # builds gridded bathymetry based on ETOPO1
grid = ImmersedBoundaryGrid(grid, GridFittedBottom(bathymetry))

# Build an ocean simulation initialized to the ECCO state estimate on Jan 1, 1993
ocean = ClimaOcean.ocean_simulation(grid)
start_date = DateTime(1993, 1, 1)
set!(ocean.model,
     T=ClimaOcean.Metadata(:temperature; dates=start_date, dataset=ClimaOcean.ECCO4Monthly()),
     S=ClimaOcean.Metadata(:salinity;    dates=start_date, dataset=ClimaOcean.ECCO4Monthly()))

# Build and run an OceanSeaIceModel (with no sea ice component) forced by JRA55 reanalysis
atmosphere = ClimaOcean.JRA55PrescribedAtmosphere(arch)
coupled_model = ClimaOcean.OceanSeaIceModel(ocean; atmosphere)
simulation = Simulation(coupled_model, Δt=5minutes, stop_time=30days)
run!(simulation)

The simulation above 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 in our scripts. For example, to plot the surface speed at the end of the simulation we write

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

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

which produces

Additional features: a utility for ocean_simulations and data wrangling

A second core abstraction in ClimaOcean is ocean_simulation. ocean_simulation configures an Oceananigans model for realistic simulations including temperature and salinity, the TEOS-10 equation of state, boundary conditions to store computed air-sea fluxes, the automatically-calibrated turbulence closure CATKEVerticalDiffusivity, and the WENOVectorInvariant advection scheme for mesoscale-turbulence-resolving simulations.

ClimaOcean also provides convenience features for wrangling datasets of bathymetry, ocean temperature, salinity, ocean velocity fields, and prescribed atmospheric states.

ClimaOcean is built on top of Oceananigans and ClimaSeaIce, so it's important that ClimaOcean users become proficient with Oceananigans. Note that though ClimaOcean is currently focused on hydrostatic modeling with Oceananigans.HydrostaticFreeSurfaceModel, realistic nonhydrostatic modeling is also within the scope of this package.

Citing

If you use ClimaOcean 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. (2025). CliMA/ClimaOcean.jl: v0.5.4 (v0.5.4). Zenodo. https://doi.org/10.5281/zenodo.15042648

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, John Marshall, Raffaele Ferrari

submitted to the Journal of Advances in Modeling Earth Systems, arXiv 2502.14148

<details><summary>bibtex</summary> <pre><code>@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 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}, }</code></pre> </details>

MixedLayerDepthField(bm, grid, tracers; difference_criterion=1e-4)

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.

• 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),
         bounding_box = 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 ECCO2Monthly(), ECCO2Daily(), ECCO4Monthly(), EN4Monthly(), RepeatYearJRA55(), and MultiYearJRA55().

• 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. 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). start_date should lie within the date range of the dataset. Default: nothing.

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

• bounding_box: Specifies the bounds of the dataset. See BoundingBox.

• dir: The directory where the dataset is stored.

Metadatum(variable_name;
          dataset,
          bounding_box = nothing,
          date = first_date(dataset, variable_name),
          dir = default_download_directory(dataset))

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

all_dates(metadata)

Extract all the dates of the given metadata formatted using the DateTime type. Note: This methods needs to be extended for any new dataset.

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.

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

An alias to construct a Metadatum of ECCO4Monthly().

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

An alias to construct a Metadatum of EN4Monthly.

JRA55FieldTimeSeries(variable_name, architecture=CPU(), FT=Float32;
                     dataset = RepeatYearJRA55(),
                     dates = all_JRA55_dates(version),
                     latitude = nothing,
                     longitude = nothing,
                     dir = download_JRA55_cache,
                     backend = InMemory(),
                     time_indexing = Cyclical())

Return a FieldTimeSeries containing atmospheric reanalysis data for variable_name, which describes one of the variables from the Japanese 55-year atmospheric reanalysis for driving ocean-sea ice models (JRA55-do). The JRA55-do dataset is described by:

Tsujino et al. (2018). JRA-55 based surface dataset for driving ocean-sea-ice models (JRA55-do), Ocean Modelling, 130(1), 79-139

The variable_names (and their shortnames used in the netCDF files) available from the JRA55-do are:

• :river_freshwater_flux ("friver")
• :rain_freshwater_flux ("prra")
• :snow_freshwater_flux ("prsn")
• :iceberg_freshwater_flux ("licalvf")
• :specific_humidity ("huss")
• :sea_level_pressure ("psl")
• :relative_humidity ("rhuss")
• :downwelling_longwave_radiation ("rlds")
• :downwelling_shortwave_radiation ("rsds")
• :temperature ("ras")
• :eastward_velocity ("uas")
• :northward_velocity ("vas")

Keyword arguments

• architecture: Architecture for the FieldTimeSeries. Default: CPU()

• dataset: The data dataset; supported datasets are: RepeatYearJRA55() and MultiYearJRA55(). MultiYearJRA55() refers to the full length of the JRA55-do dataset; RepeatYearJRA55() refers to the "repeat-year forcing" dataset derived from JRA55-do. Default: RepeatYearJRA55().

  Repeat-year forcing

  For more information about the derivation of the repeat-year forcing dataset, see:

  Stewart et al. (2020). JRA55-do-based repeat year forcing datasets for driving ocean–sea-ice models, Ocean Modelling, 147, 101557, https://doi.org/10.1016/j.ocemod.2019.101557.

  The repeat year in RepeatYearJRA55() corresponds to May 1st, 1990 - April 30th, 1991. However, the returned dataset has dates that range from January 1st to December 31st. This implies that the first 4 months of the JRA55RepeatYear() dataset correspond to year 1991 from the JRA55 reanalysis and the rest 8 months from 1990.

• start_date: The starting date to use for the dataset. Default: first_date(dataset, variable_name).

• end_date: The ending date to use for the dataset. Default: end_date(dataset, variable_name).

• dir: The directory of the data file. Default: ClimaOcean.JRA55.download_JRA55_cache.

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

• latitude: Guiding latitude bounds for the resulting grid. Used to slice the data when loading into memory. Default: nothing, which retains the latitude range of the native grid.

• longitude: Guiding longitude bounds for the resulting grid. Used to slice the data when loading into memory. Default: nothing, which retains the longitude range of the native grid.

• backend: Backend for the FieldTimeSeries. The two options are:

  • InMemory(): the whole time series is loaded into memory.
  • JRA55NetCDFBackend(total_time_instances_in_memory): only a subset of the time series

  is loaded into memory. Default: InMemory().

JRA55PrescribedAtmosphere([architecture = CPU(), FT = Float32];
                          dataset = RepeatYearJRA55(),
                          start_date = first_date(dataset, :temperature),
                          end_date = last_date(dataset, :temperature),
                          backend = JRA55NetCDFBackend(10),
                          time_indexing = Cyclical(),
                          surface_layer_height = 10,  # meters
                          include_rivers_and_icebergs = false,
                          other_kw...)

Return a PrescribedAtmosphere representing JRA55 reanalysis data. The atmospheric data will be held in JRA55FieldTimeSeries objects containing. For a detailed description of the keyword arguments, see the JRA55FieldTimeSeries constructor.

download_bathymetry(; dir = download_bathymetry_cache,
                      url = etopo_url,
                      filename = "ETOPO_2022_v1_60s_N90W180_surface.nc")

Download the bathymetry from url and saves it under filename in the directory dir and return the full filepath where the bathymetry is saved.

regrid_bathymetry(target_grid;
                  height_above_water = nothing,
                  minimum_depth = 0,
                  dir = download_bathymetry_cache,
                  url = ClimaOcean.Bathymetry.etopo_url,
                  filename = "ETOPO_2022_v1_60s_N90W180_surface.nc",
                  interpolation_passes = 1,
                  major_basins = 1)

Return bathymetry associated with the NetCDF file at path = joinpath(dir, filename) regridded onto target_grid. If path does not exist, then a download is attempted from joinpath(url, filename).

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.

• dir: directory of the bathymetry-containing file. Default: download_bathymetry_cache.

• filename: file containing bathymetric data. Must be netCDF with fields:

  1. lat vector of latitude nodes
  2. lon vector of longitude nodes
  3. z matrix of depth values

• 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.

retrieve_bathymetry(grid, filename; kw...)

Retrieve the bathymetry data from a file or generate it using a grid and save it to a file.

Arguments

• grid: The grid used to generate the bathymetry data.
• filename: The name of the file to read or save the bathymetry data.
• kw...: Additional keyword arguments.

Returns

• bottom_height: The retrieved or generated bathymetry data.

If the specified file exists, the function reads the bathymetry data from the file. Otherwise, it generates the bathymetry data using the provided grid and saves it to the file before returning it.

stretched_vertical_faces(; surface_layer_Δz = 5.0,
                           surface_layer_height = 100.0,
                           constant_bottom_spacing_depth = Inf,
                           maximum_Δz = Inf,
                           stretching = PowerLawStretching(1.02),
                           rounding_digits = 1,
                           depth = 5000)

Return an array of cell interfaces with surface_layer_Δz spacing in a surface layer of height surface_layer_height, and stretched according to the function stretching(Δz_above, z_above) down to depth. The interfaces extends from Lz = -z[1] to 0 = z[end], where Lz ≥ depth.

The grid spacing Δz is limited to be less than maximum_Δz. The grid is also uniformly-spaced below constant_bottom_spacing_depth.

rounding_digits controls the accuracy with which the grid face positions are saved.

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.

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.

ComponentInterfaces(atmosphere, ocean, sea_ice=nothing;
                    radiation = Radiation(),
                    freshwater_density = 1000,
                    atmosphere_ocean_flux_formulation = SimilarityTheoryFluxes(),
                    atmosphere_sea_ice_flux_formulation = SimilarityTheoryFluxes(eltype(ocean.model.grid)),
                    atmosphere_ocean_interface_temperature = BulkTemperature(),
                    atmosphere_ocean_interface_specific_humidity = default_ao_specific_humidity(ocean),
                    atmosphere_sea_ice_interface_temperature = default_ai_temperature(sea_ice),
                    ocean_reference_density = reference_density(ocean),
                    ocean_heat_capacity = heat_capacity(ocean),
                    ocean_temperature_units = DegreesCelsius(),
                    sea_ice_temperature_units = DegreesCelsius(),
                    sea_ice_reference_density = reference_density(sea_ice),
                    sea_ice_heat_capacity = heat_capacity(sea_ice),
                    gravitational_acceleration = g_Earth)

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.

Radiation([arch = CPU(), FT=Float64];
          ocean_emissivity = 0.97,
          sea_ice_emissivity = 1.0,
          ocean_albedo = 0.05,
          sea_ice_albedo = 0.7,
          stefan_boltzmann_constant = 5.67e-8)

Constructs a Radiation object that represents the radiation properties of the ocean and sea ice.

Arguments

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

Keyword Arguments

• ocean_emissivity: The emissivity of the ocean surface. Default: 0.97.
• sea_ice_emissivity: The emissivity of the sea ice surface. Default: 1.0.
• ocean_albedo: The albedo of the ocean surface. Default: 0.05.
• sea_ice_albedo: The albedo of the sea ice surface. Default: 0.7.
• stefan_boltzmann_constant: The Stefan-Boltzmann constant. Default: 5.67e-8.

SimilarityTheoryFluxes(FT::DataType = Float64;
                       gravitational_acceleration = 9.81,
                       von_karman_constant = 0.4,
                       turbulent_prandtl_number = 1,
                       gustiness_parameter = 1,
                       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: Increases surface fluxes in low wind conditions. Default: 1.
• 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.

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.