Update Natural GeoH2 with Yearly CF (#552)

kbrunik · web-flow · commit dda9dc817420 · 2026-03-02T20:45:51.000-06:00
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -20,6 +20,8 @@
     all tests must be marked via `@pytest.mark.<test-type>`.
   - Partial testing suite refactor to parameterize many of the common fixtures and test routines.
   - `unittest` style tests are refactored to be `pytest` style tests for test consistency.
+- Update finance models to use annual capacity factor and rated production rather than annual production.
+- Update `NaturalGeoH2PerformanceModel` to output yearly metrics.
 - Added a pre-commit hook for `yamlfix` to auto-format YAML files and `yamlfix`'d all YAML files for consistent formatting
 - Insert model names for technologies with control strategies to simplify Pyomo workflows.
 - Refactored pyomo code by splitting apart classes into separate files and removing unused properties
diff --git a/examples/04_geo_h2/run_geo_h2.py b/examples/04_geo_h2/run_geo_h2.py
@@ -12,12 +12,16 @@
 
 
 hydrogen_out = h2i_nat.prob.model.plant.geoh2_well_subsurface.NaturalGeoH2PerformanceModel.get_val(
-    "hydrogen_out", units="kg/h"
+    "annual_hydrogen_produced", units="t/year"
 )
-plt.plot(hydrogen_out)
-plt.xlabel("Time (hours)")
-plt.ylabel("Wellhead Gas Flow (kg/h)")
-plt.title("Wellhead Gas Flow Profile Over First Year")
+start_year = (
+    h2i_nat.prob.model.plant.finance_subgroup_h2.hydrogen_finance_default.params.analysis_start_year
+)
+years = range(start_year, start_year + len(hydrogen_out))
+plt.plot(years, hydrogen_out)
+plt.xlabel("Year")
+plt.ylabel("Hydrogen Production (tonne per annum)")
+plt.title("Hydrogen Production over Well Lifetime")
 plt.grid()
 plt.show()
 
diff --git a/examples/test/test_all_examples.py b/examples/test/test_all_examples.py
@@ -2115,7 +2115,7 @@ def test_natural_geoh2(subtests, temp_copy_of_example):
             pytest.approx(
                 h2i_nat.prob.get_val("finance_subgroup_h2.LCOH", units="USD/kg"), rel=1e-6
             )
-            == 1.5870496689
+            == 1.3089029
         )
     with subtests.test("subsurface Capex"):
         assert (
diff --git a/h2integrate/converters/hydrogen/geologic/simple_natural_geoh2.py b/h2integrate/converters/hydrogen/geologic/simple_natural_geoh2.py
@@ -118,8 +118,27 @@ class NaturalGeoH2PerformanceModel(GeoH2SubsurfacePerformanceBaseClass):
             Hourly wellhead gas production profile from natural accumulations,
             covering one simulated year (8760 hours), in kg/h.
 
+        wellhead_gas_out (ndarray):
+            Hourly wellhead gas production profile used for downstream modeling, in kg/h.
+
         max_wellhead_gas (float):
             Maximum wellhead gas output over the system lifetime, in kg/h.
+
+        rated_hydrogen_production (float):
+            Rated hydrogen production at the wellhead, in kg/h.
+
+        total_wellhead_gas_produced (float):
+            Total mass of gas produced at the wellhead over the simulation period, in kg/year.
+
+        total_hydrogen_produced (float):
+            Total mass of hydrogen produced at the wellhead over the simulation period, in kg/year.
+
+        annual_hydrogen_produced (list):
+            List of total hydrogen produced for each year of the simulation, in kg/year.
+
+        capacity_factor (list):
+            List of capacity factors for each year of the simulation, calculated as the ratio
+            of annual hydrogen production to the maximum hydrogen production of the well.
     """
 
     def setup(self):
@@ -229,40 +248,42 @@ def compute(self, inputs, outputs):
         h2_mw = 2.016
         x_h2 = wh_h2_conc / 100
         w_h2 = x_h2 * h2_mw / (x_h2 * h2_mw + (1 - x_h2) * balance_mw)
-        avg_h2_flow = w_h2 * wh_flow_profile
+        h2_flow = w_h2 * wh_flow_profile
+
+        # calculate the yearly capacity factors and add to dict
+        yearly_h2_cf = []
+        yearly_h2 = []
+        # for each 8760 in the flow profile, calculate the capacity factor for
+        # that year and add to array
+        for year in range(self.options["plant_config"]["plant"]["plant_life"]):
+            start_idx = year * 8760
+            end_idx = start_idx + 8760
+            max_h2_produced = ramp_up_flow * w_h2 * 8760
+            yearly_h2_produced = np.sum(h2_flow[start_idx:end_idx])
+            yearly_h2.append(yearly_h2_produced)
+            yearly_h2_cf.append(yearly_h2_produced / max_h2_produced)
 
         # Parse outputs
         outputs["wellhead_h2_concentration_mass"] = w_h2 * 100
         outputs["wellhead_h2_concentration_mol"] = wh_h2_conc
         outputs["lifetime_wellhead_flow"] = np.average(wh_flow_profile)
-        # fill "wellhead_gas_out_natural" with first year profile from wh_flow_profile
-        # updated to average value over 8760 until simulation length handling is improved
-        # commented code on lines 240 to 246 have the original intended functionality
-        outputs["wellhead_gas_out_natural"] = np.full(
-            n_timesteps, np.average(wh_flow_profile)
-        )  # wh_flow_profile[:n_timesteps]
-        outputs["wellhead_gas_out"] = np.full(
-            n_timesteps, np.average(wh_flow_profile)
-        )  # wh_flow_profile[:n_timesteps]
-        outputs["hydrogen_out"] = np.full(
-            n_timesteps, np.average(avg_h2_flow)
-        )  # avg_h2_flow[:n_timesteps]
+        # use lifetime average because H2I is unable to handle multiyear
+        # commodity_out. Noted in issue #475.
+        outputs["wellhead_gas_out_natural"] = np.full(n_timesteps, np.average(wh_flow_profile))
+        outputs["wellhead_gas_out"] = np.full(n_timesteps, np.average(wh_flow_profile))
+        outputs["hydrogen_out"] = np.full(n_timesteps, np.average(h2_flow))
+
         outputs["max_wellhead_gas"] = ramp_up_flow
-        # this is lifetime flow which decreases over time
-        outputs["total_wellhead_gas_produced"] = (
-            np.average(wh_flow_profile) * n_timesteps
-        )  # np.sum(wh_flow_profile)
-        outputs["total_hydrogen_produced"] = (
-            np.average(avg_h2_flow) * n_timesteps
-        )  # np.sum(avg_h2_flow)
-
-        outputs["annual_hydrogen_produced"] = (
-            outputs["total_hydrogen_produced"] / self.fraction_of_year_simulated
-        )
         outputs["rated_hydrogen_production"] = ramp_up_flow * w_h2
-        outputs["capacity_factor"] = outputs["total_hydrogen_produced"] / (
-            outputs["rated_hydrogen_production"] * self.n_timesteps
-        )
+        # total is amount produced over simulation, which is a single year
+        # for now so lifetime average is more accurate for model
+        outputs["total_wellhead_gas_produced"] = np.average(wh_flow_profile) * n_timesteps
+        outputs["total_hydrogen_produced"] = np.average(h2_flow) * n_timesteps
+
+        # output array of hydrogen produced and capacity factors
+        # for each year of the simulation
+        outputs["annual_hydrogen_produced"] = yearly_h2
+        outputs["capacity_factor"] = yearly_h2_cf
 
     def arps_decline_curve_fit(self, t, qi, Di, b):
         """Arps decline curve model based on Arps (1945)
diff --git a/h2integrate/converters/hydrogen/geologic/test/test_geologic_h2.py b/h2integrate/converters/hydrogen/geologic/test/test_geologic_h2.py
@@ -142,17 +142,99 @@ def test_natural_geoh2_well_performance(subtests, plant_config):
     prob.setup()
     prob.run_model()
 
+    with subtests.test("max wellhead gas"):
+        assert (
+            pytest.approx(prob.model.get_val("perf.max_wellhead_gas", units="kg/h"), rel=1e-6)
+            == 42000.0
+        )
+    # test wellhead_h2_concentration_mass wellhead_h2_concentration_mol well_head_gas_out_natural
+    with subtests.test("wellhead_h2_concentration_mass"):
+        assert (
+            pytest.approx(
+                prob.model.get_val("perf.wellhead_h2_concentration_mass", units="percent"), rel=1e-6
+            )
+            == 62.15760093
+        )
+
+    with subtests.test("wellhead_h2_concentration_mol"):
+        assert (
+            pytest.approx(
+                prob.model.get_val("perf.wellhead_h2_concentration_mol", units="percent"), rel=1e-6
+            )
+            == 95.0
+        )
+    with subtests.test("well_head_gas_out_natural"):
+        assert (
+            pytest.approx(
+                np.sum(prob.model.get_val("perf.wellhead_gas_out_natural", units="kg/h")), rel=1e-6
+            )
+            == 85426105.21216634
+        )
+    with subtests.test("wellhead_gas_out"):
+        assert (
+            pytest.approx(
+                np.sum(prob.model.get_val("perf.wellhead_gas_out", units="kg/h")), rel=1e-6
+            )
+            == 85426105.21217233
+        )
     with subtests.test("Well hydrogen production"):
         assert (
             pytest.approx(np.mean(prob.model.get_val("perf.hydrogen_out", units="kg/h")), rel=1e-6)
             == 6061.508855232839
-        ), 1e-6
+        )
 
     with subtests.test("total h2 out"):
         assert (
             pytest.approx(prob.model.get_val("perf.total_hydrogen_produced", units="kg"), rel=1e-6)
             == 53098817.57183966
-        ), 1e-6
+        )
+    with subtests.test("total_wellhead_gas_produced"):
+        assert (
+            pytest.approx(
+                prob.model.get_val("perf.total_wellhead_gas_produced", units="kg/year"), rel=1e-6
+            )
+            == 85426105.21217233
+        )
+
+    with subtests.test("cf"):
+        assert pytest.approx(
+            prob.model.get_val("perf.capacity_factor", units="unitless"), rel=1e-6
+        ) == [
+            0.8675885,
+            0.42536644,
+            0.25890337,
+            0.18514585,
+            0.14357352,
+            0.1169487,
+            0.09847076,
+            0.08491557,
+            0.07455868,
+            0.06639494,
+        ]
+
+    with subtests.test("annual_hydrogen_production"):
+        assert pytest.approx(
+            prob.model.get_val("perf.annual_hydrogen_produced", units="kg/yr"), rel=1e-6
+        ) == [
+            198409026.01277646,
+            97277155.59721786,
+            59208675.263037845,
+            42341049.07389179,
+            32833862.72029864,
+            26745026.520159103,
+            22519302.237744503,
+            19419361.600134037,
+            17050841.849094056,
+            15183874.84404233,
+        ]
+
+    with subtests.test("rated h2 production"):
+        assert (
+            pytest.approx(
+                prob.model.get_val("perf.rated_hydrogen_production", units="kg/h"), rel=1e-6
+            )
+            == 26106.19239257
+        )
 
 
 @pytest.mark.unit

Original file line number	Diff line number	Diff line change
`@@ -2115,7 +2115,7 @@ def test_natural_geoh2(subtests, temp_copy_of_example):`
`2115`	`2115`	`pytest.approx(`
`2116`	`2116`	`h2i_nat.prob.get_val("finance_subgroup_h2.LCOH", units="USD/kg"), rel=1e-6`
`2117`	`2117`	`)`
`2118`		`- == 1.5870496689`
	`2118`	`+ == 1.3089029`
`2119`	`2119`	`)`
`2120`	`2120`	`with subtests.test("subsurface Capex"):`
`2121`	`2121`	`assert (`