Calculate S-DAC carbon revenue (BICYCLE); Add SAM-EM S-DAC support [v3.15]

README.rst

@@ -392,6 +392,10 @@ Example-specific web interface deeplinks are listed in the Link column.
      - `example_SAM-single-owner-PPA-9_cooling.txt <tests/examples/example_SAM-single-owner-PPA-9_cooling.txt>`__
      - `.out <tests/examples/example_SAM-single-owner-PPA-9_cooling.out>`__
      - `link <https://gtp.scientificwebservices.com/geophires?geophires-example-id=example_SAM-single-owner-PPA-9_cooling>`__
+   * - SAM Single Owner PPA: S-DAC
+     - `S-DAC-GT-2.txt <tests/examples/S-DAC-GT-2.txt>`__
+     - `.out <tests/examples/S-DAC-GT-2.out>`__
+     - `link <https://gtp.scientificwebservices.com/geophires?geophires-example-id=S-DAC-GT-2>`__
 .. raw:: html
 
    <embed>

src/geophires_x/Economics.py

@@ -2942,6 +2942,27 @@ def Calculate(self, model: Model) -> None:
         if self.DoSDACGTCalculations.value:
             model.sdacgteconomics.Calculate(model)
 
+            # Consolidate S-DAC-GT CAPEX and OPEX into the main plant ledgers
+            max_carbon_capacity_tonnes = np.max(model.sdacgteconomics.CarbonExtractedAnnually.value)
+            sdac_overnight_capex_musd = (
+                                                    model.sdacgteconomics.CAPEX.value * model.sdacgteconomics.CAPEX_mult.value * max_carbon_capacity_tonnes) / 1_000_000.0
+            self.CCap.value += sdac_overnight_capex_musd
+
+            avg_carbon_extracted_tonnes = np.average(model.sdacgteconomics.CarbonExtractedAnnually.value)
+            sdac_annual_opex_usd = (
+                                               model.sdacgteconomics.OPEX.value + model.sdacgteconomics.storage.value + model.sdacgteconomics.transport.value) * avg_carbon_extracted_tonnes
+
+            if model.sdacgteconomics.sorbent_replacement_frequency.value > 0:
+                max_carbon_capacity_tonnes = np.max(model.sdacgteconomics.CarbonExtractedAnnually.value)
+                replacements_per_lifetime = int(
+                    model.surfaceplant.plant_lifetime.value / model.sdacgteconomics.sorbent_replacement_frequency.value)
+                annualized_replacement_usd = (
+                                                         max_carbon_capacity_tonnes * model.sdacgteconomics.sorbent_replacement_cost.value * replacements_per_lifetime) / model.surfaceplant.plant_lifetime.value
+                sdac_annual_opex_usd += annualized_replacement_usd
+
+            sdac_annual_opex_musd = sdac_annual_opex_usd / 1_000_000.0
+            self.Coam.value += sdac_annual_opex_musd
+
         self.calculate_cashflow(model)
 
         # Calculate more financial values using numpy financials
@@ -3822,6 +3843,13 @@ def calculate_cashflow(self, model: Model) -> None:
                     self.TotalRevenue.value[i] = self.TotalRevenue.value[i] + self.CarbonRevenue.value[i]
                     #self.TotalCummRevenue.value[i] = self.TotalCummRevenue.value[i] + self.CarbonCummCashFlow.value[i]
 
+            if self.DoSDACGTCalculations.value:
+                for i in range(model.surfaceplant.construction_years.value,
+                               model.surfaceplant.plant_lifetime.value + model.surfaceplant.construction_years.value,
+                               1):
+                    sdac_index = i - model.surfaceplant.construction_years.value
+                    self.TotalRevenue.value[i] += (model.sdacgteconomics.CarbonRevenue.value[sdac_index] / 1_000_000.0)
+
             # for the sake of display, insert zeros at the beginning of the pricing arrays
             for i in range(0, model.surfaceplant.construction_years.value, 1):
                 self.ElecPrice.value.insert(0, 0.0)

src/geophires_x/EconomicsSam.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import copy
 import json
 import logging
 import os
@@ -287,6 +288,11 @@ def sf(_v: float, num_sig_figs: int = 5) -> float:
     sam_economics.nominal_discount_rate.value, sam_economics.wacc.value = _calculate_nominal_discount_rate_and_wacc_pct(
         model, single_owner
     )
