Calculate S-DAC carbon revenue (addresses NatLabRockies#341)

Author: softwareengineerprogrammer

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,19 @@ def __init__(self, model: Model):
             ErrMessage="assume default Percent Energy Devoted To Process (50%)",
             ToolTipText="Percent Energy Devoted To Process (%)"
         )
+        self.carbon_credit_price = 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.ParameterDict[self.carbon_credit_price.Name] = self.carbon_credit_price
 
         # local variable initiation
         # Capital Recovery Rate or Fixed Charge Factor - set initially for definitions
@@ -335,6 +348,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}")
 
@@ -599,12 +627,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 +650,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 +667,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 +691,28 @@ 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
+        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):
+            self.CarbonRevenue.value[i] = self.CarbonExtractedAnnually.value[i] * self.carbon_credit_price.value
+            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/OutputsS_DAC_GT.py

@@ -81,21 +81,26 @@ def PrintOutputs(self, model) -> tuple:
                     model.sdacgteconomics.S_DAC_GTCummCashFlow.value
                 sdac_df[f'Cum. Cost Per Tonne ({model.sdacgteconomics.CummCostPerTonne.PreferredUnits.value})|:,.2f'] = \
                     model.sdacgteconomics.CummCostPerTonne.value
+                sdac_df[f'Carbon Revenue ({model.sdacgteconomics.CarbonRevenue.PreferredUnits.value})|:,.2f'] = \
+                    model.sdacgteconomics.CarbonRevenue.value
+                sdac_df[f'Cum. Carbon Revenue ({model.sdacgteconomics.CarbonCummCashFlow.PreferredUnits.value})|:,.2f'] = \
+                    model.sdacgteconomics.CarbonCummCashFlow.value
 
                 f.write(NL)
                 f.write("                **********************" + NL)
-                f.write("                *  S-DAC-GT PROFILE  *" + NL)
+                f.write("                * S-DAC-GT PROFILE  *" + NL)
                 f.write("                **********************" + NL)
-                f.write("Year       Carbon      Cumm. Carbon     S-DAC-GT           S-DAC-GT Cumm.      Cumm. Cost" + NL)
-                f.write("Since      Captured     Captured       Annual Cost          Cash Flow        Cost Per Tonne" + NL)
+                f.write("Year       Carbon      Cumm. Carbon     S-DAC-GT           S-DAC-GT Cumm.      Cumm. Cost        Annual Carbon" + NL)
+                f.write("Since      Captured     Captured       Annual Cost          Cash Flow        Cost Per Tonne         Revenue" + NL)
                 f.write("Start     ("+model.sdacgteconomics.CarbonExtractedAnnually.PreferredUnits.value +
                                               ")   ("+model.sdacgteconomics.S_DAC_GTCummCarbonExtracted.PreferredUnits.value +
                                                              ")          ("+model.sdacgteconomics.S_DAC_GTAnnualCost.PreferredUnits.value +
                                              ")               ("+model.sdacgteconomics.S_DAC_GTCummCashFlow.PreferredUnits.value +
-                                             ")           ("+model.sdacgteconomics.CummCostPerTonne.PreferredUnits.value + ")" +NL)
+                                             ")           ("+model.sdacgteconomics.CummCostPerTonne.PreferredUnits.value +
+                                             ")         (" + model.sdacgteconomics.CarbonRevenue.PreferredUnits.value + ")" + NL)
                 i = 0
                 for i in range(0, model.surfaceplant.plant_lifetime.value, 1):
-                    f.write(f"   {i+1:3.0f}    {model.sdacgteconomics.CarbonExtractedAnnually.value[i]:,.2f}   {model.sdacgteconomics.S_DAC_GTCummCarbonExtracted.value[i]:,.2f}    {model.sdacgteconomics.S_DAC_GTAnnualCost.value[i]:,.2f}         {model.sdacgteconomics.S_DAC_GTCummCashFlow.value[i]:,.2f}         {model.sdacgteconomics.CummCostPerTonne.value[i]:.2f}" + NL)
+                    f.write(f"   {i+1:3.0f}    {model.sdacgteconomics.CarbonExtractedAnnually.value[i]:,.2f}   {model.sdacgteconomics.S_DAC_GTCummCarbonExtracted.value[i]:,.2f}    {model.sdacgteconomics.S_DAC_GTAnnualCost.value[i]:,.2f}         {model.sdacgteconomics.S_DAC_GTCummCashFlow.value[i]:,.2f}         {model.sdacgteconomics.CummCostPerTonne.value[i]:.2f}         {model.sdacgteconomics.CarbonRevenue.value[i]:,.2f}" + NL)
                     i = i + 1
 
         except BaseException as ex:
@@ -106,9 +111,9 @@ def PrintOutputs(self, model) -> tuple:
             print(msg)
             model.logger.critical(str(ex))
             model.logger.critical(msg)
-            raise RuntimeError(msg, e)
+            raise RuntimeError(msg, ex)
 
         model.logger.info(f'Complete {str(__class__)}: {__name__}')
 
         sdac_df = sdac_df.reset_index()
-        return sdac_df, sdac_results
+        return sdac_df, sdac_results