ACCF#

Interface to the ACCF (algorithmic climate change functions) from the DLR library climaccf.

References#

  • Dietmüller, Simone. “Dlr-Pa/Climaccf: Dataset Update for GMDD.” Zenodo, September 13, 2022. https://doi.org/10.5281/zenodo.7074582.

  • Dietmüller, Simone, Sigrun Matthes, Katrin Dahlmann, Hiroshi Yamashita, Abolfazl Simorgh, Manuel Soler, Florian Linke, et al. “A Python Library for Computing Individual and Merged Non-CO2 Algorithmic Climate Change Functions: CLIMaCCF V1.0.” Preprint. Atmospheric sciences, October 17, 2022. https://doi.org/10.5194/gmd-2022-203.

[1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from pycontrails import Flight
from pycontrails.datalib.ecmwf import ERA5
from pycontrails.models.accf import ACCF
from pycontrails.physics import units

Set Domain#

time / altitude

[2]:
time_bounds = ("2022-01-01 12:00:00", "2022-01-01 15:00:00")
pressure_levels = [200, 225, 250]

Download meteorology data from ECMWF#

This demo uses ERA5 via the Copernicus Data Store (CDS) for met data. This requires account with Copernicus Data Portal and local ~/.cdsapirc file with credentials.

[3]:
# pressure level data
era5_pl = ERA5(
    time=time_bounds,
    variables=[
        "air_temperature",
        "specific_humidity",
        "potential_vorticity",
        "geopotential",
        "relative_humidity",
        "northward_wind",
        "eastward_wind",
    ],
    pressure_levels=pressure_levels,
)

# single level data (radiation)
era5_sl = ERA5(
    time=time_bounds,
    variables=["surface_solar_downward_radiation", "top_net_thermal_radiation"],
)

pl = era5_pl.open_metdataset()
sl = era5_sl.open_metdataset()

Create example flight#

[4]:
# Note: ACCF output is in K/kg of fuel so you will need to provide fuel burn of the flight to get impact in K
n = 10000
longitude = np.linspace(45, 58, n) + np.linspace(0, 1, n)
latitude = np.linspace(45, 54, n) - np.linspace(0, 1, n)
altitude = units.pl_to_m(225) * np.ones(n)

start = np.datetime64(time_bounds[0])
end = start + np.timedelta64(90, "m")  # 90 minute flight
time = pd.date_range(start, end, periods=n)
fl = Flight(
    longitude=longitude,
    latitude=latitude,
    altitude=altitude,
    time=time,
    aircraft_type="B737",
    flight_id=17,
)
fl["fuel_flow"] = np.linspace(2.1, 2.4, n)  # kg/s
[5]:
fl
[5]:
Flight [5 keys x 10000 length, 3 attributes]

Attributes
time[2022-01-01 12:00:00, 2022-01-01 13:30:00]
longitude[45.0, 59.0]
latitude[45.0, 53.0]
altitude[11037.0, 11037.0]
aircraft_typeB737
flight_id17
crsEPSG:4326
longitude latitude time altitude fuel_flow
0 45.000000 45.0000 2022-01-01 12:00:00.000000000 11037.01126 2.10000
1 45.001400 45.0008 2022-01-01 12:00:00.540054005 11037.01126 2.10003
2 45.002800 45.0016 2022-01-01 12:00:01.080108010 11037.01126 2.10006
3 45.004200 45.0024 2022-01-01 12:00:01.620162016 11037.01126 2.10009
4 45.005601 45.0032 2022-01-01 12:00:02.160216021 11037.01126 2.10012
... ... ... ... ... ...
9995 58.994399 52.9968 2022-01-01 13:29:57.839783978 11037.01126 2.39988
9996 58.995800 52.9976 2022-01-01 13:29:58.379837983 11037.01126 2.39991
9997 58.997200 52.9984 2022-01-01 13:29:58.919891989 11037.01126 2.39994
9998 58.998600 52.9992 2022-01-01 13:29:59.459945994 11037.01126 2.39997
9999 59.000000 53.0000 2022-01-01 13:30:00.000000000 11037.01126 2.40000

10000 rows × 5 columns

Initialize ACCF model and evaluate on flight#

[7]:
ac = ACCF(met=pl, surface=sl)
[8]:
# Evaluate ACCFs over this flight
fl = ac.eval(fl)
selected forecast step:  6.0 HR

 UserWarning: For this configuration formation of persistent contrails is possible, if temperatures are low enough (below 235K) and relative humidity (with respect to ice) is above or at 90.0%. However keep in mind that the threshold value for the relative humidity varies with the used forecast model and its resolution. In order to choose the appropriate threshold value, you should read the details given in Section 5.1 of the connected publication of Dietmueller et al. 2022 

Explore the output#

Raw output:

[9]:
fl
[9]:
Flight [27 keys x 10000 length, 7 attributes]

Attributes
time[2022-01-01 12:00:00, 2022-01-01 13:30:00]
longitude[45.0, 59.0]
latitude[45.0, 53.0]
altitude[11037.0, 11037.0]
aircraft_typeB737
flight_id17
crsEPSG:4326
met_source_providerECMWF
met_source_datasetERA5
met_source_productreanalysis
pycontrails_version0.47.3.dev86
longitude latitude time altitude fuel_flow flight_id air_temperature specific_humidity potential_vorticity geopotential ... aCCF_CH4 aCCF_O3 aCCF_H2O aCCF_nCont aCCF_dCont aCCF_Cont aCCF_CO2 aCCF_merged aCCF_NOx Fin
0 45.000000 45.0000 2022-01-01 12:00:00.000000000 11037.01126 2.10000 17.0 215.779129 0.000007 0.000007 105447.960938 ... -4.455085e-13 1.637062e-12 7.799343e-16 0.000000e+00 -0.000000e+00 0.000000e+00 6.940000e-16 1.696413e-14 1.191553e-12 509.002230
1 45.001400 45.0008 2022-01-01 12:00:00.540054005 11037.01126 2.10003 17.0 215.779785 0.000007 0.000007 105447.773438 ... -4.455091e-13 1.637061e-12 7.799481e-16 0.000000e+00 -0.000000e+00 0.000000e+00 6.940000e-16 1.696413e-14 1.191552e-12 508.984589
2 45.002800 45.0016 2022-01-01 12:00:01.080108010 11037.01126 2.10006 17.0 215.780441 0.000007 0.000007 105447.578125 ... -4.455098e-13 1.637061e-12 7.799622e-16 0.000000e+00 -0.000000e+00 0.000000e+00 6.940000e-16 1.696412e-14 1.191551e-12 508.966947
3 45.004200 45.0024 2022-01-01 12:00:01.620162016 11037.01126 2.10009 17.0 215.781097 0.000007 0.000007 105447.390625 ... -4.455104e-13 1.637060e-12 7.799766e-16 0.000000e+00 -0.000000e+00 0.000000e+00 6.940000e-16 1.696412e-14 1.191550e-12 508.949306
4 45.005601 45.0032 2022-01-01 12:00:02.160216021 11037.01126 2.10012 17.0 215.781754 0.000007 0.000007 105447.195312 ... -4.455110e-13 1.637060e-12 7.799912e-16 0.000000e+00 -0.000000e+00 0.000000e+00 6.940000e-16 1.696412e-14 1.191549e-12 508.931664
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
9995 58.994399 52.9968 2022-01-01 13:29:57.839783978 11037.01126 2.39988 17.0 209.364670 0.000013 0.000007 104271.609375 ... -4.495710e-13 1.523755e-12 7.788148e-16 1.717789e-28 -1.052829e-27 1.717789e-28 6.940000e-16 1.543721e-14 1.074184e-12 328.603173
9996 58.995800 52.9976 2022-01-01 13:29:58.379837983 11037.01126 2.39991 17.0 209.364761 0.000013 0.000007 104271.421875 ... -4.495715e-13 1.523748e-12 7.789030e-16 1.717655e-28 -1.052830e-27 1.717655e-28 6.940000e-16 1.543720e-14 1.074177e-12 328.584764
9997 58.997200 52.9984 2022-01-01 13:29:58.919891989 11037.01126 2.39994 17.0 209.364868 0.000013 0.000007 104271.234375 ... -4.495719e-13 1.523742e-12 7.789914e-16 1.717522e-28 -1.052832e-27 1.717522e-28 6.940000e-16 1.543720e-14 1.074170e-12 328.566356
9998 58.998600 52.9992 2022-01-01 13:29:59.459945994 11037.01126 2.39997 17.0 209.364960 0.000013 0.000007 104271.046875 ... -4.495724e-13 1.523735e-12 7.790802e-16 1.717388e-28 -1.052834e-27 1.717388e-28 6.940000e-16 1.543720e-14 1.074163e-12 328.547947
9999 59.000000 53.0000 2022-01-01 13:30:00.000000000 11037.01126 2.40000 17.0 209.365051 0.000013 0.000007 104270.859375 ... -4.495728e-13 1.523728e-12 7.791693e-16 1.717254e-28 -1.052836e-27 1.717254e-28 6.940000e-16 1.543719e-14 1.074156e-12 328.529539

10000 rows × 27 columns

Now convert fuel flow to total fuel burned by waypoint, and get the impact of each flight waypoint and plot

[11]:
# Waypoint duration in seconds
dt_sec = fl.segment_duration()

# kg fuel per contrail
fuel_burn = fl["fuel_flow"] * dt_sec

# Get impacts in degrees K per waypoint
warming_contrails = fuel_burn * fl["aCCF_Cont"]
warming_merged = fuel_burn * fl["aCCF_merged"]
[12]:
f, (ax0, ax1) = plt.subplots(2, 1, sharex=True, figsize=(12, 8))
ax0.plot(fl["time"], warming_contrails, label="Contrails")
ax0.plot(fl["time"], warming_merged, label="Combined ACCFs")
ax0.set_ylabel("Degrees K")
ax0.set_title("Warming impact by waypoint")
ax0.legend()

ax1.plot(fl["time"], np.cumsum(warming_contrails), label="Contrails")
ax1.plot(fl["time"], np.cumsum(warming_merged), label="Combined ACCFs")
ax1.legend()
ax1.set_xlabel("Waypoint Time")
ax1.set_ylabel("Degrees K")
ax1.set_title("Cumulative warming impact");
../_images/integrations_ACCF_16_0.png

Grid ACCF with Parameters#

[13]:
ac = ACCF(pl, sl, {"emission_scenario": "future_scenario", "accf_v": "V1.1"})
ds = ac.eval()  # This will evaluate over the met grid
selected forecast step:  6.0 HR

 UserWarning: For this configuration formation of persistent contrails is possible, if temperatures are low enough (below 235K) and relative humidity (with respect to ice) is above or at 90.0%. However keep in mind that the threshold value for the relative humidity varies with the used forecast model and its resolution. In order to choose the appropriate threshold value, you should read the details given in Section 5.1 of the connected publication of Dietmueller et al. 2022 

[14]:
fig, ax = plt.subplots()
p = ax.pcolor(
    ds["longitude"].data,
    ds["latitude"].data,
    ds["aCCF_merged"].data[:, :, 0, 0].T,
    cmap="coolwarm",
    vmin=-6e-13,
    vmax=6e-13,
)
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.set_title("Merged ACCF")
cbar = plt.colorbar(p)
cbar.set_label("Degrees K/km");
../_images/integrations_ACCF_19_0.png