diff --git a/CHANGELOG.md b/CHANGELOG.md index 486553cc9..137844dd3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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] diff --git a/docs/technology_models/geologic_hydrogen.md b/docs/technology_models/geologic_hydrogen.md index c9f8ddd2d..99b52a847 100644 --- a/docs/technology_models/geologic_hydrogen.md +++ b/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 diff --git a/examples/04_geo_h2/run_geo_h2.py b/examples/04_geo_h2/run_geo_h2.py index 5b0938936..19215c6b4 100644 --- a/examples/04_geo_h2/run_geo_h2.py +++ b/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") diff --git a/examples/04_geo_h2/tech_config_natural.yaml b/examples/04_geo_h2/tech_config_natural.yaml index d23814653..875931efc 100644 --- a/examples/04_geo_h2/tech_config_natural.yaml +++ b/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 diff --git a/examples/test/test_all_examples.py b/examples/test/test_all_examples.py index 8141b5497..91e2d74f4 100644 --- a/examples/test/test_all_examples.py +++ b/examples/test/test_all_examples.py @@ -1689,7 +1689,7 @@ 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"): @@ -1697,7 +1697,7 @@ def test_natural_geoh2(subtests): pytest.approx( h2i_nat.prob.get_val("finance_subgroup_h2.LCOH", units="USD/kg"), rel=1e-6 ) - == 1.59307314 + == 0.6026907731 ) with subtests.test("subsurface Capex"): assert ( @@ -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( diff --git a/h2integrate/converters/hydrogen/geologic/simple_natural_geoh2.py b/h2integrate/converters/hydrogen/geologic/simple_natural_geoh2.py index 0710244d2..c23f065e4 100644 --- a/h2integrate/converters/hydrogen/geologic/simple_natural_geoh2.py +++ b/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,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 + 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) diff --git a/h2integrate/converters/hydrogen/geologic/test/test_geologic_h2.py b/h2integrate/converters/hydrogen/geologic/test/test_geologic_h2.py index 417d5c1ab..477599833 100644 --- a/h2integrate/converters/hydrogen/geologic/test/test_geologic_h2.py +++ 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.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} @@ -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( @@ -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( @@ -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"):