Oxygen output from PEM electrolyzer (#642)

* added oxygen output from electrolyzer

* updated example 14 to calc LCOO

---------

Co-authored-by: John Jasa <johnjasa11@gmail.com>

Author: elenya-grant

CHANGELOG.md

@@ -32,6 +32,7 @@
 - Made generating an XDSM diagram from connections in a model optional and added documentation on model visualization. [PR 629](https://github.com/NatLabRockies/H2Integrate/pull/629)
 - Added a storage performance baseclass model `StoragePerformanceBase` and updated the other storage performance models to inherit it [PR 624](https://github.com/NatLabRockies/H2Integrate/pull/624)
 - Modified the calc tilt angle function for pysam solar to support latitudes in the southern hemisphere [PR 646](https://github.com/NatLabRockies/H2Integrate/pull/646)
+- Added oxygen production metrics and as outputs to `ECOElectrolyzerPerformanceModel` [PR 642](https://github.com/NatLabRockies/H2Integrate/pull/642)
 
 ## 0.7.1 [March 13, 2026]
 

examples/14_wind_hydrogen_dispatch/inputs/plant_config.yaml

@@ -59,3 +59,7 @@ finance_parameters:
       commodity: hydrogen
       commodity_stream: h2_storage  # use only dispatched hydrogen from h2_storage controller in finance calc
       technologies: [wind, electrolyzer, h2_storage]
+    oxygen:
+      commodity: oxygen
+      commodity_stream: electrolyzer
+      technologies: [wind, electrolyzer]

examples/test/test_all_examples.py

@@ -827,6 +827,14 @@ def test_hydrogen_dispatch_example(subtests, temp_copy_of_example):
             )
             == 7.564000289456695
         )
+    with subtests.test("Check LCOO"):
+        assert (
+            pytest.approx(
+                model.prob.get_val("finance_subgroup_oxygen.LCOO", units="USD/kg")[0],
+                rel=1e-5,
+            )
+            == 0.666523050
+        )
 
 
 @pytest.mark.integration

h2integrate/converters/hydrogen/pem_electrolyzer.py

@@ -90,6 +90,11 @@ def setup(self):
             units="MW",
             desc="Size of the electrolyzer in MW",
         )
+
+        self.add_output("oxygen_out", val=0, shape=self.n_timesteps, units="kg/h")
+        self.add_output("rated_oxygen_production", val=0, shape=1, units="kg/h")
+        self.add_output("annual_oxygen_produced", val=0, shape=self.plant_life, units="kg/yr")
+
         self.add_input("cluster_size", val=-1.0, units="MW")
         self.add_input("max_hydrogen_capacity", val=1000.0, units="kg/h")
         # TODO: add feedstock inputs and consumption outputs
@@ -188,8 +193,13 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
             "Capacity Factor [-]"
         ].values
         outputs["annual_hydrogen_produced"] = H2_Results["Life: Annual H2 production [kg/year]"]
-
         # TODO: replace above line w below
         # outputs["annual_hydrogen_produced"] = H2_Results["Performance Schedules"][
         #     "Annual H2 Production [kg/year]"
         # ].values
+
+        outputs["oxygen_out"] = H2_Results["Oxygen Hourly Production [kg/hr]"]
+        outputs["rated_oxygen_production"] = H2_Results["Rated BOL: O2 Production [kg/hr]"]
+        outputs["annual_oxygen_produced"] = H2_Results["Performance Schedules"][
+            "Annual O2 Production [kg/year]"
+        ]

h2integrate/converters/hydrogen/pem_model/PEM_H2_LT_electrolyzer_Clusters.py

@@ -29,6 +29,8 @@
 import rainflow
 from scipy import interpolate
 
+from h2integrate.tools.constants import O2_MW
+
 
 np.set_printoptions(threshold=sys.maxsize)
 
@@ -174,7 +176,7 @@ def run(self, input_external_power_kw):
         cluster_cycling = np.array(cluster_cycling)
 
         # how much to reduce h2 by based on cycling status
-        h2_multiplier = np.where(cluster_cycling > 0, startup_ratio, 1)
+        production_multiplier = np.where(cluster_cycling > 0, startup_ratio, 1)
         # number of "stacks" on at a single time
         self.n_stacks_op = self.max_stacks * self.cluster_status
         # n_stacks_op is now either number of pem per cluster or 0 if cluster is off!
@@ -218,10 +220,15 @@ def run(self, input_external_power_kw):
         # h20_gal_used_system=self.water_supply(h2_kg_hr_system_init)
         p_consumed_max, rated_h2_hr = self.rated_h2_prod()
         # scales h2 production to account for start-up time if going from off->on
-        h2_kg_hr_system = h2_kg_hr_system_init * h2_multiplier
+        h2_kg_hr_system = h2_kg_hr_system_init * production_multiplier
 
         h20_gal_used_system = self.water_supply(h2_kg_hr_system)
 
+        # Get oxygen production and rated used
+        rated_o2_hr = self.rated_o2_prod()
+        o2_kg_hr_system_init = self.o2_production_rate(stack_current, self.n_stacks_op)
+        o2_kg_hr_system = o2_kg_hr_system_init * production_multiplier
+
         pem_cf = np.sum(h2_kg_hr_system) / (rated_h2_hr * len(input_power_kw) * self.max_stacks)
         efficiency = self.system_efficiency(input_power_kw, stack_current)  # Efficiency as %-HHV
 
@@ -231,6 +238,7 @@ def run(self, input_external_power_kw):
         h2_results["hydrogen production no start-up time"] = h2_kg_hr_system_init
         h2_results["hydrogen_hourly_production"] = h2_kg_hr_system
         h2_results["water_hourly_usage_kg"] = h20_gal_used_system * 3.79
+        h2_results["oxygen_hourly_production"] = o2_kg_hr_system
         h2_results["electrolyzer_total_efficiency_perc"] = efficiency
         h2_results["kwh_per_kgH2"] = input_power_kw / h2_kg_hr_system
         h2_results["Power Consumed [kWh]"] = system_power_consumed
@@ -246,13 +254,17 @@ def run(self, input_external_power_kw):
             p_consumed_max * self.max_stacks
         )
         h2_results_aggregates["Cluster Rated H2 Production [kg/hr]"] = rated_h2_hr * self.max_stacks
