hip-ra-x-result.json [WIP]

.gitignore

@@ -88,6 +88,7 @@ nosetests.xml
 *.iml
 *.komodoproject
 .idea
+.junie
 .project
 .pydevproject
 .vscode
@@ -124,6 +125,9 @@ src/geophires_docs/fervo_project_red-2026_graph-data-extraction.xcf
 /Useful sites for Sphinx docstrings.txt
 /.github/workflows/workflows.7z
 tmp.patch
+project-structure.txt
+geophires-aliases.sh
+*message.txt
 
 # Mypy Cache
 .mypy_cache/

README.rst

@@ -384,6 +384,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

@@ -2901,6 +2901,18 @@ 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_musd = ((
+                                                 model.sdacgteconomics.OPEX.value + model.sdacgteconomics.storage.value + model.sdacgteconomics.transport.value)
+                                     * avg_carbon_extracted_tonnes) / 1_000_000.0
+            self.Coam.value += sdac_annual_opex_musd
+
         self.calculate_cashflow(model)
 
         # Calculate more financial values using numpy financials
@@ -3723,6 +3735,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/EconomicsS_DAC_GT.py

@@ -219,6 +219,31 @@ def __init__(self, model: Model):
             ErrMessage="assume default Percent Energy Devoted To Process (50%)",
             ToolTipText="Percent Energy Devoted To Process (%)"
         )
+        self.carbon_credit_price = self.ParameterDict[self.carbon_credit_price.Name] = floatParameter(
+            "S-DAC-GT Carbon Credit Price",
+            value=180.0,
+            DefaultValue=180.0,
+            Min=0.0,
+            Max=1000.0,
+            UnitType=Units.COSTPERMASS,
+            PreferredUnits=CostPerMassUnit.DOLLARSPERTONNE,
+            CurrentUnits=CostPerMassUnit.DOLLARSPERTONNE,
+            ErrMessage="assume default Carbon Credit Price (180 USD per tonne CO2)",
+            ToolTipText="Carbon Credit or Market Price (USD per tonne CO2)"
+        )
+
+        self.carbon_credit_duration = self.ParameterDict[self.carbon_credit_duration.Name] = floatParameter(
+            "S-DAC-GT Carbon Credit Duration",
+            value=12.0,
+            DefaultValue=12.0,
+            Min=0.0,
+            Max=100.0,
+            UnitType=Units.TIME,
+            PreferredUnits=TimeUnit.YEAR,
+            CurrentUnits=TimeUnit.YEAR,
+            ErrMessage="assume default Carbon Credit Duration (12 years)",
+            ToolTipText="Duration for which the carbon credit can be claimed (e.g., 12 years for US 45Q)"
+        )
 
         # local variable initiation
         # Capital Recovery Rate or Fixed Charge Factor - set initially for definitions
@@ -335,6 +360,21 @@ def __init__(self, model: Model):
             PreferredUnits=CostPerMassUnit.DOLLARSPERTONNE,
             CurrentUnits=CostPerMassUnit.DOLLARSPERTONNE
         )
+        self.CarbonRevenue = OutputParameter(
+            Name="Annual Carbon Revenue",
+            UnitType=Units.CURRENCYFREQUENCY,
+            PreferredUnits=CurrencyFrequencyUnit.DOLLARSPERYEAR,
+            CurrentUnits=CurrencyFrequencyUnit.DOLLARSPERYEAR
+        )
+        self.OutputParameterDict[self.CarbonRevenue.Name] = self.CarbonRevenue
+
+        self.CarbonCummCashFlow = OutputParameter(
+            Name="Cumulative Carbon Revenue",
+            UnitType=Units.CURRENCY,
+            PreferredUnits=CurrencyUnit.DOLLARS,
+            CurrentUnits=CurrencyUnit.DOLLARS
+        )
+        self.OutputParameterDict[self.CarbonCummCashFlow.Name] = self.CarbonCummCashFlow
 
         model.logger.info(f"Complete {str(__class__)}: {sys._getframe().f_code.co_name}")
 
@@ -483,6 +523,13 @@ def range_check(self) -> tuple:
                                                                                                   storage_max)
             return True, error_message
 
