Standard configurations for CoCiP simulations

Contrails.org maintains an implementation of CoCiP in pycontrails and uses it to produce two products: an Impact Inventory that backs our Impact Explorer and research-grade simulations used in scientific publications.

CoCiP is an evolving model with a relatively large number of user-configurable inputs and parameters. Our Impact Inventory and recent research-grade simulations use different CoCiP versions in different configurations and–even for otherwise-identical flights–often produce slightly different results.

This notebook describes the rationale behind the CoCiP configurations used in the Impact Inventory and for research, shows how to access results using both configurations via the Contrails.org API, demonstrates how to reproduce results using pycontrails, and provides recommendations to help users choose between the two configurations.

Sample flight

The examples in this notebook use a sample flight from the Canary Islands to London Luton. The flight trajectory (longitude, latitude, altitude, and time) was produced from raw ADS-B telemetry processed by our open-source flights-pipeline.

# import dependencies for loading and plotting sample flight
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import pandas as pd

Open the sample trajectory.

df = pd.read_csv("trajectory.csv")

ax = plt.subplot(111, projection=ccrs.PlateCarree())
ax.set_adjustable("datalim")
ax.plot(df["longitude"], df["latitude"], "k-", transform=ccrs.Geodetic())
ax.add_feature(cfeature.LAND, facecolor="whitesmoke")
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor="gray")
ax.add_feature(cfeature.STATES, linewidth=0.2, edgecolor="gray")
ax.set_extent(ax.get_extent())

Impact Inventory

The Impact Inventory is the dataset backing our Impact Explorer. It contains contrail impact modeling for millions of flights and many years. Each of the flights in the inventory was processed from raw ADS-B data through trajectory QA/QC processes, then run through CoCiP to model contrail formation and evolution. The process for ingesting and processing the ADS-B data, cleaning and sanitizing the trajectories, then modeling the flight contrail impacts with CoCiP are in our flights-pipeline.

The design of our flights-pipeline emphasizes scalability and stability. The emphasis on scalability permits efficient queries of inventory results, and the emphasis on stability supports analyses of temporal trends. Users should not expect frequent changes to model configurations. When changes are made, they will be accompanied by re-publication of a new version of the full impact inventory.

We recommend use of this CoCiP configuration for production-grade pipelines that require a foundation in published research but prioritize long-term stability over rapid incorporation of the latest scientific advances.

Contrails API

Results from CoCiP simulations archived in the Impact Inventory can be queried using the v1/inventory/segments endpoint.

Because results returned by v1/inventory endpoints are taken from a precomputed impact inventory, queries are cheap but updates to methodology are expensive, requiring reprocessing of many millions of flights.

# import dependencies for making API requests
import io
import os
import requests

url = "https://api.contrails.org/v1/inventory/segments"
headers = {"x-api-key": os.environ["CONTRAILS_API_KEY"]}
times = pd.date_range(
    pd.to_datetime(df["time"].min()).floor("1h"),
    pd.to_datetime(df["time"].max()).floor("1h"),
    freq="1h",
)
flight_id = df["flight_id"][0]

inventory_list = []

for time in times:
    params = {"date": time.strftime("%Y-%m-%dT%H")}
    resp = requests.get(url, headers=headers, params=params)

    with io.BytesIO() as f:
        f.write(resp.content)
        f.seek(0)
        inv = pd.read_parquet(f)

    inv = inv.loc[inv["flight_id"] == flight_id]
    inventory_list.append(inv)

inventory = pd.concat(inventory_list)

print(f"Contrail EF (impact inventory): {inventory["sum_ef_mj"].sum():.6e} MJ")

Contrail EF (impact inventory): 1.628445e+09 MJ

Equivalent pycontrails code

Impact Inventory results are generated by the flights-pipeline using pycontrails v0.60.3. Individual CoCiP simulations are handled by the CocipTrajectoryHandler class.

import pycontrails

assert pycontrails.__version__ == "0.60.3"

Meteorology

The Impact Inventory uses ERA5 reanalysis meteorological data read from a zarr store. This zarr store contains fields on a quarter-degree grid at hourly temporal resolution, with data interpolated from model levels 66-93 to pressure levels at 134, 141, 148, 155, 163, 171, 180, 188, 197, 207, 217, 227, 237, 248, 260, 272, 284, 297, 310, 323, 337, 352, 367, 383, 399, and 416 hPa. (Sufficiently high vertical resolution is important to avoid missing relatively thin contrail-forming regions.)