+        h2_results_aggregates["Cluster Rated O2 Production [kg/hr]"] = rated_o2_hr * self.max_stacks
         h2_results_aggregates["gal H20 per kg H2"] = np.sum(h20_gal_used_system) / np.sum(
             h2_kg_hr_system
         )
         h2_results_aggregates["Stack Rated Efficiency [kWh/kg]"] = p_consumed_max / rated_h2_hr
         h2_results_aggregates["Cluster Rated H2 Production [kg/yr]"] = (
             rated_h2_hr * len(input_power_kw) * self.max_stacks
         )
+        h2_results_aggregates["Cluster Rated O2 Production [kg/yr]"] = (
+            rated_o2_hr * len(input_power_kw) * self.max_stacks
+        )
         h2_results_aggregates["Operational Time / Simulation Time (ratio)"] = (
             self.percent_of_sim_operating
         )  # added
@@ -367,13 +379,17 @@ def make_yearly_performance_dict(self, power_in_kW, V_deg, V_cell, I_op, grid_co
         cluster_cycling = [0, *list(np.diff(self.cluster_status))]  # no delay at beginning of sim
         cluster_cycling = np.array(cluster_cycling)
         startup_ratio = 1 - (600 / 3600)  # TODO: don't have this hard-coded
-        h2_multiplier = np.where(cluster_cycling > 0, startup_ratio, 1)
+        production_multiplier = np.where(cluster_cycling > 0, startup_ratio, 1)
 
         _, rated_h2_pr_stack_BOL = self.rated_h2_prod()
+        rated_o2_pr_stack_BOL = self.rated_o2_prod()
         rated_h2_pr_sim = rated_h2_pr_stack_BOL * self.max_stacks * sim_length
+        rated_o2_pr_sim = rated_o2_pr_stack_BOL * self.max_stacks * sim_length
 
         kg_h2_pr_sim = np.zeros(int(self.plant_life_years))
+        kg_o2_pr_sim = np.zeros(int(self.plant_life_years))
         capfac_per_sim = np.zeros(int(self.plant_life_years))
+        o2_capfac_per_sim = np.zeros(int(self.plant_life_years))
         d_sim = np.zeros(int(self.plant_life_years))
         power_pr_yr_kWh = np.zeros(int(self.plant_life_years))
         Vdeg0 = 0
@@ -396,20 +412,28 @@ def make_yearly_performance_dict(self, power_in_kW, V_deg, V_cell, I_op, grid_co
                     power_in_kW, V_cell, V_deg_pr_sim
                 )
                 h2_kg_hr_system_init = self.h2_production_rate(stack_current, self.n_stacks_op)
+                o2_kg_hr_system_init = self.o2_production_rate(stack_current, self.n_stacks_op)
                 # total_sim_input_power = self.max_stacks*np.sum(power_in_kW)
                 power_pr_yr_kWh[i] = self.max_stacks * np.sum(power_in_kW)
             else:
                 h2_kg_hr_system_init = self.h2_production_rate(I_op, self.n_stacks_op)
                 h2_kg_hr_system_init = h2_kg_hr_system_init * np.ones(len(power_in_kW))
+
+                o2_kg_hr_system_init = self.o2_production_rate(I_op, self.n_stacks_op)
+                o2_kg_hr_system_init = o2_kg_hr_system_init * np.ones(len(power_in_kW))
+
                 annual_power_consumed_kWh = (
                     self.max_stacks * I_op * (V_cell + V_deg_pr_sim) * self.N_cells / 1000
                 )
                 # total_sim_input_power = np.sum(annual_power_consumed_kWh)
                 power_pr_yr_kWh[i] = np.sum(annual_power_consumed_kWh)
 
-            h2_kg_hr_system = h2_kg_hr_system_init * h2_multiplier
+            h2_kg_hr_system = h2_kg_hr_system_init * production_multiplier
+            o2_kg_hr_system = o2_kg_hr_system_init * production_multiplier
+            kg_o2_pr_sim[i] = np.sum(o2_kg_hr_system)
             kg_h2_pr_sim[i] = np.sum(h2_kg_hr_system)
             capfac_per_sim[i] = np.sum(h2_kg_hr_system) / rated_h2_pr_sim
+            o2_capfac_per_sim[i] = np.sum(o2_kg_hr_system) / rated_o2_pr_sim
             d_sim[i] = V_deg_pr_sim[sim_length - 1]
             Vdeg0 = V_deg_pr_sim[sim_length - 1]
         performance_by_year = {}
@@ -427,7 +451,8 @@ def make_yearly_performance_dict(self, power_in_kW, V_deg, V_cell, I_op, grid_co
             zip(year, self.eta_h2_hhv / (power_pr_yr_kWh / kg_h2_pr_sim))
         )
         performance_by_year["Annual Energy Used [kWh/year]"] = dict(zip(year, power_pr_yr_kWh))
-
+        performance_by_year["Annual O2 Production [kg/year]"] = dict(zip(year, kg_o2_pr_sim))
+        performance_by_year["O2 Capacity Factor [-]"] = dict(zip(year, o2_capfac_per_sim))
         return performance_by_year
 
     def reset_uptime_degradation_rate(self, uptime_hours_until_eol):
@@ -606,6 +631,22 @@ def rated_h2_prod(self):
         max_h2_stack_kg = self.h2_production_rate(I_max, 1)
         return P_consumed_stack_kw, max_h2_stack_kg
 
+    def rated_o2_prod(self):
+        """Calculates the oxygen production per stack at the rated current in kg
+
+        Returns:
+            float: rated oxygen production in kg/dt
+        """
+        i = self.output_dict["BOL Efficiency Curve Info"].index[
+            self.output_dict["BOL Efficiency Curve Info"]["Power Sent [kWh]"]
+            == self.stack_rating_kW
+        ]
+        I_max = self.output_dict["BOL Efficiency Curve Info"]["Current"].iloc[i].values[0]
+        # I_max = calc_current((self.stack_rating_kW,self.T_C),*self.curve_coeff)
+
+        max_o2_stack_kg = self.o2_production_rate(I_max, 1)
+        return max_o2_stack_kg
+
     def external_power_supply(self, input_external_power_kw):
         """
         External power source (grid or REG) which will need to be stepped
@@ -892,6 +933,33 @@ def h2_production_rate(self, stack_current, n_stacks_op):
 
         return h2_produced_kg_hr_system
 
+    def o2_production_rate(self, stack_current, n_stacks_op):
+        """Calculate the oxygen production rate of the cluster without warm-up losses.
+        These equations are based on Faraday's law:
+
+        O2 [mol/sec] = eta_F*(n_cells*I)/(4F)
+
+        where eta_F is the faradaic efficiency. This can be found in Equation 12 of
+        `Simple PEM water electrolyser model and experimental validation <http://dx.doi.org/10.1016/j.ijhydene.2011.09.027>`_.
+
+        Args:
+            stack_current (float | np.array): current of the stack in A
+            n_stacks_op (int | np.array[int]): number of stacks that are operating in the stack.
+
+        Returns:
+            float | np.array: oxygen production profile of the cluster in kg/dt
+        """
+        # Calculate the faradaic efficiency
+        n_Tot = self.faradaic_efficiency(stack_current)
+        # Faraday's law
+        o2_production_rate = n_Tot * ((self.N_cells * stack_current) / (4 * self.F))  # mol/s
+        # O2_MW is in g/mol
+        # mol/s * g/mol = g/s
+        o2_production_rate_g_s = o2_production_rate * O2_MW
+        o2_produced_kg_dt = o2_production_rate_g_s * self.dt / 1000
+        o2_produced_kg_dt_system = o2_produced_kg_dt * n_stacks_op
+        return o2_produced_kg_dt_system
+
     def water_supply(self, h2_kg_hr):
         """
         Calculate water supply rate based system efficiency and H2 production
@@ -920,7 +988,7 @@ def run_grid_connected_workaround(self, power_input_signal, current_signal):
         cluster_cycling = np.array(cluster_cycling)
         power_per_stack = np.where(self.n_stacks_op > 0, power_input_signal / self.n_stacks_op, 0)
 
-        h2_multiplier = np.where(cluster_cycling > 0, startup_ratio, 1)
+        production_multiplier = np.where(cluster_cycling > 0, startup_ratio, 1)
         self.n_stacks_op = self.max_stacks * self.cluster_status
 
         V_init = self.cell_design(self.T_C, current_signal)
@@ -936,7 +1004,9 @@ def run_grid_connected_workaround(self, power_input_signal, current_signal):
 
         h2_kg_hr_system_init = self.h2_production_rate(current_signal, self.n_stacks_op)
         p_consumed_max, rated_h2_hr = self.rated_h2_prod()
-        h2_kg_hr_system = h2_kg_hr_system_init * h2_multiplier  # scales h2 production to account
+        h2_kg_hr_system = (
+            h2_kg_hr_system_init * production_multiplier
+        )  # scales h2 production to account
         # for start-up time if going from off->on
         h20_gal_used_system = self.water_supply(h2_kg_hr_system)
 

h2integrate/converters/hydrogen/pem_model/run_h2_PEM.py

@@ -13,6 +13,8 @@ def clean_up_final_outputs(h2_tot, h2_ts):
             "Cluster Rated H2 Production [kg/yr]",
             "Stack Rated H2 Production [kg/hr]",
             "Stack Rated Power Consumed [kWh]",
+            "Cluster Rated O2 Production [kg/hr]",
+            "Cluster Rated O2 Production [kg/yr]",
         ]
     )
     h2_ts.sum(axis=1)
