[BUG] Soil sink in EmisCH4_Total and sectoral breakdown

[djvaron]

Name: Daniel Varon Institution: Harvard University

Soil sink bug when computing total emissions, sectoral breakdown

The visualization notebook creates the posterior emissions from the prior emissions and scale factors:

# Prior emissions
prior = xr.load_dataset(prior_pth)["EmisCH4_Total"].isel(time=0)
if config["KalmanMode"]:
    prior_sf = xr.load_dataset(prior_sf_pth)
    prior = prior * prior_sf["ScaleFactor"]

# Optimized scale factors
scale = xr.load_dataset(results_pth)["ScaleFactor"]

# Posterior emissions
posterior = prior * scale

The soil absorption sink should be removed from EmisCH4_Total before applying posterior scale factors, because soil absorption is not currently optimized in the inversion. Not doing so leads to incorrect posterior emission output and sectoral breakdown.

This issue may affect other parts of the codebase where EmisCH4_Total is called, including the Kalman filter feature (e.g., prepare_sf.py)

The posterior simulation should NOT be affected because soil absorption is correctly removed before the scale factors are applied.

[djvaron]

Here are potential fixes for the visualization notebook.

Defining the prior and posterior emissions:

# Prior emissions
hemco_emissions = xr.load_dataset(prior_pth).isel(time=0)
prior = hemco_emissions["EmisCH4_Total"].copy()
prior_noSoilAbsorb = prior - hemco_emissions["EmisCH4_SoilAbsorb"]

# Optimized scale factors
scale = xr.load_dataset(results_pth)["ScaleFactor"]

# Posterior emissions
posterior = prior_noSoilAbsorb * scale + hemco_emissions["EmisCH4_SoilAbsorb"]

Sectoral breakdown:

def my_autopct(
    pct,
):  # Only plots the data labels for sectoral emissions that are larger than 15%
    return ("%.1f%%" % pct) if pct > 15 else ""

# load in full sectors based on prior distribution
posterior_sectors = xr.load_dataset(prior_pth) * scale

# overwrite scaled soil absorption with original soil absorption
soil_absorb = xr.load_dataset(prior_pth)["EmisCH4_SoilAbsorb"]
posterior_sectors["EmisCH4_SoilAbsorb"] = soil_absorb

# extract sector names
sector_list = [var for var in list(posterior_sectors.keys()) if "EmisCH4" in var]
sector_list.remove("EmisCH4_Total")

# calculate total emissions for each sector and print
emission_types = {}
for sector in sector_list:
    emission = sum_total_emissions(
        posterior_sectors[sector].isel(time=0), areas, mask
    )
    print(f"{sector} = {emission} Tg/yr")
    if emission > 0:
        emission_types[sector] = emission

title = plt.title("CH4 emissions by sector")
title.set_ha("center")
plt.gca().axis("equal")

pie = plt.pie(np.array(list(emission_types.values()))/total_posterior_emissions, startangle=0, autopct=my_autopct)
plt.legend(
    pie[0],
    emission_types.keys(),
    bbox_to_anchor=(1, 0.5),
    loc="center right",
    fontsize=10,
    bbox_transform=plt.gcf().transFigure,
)
plt.subplots_adjust(left=0.0, bottom=0.1, right=0.65)

[laestrada]

Closing with #213