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 matplotlib.pyplot as plt
import numpy as np
import pandas as pd
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]:
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_type | B737 |
flight_id | 17 |
crs | EPSG: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]:
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_type | B737 |
flight_id | 17 |
crs | EPSG:4326 |
met_source_provider | ECMWF |
met_source_dataset | ERA5 |
met_source_product | reanalysis |
pycontrails_version | 0.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");
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");