GeoH2: Arps decline curve

CHANGELOG.md

@@ -37,6 +37,7 @@
   enforced by a newly added test to ensure adherence.
 - Remove `pytest-subtests` as it's incorporated into pytest as of v9, and is an archived project.
 - Added `PerformanceModelBaseClass` and standardized outputs of converter performance models
+- Add Arps decline rate to natural geologic hydrogen model
 
 ## 0.5.1 [December 18, 2025]
 

docs/technology_models/geologic_hydrogen.md

@@ -24,6 +24,12 @@ The modeling approach in this performance model is informed by:
 - [Gelman et al. (USGS)](https://doi.org/10.3133/pp1900)
 - [Tang et al. (Southwest Petroleum University)](https://doi.org/10.1016/j.petsci.2024.07.029)
 
+The natural geologic hydrogen model is able to model well decline over time and there are two methods of decline.
+1. If not specified or `use_arps_decline_curve` is `False`: the decline rate will be linear over the lifetime of the well as defined in the attribute `plant_config["plant"]["plant_life"]`
+2. If `use_arps_decline_curve` is `True`: The well production will decline according to the Arps model as defined in [Tang et al. (Southwest Petroleum University)](https://doi.org/10.1016/j.petsci.2024.07.029). There are several options for using the decline curve to model well production.
+    1. `decline_fit_params` is a dictionary where a user can specify the decline rate and loss rate. It should be noted that typically the modeling of these decline rates are monthly. This model uses hourly timesteps so the decline and loss rate should be modified accordingly.
+    2. `decline_fit_params["fit_name"]` is an optional string that can be specified within the `decline_fit_params` dictionary that will model a decline rate similar to those noted at either the "Bakken", "Eagle Ford" or "Permian" shale wells documented in [Tang et al. (Southwest Petroleum University)](https://doi.org/10.1016/j.petsci.2024.07.029) Figure 7.
+
 (templeton-serpentinization-geoh2-performance)=
 ### Templeton Serpentinization GeoH2 Performance
 

examples/04_geo_h2/run_geo_h2.py

@@ -8,6 +8,21 @@
 h2i_nat.run()
 h2i_nat.post_process()
 
+import matplotlib.pyplot as plt
+
+
+hydrogen_out = (
+    h2i_nat.prob.model.plant.geoh2_well_subsurface.simple_natural_geoh2_performance.get_val(
+        "hydrogen_out"
+    )
+)
+plt.plot(hydrogen_out)
+plt.xlabel("Time (hours)")
+plt.ylabel("Wellhead Gas Flow (kg/h)")
+plt.title("Wellhead Gas Flow Profile Over First Year")
+plt.grid()
+plt.show()
+
 
 # Create an H2I model for stimulated geologic hydrogen production
 h2i_stim = H2IntegrateModel("04_geo_h2_stimulated.yaml")

examples/04_geo_h2/tech_config_natural.yaml

@@ -19,7 +19,13 @@ technologies:
         site_prospectivity: 0.7
         wellhead_h2_concentration: 95
         initial_wellhead_flow: 4000
+        gas_flow_density: 0.11741 # for 95% h2
+        ramp_up_time_months: 6 #months
+        percent_increase_during_rampup: 5 # percent
         gas_reservoir_size: 1000000
+        use_arps_decline_curve: True
+        decline_fit_params:
+          fit_name: Eagle_Ford
       cost_parameters:
         use_cost_curve: True
         test_drill_cost: 500000

examples/test/test_all_examples.py

@@ -1693,15 +1693,15 @@ def test_natural_geoh2(subtests):
                 np.mean(h2i_nat.model.get_val("geoh2_well_subsurface.hydrogen_out", units="kg/h")),
                 rel=1e-6,
             )
-            == 603.4286677531819
+            == 606.1508855232839
         )
 
     with subtests.test("integrated LCOH"):
         assert (
             pytest.approx(
                 h2i_nat.prob.get_val("finance_subgroup_h2.LCOH", units="USD/kg"), rel=1e-6
             )
-            == 1.59307314
+            == 1.5870496689
         )
     with subtests.test("subsurface Capex"):
         assert (
@@ -1725,14 +1725,15 @@ def test_natural_geoh2(subtests):
 
     with subtests.test("surface Capex"):
         assert (
-            pytest.approx(h2i_nat.model.get_val("geoh2_well_surface.CapEx"), rel=1e-6) == 1795733.55
+            pytest.approx(h2i_nat.model.get_val("geoh2_well_surface.CapEx"), rel=1e-6)
+            == 1800711.83796
         )
     with subtests.test("surface fixed Opex"):
         assert pytest.approx(h2i_nat.model.get_val("geoh2_well_surface.OpEx"), rel=1e-6) == 4567464
     with subtests.test("surface variable Opex"):
         assert (
             pytest.approx(h2i_nat.model.get_val("geoh2_well_surface.VarOpEx"), rel=1e-6)
-            == 984842.53
+            == 989213.8787
         )
     with subtests.test("surface adjusted opex"):
         surface_adjusted_opex = h2i_nat.prob.get_val(

h2integrate/converters/hydrogen/geologic/simple_natural_geoh2.py

@@ -2,6 +2,7 @@
 from attrs import field, define
 
 from h2integrate.core.utilities import merge_shared_inputs
+from h2integrate.core.validators import range_val
 from h2integrate.converters.hydrogen.geologic.h2_well_subsurface_baseclass import (
     GeoH2SubsurfacePerformanceConfig,
     GeoH2SubsurfacePerformanceBaseClass,
@@ -32,15 +33,39 @@ class NaturalGeoH2PerformanceConfig(GeoH2SubsurfacePerformanceConfig):
             Hydrogen flow rate measured immediately after well completion, in kilograms
             per hour (kg/h).
 
+        gas_flow_density (float):
+            Density of the wellhead gas flow, in kilograms per cubic meter (kg/m^3).
+
+        ramp_up_time_months (float):
+            Number of months after initial flow from the well before full utilization.
+
+        percent_increase_during_rampup (float):
+            Percent increase in wellhead flow during ramp-up period in percent (%).
+
         gas_reservoir_size (float):
             Total amount of hydrogen stored in the geologic accumulation, in tonnes (t).
+
+        use_arps_decline_curve (bool):
+            Whether to use the Arps decline curve model for well production decline.
+
+        decline_fit_params (dict):
+            (Optional) Parameters for the Arps decline curve model, including:
+                - 'Di' (float): Decline rate.
+                - 'b' (float): Loss rate.
+                - 'fit_name' (str): Name of the well fit to use. If provided, overrides Di and b.
+                    Options are "Eagle_Ford" or "Permian" or "Bakken".
     """
 
     use_prospectivity: bool = field()
     site_prospectivity: float = field()
     wellhead_h2_concentration: float = field()
     initial_wellhead_flow: float = field()
+    gas_flow_density: float = field()
+    ramp_up_time_months: float = field()
+    percent_increase_during_rampup: float = field(validator=range_val(0, 100))
     gas_reservoir_size: float = field()
+    use_arps_decline_curve: bool = field()
+    decline_fit_params: dict = field(default=None)
 
 
 class NaturalGeoH2PerformanceModel(GeoH2SubsurfacePerformanceBaseClass):
@@ -110,7 +135,15 @@ def setup(self):
             "wellhead_h2_concentration", units="percent", val=self.config.wellhead_h2_concentration
         )
         self.add_input("initial_wellhead_flow", units="kg/h", val=self.config.initial_wellhead_flow)
+        self.add_input("gas_flow_density", units="kg/m**3", val=self.config.gas_flow_density)
         self.add_input("gas_reservoir_size", units="t", val=self.config.gas_reservoir_size)
+        self.add_input("ramp_up_time", units="yr/12", val=self.config.ramp_up_time_months)
+        self.add_input(
+            "percent_increase_during_rampup",
+            units="percent",
+            val=self.config.percent_increase_during_rampup,
+            desc="Percent increase in wellhead flow during ramp-up period in percent (%)",
+        )
 
         self.add_output("wellhead_h2_concentration_mass", units="percent")
         self.add_output("wellhead_h2_concentration_mol", units="percent")
@@ -119,6 +152,13 @@ def setup(self):
         self.add_output("max_wellhead_gas", units="kg/h")
 
     def compute(self, inputs, outputs):
+        n_timesteps = self.options["plant_config"]["plant"]["simulation"]["n_timesteps"]
+
+        # Coerce scalar inputs to Python scalars (handles 0-d and 1-d arrays)
+        ramp_up_time = float(np.asarray(inputs["ramp_up_time"]).item())
+        percent_increase = float(np.asarray(inputs["percent_increase_during_rampup"]).item())
+        init_wh_flow = float(np.asarray(inputs["initial_wellhead_flow"]).item())
+
         if self.config.rock_type == "peridotite":  # TODO: sub-models for different rock types
             # Calculate expected wellhead h2 concentration from prospectivity
             prospectivity = inputs["site_prospectivity"]
@@ -129,31 +169,114 @@ def compute(self, inputs, outputs):
 
         # Calculated average wellhead gas flow over well lifetime
         init_wh_flow = inputs["initial_wellhead_flow"]
-        lifetime = self.options["plant_config"]["plant"]["plant_life"]
-        n_timesteps = self.options["plant_config"]["plant"]["simulation"]["n_timesteps"]
-        avg_wh_flow = (-0.193 * np.log(lifetime) + 0.6871) * init_wh_flow  # temp. fit to Arps data
+
+        # Coerce scalar inputs to Python scalars (handles 0-d and 1-d arrays)
+        ramp_up_time = float(np.asarray(inputs["ramp_up_time"]).item())
+        percent_increase = float(np.asarray(inputs["percent_increase_during_rampup"]).item())
+        init_wh_flow = float(np.asarray(inputs["initial_wellhead_flow"]).item())
+
+        # Apply ramp-up assumed linear increase
+        ramp_up_steps = int(ramp_up_time * (n_timesteps / 12))  # hrs
+        if ramp_up_steps > 0:
+            ramp_up_flow = init_wh_flow * ((100 + percent_increase) / 100)
+            ramp_up_profile = np.linspace(init_wh_flow, ramp_up_flow, ramp_up_steps)
+        else:
+            ramp_up_flow = init_wh_flow
+        remaining_steps = (
+            n_timesteps * self.options["plant_config"]["plant"]["plant_life"] - ramp_up_steps
+        )  # remaining time steps in lifetime
+
+        # Use decline curve modeling if selected
+        if self.config.use_arps_decline_curve:
+            t = np.arange(remaining_steps)  # hrs
+            if self.config.decline_fit_params and "fit_name" in self.config.decline_fit_params:
+                # decline curves from literature is in million standard cubic feet per hour
+                ramp_up_flow_m3 = ramp_up_flow / inputs["gas_flow_density"]  # m3/h
+                # convert from m3/h to million standard cubic feet per hour (MMSCF/h)
+                ramp_up_flow_mmscf = ramp_up_flow_m3 / 28316.846592  # 1 MMSCF = 28316.846592 m3
+
+                # fits for MMSCF/h based on flow rates Figure 7 in Tang et al. (2024)
+                fit_name = self.config.decline_fit_params["fit_name"]
+                if fit_name == "Eagle_Ford":
+                    Di = 0.000157
+                    b = 0.932
+                elif fit_name == "Permian":
+                    Di = 0.000087
+                    b = 0.708
+                elif fit_name == "Bakken":
+                    Di = 0.000076
+                    b = 0.784
+                else:
+                    msg = f"Unknown fit_name '{fit_name}' \
+                        for Arps decline curve. Valid options are \
+                        'Eagle_Ford', 'Permian', or 'Bakken'."
+                    raise ValueError(msg)
+                decline_profile = self.arps_decline_curve_fit(t, ramp_up_flow_mmscf, Di, b)
+                # convert back to kg/h from MMSCF/h
+                decline_profile = decline_profile * 28316.846592 * inputs["gas_flow_density"]
+            else:
+                Di = self.config.decline_fit_params.get("Di")
+                b = self.config.decline_fit_params.get("b")
+                decline_profile = self.arps_decline_curve_fit(t, ramp_up_flow, Di, b)
+        else:
+            # linear decline for rest of lifetime
+            decline_profile = np.linspace(ramp_up_flow, 0, remaining_steps)
+
+        wh_flow_profile = np.concatenate((ramp_up_profile, decline_profile))
 
         # Calculated hydrogen flow out
         balance_mw = 23.32  # Note: this is based on Aspen models in aspen_surface_processing.py
         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 * avg_wh_flow
+        avg_h2_flow = w_h2 * wh_flow_profile
 
         # Parse outputs
         outputs["wellhead_h2_concentration_mass"] = w_h2 * 100
         outputs["wellhead_h2_concentration_mol"] = wh_h2_conc
-        outputs["lifetime_wellhead_flow"] = avg_wh_flow
-        outputs["wellhead_gas_out_natural"] = np.full(n_timesteps, avg_wh_flow)
-        outputs["wellhead_gas_out"] = np.full(n_timesteps, avg_wh_flow)
-        outputs["hydrogen_out"] = np.full(n_timesteps, avg_h2_flow)
-        outputs["max_wellhead_gas"] = init_wh_flow
-        outputs["total_wellhead_gas_produced"] = np.sum(outputs["wellhead_gas_out"])
-        outputs["total_hydrogen_produced"] = np.sum(outputs["hydrogen_out"])
-        outputs["annual_hydrogen_produced"] = outputs["total_hydrogen_produced"] * (
-            1 / self.fraction_of_year_simulated
-        )
-        outputs["rated_hydrogen_production"] = init_wh_flow  # TODO: double check
-        outputs["capacity_factor"] = outputs["total_hydrogen_produced"] / (
-            outputs["rated_hydrogen_production"] * self.n_timesteps
-        )
+        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]
+        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)
+
+    def arps_decline_curve_fit(self, t, qi, Di, b):
+        """Arps decline curve model based on Arps (1945)
+            https://doi.org/10.2118/945228-G.
+
+        Other Relevant literature:
+            Tang et al. (2024) https://doi.org/10.1016/j.jngse.2021.103818
+            Adapted the Arps model from Table 2 to fit the
+            monthly gas rates from Figure 7 to characterize natural hydrogen
+            well production decline for the three oil shale wells
+            (Bakken, Eagle Ford and Permian).
+
+        Args:
+            t (np.array): Well production duration from max production.
+            qi (float): Maximum initial production rate.
+            Di (float): Decline rate.
+            b (float): Loss rate.
+
+        Returns:
+            (np.array): Production rate at time t.
+        """
+        if np.isclose(b, 0):
+            return qi * np.exp(-Di * t)
+        else:
+            return qi / (1 + b * Di * t) ** (1 / b)