Add optional CH4 effects

dynamic_characterization/dynamic_characterization.py

@@ -1,15 +1,12 @@
 import json
 import os
-import warnings
 from collections.abc import Collection
 from datetime import datetime
 from typing import Callable, Dict, Tuple
 from loguru import logger
 
-import bw2data as bd
 import numpy as np
 import pandas as pd
-from bw2data.utils import UnknownObject
 
 from dynamic_characterization.classes import CharacterizedRow
 from dynamic_characterization.ipcc_ar6.radiative_forcing import (
@@ -40,6 +37,7 @@ def characterize(
     characterization_function_co2: Callable = None,
     time_varying_re: bool = False,
     fallback_to_ipcc: bool = True,
+    indirect_effects: str = "all",
 ) -> pd.DataFrame:
     """
     Characterizes the dynamic inventory, formatted as a Dataframe, by evaluating each emission (row in DataFrame) using given dynamic characterization functions.
@@ -88,6 +86,11 @@ def characterize(
     fallback_to_ipcc: bool, optional
         Only for prospective metrics. If True (default), use IPCC AR6 characterization functions for GHGs not available
         in the Watanabe module (e.g., CO and other GHGs). If False, only GHGs available in Watanabe (CO2, CH4, N2O) are characterized.
+    indirect_effects : str
+        Indirect effect treatment. One of:
+        - all: all indirect effects. Matches Watanabe et al. Default.
+        - no_carbon_cyle: all indirect effects except carbon cycle feedbacks.
+        - none: only direct effects of the gasae.
 
     Returns
     -------
@@ -163,6 +166,7 @@ def characterize(
                     original_time_horizon=time_horizon,
                     dynamic_time_horizon=dynamic_time_horizon,
                     time_varying_re=time_varying_re,
+                    indirect_effects=indirect_effects,
                 )
             )
 
@@ -184,6 +188,7 @@ def characterize(
                     row=row,
                     time_horizon=dynamic_time_horizon,
                     time_varying_re=time_varying_re,
+                    indirect_effects=indirect_effects,
                 )
             )
 
@@ -244,6 +249,16 @@ def create_characterization_functions_from_method(
 
     characterization_functions = dict()
 
+    try:
+        import bw2data as bd
+        from bw2data.utils import UnknownObject
+    except ModuleNotFoundError as exc:
+        raise ModuleNotFoundError(
+            "bw2data is required to build default characterization functions from an LCIA method. "
+            "Install bw2data with its dependencies, or provide custom characterization_functions "
+            "when calling characterize()."
+        ) from exc
+
     filepath = os.path.join(
         os.path.dirname(os.path.abspath(__file__)),
         "ipcc_ar6",
@@ -421,23 +436,31 @@ def _characterize_pgwp(
     original_time_horizon,
     dynamic_time_horizon,
     time_varying_re: bool = False,
+    indirect_effects: str = "all",
 ) -> CharacterizedRow:
     """
     Calculate prospective GWP using Watanabe et al. (2026) characterization.
 
     Uses AGWP_gas / AGWP_CO2 to calculate kg CO2 equivalent.
     Emission year is extracted from the row's date.
     """
-    # Get emission year from the row date
-    emission_date = row.date
-    emission_year = int(str(emission_date.to_numpy())[:4])
 
     # Calculate AGWP for the gas using its characterization function
-    radiative_forcing_ghg = characterization_functions[row.flow](
-        row,
-        dynamic_time_horizon,
-        time_varying_re=time_varying_re,
-    )
+    char_func = characterization_functions[row.flow]
+    if char_func == prospective_characterize_ch4:
+        # Only CH4 supports indirect_effects
+        radiative_forcing_ghg = char_func(
+            row,
+            dynamic_time_horizon,
+            time_varying_re=time_varying_re,
+            indirect_effects=indirect_effects,
+        )
+    else:
+        radiative_forcing_ghg = char_func(
+            row,
+            dynamic_time_horizon,
+            time_varying_re=time_varying_re,
+        )
 
     # Calculate reference AGWP for 1 kg of CO2
     radiative_forcing_co2 = prospective_characterize_co2(
@@ -533,6 +556,7 @@ def _characterize_prospective_radiative_forcing(
     row,
     time_horizon,
     time_varying_re: bool = False,
+    indirect_effects: str = "all",
 ) -> CharacterizedRow:
     """
     Calculate prospective radiative forcing using Watanabe et al. (2026) characterization.
@@ -550,6 +574,14 @@ def _characterize_prospective_radiative_forcing(
         prospective_characterize_n2o,
     )
 
+    if char_func == prospective_characterize_ch4:
+        # CH4 supports multiple treatments of indirect effects
+        return char_func(
+            row,
+            time_horizon,
+            time_varying_re=time_varying_re,
+            indirect_effects=indirect_effects,
+        )
     if char_func in prospective_functions:
         # Watanabe functions support time_varying_re
         return char_func(row, time_horizon, time_varying_re=time_varying_re)

dynamic_characterization/prospective/agwp.py

@@ -13,8 +13,8 @@
 
 import numpy as np
 
-from .config import get_scenario
-from .data_loader import (
+from dynamic_characterization.prospective.config import get_scenario
+from dynamic_characterization.prospective.data_loader import (
     load_irf_ch4,
     load_irf_co2,
     load_irf_n2o,
@@ -46,6 +46,134 @@
 CONST_CH4 = 1.0 / PPB_TO_KG_CH4  # ppb/kg
 CONST_N2O = 1.0 / PPB_TO_KG_N2O  # ppb/kg
 
+CH4_INDIRECT_RE = 0.00014 + 0.00004
+
+
+def _check_indirect_effects_mode(mode: str) -> str:
+    """Validate the selected mode for indirect effects."""
+    valid_modes = {"none", "no_carbon_cycle", "all"}
+    if mode not in valid_modes:
+        raise ValueError(
+            f"Invalid indirect effect mode '{mode}'. Valid modes: {sorted(valid_modes)}"
+        )
+    return mode
+
+
+def _get_re_value(
+    re_series: np.ndarray,
+    year_idx: int,
+    step: int,
+    time_varying_re: bool,
+) -> float:
+    """Configuration-aware radiate efficiency getter."""
+    if time_varying_re:
+        # RE evolves: use RE at emission_year + t
+        re_idx = min(year_idx + step, len(re_series) - 1)
+        re_t = re_series[re_idx]
+    else:
+        re_t = re_series[year_idx]
+    return re_t
+
+
+def _agwp_ch4_cumulative(
+    emission_year: int,
+    time_horizon: int = 100,
+    time_varying_re: bool = False,
+    indirect_effects: str = "all",
+) -> np.ndarray:
+    """
+    Return cumulative CH4 AGWP values for each year of the selected horizon.
+
+    Based on Watanabe's Code #4 CH4 AGWP function.
+
+    Parameters
+    ----------
+    emission_year : int
+        Year of emission (2030-2100, clamped if outside)
+    time_horizon : int
+        Integration period in years
+    time_varying_re : bool
+        If True, use RE that evolves over the decay period.
+    indirect_effects : str
+        Indirect effects to include. Defaults to "all".
+
+    Returns
+    -------
+    np.ndarray
+        Cumulative AGWP in W*yr/m^2/kg
+    """
+    indirect_effects = _check_indirect_effects_mode(indirect_effects)
+
+    scenario = get_scenario()
+    iam, ssp, rcp = scenario["iam"], scenario["ssp"], scenario["rcp"]
+
+    # Load data
+    re_data_ch4 = load_re_ch4(iam)
+    re_data_co2 = load_re_co2(iam)
+    irf_series = load_irf_ch4()
+
+    years = re_data_ch4["_years"]
+    re_series_ch4 = re_data_ch4[ssp][rcp]
+    re_series_co2 = re_data_co2[ssp][rcp]
+
+    year_idx = _get_year_index(emission_year, years)
+    max_years = min(time_horizon, len(irf_series))
+
+    if max_years <= 0:
+        # Year 0 -> no effect
+        return np.array([0], dtype="float64")
+
+    total_rf = np.zeros(max_years, dtype="float64")
+
+    temp_s = np.zeros(max_years, dtype="float64")
+    temp_d = np.zeros(max_years, dtype="float64")
+    c1 = np.zeros(max_years, dtype="float64")
+    c2 = np.zeros(max_years, dtype="float64")
+    c3 = np.zeros(max_years, dtype="float64")
+
+    for t in range(max_years):
+        re_ch4 = _get_re_value(re_series_ch4, year_idx, t, time_varying_re)
+        re_co2 = _get_re_value(re_series_co2, year_idx, t, time_varying_re)
+
+        # Apply indirect effects if requested
+        re_ch4 = (
+            re_ch4 + CH4_INDIRECT_RE
+            if indirect_effects in {"no_carbon_cycle", "all"}
+            else re_ch4
+        )
+        factor1 = re_ch4 * irf_series[t] * CONST_CH4
+        factor1_model = factor1 * 1e12
+
+        factor2 = 0.0
+        if indirect_effects == "all" and t > 0:
+            temp_s_prev = temp_s[t - 1] if t > 0 else 0
+            temp_d_prev = temp_d[t - 1] if t > 0 else 0
+            c1_prev = c1[t - 1] if t > 0 else 0
+
+            temp_s[t] = (
+                temp_s_prev
+                + factor1_model / 7.7
+                - 1.31 * temp_s_prev / 7.7
+                - 0.88 * (temp_s_prev - temp_d_prev) / 7.7
+            )
+            temp_d[t] = temp_d_prev + (0.88 / 1.03) * (temp_s_prev - temp_d_prev) / 147
+
+            c1[t] = temp_s[t] * 11.06 * 0.6368 + c1_prev * np.exp(-1 / 2.376)
+            c2[t] = temp_s[t] * 11.06 * 0.3322 + c1_prev * np.exp(-1 / 30.14)
+            c3[t] = temp_s[t] * 11.06 * 0.031 + c1_prev * np.exp(-1 / 490.1)
+
+            factor2 = (c1[t] + c2[t] + c3[t]) * re_co2 * CONST_CO2
+
+        total_rf[t] = factor1 + factor2
+
+    cumulative_agwp = np.zeros(max_years, dtype="float64")
+    for t in range(1, max_years):
+        cumulative_agwp[t] = (
+            cumulative_agwp[t - 1] + (total_rf[t - 1] + total_rf[t]) / 2
+        )
+
+    return cumulative_agwp
+
 
 def _get_year_index(emission_year: int, years: np.ndarray) -> int:
     """
@@ -119,14 +247,7 @@ def agwp_co2(
     for t in range(max_years):
         irf_t = irf_series[t]
 
-        if time_varying_re:
-            # RE evolves: use RE at emission_year + t
-            re_idx = min(year_idx + t, len(re_series) - 1)
-            re_t = re_series[re_idx]
-        else:
-            # Fixed RE from emission year
-            re_t = re_series[year_idx]
-
+        re_t = _get_re_value(re_series, year_idx, t, time_varying_re)
         # AGWP contribution: RE * IRF * conversion * dt (dt=1 year)
         agwp += re_t * irf_t * CONST_CO2
 
@@ -137,6 +258,7 @@ def agwp_ch4(
     emission_year: int,
     time_horizon: int = 100,
     time_varying_re: bool = False,
+    indirect_effects: str = "all",
 ) -> float:
     """
     Calculate AGWP for 1 kg CH4.
@@ -149,39 +271,21 @@ def agwp_ch4(
         Integration period in years
     time_varying_re : bool
         If True, use RE that evolves over the decay period.
+    indirect_effects : str
+        Indirect effects to include. Defaults to "all".
 
     Returns
     -------
     float
         AGWP in W*yr/m^2/kg
     """
-    scenario = get_scenario()
-    iam, ssp, rcp = scenario["iam"], scenario["ssp"], scenario["rcp"]
-
-    # Load data
-    re_data = load_re_ch4(iam)
-    irf_series = load_irf_ch4()
-
-    years = re_data["_years"]
-    re_series = re_data[ssp][rcp]  # W/m^2/ppb
-
-    year_idx = _get_year_index(emission_year, years)
-
-    max_years = min(time_horizon, len(irf_series))
-
-    agwp = 0.0
-    for t in range(max_years):
-        irf_t = irf_series[t]
-
-        if time_varying_re:
-            re_idx = min(year_idx + t, len(re_series) - 1)
-            re_t = re_series[re_idx]
-        else:
-            re_t = re_series[year_idx]
-
-        agwp += re_t * irf_t * CONST_CH4
-
-    return agwp
+    cumulative_agwp = _agwp_ch4_cumulative(
+        emission_year=emission_year,
+        time_horizon=time_horizon,
+        time_varying_re=time_varying_re,
+        indirect_effects=indirect_effects,
+    )
+    return float(cumulative_agwp[-1])
 
 
 def agwp_n2o(
@@ -224,12 +328,7 @@ def agwp_n2o(
     for t in range(max_years):
         irf_t = irf_series[t]
 
-        if time_varying_re:
-            re_idx = min(year_idx + t, len(re_series) - 1)
-            re_t = re_series[re_idx]
-        else:
-            re_t = re_series[year_idx]
-
+        re_t = _get_re_value(re_series, year_idx, t, time_varying_re)
         agwp += re_t * irf_t * CONST_N2O
 
     return agwp