Our internal zarr store is not publicly accessible, but the fields in it can be approximately reproduced (with identical resolution) using the ERA5ModelLevel pycontrails datalib.

import xarray as xr

from pycontrails.core import MetDataset
from pycontrails.models.cocip import Cocip

dates = pd.date_range(
    (pd.Timestamp(df["time"].min()) - pd.Timedelta(hours=1)).floor("1D"),
    (pd.Timestamp(df["time"].min()) + pd.Timedelta(hours=13)).floor("1D"),
    freq="1D",
)

pl_list = []
sl_list = []
for date in dates:
    ds = xr.open_zarr(f"gs://contrails-301217-ecmwf-era5-zarr-v2/{date:%Y%m%d}_pl.zarr")
    pl_list.append(ds)

    ds = xr.open_zarr(f"gs://contrails-301217-ecmwf-era5-zarr-v2/{date:%Y%m%d}_sl.zarr")
    sl_list.append(ds)

pl = xr.concat(pl_list, dim="time")
sl = xr.concat(sl_list, dim="time")

met = MetDataset(pl).standardize_variables(Cocip.ecmwf_met_variables())
rad = MetDataset(sl).standardize_variables(Cocip.ecmwf_rad_variables())

Aircraft performance modeling

The flights-pipeline uses the Poll-Schumann aircraft performance model for supported aircraft types with BADA as a fallback. Flights using aircraft types not supported by either aircraft performance model are ejected. Our sample flight uses an Airbus A-320, which is supported by the Poll-Schumann model.

from pycontrails.models.ps_model import PSFlight

ap_model = PSFlight(
    fill_low_altitude_with_isa_temperature=True,
    fill_low_altitude_with_zero_wind=True,
)
aircraft_type = df["aircraft_type_icao"][0]
assert ap_model.check_aircraft_type_availability(aircraft_type)

CoCiP configuration

The CoCiP configuration used for the Impact Inventory is based on Teoh et. al. 2024 with minor modification (e.g., the removal of a 13 km maximum altitude).

import numpy as np

from pycontrails.core import Flight
from pycontrails.models.humidity_scaling import ExponentialBoostLatitudeCorrectionHumidityScaling

inventory_params = {
    # Use short model time step to resolve lifecycle of short-lived contrails
    "dt_integration": np.timedelta64(5, "m"),
    # Use humidity scaling from Teoh et al 2024
    # https://doi.org/10.5194/acp-24-6071-2024
    "humidity_scaling": ExponentialBoostLatitudeCorrectionHumidityScaling(),
    # Use maximum contrail age of 12 h
    # Contrails older than 12 h contribute relatively little to simulated climate impact
    "max_age": np.timedelta64(12, "h"),
    # Remaining parameters are primarily for performance but may affect reproducibility
    "interpolation_use_indices": True,
    "interpolation_bounds_error": False,
    "filter_sac": True,
    "copy_source": True,
    "met_longitude_buffer": (10.0, 10.0),
    "met_latitude_buffer": (10.0, 10.0),
    "met_level_buffer": (20, 20),
}

inventory_fl = Flight(
    longitude=df["longitude"],
    latitude=df["latitude"],
    altitude_ft=df["altitude_ft"],
    time=pd.to_datetime(df["time"]),
    attrs={
        "flight_id": flight_id,
        "aircraft_type": aircraft_type,
        "engine_uid": df["engine_uid"][0],
    },
)

cocip = Cocip(met, rad, **inventory_params)
inventory_result = cocip.eval(inventory_fl, aircraft_performance=ap_model)

The energy forcing simulated by CoCiP matches the energy forcing from the impact inventory to about 1 part in 1 million.

print(f"Contrail EF (impact inventory):\t\t\t{inventory["sum_ef_mj"].sum():.6e} MJ")
print(
    f"Contrail EF (CoCiP, inventory configuration):\t{1e-6 * inventory_result["ef"].sum():.6e} MJ"
)

