Skip to content

GeoH2: Arps decline curve #454

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
kbrunik wants to merge 15 commits into
base: develop
Choose a base branch
from
+213 −20
Open

GeoH2: Arps decline curve #454

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
3e49221
WIP arps decline curve
kbrunik Jan 14, 2026
fbee5ba
update decline curve
kbrunik Jan 20, 2026
db2aa10
update geoh2 and docstring
kbrunik Jan 20, 2026
18c0dc9
update test
kbrunik Jan 26, 2026
0567f10
Merge branch 'develop' into decline_curves
kbrunik Jan 26, 2026
a7573f4
update to be percent
kbrunik Jan 27, 2026
02e5145
update documentation
kbrunik Jan 27, 2026
9b590cd
Merge branch 'decline_curves' of https://github.com/kbrunik/H2Integra…
kbrunik Jan 27, 2026
ee1b17a
Merge branch 'develop' into decline_curves
kbrunik Jan 27, 2026
9f4bdac
fix changelog
kbrunik Jan 27, 2026
0eb034d
Merge branch 'develop' into decline_curves
kbrunik Jan 27, 2026
ee5fddf
Merge branch 'develop' into decline_curves
kbrunik Jan 27, 2026
0a6b686
Update h2integrate/converters/hydrogen/geologic/simple_natural_geoh2.py
kbrunik Jan 29, 2026
b77ba18
Update h2integrate/converters/hydrogen/geologic/simple_natural_geoh2.py
kbrunik Jan 29, 2026
259cbbf
precommit fix
kbrunik Jan 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@
- Remove unused dependencies.
- Fixes typos for skipped folders.
- Fixes missing dependencies for `gis` modifier used in new iron mapping tests.
- Add Arps decline rate to natural geologic hydrogen model

## 0.5.1 [December 18, 2025]

Expand Down
6 changes: 6 additions & 0 deletions docs/technology_models/geologic_hydrogen.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
15 changes: 15 additions & 0 deletions examples/04_geo_h2/run_geo_h2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
6 changes: 6 additions & 0 deletions examples/04_geo_h2/tech_config_natural.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
8 changes: 4 additions & 4 deletions examples/test/test_all_examples.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1689,15 +1689,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
== 2264.943
)

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
== 0.6026907731
Copy link
Collaborator

@jmartin4u jmartin4u Jan 28, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LCOH is too low, because "total_hydrogen_produced" is only based off the first year of production

)
with subtests.test("subsurface Capex"):
assert (
Expand All @@ -1721,14 +1721,14 @@ 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) == 4419318.42
)
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
== 3613663.44
)
with subtests.test("surface adjusted opex"):
surface_adjusted_opex = h2i_nat.prob.get_val(
Expand Down
140 changes: 129 additions & 11 deletions h2integrate/converters/hydrogen/geologic/simple_natural_geoh2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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")
Expand All @@ -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"]
Expand All @@ -128,24 +168,102 @@ 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["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"] = ramp_up_flow
# this is lifetime flow which decreases over time
Copy link
Collaborator

@jmartin4u jmartin4u Jan 28, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Going to suggest a couple changes here to temporarily fix the problem of not having declining production over the years. This will make total_hydrogen_produced be based off of the average lifetime flow, rather than the production over the first year

outputs["total_wellhead_gas_produced"] = np.average(wh_flow_profile) * n_timesteps
outputs["total_hydrogen_produced"] = np.average(avg_h2_flow) * n_timesteps

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)
57 changes: 52 additions & 5 deletions h2integrate/converters/hydrogen/geologic/test/test_geologic_h2.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,7 +55,12 @@ def geoh2_subsurface_well():
"site_prospectivity": 0.7,
"wellhead_h2_concentration": 95,
"initial_wellhead_flow": 4000,
"gas_flow_density": 0.11741,
"ramp_up_time_months": 6,
"percent_increase_during_rampup": 5,
"gas_reservoir_size": 1000000,
"use_arps_decline_curve": True,
"decline_fit_params": {"fit_name": "Eagle_Ford"},
},
}
tech_config_dict = {"model_inputs": subsurface_perf_config}
Expand Down Expand Up @@ -91,6 +96,48 @@ def aspen_geoh2_config():
return {"model_inputs": model_inputs}


def test_natural_geoh2_well_performance(subtests, plant_config):
subsurface_perf_config = {
"shared_parameters": {
"borehole_depth": 300,
"well_diameter": "small",
"well_geometry": "vertical",
},
"performance_parameters": {
"rock_type": "peridotite",
"grain_size": 0.01,
"use_prospectivity": False,
"site_prospectivity": 0.7,
"wellhead_h2_concentration": 95,
"initial_wellhead_flow": 40000,
"gas_flow_density": 0.11741,
"ramp_up_time_months": 6,
"percent_increase_during_rampup": 5,
"gas_reservoir_size": 10000000,
"use_arps_decline_curve": True,
"decline_fit_params": {"fit_name": "Eagle_Ford"},
},
}
tech_config_dict = {"model_inputs": subsurface_perf_config}
subsurface_comp = NaturalGeoH2PerformanceModel(
plant_config=plant_config,
tech_config=tech_config_dict,
driver_config={},
)

prob = om.Problem()
prob.model.add_subsystem("perf", subsurface_comp, promotes=["*"])

prob.setup()
prob.run_model()

with subtests.test("Well hydrogen production"):
assert (
pytest.approx(np.mean(prob.model.get_val("perf.hydrogen_out", units="kg/h")), rel=1e-6)
== 22649.432193239327
), 1e-6


def test_aspen_geoh2_performance(subtests, plant_config, geoh2_subsurface_well, aspen_geoh2_config):
prob = om.Problem()
perf_comp = AspenGeoH2SurfacePerformanceModel(
Expand All @@ -114,16 +161,16 @@ def test_aspen_geoh2_performance(subtests, plant_config, geoh2_subsurface_well,
with subtests.test("Well hydrogen production"):
assert (
pytest.approx(np.mean(prob.model.get_val("well.hydrogen_out", units="kg/h")), rel=1e-6)
== 603.4286677531819
== 2264.9432193239327
), 1e-6


def test_aspen_geoh2_performance_cost(
subtests, plant_config, geoh2_subsurface_well, aspen_geoh2_config
):
expected_capex = 1795733.55
expected_capex = 4419318.42
expected_opex = 4567464
expected_varopex = 984842.53
expected_varopex = 3613663.44

prob = om.Problem()
perf_comp = AspenGeoH2SurfacePerformanceModel(
Expand Down Expand Up @@ -198,8 +245,8 @@ def test_aspen_geoh2_refit_coeffs(

with subtests.test("Well hydrogen production"):
assert (
pytest.approx(prob.model.get_val("well.hydrogen_out", units="kg/h"), rel=1e-6)
== 603.4286677531819
pytest.approx(np.mean(prob.model.get_val("well.hydrogen_out", units="kg/h")), rel=1e-6)
== 2264.9432193239327
), 1e-6

with subtests.test("Refit Performance Coeff File"):
Expand Down