+
+    # noinspection SpellCheckingInspection
+    if hasattr(model.economics, 'DoSDACGTCalculations') and model.economics.DoSDACGTCalculations.value:
+        sam_economics.s_dac_carbon_extracted_annually = copy.deepcopy(model.sdacgteconomics.CarbonExtractedAnnually)
+
     sam_economics.moic.value = _calculate_moic(sam_economics.sam_cash_flow_profile, model)
     sam_economics.project_vir.value = _calculate_project_vir(sam_economics.sam_cash_flow_profile, model)
     sam_economics.project_payback_period.value = _calculate_project_payback_period(
@@ -598,6 +604,13 @@ def _get_utility_rate_parameters(m: Model) -> dict[str, Any]:
 
     max_total_kWh_produced = np.max(m.surfaceplant.TotalkWhProduced.quantity().to(convertible_unit('kWh')).magnitude)
 
+    if econ.DoSDACGTCalculations.value:
+        # Restore the true gross maximum before S-DAC in-place mutation decremented it
+        sdac_elec_consumption_kwh = (
+            np.max(m.sdacgteconomics.CarbonExtractedAnnually.value) * m.sdacgteconomics.elec.value
+        )
+        max_total_kWh_produced += sdac_elec_consumption_kwh
+
     net_kwh_produced_series: Iterable | float | int = (
         m.surfaceplant.NetkWhProduced.quantity().to(convertible_unit('kWh')).magnitude
     )
@@ -661,8 +674,36 @@ def _get_single_owner_parameters(model: Model) -> dict[str, Any]:
     royalty_supplemental_payments_by_year_usd = econ.get_royalty_supplemental_payments_schedule_usd(model)[
         _pre_revenue_years_count(model) :
     ]
+
+    sdac_average_opex_usd = 0.0
+    sdac_opex_by_year_usd = [0.0] * model.surfaceplant.plant_lifetime.value
+    if econ.DoSDACGTCalculations.value:
+        sdac_opex_by_year_usd = model.sdacgteconomics.AnnualOPEX_USD
+        avg_carbon_extracted_tonnes = np.average(model.sdacgteconomics.CarbonExtractedAnnually.value)
+        sdac_average_opex_usd = (
+            model.sdacgteconomics.OPEX.value
+            + model.sdacgteconomics.storage.value
+            + model.sdacgteconomics.transport.value
+        ) * avg_carbon_extracted_tonnes
+        if model.sdacgteconomics.sorbent_replacement_frequency.value > 0:
+            max_carbon_capacity_tonnes = np.max(model.sdacgteconomics.CarbonExtractedAnnually.value)
+            replacements_per_lifetime = int(
+                model.surfaceplant.plant_lifetime.value / model.sdacgteconomics.sorbent_replacement_frequency.value
+            )
+            annualized_replacement_usd = (
+                max_carbon_capacity_tonnes
+                * model.sdacgteconomics.sorbent_replacement_cost.value
+                * replacements_per_lifetime
+            ) / model.surfaceplant.plant_lifetime.value
+            sdac_average_opex_usd += annualized_replacement_usd
+
+    opex_base_without_sdac_usd = opex_base_usd - sdac_average_opex_usd
     for year_index in range(model.surfaceplant.plant_lifetime.value):
-        opex_by_year_usd.append(opex_base_usd + royalty_supplemental_payments_by_year_usd[year_index])
+        opex_by_year_usd.append(
+            opex_base_without_sdac_usd
+            + royalty_supplemental_payments_by_year_usd[year_index]
+            + sdac_opex_by_year_usd[year_index]
+        )
 
     if model.surfaceplant.enduse_option.value == EndUseOptions.HEAT:
         # For pure direct-use, pumping is a grid purchase.
@@ -688,18 +729,62 @@ def _get_single_owner_parameters(model: Model) -> dict[str, Any]:
     ret['federal_tax_rate'], ret['state_tax_rate'] = _get_fed_and_state_tax_rates(econ.CTR.value)
 
     geophires_itc_tenths = Decimal(econ.RITC.value)
-    ret['itc_fed_percent'] = [float(geophires_itc_tenths * Decimal(100))]
 
     geophires_state_itc_usd = Decimal(econ.ritc_state_amount.quantity().to(convertible_unit('USD')).magnitude)
     ret['itc_sta_amount'] = [float(geophires_state_itc_usd)]
 
-    if econ.PTCElec.Provided:
-        ret['ptc_fed_amount'] = [econ.PTCElec.quantity().to(convertible_unit('USD/kWh')).magnitude]
+    if econ.DoSDACGTCalculations.value and geophires_itc_tenths > 0:
+        # Shield DAC CAPEX from ITC to prevent unlawful MACRS basis reduction
+        max_carbon_capacity_tonnes = max(model.sdacgteconomics.CarbonExtractedAnnually.value)
+        sdac_capex_usd = (
+            model.sdacgteconomics.CAPEX.value * model.sdacgteconomics.CAPEX_mult.value * max_carbon_capacity_tonnes
+        )
 
-        ret['ptc_fed_term'] = econ.PTCDuration.quantity().to(convertible_unit('yr')).magnitude
+        # Geothermal eligible basis = Total Installed Cost - DAC CAPEX
+        eligible_geothermal_basis_usd = ret['total_installed_cost'] - sdac_capex_usd
+        itc_fed_amount_usd = float(Decimal(eligible_geothermal_basis_usd) * geophires_itc_tenths)
 
-        if econ.PTCInflationAdjusted.value:
-            ret['ptc_fed_escal'] = _pct(econ.RINFL)
+        ret['itc_fed_percent'] = [0.0]  # Disable percentage-based ITC on the whole project
+        ret['itc_fed_amount'] = [itc_fed_amount_usd]  # Inject fixed ITC amount for geothermal only
+    else:
+        ret['itc_fed_percent'] = [float(geophires_itc_tenths * Decimal(100))]
+
+    # Build a year-by-year schedule for the Federal Production Tax Credit
+    ptc_fed_amount_schedule = [0.0] * model.surfaceplant.plant_lifetime.value
+
+    # 1. Base Electricity PTC
+    if econ.PTCElec.Provided:
+        base_ptc_rate = econ.PTCElec.quantity().to(convertible_unit('USD/kWh')).magnitude
+        ptc_term = int(econ.PTCDuration.quantity().to(convertible_unit('yr')).magnitude)
+        for i in range(min(ptc_term, model.surfaceplant.plant_lifetime.value)):
+            escalation = (1.0 + _pct(econ.RINFL) / 100.0) ** i if econ.PTCInflationAdjusted.value else 1.0
+            ptc_fed_amount_schedule[i] = base_ptc_rate * escalation
+
+    # 2. S-DAC 45Q Equivalent (Mapped as Non-Taxable PBI to avoid MACRS/Exclusivity issues)
+    if econ.DoSDACGTCalculations.value:
+        pbi_oth_amount_schedule = [0.0] * model.surfaceplant.plant_lifetime.value
+
+        for i in range(model.surfaceplant.plant_lifetime.value):
+            # The statutory duration cutoff (e.g., 12 years) is already handled inside EconomicsS_DAC_GT.py
+            # so CarbonRevenue will naturally be 0.0 for years > 12.
+            sdac_revenue_usd = model.sdacgteconomics.CarbonRevenue.value[i]
+            net_kwh_produced = model.surfaceplant.NetkWhProduced.value[i]
+
+            # Convert absolute S-DAC revenue (USD) into an equivalent USD/kWh PBI rate for SAM
+            if net_kwh_produced > 0:
+                pbi_oth_amount_schedule[i] = sdac_revenue_usd / net_kwh_produced
+
+        if any(rate > 0.0 for rate in pbi_oth_amount_schedule):
+            ret['pbi_oth_amount'] = pbi_oth_amount_schedule
+            ret['pbi_oth_term'] = model.surfaceplant.plant_lifetime.value
+            ret['pbi_oth_tax_fed'] = 0.0  # Strictly exclude from Federal taxable gross income
+            ret['pbi_oth_tax_sta'] = 0.0  # Strictly exclude from State taxable gross income
+
+    # Inject the combined array into SAM
+    if any(rate > 0.0 for rate in ptc_fed_amount_schedule):
+        ret['ptc_fed_amount'] = ptc_fed_amount_schedule
+        ret['ptc_fed_term'] = model.surfaceplant.plant_lifetime.value
+        ret['ptc_fed_escal'] = 0.0  # Escalation is already manually calculated in the array above
 
     # 'Property Tax Rate'
     geophires_ptr_tenths = Decimal(econ.PTR.value)

src/geophires_x/EconomicsSamCalculations.py

@@ -22,6 +22,7 @@
     investment_tax_credit_output_parameter,
     lcoh_output_parameter,
     lcoc_output_parameter,
+    carbon_extracted_annually_output_parameter,
 )
 from geophires_x.GeoPHIRESUtils import is_float, quantity, is_int
 from geophires_x.Parameter import OutputParameter
