pycontrails.models.cocip.Cocip

class pycontrails.models.cocip.Cocip(met, rad, params=None, **params_kwargs)

Bases: Model

Contrail Cirrus Prediction Model (CoCiP).

Published by Ulrich Schumann et. al. (DLR Institute of Atmospheric Physics) in [Schumann, 2012], [Schumann et al., 2012].

Parameters:

• met (MetDataset) – Pressure level dataset containing met_variables variables. See Notes for variable names by data source.

• rad (MetDataset) – Single level dataset containing top of atmosphere radiation fluxes. See Notes for variable names by data source.

• params (dict[str, Any], optional) – Override Cocip model parameters with dictionary. See CocipFlightParams for model parameters.

• **params_kwargs (Any) – Override Cocip model parameters with keyword arguments. See CocipFlightParams for model parameters.

Notes

Inputs

The required meteorology variables depend on the data source (e.g. ECMWF, GFS).

See met_variables and rad_variables for the list of required variables to the met and rad parameters, respectively. When an item in one of these arrays is a tuple, variable keys depend on data source.

The current list of required variables (labelled by "standard_name"):

Variable keys for pressure level data

Parameter	ECMWF	GFS
Air Temperature	air_temperature	air_temperature
Specific Humidity	specific_humidity	specific_humidity
Eastward wind	eastward_wind	eastward_wind
Northward wind	northward_wind	northward_wind
Vertical velocity	lagrangian_tendency_of_air_pressure	lagrangian_tendency_of_air_pressure
Ice water content	specific_cloud_ice_water_content	ice_water_mixing_ratio

Variable keys for single-level radiation data

Parameter	ECMWF	GFS
Top solar radiation	top_net_solar_radiation	toa_upward_shortwave_flux
Top thermal radiation	top_net_thermal_radiation	toa_upward_longwave_flux

Modifications

This implementation differs from original CoCiP (Fortran) implementation in a few places:

• This model uses aircraft performance and emissions models to calculate nvPM, fuel flow, and overall propulsion efficiency, if not already provided.

• As described in [Teoh et al., 2022], this implementation sets the initial ice particle activation rate to be a function of the difference between the ambient temperature and the critical SAC threshold temperature. See pycontrails.models.sac.T_critical_sac().

• Isobaric heat capacity calculation. The original model uses a constant value of 1004 \(J \ kg^{-1} \ K^{-1}\), whereas this model calculates isobaric heat capacity as a function of specific humidity. See pycontrails.physics.thermo.c_pm().

• Solar direct radiation. The original algorithm uses ECMWF radiation variable tisr (top incident solar radiation) as solar direct radiation value. This implementation calculates the theoretical solar direct radiation at any arbitrary point in the atmosphere. See pycontrails.physics.geo.solar_direct_radiation().

• Segment angle. The segment angle calculations for flights and contrail segments have been updated to use more precise spherical geometry instead of a triangular approximation. As the triangle approaches zero, the two calculations agree. See pycontrails.physics.geo.segment_angle().

• Integration. This implementation consistently uses left-Riemann sums in the time integration of contrail segments.

• Segment length ratio. Instead of taking a geometric mean between contrail segments before/after advection, a simple ratio is computed. See contrail_properties.segment_length_ratio().

• Segment energy flux. This implementation does not average spatially contiguous contrail segments when calculating the mean energy flux for the segment of interest. See contrail_properties.mean_energy_flux_per_m().

This implementation is regression tested against results from [Teoh et al., 2022]. See tests/benchmark/north-atlantic-study/validate.py.

Outputs

NaN values may appear in model output. Specifically, np.nan values are used to indicate:

• Flight waypoint or contrail waypoint is not contained with the met domain.

• The variable was NOT computed during the model evaluation. For example, at flight waypoints not producing any persistent contrails, “radiative” variables (rsr, olr, rf_sw, rf_lw, rf_net) are not computed. Consequently, the corresponding values in the output of eval() are NaN. One exception to this rule is found on ef (energy forcing) contrail_age predictions. For these two “cumulative” variables, waypoints not producing any persistent contrails are assigned 0 values.

References

• [Schumann and Wendling, 1990]

• [Schumann, 2010]

• [Voigt et al., 2010]

• [Schumann et al., 2011]

• [Schumann, 2012]

• [Schumann et al., 2012]

• [Schumann et al., 2012]

• [Schumann et al., 2011]

• [Schumann et al., 2015]

• [Teoh et al., 2020]

• [Schumann et al., 2021]

• [Schumann et al., 2021]

• [Teoh et al., 2022]

See also

CocipFlightParams, wake_vortex, contrail_properties, radiative_forcing, humidity_scaling, Emissions, sac, tau_cirrus

__init__(met, rad, params=None, **params_kwargs)