+        if not (self.carbon_credit_duration.Min
+                <= self.carbon_credit_duration.value
+                <= self.carbon_credit_duration.Max):
+            error_message = "S-DAC-GT ERROR: Carbon Credit Duration should be between {} and {}".format(
+                self.carbon_credit_duration.Min, self.carbon_credit_duration.Max)
+            return True, error_message
+
         return False, ""
 
     def geo_therm_cost(self, power_cost: float, CAPEX_mult: float, OPEX_mult: float, depth: float,
@@ -599,12 +646,13 @@ def Calculate(self, model: Model) -> None:
         # Convert from $/McF to $/kWh_th, but don't change any parameters value directly - it will throw off the rehydration
         NG_price = self.NG_price.value / self.NG_EnergyDensity.value
         NG_totalcost = self.therm.value * NG_price
-        self.LCOH.value, self.kWh_e_per_kWh_th.value = self.geo_therm_cost(model.surfaceplant.electricity_cost_to_buy.value,
-                                                                           self.CAPEX_mult.value, self.OPEX_mult.value,
-                                                                           model.reserv.depth.value * 3280.84,
-                                                                           np.average(model.wellbores.ProducedTemperature.value),
-                                                                           model.wellbores.Tinj.value,
-                                                                           model.wellbores.nprod.value * model.wellbores.prodwellflowrate.value)
+        self.LCOH.value, self.kWh_e_per_kWh_th.value = self.geo_therm_cost(
+            model.surfaceplant.electricity_cost_to_buy.value,
+            self.CAPEX_mult.value, self.OPEX_mult.value,
+            model.reserv.depth.value * 3280.84,
+            np.average(model.wellbores.ProducedTemperature.value),
+            model.wellbores.Tinj.value,
+            model.wellbores.nprod.value * model.wellbores.prodwellflowrate.value)
         geothermal_totalcost = self.LCOH.value * self.therm.value
         co2_power = self.elec.value / 1000 * self.power_co2intensity.value
         co2_elec_heat = self.therm.value / 1000 * self.power_co2intensity.value
@@ -621,7 +669,8 @@ def Calculate(self, model: Model) -> None:
 
         # calculate the net impact of S-DAC-GT on the annual production of the model
         avg_first_law_eff = np.average(model.surfaceplant.FirstLawEfficiency.value)
-        self.tot_heat_energy_consumed_per_tonne.value = (self.elec.value / avg_first_law_eff) + self.therm.value  # kWh_th/tonne
+        self.tot_heat_energy_consumed_per_tonne.value = (
+                                                                self.elec.value / avg_first_law_eff) + self.therm.value  # kWh_th/tonne
         self.tot_cost_per_tonne.value = CAPEX + self.OPEX.value + self.storage.value + self.transport.value  # USD/tonne
         self.percent_thermal_energy_going_to_heat.value = self.therm.value / self.tot_heat_energy_consumed_per_tonne.value
 
@@ -637,18 +686,22 @@ def Calculate(self, model: Model) -> None:
         # That then gives us the revenue, since we have a carbon price model
         # We can also get annual cash flow from it.
         for i in range(0, model.surfaceplant.plant_lifetime.value, 1):
-            self.CarbonExtractedAnnually.value[i] = (self.EnergySplit.value * model.surfaceplant.HeatkWhExtracted.value[i]) / self.tot_heat_energy_consumed_per_tonne.value
+            self.CarbonExtractedAnnually.value[i] = (self.EnergySplit.value * model.surfaceplant.HeatkWhExtracted.value[
+                i]) / self.tot_heat_energy_consumed_per_tonne.value
             if i == 0:
                 self.S_DAC_GTCummCarbonExtracted.value[i] = self.CarbonExtractedAnnually.value[i]
             else:
-                self.S_DAC_GTCummCarbonExtracted.value[i] = self.S_DAC_GTCummCarbonExtracted.value[i - 1] + self.CarbonExtractedAnnually.value[i]
+                self.S_DAC_GTCummCarbonExtracted.value[i] = self.S_DAC_GTCummCarbonExtracted.value[i - 1] + \
+                                                            self.CarbonExtractedAnnually.value[i]
             self.CarbonExtractedTotal.value = self.CarbonExtractedTotal.value + self.CarbonExtractedAnnually.value[i]
             self.S_DAC_GTAnnualCost.value[i] = self.CarbonExtractedAnnually.value[i] * self.tot_cost_per_tonne.value
             if i == 0:
                 self.S_DAC_GTCummCashFlow.value[i] = self.S_DAC_GTAnnualCost.value[i]
             else:
-                self.S_DAC_GTCummCashFlow.value[i] = self.S_DAC_GTCummCashFlow.value[i - 1] + self.S_DAC_GTAnnualCost.value[i]
-            self.CummCostPerTonne.value[i] = self.S_DAC_GTCummCashFlow.value[i] / self.S_DAC_GTCummCarbonExtracted.value[i]
+                self.S_DAC_GTCummCashFlow.value[i] = self.S_DAC_GTCummCashFlow.value[i - 1] + \
+                                                     self.S_DAC_GTAnnualCost.value[i]
+            self.CummCostPerTonne.value[i] = self.S_DAC_GTCummCashFlow.value[i] / \
+                                             self.S_DAC_GTCummCarbonExtracted.value[i]
 
         # We need to update the heat and electricity generated because we have consumed
         # some (all) of it to do the capture, so when they get used in the final economic calculation (below),
@@ -657,28 +710,31 @@ def Calculate(self, model: Model) -> None:
             if model.surfaceplant.enduse_option.value is not EndUseOptions.HEAT:
                 # all these end-use options have an electricity generation component
                 model.surfaceplant.TotalkWhProduced.value[i] = model.surfaceplant.TotalkWhProduced.value[i] - (
-                    self.CarbonExtractedAnnually.value[i] * self.elec.value)
+                        self.CarbonExtractedAnnually.value[i] * self.elec.value)
                 model.surfaceplant.NetkWhProduced.value[i] = model.surfaceplant.NetkWhProduced.value[i] - (
-                    self.CarbonExtractedAnnually.value[i] * self.elec.value)
+                        self.CarbonExtractedAnnually.value[i] * self.elec.value)
                 if model.surfaceplant.enduse_option.value is not EndUseOptions.ELECTRICITY:
                     model.surfaceplant.HeatkWhProduced.value[i] = model.surfaceplant.HeatkWhProduced.value[i] - (
-                        self.CarbonExtractedAnnually.value[i] * self.therm.value)
+                            self.CarbonExtractedAnnually.value[i] * self.therm.value)
             else:
                 # all the end-use option of direct-use only component
                 model.surfaceplant.HeatkWhProduced.value[i] = (model.surfaceplant.HeatkWhProduced.value[i] -
-                                                               (self.CarbonExtractedAnnually.value[i] * self.therm.value))
-
-        # FIXME TODO https://github.com/NREL/GEOPHIRES-X/issues/341?title=S-DAC+does+not+calculate+carbon+revenue
-        # Build a revenue generation model for the carbon capture, assuming the capture is being sequestered and that
-        # there is some sort of credit involved for doing that sequestering
-        # note that there may already be values in the CarbonRevenue array, so we need to
-        # add to them, not just set them. If there isn't values, there, the array will be filed with zeros, so adding won't be a problem
-        #total_duration = model.surfaceplant.plant_lifetime.value
-        #for i in range(0, total_duration, 1):
-        #    model.sdacgteconomics.CarbonRevenue.value[i] = (model.sdacgteconomics.CarbonRevenue.value[i] +
-        #                                              (self.CarbonExtractedAnnually.value[i] * model.economics.CarbonPrice.value[i]))
-        #   if i > 0:
-        #       model.economics.CarbonCummCashFlow.value[i] = model.economics.CarbonCummCashFlow.value[i - 1] + model.economics.CarbonRevenue.value[i]
+                                                               (self.CarbonExtractedAnnually.value[
+                                                                    i] * self.therm.value))
+
+        # Calculate Carbon Revenue based on S-DAC-GT specific credit price and duration
+        self.CarbonRevenue.value = [0.0] * model.surfaceplant.plant_lifetime.value
+        self.CarbonCummCashFlow.value = [0.0] * model.surfaceplant.plant_lifetime.value
+
+        for i in range(0, model.surfaceplant.plant_lifetime.value, 1):
+            # Enforce the parameterized statutory limit for tax credits
+            applicable_price = self.carbon_credit_price.value if i < self.carbon_credit_duration.value else 0.0
 
+            self.CarbonRevenue.value[i] = self.CarbonExtractedAnnually.value[i] * applicable_price
+            if i == 0:
+                self.CarbonCummCashFlow.value[i] = self.CarbonRevenue.value[i]
+            else:
+                self.CarbonCummCashFlow.value[i] = self.CarbonCummCashFlow.value[i - 1] + \
+                                                   self.CarbonRevenue.value[i]
         self._calculate_derived_outputs(model)
-        model.logger.info(f'Complete {str(__class__)}: {sys._getframe().f_code.co_name}')
+        model.logger.info(f'Complete {str(__class__)}: {sys._getframe().f_code.co_name}')

src/geophires_x/EconomicsSam.py

@@ -678,15 +678,59 @@ 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))]
 
-    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
+        )
+
+        # 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)
 
-        ret['ptc_fed_term'] = econ.PTCDuration.quantity().to(convertible_unit('yr')).magnitude
+        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
 
-        if econ.PTCInflationAdjusted.value:
-            ret['ptc_fed_escal'] = _pct(econ.RINFL)
+    # 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)