@@ -83,6 +84,9 @@ class SamEconomicsCalculations:
     investment_tax_credit: OutputParameter = field(default_factory=investment_tax_credit_output_parameter)
 
     capacity_payment_revenue_sources: list[CapacityPaymentRevenueSource] = field(default_factory=list)
+    s_dac_carbon_extracted_annually: OutputParameter | None = (
+        None  # field(default_factory=carbon_extracted_annually_output_parameter)
+    )
 
     @property
     def _pre_revenue_years_count(self) -> int:
@@ -169,6 +173,8 @@ def _get_row(row_name__: str) -> list[Any]:
 
         ret = self._insert_capacity_payment_line_items(ret)
 
+        ret = self._insert_s_dac_line_items(ret)
+
         ret = self._insert_calculated_levelized_metrics_line_items(ret)
 
         if self._may_consume_grid_electricity:
@@ -261,6 +267,42 @@ def _for_operational_years(_row: list[Any]) -> list[Any]:
 
         return ret
 
+    def _insert_s_dac_line_items(self, cf_ret: list[list[Any]]) -> list[list[Any]]:
+        ret: list[list[Any]] = cf_ret.copy()
+
+        if self.s_dac_carbon_extracted_annually is None:
+            return ret
+
+        def _get_row_index(row_name_: str) -> int:
+            return [it[0] for it in ret].index(row_name_)
+
+        def _insert_row_before(before_row_name: str, row_name: str, row_content: list[Any]) -> None:
+            ret.insert(
+                _get_row_index(before_row_name),
+                [row_name, *row_content],
+            )
+
+        def _insert_blank_line_before(before_row_name: str) -> None:
+            _insert_row_before(before_row_name, '', ['' for _it in ret[_get_row_index(before_row_name)]][1:])
+
+        REVENUE_CATEGORY_ROW_NAME = 'REVENUE'
+        # ENERGY_CATEGORY_ROW_NAME = 'ENERGY'
+
+        def _for_operational_years(_row: list[Any]) -> list[Any]:
+            return [*([''] * (self._pre_revenue_years_count - 1)), 0, *_row]
+
+        line_item_display_name = self.s_dac_carbon_extracted_annually.Name.replace(
+            'Tonnes per Year CO2 extracted', 'S-DAC CO2 extracted'
+        )
+        _insert_row_before(
+            REVENUE_CATEGORY_ROW_NAME,
+            f'{line_item_display_name} ({self.s_dac_carbon_extracted_annually.CurrentUnits.value})',
+            _for_operational_years(self.s_dac_carbon_extracted_annually.value),
+        )
+        _insert_blank_line_before(REVENUE_CATEGORY_ROW_NAME)
+
+        return ret
+
     def _insert_calculated_levelized_metrics_line_items(self, cf_ret: list[list[Any]]) -> list[list[Any]]:
         ret = cf_ret.copy()
 