Methods

__init__(met, rad[, params])
downselect_met()	Downselect met domain to the max/min bounds of source.
eval([source])	Run CoCiP simulation on flight.
get_source_param(key[, default, set_attr])	Get source data with default set by parameter key.
require_met()	Ensure that met is a MetDataset.
require_source_type(type_)	Ensure that source is type_.
set_source([source])	Attach original or copy of input source to source.
set_source_met([optional, variable])	Ensure or interpolate each required met_variables on source .
transfer_met_source_attrs([source])	Transfer met source metadata from met to source.
update_params([params])	Update model parameters on params.

Attributes

rad	Radiation data formatted as a MetDataset at a single pressure level [-1]
contrail	Contrail evolution output from model.
contrail_dataset	xr.Dataset representation of contrail evolution.
contrail_list	List of GeoVectorDataset contrail objects - one for each timestep
timesteps	Array of numpy.datetime64 time steps for contrail evolution
hash	Generate a unique hash for model instance.
interp_kwargs	Shortcut to create interpolation arguments from params.
long_name
met	Met data is not optional
met_required	Require meteorology is not None on __init__()
met_variables	Required meteorology pressure level variables.
name
optional_met_variables	Additional met variables used to support outputs
params	Instantiated model parameters, in dictionary form
processed_met_variables	Minimal set of met variables needed to run the model after pre-processing.
rad_variables	Required single-level top of atmosphere radiation variables.
source	Last Flight modeled in eval()

contrail

Contrail evolution output from model.

Set to None when no contrails are formed. Otherwise, this is a pandas.DataFrame describing the evolution of the contrail. Columns include:

• waypoint: The index of the waypoint in the original flight creating the contrail. This can be used to join the contrail DataFrame to the source.

• formation_time: Time of contrail formation. Agrees with the time column in source.

• continuous: Boolean indicating whether the contrail is continuous or not.

• persistent: Boolean indicating whether the contrail is persistent or not. A contrail segment is considered continuous if both the current and the next contrail waypoint at the same time step persist.

• segment_length: Length of the contrail segment, [\(m\)].

• sin_a, cos_a: Sine and cosine of the segment angle.

• width, depth: Contrail width and depth, [\(m\)].

• sigma_yz: The yz component of the covariance matrix, [\(m^{2}\)]. See contrail_properties.plume_temporal_evolution().

• q_sat: Saturation specific humidity over ice, [\(kg \ kg^{-1}\)].

• n_ice_per_m: Number of ice particles per distance, [\(m^{-1}\)].

• iwc: Ice water content, [\(kg_{ice} kg_{air}^{-1}\)].

• tau_contrail: Optical depth of the contrail. See contrail_properties.contrail_optical_depth().

• rf_sw, rf_lw, rf_net: Shortwave, longwave, and net instantaneous radiative forcing, [\(W \ m^{-2}\)] at the contrail waypoint.

• ef: Energy forcing, [\(J\)] at the contrail waypoint. See contrail_properties.energy_forcing().

contrail_dataset

xr.Dataset representation of contrail evolution.

contrail_list

List of GeoVectorDataset contrail objects - one for each timestep

default_params

alias of CocipFlightParams

eval(source=None, **params)

Run CoCiP simulation on flight.

Simulates the formation and evolution of contrails from a Flight using the contrail cirrus prediction model (CoCiP) from Schumann (2012) [Schumann, 2012].

Changed in version 0.25.11: Previously, any waypoint not surviving the wake vortex downwash phase of CoCiP was assigned a nan-value in the ef array within the model output. This is no longer the case. Instead, energy forcing is set to 0.0 for all waypoints which fail to produce persistent contrails. In particular, nan values in the ef array are only used to indicate an out-of-met-domain waypoint. The same convention is now used for output variables contrail_age and cocip as well.

Parameters:

• source (Flight | Sequence[Flight] | None) – Input Flight(s) to model.

• **params (Any) – Overwrite model parameters before eval.

Returns:

Flight | list[Flight] | NoReturn – Flight(s) with updated Contrail data. The model parameter “verbose_outputs” determines the variables on the return flight object.

References

• [Schumann, 2012]

long_name = 'Contrail Cirrus Prediction Model'

met_required = True

Require meteorology is not None on __init__()

