Added additional estimates from GR website

src/geophires_x/CylindricalReservoir.py

@@ -6,7 +6,7 @@
 import numpy as np
 from pint.facets.plain import PlainQuantity
 
-from geophires_x.GeoPHIRESUtils import density_water_kg_per_m3, static_pressure_MPa, quantity
+from geophires_x.GeoPHIRESUtils import density_water_kg_per_m3, lithostatic_pressure_MPa, quantity
 
 from geophires_x.GeoPHIRESUtils import heat_capacity_water_J_per_kg_per_K
 import geophires_x.Model as Model
@@ -43,6 +43,7 @@ def __init__(self, model: Model):
 
         self.InputDepth = self.ParameterDict[self.InputDepth.Name] = floatParameter(
             "Cylindrical Reservoir Input Depth",
+            value=3.0,
             DefaultValue=3.0,
             Min=0.1,
             Max=15,
@@ -56,6 +57,7 @@ def __init__(self, model: Model):
 
         self.OutputDepth = self.ParameterDict[self.OutputDepth.Name] = floatParameter(
             "Cylindrical Reservoir Output Depth",
+            value=self.InputDepth.value,
             DefaultValue=self.InputDepth.value,
             Min=0.1,
             Max=15,
@@ -68,6 +70,7 @@ def __init__(self, model: Model):
         )
         self.Length = self.ParameterDict[self.Length.Name] = floatParameter(
             "Cylindrical Reservoir Length",
+            value=4.0,
             DefaultValue=4.0,
             Min=0.1,
             Max=10.0,
@@ -80,6 +83,7 @@ def __init__(self, model: Model):
         )
         self.RadiusOfEffect = self.ParameterDict[self.RadiusOfEffect.Name] = floatParameter(
             "Cylindrical Reservoir Radius of Effect",
+            value=30.0,
             DefaultValue=30.0,
             Min=0,
             Max=1000.0,
@@ -92,13 +96,14 @@ def __init__(self, model: Model):
         )
         self.RadiusOfEffectFactor = self.ParameterDict[self.RadiusOfEffectFactor.Name] = floatParameter(
             "Cylindrical Reservoir Radius of Effect Factor",
+            value=1.0,
             DefaultValue=1.0,
             Min=0.0,
             Max=10.0,
             UnitType=Units.PERCENT,
             PreferredUnits=PercentUnit.TENTH,
             CurrentUnits=PercentUnit.TENTH,
-            ErrMessage="assume default cylindrical reservoir radius of effect reduction factor (0.1)",
+            ErrMessage="assume default cyclindrical reservoir radius of effect reduction factor (0.1)",
             ToolTipText="The radius of effect reduction factor - to account for the fact that we cannot extract 100%"
                         + " of the heat in the cylinder.",
         )
@@ -109,6 +114,7 @@ def __init__(self, model: Model):
         # internal values required for calculations
         self.depth = self.ParameterDict[self.depth.Name] = floatParameter(
             "Drilled length",
+            value=0.0,
             DefaultValue=0.0,
             Min=0.0,
             Max=150,
@@ -121,6 +127,7 @@ def __init__(self, model: Model):
         )
         self.waterloss = self.ParameterDict[self.waterloss.Name] = floatParameter(
             "Water Loss Fraction",
+            value=0.0,
             DefaultValue=0.0,
             Min=0.0,
             Max=0.99,
@@ -222,8 +229,8 @@ def Calculate(self, model: Model) -> None:
         # initialize with the Initial reservoir temperature
         self.Tresoutput.value = np.array(len(self.timevector.value) * [self.Trock.value])
         # depth in this case is actually the total length of the drilled assembly
-        # as the average of the InputDepth and the OutputDepth
-        self.depth.value = (self.InputDepth.value + self.OutputDepth.value) / 2.0
+        self.depth.value = ((self.InputDepth.quantity() + self.OutputDepth.quantity() + self.Length.quantity()
+             ).to(self.depth.CurrentUnits).magnitude)
 
         # Total volume of all laterals but hollow cylinder - doesn't include drilled-out area, units = m3
         self.resvolcalc.value = (
@@ -244,17 +251,19 @@ def Calculate(self, model: Model) -> None:
                                                  ) / 1e15  # 10^15 J
         self.cpwater.value = heat_capacity_water_J_per_kg_per_K(
             model.wellbores.Tinj.value * 0.5 + (self.Trock.value * 0.9 + model.wellbores.Tinj.value * 0.1) * 0.5,
-            pressure=self.hydrostatic_pressure()
+            pressure=model.reserv.lithostatic_pressure()
         )
         self.rhowater.value = density_water_kg_per_m3(
             model.wellbores.Tinj.value * 0.5 + (self.Trock.value * 0.9 + model.wellbores.Tinj.value * 0.1) * 0.5,
-            pressure=self.hydrostatic_pressure()
+            pressure=model.reserv.lithostatic_pressure()
         )
 
         model.logger.info(f'complete {str(__class__)}: {sys._getframe().f_code.co_name}')
 
     def lithostatic_pressure(self) -> PlainQuantity:
         """@override"""
 
+        Standard reservoir implementation uses depth but CylindricalReservoir sets depth to total drilled length
+        """
         return quantity(static_pressure_MPa(self.rhorock.quantity().to('kg/m**3').magnitude,
                                             self.InputDepth.quantity().to('m').magnitude), 'MPa')

src/geophires_x/Economics.py

@@ -556,7 +556,7 @@ def CalculateLCOELCOHLCOC(self, model: Model) -> tuple:
             NPVgrt = self.GTR.value / (1 - self.GTR.value) * (NPVcap + NPVoandm + NPVfc + NPVit - NPVitc)
             LCOH = (NPVcap + NPVoandm + NPVfc + NPVit + NPVgrt - NPVitc) / np.sum(
                 model.surfaceplant.HeatkWhProduced.value * inflationvector * discountvector) * 1E8
-            LCOH = self.LCOH.value * 2.931  # $/MMBTU
+            LCOH = LCOH * 2.931  # $/MMBTU
 
         elif model.surfaceplant.enduse_option.value == EndUseOptions.HEAT and model.surfaceplant.plant_type.value == PlantType.DISTRICT_HEATING:
             PumpingCosts = model.surfaceplant.PumpingkWh.value * model.surfaceplant.electricity_cost_to_buy.value / 1E6
@@ -1497,11 +1497,46 @@ def __init__(self, model: Model):
             self.jobs_created_per_MW_electricity.Name] = floatParameter(
             "Estimated Jobs Created per MW of Electricity Produced",
             DefaultValue=2.13,
-            UnitType=Units.NONE,
+            UnitType=Units.JOBS_PER_ENERGY,
+            PreferredUnits=JobsPerEnergyUnit.JOBSPERMW,
+            CurrentUnits=JobsPerEnergyUnit.JOBSPERMW,
             Required=False,
             ToolTipText="Estimated jobs created per MW of electricity produced, per https://geothermal.org/resources/geothermal-basics"
         )
 
+        self.property_tax_per_MW_electricity = self.ParameterDict[
+            self.property_tax_per_MW_electricity.Name] = floatParameter(
+            "Estimated Property Tax per MW of Electricity Produced",
+            DefaultValue=0.210,
+            UnitType=Units.ROYALTY_PER_ENERGY,
+            PreferredUnits=RoyaltyPerEnergyUnit.ROYALTYPERMW,
+            CurrentUnits=RoyaltyPerEnergyUnit.ROYALTYPERMW,
+            Required=False,
+            ToolTipText="Estimated property tax per MW of electricity produced, per https://geothermal.org/resources/geothermal-basics"
+        )
+
+        self.gov_royalty_per_MW_electricity = self.ParameterDict[
+            self.gov_royalty_per_MW_electricity.Name] = floatParameter(
+            "Estimated Governmental Royalty per MW of Electricity Produced",
+            DefaultValue=0.315,
+            UnitType=Units.ROYALTY_PER_ENERGY,
+            PreferredUnits=RoyaltyPerEnergyUnit.ROYALTYPERMW,
+            CurrentUnits=RoyaltyPerEnergyUnit.ROYALTYPERMW,
+            Required=False,
+            ToolTipText="Estimated Estimated Governmental Royalty per MW of electricity produced, per https://geothermal.org/resources/geothermal-basics"
+        )
+
+        self.total_royalty_per_MW_electricity = self.ParameterDict[
+            self.total_royalty_per_MW_electricity.Name] = floatParameter(
+            "Estimated Jobs Created per MW of Electricity Produced",
+            DefaultValue=0.420,
+            UnitType=Units.ROYALTY_PER_ENERGY,
+            PreferredUnits=RoyaltyPerEnergyUnit.ROYALTYPERMW,
+            CurrentUnits=RoyaltyPerEnergyUnit.ROYALTYPERMW,
+            Required=False,
+            ToolTipText="Estimated total royalty per MW of electricity produced, per https://geothermal.org/resources/geothermal-basics"
+        )
+
         # local variable initialization
         self.CAPEX_cost_electricity_plant = 0.0
         self.CAPEX_cost_heat_plant = 0.0
@@ -1813,7 +1848,27 @@ def __init__(self, model: Model):
         )
         self.jobs_created = self.OutputParameterDict[self.jobs_created.Name] = OutputParameter(
             Name="Estimated Jobs Created",
-            UnitType=Units.NONE,
+            UnitType=Units.JOBS,
+            CurrentUnits=JobsUnit.JOBS,
+            PreferredUnits=JobsUnit.JOBS
+        )
+        self.property_tax_created = self.OutputParameterDict[self.property_tax_created.Name] = OutputParameter(
+            Name="Estimated amount of property tax that will be paid",
+            UnitType=Units.CURRENCY,
+            CurrentUnits=CurrencyUnit.MDOLLARS,
+            PreferredUnits=CurrencyUnit.MDOLLARS
+        )
+        self.total_royalties_created = self.OutputParameterDict[self.total_royalties_created.Name] = OutputParameter(
+            Name="Estimated total royalties that will have to be paid",
+            UnitType=Units.CURRENCY,
+            CurrentUnits=CurrencyUnit.MDOLLARS,
+            PreferredUnits=CurrencyUnit.MDOLLARS
+        )
+        self.gov_royalties_created = self.OutputParameterDict[self.gov_royalties_created.Name] = OutputParameter(
+            Name="Estimated governmental royalties to be paid",
+            UnitType=Units.CURRENCY,
+            CurrentUnits=CurrencyUnit.MDOLLARS,
+            PreferredUnits=CurrencyUnit.MDOLLARS
         )
 
         model.logger.info(f'Complete {__class__!s}: {sys._getframe().f_code.co_name}')
@@ -2892,6 +2947,21 @@ def Calculate(self, model: Model) -> None:
             np.average(model.surfaceplant.ElectricityProduced.quantity().to(
                 'MW').magnitude * self.jobs_created_per_MW_electricity.value))
 
+        # https://github.com/NREL/GEOPHIRES-X/issues/232
+        self.property_tax_created.value = (
+            np.average(model.surfaceplant.ElectricityProduced.quantity().to(
+                'MW').magnitude * self.property_tax_per_MW_electricity.value))
+
+        # https://github.com/NREL/GEOPHIRES-X/issues/232
+        self.total_royalties_created.value = (
+            np.average(model.surfaceplant.ElectricityProduced.quantity().to(
+                'MW').magnitude * self.total_royalty_per_MW_electricity.value))
+
+        # https://github.com/NREL/GEOPHIRES-X/issues/232
+        self.gov_royalties_created.value = (
+            np.average(model.surfaceplant.ElectricityProduced.quantity().to(
+                'MW').magnitude * self.gov_royalty_per_MW_electricity.value))
+
         model.logger.info(f'complete {__class__!s}: {sys._getframe().f_code.co_name}')
 
     def __str__(self):

src/geophires_x/Outputs.py

@@ -1641,7 +1641,11 @@ def PrintOutputs(self, model: Model):
                     f.write(f'      CHP: Percent cost allocation for electrical plant: {model.economics.CAPEX_heat_electricity_plant_ratio.value*100.0:10.2f} %\n')
 
                 if model.surfaceplant.enduse_option.value in [EndUseOptions.ELECTRICITY]:
-                    f.write(f'      Estimated Jobs Created:                                 {model.economics.jobs_created.value}\n')
+                    f.write(f'      Estimates from https://geothermal.org/resources/geothermal-basics:' + NL)
+                    f.write(f'         Estimated Jobs Created:                       {model.economics.jobs_created.value}' + NL)
+                    f.write(f'         Estimated Property tax that will be paid:     {model.economics.property_tax_created.value:10.2f }' + model.economics.property_tax_created.PreferredUnits.value + NL)
+                    f.write(f'         Estimated total royalties that will be paid:  {model.economics.total_royalties_created.value:10.2f }' + model.economics.total_royalties_created.PreferredUnits.value + NL)
+                    f.write(f'         Estimated government royalties to be paid:    {model.economics.gov_royalties_created.value:10.2f }'  + model.economics.gov_royalties_created.PreferredUnits.value+ NL)
 
 
                 f.write(NL)

src/geophires_x/Units.py

@@ -60,6 +60,9 @@ class Units(IntEnum):
     POWERPERUNITVOLUME = auto()
     DECAY_RATE=auto()
     INFLATION_RATE=auto()
+    JOBS_PER_ENERGY=auto()
+    ROYALTY_PER_ENERGY=auto()
+    JOBS = auto()
 
 
 class TemperatureUnit(str, Enum):
@@ -350,3 +353,19 @@ class Decay_RateUnit(str,Enum):
 class Inflation_RateUnit(str,Enum):
     """Decay rate Units"""
     KPASCALPERYEAR = "kPa/yr"
+
+
+class JobsPerEnergyUnit(str,Enum):
+    """Jobs per energy Units"""
+    JOBSPERKW = "jobs/kW"
+    JOBSPERMW = "jobs/MW"
+    JOBSPERGW = "jobs/GW"
+
+class RoyaltyPerEnergyUnit(str,Enum):
+    """Royalty per energy Units"""
+    ROYALTYPERMW = "MUSD/MW"
+    ROYALTYPERKW = "MUSD/kW"
+
+class JobsUnit(str,Enum):
+    """Jobs Units"""
+    JOBS = "jobs"