Contrail EF (impact inventory):                 1.628445e+09 MJ
Contrail EF (CoCiP, inventory configuration):   1.628444e+09 MJ

Reproducibility using v1/inventory/segments trajectories

Impact Inventory results served from the v1/inventory/segments endpoint include full trajectories. Each 1-minute segment is an output from CoCiP and has a lat_start, lat_end, lon_start, lon_end, time_start, time_end, and mean_altitude_ft. CoCiP breaks segments up in a forward-looking way. That is, the 0th element represents the segment from the 0th to 1st elements of the input. The final element is dropped, so CoCiP outputs and outputs of the v1/inventory/segments endpoint are one element shorter than the trajectory inputs used as inputs to the CoCiP simulations.

Trajectories from the v1/inventory/segments endpoint can be used as input to CoCiP simulations, but results will not match CoCiP output from the v1/inventory/segments endpoint exactly. The key difference is that Impact Inventory results provide the mean altitude in feet for each segment, whereas input trajectories used to generate the Impact Inventory include altitudes at segment endpoints. This means the altitude from the v1/inventory/segments endpoint does not exactly match the altitudes in the input trajectory used to generate the Impact Inventory. Using the trajectory from the v1/inventory/segments to try to reproduce the CoCiP results from the same endpoint output produces small differences due to accumulation of small differences between the trajectories, particularly during ascents and descents.

Additionally, the results stored in the Impact Inventory are truncated to whole megajoules per minute-segment, which can also accumulate some small differences in EF compared with raw CoCiP outputs.

Research-grade simulations

Although the CoCiP configuration used for the Impact Inventory is based on peer-reviewed scientific research, the Inventory’s emphasis on stability precludes frequent updates as the science of contrail modeling improves.

As a result, we maintain a separate recommended model configuration that we use for research-grade simulations with recent versions of pycontrails. This configuration emphasizes rapid development at the expense of stability, and users should expect frequent changes as scientific understanding evolves.

We recommend that this configuration be used primarily for scientific research. It may also be suitable for production-grade pipelines without a model stability requirement, though we caution users that–while we make our best efforts to recommend updates only when they improve model fidelity–we cannot guarantee that the current recommended configuration is the “best” configuration for all applications.

Contrails API

Output from research-grade simulations can be generated using the v0/trajectory/cocip endpoint.

The v0/trajectory endpoints compute results on-the-fly, so queries are expensive (slow) but updates to methodology are cheap (updated more frequently to keep pace with the models’ cutting edge).

url = "https://api.contrails.org/v0/trajectory/cocip"
headers = {"x-api-key": os.environ["CONTRAILS_API_KEY"]}
data = {
    "longitude": df["longitude"].tolist(),
    "latitude": df["latitude"].tolist(),
    "altitude": df["altitude_ft"].tolist(),
    "time": df["time"].tolist(),
    "aircraft_type": df["aircraft_type_icao"][0],
    "engine_uid": df["engine_uid"][0],
}

resp = requests.post(url, headers=headers, json=data)
trajectory = pd.DataFrame.from_dict(resp.json())
print(f"Contrail EF (trajectory endpoint): {1e-6 * trajectory["energy_forcing"].sum():.6e} MJ")

Contrail EF (trajectory endpoint): 1.253250e+09 MJ

Equivalent pycontrails code

Results from the v0/trajectory/cocip endpoint use pycontrails v0.63.4 (required for some recommended research-grade parameters).

import pycontrails

assert pycontrails.__version__ == "0.63.4"

Meteorology

For flights far enough in the past for ERA5 data to be available, the v0/trajectory/cocip endpoint uses the same internal zarr store as the Impact Inventory. For recent flights, the v0/trajectory/cocip endpoint may use IFS forecast data instead. The meteorology data source is noted in the response body.

trajectory[["met_source_provider", "met_source_dataset", "met_source_product"]].head()

	met_source_provider	met_source_dataset	met_source_product
0	ECMWF	ERA5	reanalysis
1	ECMWF	ERA5	reanalysis
2	ECMWF	ERA5	reanalysis
3	ECMWF	ERA5	reanalysis
4	ECMWF	ERA5	reanalysis

Aircraft performance modeling

Our research-grade recommendation for aircraft performance modeling also matches the Impact Inventory: use the Poll-Schumann model by default with BADA (if available) as a fallback.