met_variables = (MetVariable(short_name='t', standard_name='air_temperature', long_name='Air Temperature', level_type='isobaricInhPa', ecmwf_id=130, grib1_id=11, grib2_id=(0, 0, 0), units='K', amip='ta', description='Air temperature is the bulk temperature of the air, not the surface (skin) temperature.'), MetVariable(short_name='q', standard_name='specific_humidity', long_name='Specific Humidity', level_type='isobaricInhPa', ecmwf_id=133, grib1_id=51, grib2_id=(0, 1, 0), units='kg kg**-1', amip='hus', description='Specific means per unit mass. Specific humidity is the mass fraction of water vapor in (moist) air.'), MetVariable(short_name='u', standard_name='eastward_wind', long_name='Eastward Wind', level_type='isobaricInhPa', ecmwf_id=131, grib1_id=33, grib2_id=(0, 2, 2), units='m s**-1', amip='ua', description='"Eastward" indicates a vector component which is positive when directed eastward (negative westward). Wind is defined as a two-dimensional (horizontal) air velocity vector, with no vertical component.'), MetVariable(short_name='v', standard_name='northward_wind', long_name='Northward Wind', level_type='isobaricInhPa', ecmwf_id=132, grib1_id=34, grib2_id=(0, 2, 3), units='m s**-1', amip='va', description='"Northward" indicates a vector component which is positive when directed northward (negative southward). Wind is defined as a two-dimensional (horizontal) air velocity vector, with no vertical component.'), MetVariable(short_name='w', standard_name='lagrangian_tendency_of_air_pressure', long_name='Vertical Velocity (omega)', level_type='isobaricInhPa', ecmwf_id=135, grib1_id=39, grib2_id=(0, 2, 8), units='Pa s**-1', amip='wap', description='The Lagrangian tendency of air pressure, often called "omega", plays the role of the upward component of air velocity when air pressure is being used as the vertical coordinate. If the vertical air velocity is upwards, it is negative when expressed as a tendency of air pressure; downwards is positive. Air pressure is the force per unit area which would be exerted when the moving gas molecules of which the air is composed strike a theoretical surface of any orientation.'), (MetVariable(short_name='ciwc', standard_name='specific_cloud_ice_water_content', long_name='Specific cloud ice water content', level_type='isobaricInhPa', ecmwf_id=247, grib1_id=None, grib2_id=(0, 1, 84), units='kg kg**-1', amip=None, description="This parameter is the mass of cloud ice particles per kilogram of the total mass of moist air. The 'total mass of moist air' is the sum of the dry air, water vapour, cloud liquid, cloud ice, rain and falling snow. This parameter represents the average value for a grid box."), MetVariable(short_name='icmr', standard_name='ice_water_mixing_ratio', long_name='Cloud ice water mixing ratio', level_type='isobaricInhPa', ecmwf_id=260019, grib1_id=None, grib2_id=(0, 1, 23), units='kg kg**-1', amip=None, description='This parameter is the mass of cloud ice particles per kilogram of the total mass of dry air. ')))

Required meteorology pressure level variables. Each element in the list is a MetVariable or a tuple[MetVariable]. If element is a tuple[MetVariable], the variable depends on the data source. Only one variable in the tuple is required.

name = 'cocip'

optional_met_variables = ((MetVariable(short_name='z', standard_name='geopotential', long_name='Geopotential', level_type='isobaricInhPa', ecmwf_id=129, grib1_id=6, grib2_id=(0, 3, 4), units='m**2 s**-2', amip=None, description='Geopotential is the sum of the specific gravitational potential energy relative to the geoid and the specific centripetal potential energy.'), MetVariable(short_name='gh', standard_name='geopotential_height', long_name='Geopotential Height', level_type='isobaricInhPa', ecmwf_id=156, grib1_id=7, grib2_id=(0, 3, 5), units='m', amip='zg', description='Geopotential is the sum of the specific gravitational potential energy relative to the geoid and the specific centripetal potential energy. Geopotential height is the geopotential divided by the standard acceleration due to gravity. It is numerically similar to the altitude (or geometric height) and not to the quantity with standard name height, which is relative to the surface.')), (MetVariable(short_name='cc', standard_name='fraction_of_cloud_cover', long_name='Cloud area fraction in atmosphere layer', level_type='isobaricInhPa', ecmwf_id=248, grib1_id=None, grib2_id=(0, 6, 32), units='[0 - 1]', amip='cl', description='This parameter is the proportion of a grid box covered by cloud (liquid or ice) at a specific pressure level.'), MetVariable(short_name='tcc', standard_name='total_cloud_cover_isobaric', long_name='Total cloud cover at isobaric surface', level_type='isobaricInhPa', ecmwf_id=228164, grib1_id=None, grib2_id=(0, 6, 1), units='%', amip=None, description='This parameter is the percentage of a grid box covered by cloud (liquid or ice) at a specific pressure level.')))

Additional met variables used to support outputs

Changed in version 0.48.0: Moved Geopotential from met_variables to optional_met_variables

