GeoH2: Arps decline curve

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.1
+        ramp_up_time_months: 6 #months
+        percent_increase_during_rampup: 0.05
         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

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):
@@ -109,7 +134,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")
@@ -118,6 +151,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"]
@@ -128,24 +168,93 @@ 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 * (1 + 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 / 28.3168 / 1e6
+
+                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)
+            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["lifetime_wellhead_flow"] = np.average(wh_flow_profile)
+        # fill "wellhead_gas_out_natural" with first year profile from wh_flow_profile
+        outputs["wellhead_gas_out_natural"] = wh_flow_profile[:n_timesteps]
+        outputs["wellhead_gas_out"] = wh_flow_profile[:n_timesteps]
+        outputs["hydrogen_out"] = avg_h2_flow[:n_timesteps]
         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"])
+
+    def arps_decline_curve_fit(self, t, qi, Di, b):
+        """Arps decline model.
+
+        Relevant literature: https://doi.org/10.1016/j.jngse.2021.103818
+
+        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.
+        """
+        qi = 0
+        if np.isclose(b, 0):
+            return qi * np.exp(-Di * t)
+        else:
+            return qi / (1 + b * Di * t) ** (1 / b)

h2integrate/converters/hydrogen/geologic/test/test_geologic_h2.py

@@ -55,7 +55,12 @@ def geoh2_subsurface_well():
             "site_prospectivity": 0.7,
             "wellhead_h2_concentration": 95,
             "initial_wellhead_flow": 4000,
+            "gas_flow_density": 0.1,
+            "ramp_up_time": 6,
+            "percent_increase_during_rampup": 0.05,
             "gas_reservoir_size": 1000000,
+            "use_arps_decline_curve": True,
+            "decline_fit_params": {"fit_name": "Eagle_Ford"},
         },
     }
     tech_config_dict = {"model_inputs": subsurface_perf_config}