src/geophires_x/EconomicsUtils.py

@@ -1,7 +1,15 @@
 from __future__ import annotations
 
 from geophires_x.Parameter import OutputParameter
-from geophires_x.Units import Units, PercentUnit, TimeUnit, CurrencyUnit, CurrencyFrequencyUnit, EnergyCostUnit
+from geophires_x.Units import (
+    Units,
+    PercentUnit,
+    TimeUnit,
+    CurrencyUnit,
+    CurrencyFrequencyUnit,
+    EnergyCostUnit,
+    MassPerTimeUnit,
+)
 
 CONSTRUCTION_CAPEX_SCHEDULE_PARAMETER_NAME = 'Construction CAPEX Schedule'
 
@@ -245,6 +253,15 @@ def investment_tax_credit_output_parameter() -> OutputParameter:
     )
 
 
+def carbon_extracted_annually_output_parameter() -> OutputParameter:
+    return OutputParameter(
+        Name="Tonnes per Year CO2 extracted",
+        UnitType=Units.MASSPERTIME,
+        PreferredUnits=MassPerTimeUnit.TONNEPERYEAR,
+        CurrentUnits=MassPerTimeUnit.TONNEPERYEAR,
+    )
+
+
 def expand_schedule_dsl(schedule_strings: list[str | float], total_years: int) -> list[float]:
     """
     Deprecated, call ParameterUtils.expand_schedule_dsl

src/geophires_x/Outputs.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import copy
 import datetime
 import math
 import time
@@ -933,6 +934,32 @@ def o(output_param: OutputParameter):
             else:
                 return output_param
 
+        econ_CarbonPrice = copy.deepcopy(econ.CarbonPrice)
+        econ_CarbonRevenue = copy.deepcopy(econ.CarbonRevenue)
+        econ_CarbonCummCashFlow = copy.deepcopy(econ.CarbonCummCashFlow)
+
+        if econ.DoSDACGTCalculations.value:
+            econ_CarbonRevenue = copy.deepcopy(model.sdacgteconomics.CarbonRevenue)
+            econ_CarbonRevenue.value = [
+                *([0]*model.surfaceplant.construction_years.value),
+                *model.sdacgteconomics.CarbonRevenue.value
+            ]
+
+            def _convert(gt_param,  econ_param) -> None:
+                gt_param.value = gt_param.quantity().to(econ_param.CurrentUnits).magnitude
+                gt_param.CurrentUnits = econ_param.CurrentUnits
+
+            _convert(econ_CarbonRevenue, econ.CarbonRevenue)
+
+            econ_CarbonCummCashFlow = copy.deepcopy(model.sdacgteconomics.CarbonCummCashFlow)
+            econ_CarbonCummCashFlow.value = [
+                *([0] * model.surfaceplant.construction_years.value),
+                *model.sdacgteconomics.CarbonCummCashFlow.value
+            ]
+
+            _convert(econ_CarbonCummCashFlow, econ.CarbonCummCashFlow)
+
+
         f.write('Start    ('
                 + o(econ.ElecPrice).CurrentUnits.value +
                 ')(' + o(econ.ElecRevenue).CurrentUnits.value +
@@ -943,9 +970,9 @@ def o(output_param: OutputParameter):
                 ')   |(' + o(econ.CoolingPrice).CurrentUnits.value +
                 ') (' + o(econ.CoolingRevenue).CurrentUnits.value +
                 ')    (' + o(econ.CoolingCummRevenue).CurrentUnits.value +
-                ')    |(' + o(econ.CarbonPrice).CurrentUnits.value +
-                ')    (' + o(econ.CarbonRevenue).CurrentUnits.value +
-                ')    (' + o(econ.CarbonCummCashFlow).CurrentUnits.value +
+                ')    |(' + o(econ_CarbonPrice).CurrentUnits.value +
+                ')    (' + o(econ_CarbonRevenue).CurrentUnits.value +
+                ')    (' + o(econ_CarbonCummCashFlow).CurrentUnits.value +
                 ')    |(' + o(econ.Coam).CurrentUnits.value +
                 ') (' + o(econ.TotalRevenue).CurrentUnits.value +
                 ')    (' + o(econ.TotalCummRevenue).CurrentUnits.value + ')\n')
@@ -959,7 +986,7 @@ def o(output_param: OutputParameter):
             else:
                 opex = o(econ.Coam).value
             f.write(
-                f'{ii:3.0f}     {o(econ.ElecPrice).value[ii]:5.2f}          {o(econ.ElecRevenue).value[ii]:5.2f}  {o(econ.ElecCummRevenue).value[ii]:5.2f}     |   {o(econ.HeatPrice).value[ii]:5.2f}    {o(econ.HeatRevenue).value[ii]:5.2f}        {o(econ.HeatCummRevenue).value[ii]:5.2f}    |   {o(econ.CoolingPrice).value[ii]:5.2f}    {o(econ.CoolingRevenue).value[ii]:5.2f}        {o(econ.CoolingCummRevenue).value[ii]:5.2f}     |   {o(econ.CarbonPrice).value[ii]:5.2f}    {o(econ.CarbonRevenue).value[ii]:5.2f}        {o(econ.CarbonCummCashFlow).value[ii]:5.2f}     | {opex:5.2f}     {o(econ.TotalRevenue).value[ii]:5.2f}     {o(econ.TotalCummRevenue).value[ii]:5.2f}\n')
+                f'{ii:3.0f}     {o(econ.ElecPrice).value[ii]:5.2f}          {o(econ.ElecRevenue).value[ii]:5.2f}  {o(econ.ElecCummRevenue).value[ii]:5.2f}     |   {o(econ.HeatPrice).value[ii]:5.2f}    {o(econ.HeatRevenue).value[ii]:5.2f}        {o(econ.HeatCummRevenue).value[ii]:5.2f}    |   {o(econ.CoolingPrice).value[ii]:5.2f}    {o(econ.CoolingRevenue).value[ii]:5.2f}        {o(econ.CoolingCummRevenue).value[ii]:5.2f}     |   {o(econ_CarbonPrice).value[ii]:5.2f}    {o(econ_CarbonRevenue).value[ii]:5.2f}        {o(econ_CarbonCummCashFlow).value[ii]:5.2f}     | {opex:5.2f}     {o(econ.TotalRevenue).value[ii]:5.2f}     {o(econ.TotalCummRevenue).value[ii]:5.2f}\n')
         f.write(NL)
 
     # noinspection PyMethodMayBeStatic