@@ -21,6 +23,7 @@ def clean_up_final_outputs(h2_tot, h2_ts):
         "Power Consumed [kWh]",
         "hydrogen production no start-up time",
         "hydrogen_hourly_production",
+        "oxygen_hourly_production",
         "water_hourly_usage_kg",
     ]
 
@@ -89,6 +92,7 @@ def run_h2_PEM(
     # time-series info (unchanged)
     energy_input_to_electrolyzer = h2_ts.loc["Input Power [kWh]"].sum()
     hydrogen_hourly_production = h2_ts.loc["hydrogen_hourly_production"].sum()
+    oxygen_hourly_production = h2_ts.loc["oxygen_hourly_production"].sum()
     hourly_system_electrical_usage = h2_ts.loc["Power Consumed [kWh]"].sum()
     water_hourly_usage = h2_ts.loc["water_hourly_usage_kg"].sum()
     avg_eff_perc = eta_h2_hhv * hydrogen_hourly_production / hourly_system_electrical_usage
@@ -102,6 +106,7 @@ def run_h2_PEM(
 
     # Beginning of Life (BOL) Rated Specs (attributes/system design)
     max_h2_pr_hr = h2_tot.loc["Cluster Rated H2 Production [kg/hr]"].sum()
+    max_o2_pr_hr = h2_tot.loc["Cluster Rated O2 Production [kg/hr]"].sum()
     max_pwr_pr_hr = h2_tot.loc["Cluster Rated Power Consumed [kWh]"].sum()
     rated_kWh_pr_kg = h2_tot.loc["Stack Rated Efficiency [kWh/kg]"].mean()
     elec_rated_h2_capacity_kgpy = h2_tot.loc["Cluster Rated H2 Production [kg/yr]"].sum()
@@ -110,6 +115,7 @@ def run_h2_PEM(
     atrribute_desc = [
         "Efficiency [kWh/kg]",
         "H2 Production [kg/hr]",
+        "O2 Production [kg/hr]",
         "Power Consumed [kWh]",
         "Annual H2 Production [kg/year]",
         "Gal H2O per kg-H2",
@@ -118,6 +124,7 @@ def run_h2_PEM(
     attributes = [
         rated_kWh_pr_kg,
         max_h2_pr_hr,
+        max_o2_pr_hr,
         max_pwr_pr_hr,
         elec_rated_h2_capacity_kgpy,
         gal_h20_pr_kg_h2,
@@ -188,6 +195,7 @@ def run_h2_PEM(
     H2_Results.update(dict(zip(life_desc, life_vals)))
     H2_Results.update({"Performance Schedules": pd.DataFrame(annual_avg_performance)})
     H2_Results.update({"Hydrogen Hourly Production [kg/hr]": hydrogen_hourly_production})
+    H2_Results.update({"Oxygen Hourly Production [kg/hr]": oxygen_hourly_production})
     H2_Results.update({"Water Hourly Consumption [kg/hr]": water_hourly_usage})
 
     if not debug_mode:

h2integrate/converters/hydrogen/test/test_pem_electrolyzer_performance.py

@@ -129,3 +129,62 @@ def test_electrolyzer_outputs(tech_config, plant_config, subtests):
         assert prob.get_val("comp.operational_life", units="yr") == plant_life
     with subtests.test("replacement_schedule value"):
         assert np.any(prob.get_val("comp.replacement_schedule", units="unitless") == 0)
+
+
+@pytest.mark.regression
+def test_electrolyzer_results(tech_config, plant_config, subtests):
+    prob = om.Problem()
+    comp = ECOElectrolyzerPerformanceModel(
+        plant_config=plant_config, tech_config=tech_config, driver_config={}
+    )
+    prob.model.add_subsystem("comp", comp, promotes=["*"])
+    prob.setup()
+    power_profile = np.full(8760, 32.0)
+    prob.set_val("comp.electricity_in", power_profile, units="MW")
+
+    prob.run_model()
+
+    with subtests.test("Total hydrogen produced"):
+        assert (
+            pytest.approx(5297814.89908964, rel=1e-6)
+            == prob.get_val("comp.hydrogen_out", units="kg/h").sum()
+        )
+
+    with subtests.test("Total oxygen produced"):
+        assert (
+            pytest.approx(42045916.90741967, rel=1e-6)
+            == prob.get_val("comp.oxygen_out", units="kg/h").sum()
+        )
+
+    with subtests.test("Year 0 capacity factor"):
+        assert (
+            pytest.approx(77.10460139, rel=1e-6)
+            == prob.get_val("comp.capacity_factor", units="percent")[0]
+        )
+
+    with subtests.test("Rated H2 production"):
+        assert pytest.approx(784.3544735823235, rel=1e-6) == prob.get_val(
+            "comp.rated_hydrogen_production", units="kg/h"
+        )
+
+    with subtests.test("Rated O2 production"):
+        assert pytest.approx(6225.00099576, rel=1e-6) == prob.get_val(
+            "comp.rated_oxygen_production", units="kg/h"
+        )
+
+    with subtests.test("H2: CF*Rated = Annual"):
+        np.testing.assert_allclose(
+            prob.get_val("comp.rated_hydrogen_production", units="kg/h")
+            * prob.get_val("comp.capacity_factor", units="unitless")
+            * 8760,
+            prob.get_val("comp.annual_hydrogen_produced", units="kg/yr"),
+            rtol=1e-6,
+        )
+    with subtests.test("O2: CF*Rated = Annual"):
+        np.testing.assert_allclose(
+            prob.get_val("comp.rated_oxygen_production", units="kg/h")
+            * prob.get_val("comp.capacity_factor", units="unitless")
+            * 8760,
+            prob.get_val("comp.annual_oxygen_produced", units="kg/yr"),
+            rtol=1e-6,
+        )

h2integrate/postprocess/test/test_sql_timeseries_to_csv.py

@@ -57,7 +57,7 @@ def test_save_csv_all_results(subtests, configuration, run_example_02_sql_fpath)
     res = save_case_timeseries_as_csv(run_example_02_sql_fpath, save_to_file=True)
 
     with subtests.test("Check number of columns"):
-        assert len(res.columns.to_list()) == 40
+        assert len(res.columns.to_list()) == 41
 
     with subtests.test("Check number of rows"):
         assert len(res) == 8760