Model parameters

However, our research-grade recommendations for CoCiP model parameters differ significantly from those used in `Teoh et. al. 2024 <>`__ and the Impact Inventory.

from pycontrails.models.emissions import Emissions

research_params = {
    # Use short model time step to resolve lifecycle of short-lived contrails
    # Matches inventory parmeters
    "dt_integration": np.timedelta64(5, "m"),
    # Use humidity scaling from Teoh et al 2024
    # https://doi.org/10.5194/acp-24-6071-2024
    # Matches inventory parameters
    "humidity_scaling": ExponentialBoostLatitudeCorrectionHumidityScaling(),
    # Use maximum contrail age of 12 h
    # Contrails older than 12 h contribute relatively little to simulated climate impact
    # Matches inventory parameters
    "max_age": np.timedelta64(12, "h"),
    # Use MEEM2/SCOPE11 methodology for emissions calculations
    # https://doi.org/10.4271/2025-01-6000
    # https://doi.org/10.1021/acs.est.8b04060
    "emissions": Emissions(use_meem=True),
    # Model VPM activation based on extended K15 model
    # https://doi.org/10.5194/acp-25-18617-2025
    "vpm_activation": True,
    # Model crystal depletion in vortex downwash based on LES results
    # https://doi.org/10.5194/acp-16-2059-2016
    "unterstrasser_ice_survival_fraction": True,
    # Allow differential radiative heating to enhance vertical mixing
    "radiative_heating_effects": True,
    # Constrain contrail properties.
    # In most cases contrails will dissipate before
    # exceeding any of these limits, but we use them
    # in global simulations for rare cases where contrail
    # properties exceed reasonable bounds.
    "min_altitude_m": 6500.0,
    "max_altitude_m": 15000.0,
    "max_depth": None,
    "max_horizontal_diffusivity": 100.0,
    "max_vertical_diffusivity": 10.0,
}

v0 endpoints convert flight coordinates to single-precision floating point values before running CoCiP, and reproducing results as precisely as possible requires using the same floating-point precision when running CoCiP locally.

research_fl = Flight(
    longitude=df["longitude"].astype(np.float32),
    latitude=df["latitude"].astype(np.float32),
    altitude_ft=df["altitude_ft"].astype(np.float32),
    time=pd.to_datetime(df["time"]),
    attrs={
        "flight_id": df["flight_id"][0],
        "aircraft_type": df["aircraft_type_icao"][0],
        "engine_uid": df["engine_uid"][0],
    },
)

cocip = Cocip(met, rad, **research_params)
research_result = cocip.eval(research_fl, aircraft_performance=ap_model)

/home/ext_tristan_abbott_breakthroughe/.local/share/virtualenvs/research/lib/python3.12/site-packages/pycontrails/models/extended_k15.py:535: UserWarning: This model is a minimal framework used to approximate the apparent
            emission index of contrail ice crystals in the jet regime. It does not fully
            represent the complexity of microphysical plume processes, including the
            formation and growth of vPM. Instead, vPM properties are prescribed as model
            inputs, which strongly impact model outputs. Therefore, the model should
            only be used for research purposes, together with thorough sensitivity
            analyses or explicit reference to the limitations outlined above.

  warnings.warn(

The energy forcing simulated by CoCiP matches the energy forcing from the v0/trajectory/cocip endpoint to about 1 part in 10,000.

print(f"Contrail EF (trajectory endpoint):\t\t{1e-6 * trajectory["energy_forcing"].sum():.6e} MJ")
print(f"Contrail EF (CoCiP, research configuration):\t{1e-6 * research_result["ef"].sum():.6e} MJ")

Contrail EF (trajectory endpoint):              1.253250e+09 MJ
Contrail EF (CoCiP, research configuration):    1.625606e+09 MJ

This small difference comes from rounding energy forcing values to multiples of 10 GJ in the API response. The API result can be reproduced to within 10 GJ by applying the same rounding logic:

print(
    f"Contrail EF (CoCiP, research configuration):\t{1e-6 * research_result["ef"].astype(float).round(-10).sum():.6e} MJ"
)

Contrail EF (CoCiP, research configuration):    1.625620e+09 MJ