processed_met_variables = (MetVariable(short_name='t', standard_name='air_temperature', long_name='Air Temperature', level_type='isobaricInhPa', ecmwf_id=130, grib1_id=11, grib2_id=(0, 0, 0), units='K', amip='ta', description='Air temperature is the bulk temperature of the air, not the surface (skin) temperature.'), MetVariable(short_name='q', standard_name='specific_humidity', long_name='Specific Humidity', level_type='isobaricInhPa', ecmwf_id=133, grib1_id=51, grib2_id=(0, 1, 0), units='kg kg**-1', amip='hus', description='Specific means per unit mass. Specific humidity is the mass fraction of water vapor in (moist) air.'), MetVariable(short_name='u', standard_name='eastward_wind', long_name='Eastward Wind', level_type='isobaricInhPa', ecmwf_id=131, grib1_id=33, grib2_id=(0, 2, 2), units='m s**-1', amip='ua', description='"Eastward" indicates a vector component which is positive when directed eastward (negative westward). Wind is defined as a two-dimensional (horizontal) air velocity vector, with no vertical component.'), MetVariable(short_name='v', standard_name='northward_wind', long_name='Northward Wind', level_type='isobaricInhPa', ecmwf_id=132, grib1_id=34, grib2_id=(0, 2, 3), units='m s**-1', amip='va', description='"Northward" indicates a vector component which is positive when directed northward (negative southward). Wind is defined as a two-dimensional (horizontal) air velocity vector, with no vertical component.'), MetVariable(short_name='w', standard_name='lagrangian_tendency_of_air_pressure', long_name='Vertical Velocity (omega)', level_type='isobaricInhPa', ecmwf_id=135, grib1_id=39, grib2_id=(0, 2, 8), units='Pa s**-1', amip='wap', description='The Lagrangian tendency of air pressure, often called "omega", plays the role of the upward component of air velocity when air pressure is being used as the vertical coordinate. If the vertical air velocity is upwards, it is negative when expressed as a tendency of air pressure; downwards is positive. Air pressure is the force per unit area which would be exerted when the moving gas molecules of which the air is composed strike a theoretical surface of any orientation.'), MetVariable(short_name='tau_cirrus', standard_name='tau_cirrus', long_name='Cirrus optical depth', level_type=None, ecmwf_id=None, grib1_id=None, grib2_id=None, units='dimensionless', amip=None, description=None))

Minimal set of met variables needed to run the model after pre-processing. The intention here is that ciwc is unnecessary after tau_cirrus has already been calculated.

rad

Radiation data formatted as a MetDataset at a single pressure level [-1]

rad_variables = ((MetVariable(short_name='tsr', standard_name='top_net_solar_radiation', long_name='Top of atmosphere net solar (shortwave) radiation', level_type='nominalTop', ecmwf_id=178, grib1_id=None, grib2_id=(0, 4, 1), units='J m**-2', amip=None, description="This parameter is the incoming solar radiation (also known as shortwave radiation) minus the outgoing solar radiation at the top of the atmosphere. It is the amount of radiation passing through a horizontal plane. The incoming solar radiation is the amount received from the Sun. The outgoing solar radiation is the amount reflected and scattered by the Earth's atmosphere and surfaceSee https://www.ecmwf.int/sites/default/files/elibrary/2015/18490-radiation-quantities-ecmwf-model-and-mars.pdf"), MetVariable(short_name='uswrf', standard_name='toa_upward_shortwave_flux', long_name='Top of atmosphere upward shortwave radiation', level_type='nominalTop', ecmwf_id=None, grib1_id=None, grib2_id=(0, 4, 193), units='W m**-2', amip=None, description='This parameter is the outgoing shortwave (solar) radiation at the nominal top of the atmosphere.')), (MetVariable(short_name='ttr', standard_name='top_net_thermal_radiation', long_name='Top of atmosphere net thermal (longwave) radiation', level_type='nominalTop', ecmwf_id=179, grib1_id=None, grib2_id=(0, 5, 5), units='J m**-2', amip=None, description='The thermal (also known as terrestrial or longwave) radiation emitted to space at the top of the atmosphere is commonly known as the Outgoing Longwave Radiation (OLR). The top net thermal radiation (this parameter) is equal to the negative of OLR.See https://www.ecmwf.int/sites/default/files/elibrary/2015/18490-radiation-quantities-ecmwf-model-and-mars.pdf'), MetVariable(short_name='ulwrf', standard_name='toa_upward_longwave_flux', long_name='Top of atmosphere upward longwave radiation', level_type='nominalTop', ecmwf_id=None, grib1_id=None, grib2_id=(0, 5, 193), units='W m**-2', amip=None, description='This parameter is the outgoing longwave (thermal) radiation at the nominal top of the atmosphere.')))

Required single-level top of atmosphere radiation variables. Variable keys depend on data source (e.g. ECMWF, GFS).

timesteps

Array of numpy.datetime64 time steps for contrail evolution