The problem is formulated using an exchange coordinate κ which is independent of the reactor size, as shown in Equation 4: κ x = ∫ 0 x j A d x n ̇ 1 , (4) with A the species being exchanged, j A as the molar flux of species A from flow 1 to flow 2, n ̇ 1 as the molar flow rate of flow 1, and x the position along the reactor. In practical terms for redox cycles, κ would denote the change in the mole fraction of any reaction species along the reactor (an example provided in the Supplementary Material). The condition for a spontaneous transfer process of species A from flow 1 to flow 2 must fulfill Equation 5: μ A, 1 κ ≥ μ A, 2 κ ∀ κ ∈ 0 , κ total , (5) with μ being the chemical potential. Applying the conservation of mass, we note that the number of moles of A that have left flow 1 must equal to the amount of moles that have entered flow 2, at any point along the reactor. For the parallel flow (co-current) case this yields in Equation 6: κ 1 = κ 2 ≡ κ (6) while for counter-current flow the exchange coordinate of flow 2 is reversed, as shown in Equation 7 yielding κ 1 = κ 2 − κ total ≡ κ . (7)

The thermodynamic upper limit for the exchange of A in a parallel flow (co-current) system is given by Equation 8: μ A, 1 κ total = μ A, 2 κ total (8) with κ total equal to κ at the reactor outlet. For the counter-current case, the chemical potential can meet at each boundary, share a common tangent somewhere within the reactor, or have all of species A transferred from flow 1 to flow 2. These are given by μ A, 1 0 = μ A, 2 0  or μ A, 1 κ total = μ A, 2 κ total (9) and ∂ μ A, 1 ∂ κ = ∂ μ A, 2 ∂ κ  and μ A, 1 κ = μ A, 2 κ (10)

In the case of a complete transfer of A from flow 1 to flow 2, neither Equation 9 nor Equation 10 are binding, but the fundamental thermodynamic condition of Equation 5 is met within the entirety of the reactor length. Two system parameters that are needed to calculate the O 2 exchange during reduction are the ratio between the molar flow rates of both streams Equation 11: ω red = n ̇ 2 n ̇ 1 = n ̇ sg,in n ̇ MO , (11) and the O 2 mole fraction in the inert sweep gas x O 2 , in . n ̇ sg,in is the molar flow rate of the inert sweep gas at the reactor inlet and n ̇ MO is the redox-active material molar flow rate. Assuming an isothermal and isobaric reactor, the chemical potential becomes a function of the O 2 partial pressure p O 2 . The O 2 partial pressure at the reactor outlet is then calculated as a function of κ . This can be obtained from the relation between the O 2 partial pressure and the O 2 mole fraction at the reactor outlet, after the oxygen exchange occurred Equation 12: p O 2 , sg red κ = x O 2 , in ω red + κ ω red + κ p red (12) with x O 2 , in as the O 2 mole fraction at the reactor inlet (function of the sweep gas purity), κ as the total O 2 exchange, and p red as the total reactor pressure during reduction. Since each mole of O 2 released from the oxide will create two O vacancies in the oxide material, we can calculate the O 2 partial pressure of the oxide using Equation 13: p O 2 , MO red κ = f T , δ ox + 2 κ , (13) with δ ox as the initial non-stoichiometry extent at the end of oxidation and κ as the total O 2 exchange, same as in Equation 12. This formulation allows us to numerically solve for increasing κ from 0 to κ max , while checking the criteria in Equations 8–10 for each κ . The maximum O 2 transfer that is calculated without violating the constraints in Equations 8–10 is denoted as κ red , being the maximum species exchange during reduction.

During the oxidation, a flow of H 2 O or CO 2 enters the reactor and re-oxidizes the redox-active material, splitting the gaseous reactant into H 2 or CO , respectively. The conversion is usually incomplete, so the reactor outlet stream consists of a mixture of unconverted reactants and products. Following the same methodology as in the reduction reaction, we calculate the chemical potential of both gas and solid phases, to ensure the solution would adhere to requirements of a spontaneous process per Equation 5. The equilibrium constant for the thermolysis reaction ( H 2 O or CO 2 ) in the gas phase can be written as shown in Equation 14: K = p prod p O 2 p ° 1 2 p reac = x prod x reac p O 2 p ° 1 2 (14) with the reactant and product mole fractions at the reactor effluent given by x reac = ω ox x reac,in − 2 κ and x prod = ω ox x prod,in + 2 κ , respectively (and using p i = x i p ). K is the relevant thermolysis reaction equilibrium constant (i.e., H 2 O or CO 2 thermolysis) and ω ox is the ratio between the molar flow rates of both streams during oxidation ω ox = n ̇ feed / n ̇ MO with n ̇ feed as the total feed molar flow rate which is given by n ̇ feed = n ̇ reac,in + n ̇ prod,in . The equilibrium constant is calculated from K = exp − Δ G ° R T . Substituting the partial pressure equations into Equation 14, assuming that p O 2 ≪ p reac and p O 2 ≪ p prod , we can calculate the O 2 partial pressure in the oxidizer stream using Equation 15: p O 2 , feed ox κ = p ° K ω ox x reac,in − 2 κ ω ox x prod,in + 2 κ 2 . (15)

The O 2 partial pressure of the oxide material is calculated via Equation 16: p O 2 , MO ox κ = f T , δ red − 2 κ (16) considering that this time the oxide is absorbing O 2 molecules from the oxidizer gaseous stream. By solving in the same manner as for the reduction, increasing κ from 0 to κ max , we can find the maximum O 2 exchange during the oxidation step κ ox that satisfies the constraints in Equation 8, Equation 9, and Equation 10. Following the analysis of both steps, the maximum extent of re-oxidation (i.e., the lowest δ possible during oxidation at the reactor outlet) is calculated from Equation 17: δ ox = δ red − 2 κ ox (17) and the extent of reduction can be calculated as shown in Equation 18: Δ δ = δ red − δ ox = 2 κ ox . (18)

Assuming full selectivity, the conversion of H 2 O during water splitting oxidation can be calculated from Equation 19: X = n ̇ prod,out n ̇ reac,in = 2 κ ox ω ox (19) with n ̇ prod,out as the product molar flow rate at the outlet ( H 2 or CO ) and n ̇ reac,in as the reactant molar flow rate into the reactor ( H 2 O or CO 2 ). Since conversion is limited to unity, the maximum value of κ during oxidation is κ max = 0.5 ω ox .

To allow for a universal analysis of different redox materials, with different defect and vacancy formation mechanisms Zinkevich et al. (2006); Warren et al. (2022); Vieten et al. (2019); Wexler et al. (2023), we define the molar oxygen content fraction in the solid in Equation 20 as ϕ = n O n MO (20) with n O as the moles of O atoms in the solid and n MO the moles of redox-active material solid. The amount of O atoms released during oxidation is then Δ ϕ = ϕ ox − ϕ red . In such a manner, we can use the same formulation for materials such as ceria and its solid solutions, perovskites, iron aluminates, and others. The relation between κ and ϕ is then given per Equation 21: ϕ f = ϕ i ± 2 κ (21) with ϕ f as the final state ϕ , ϕ i as the initial state ϕ , and the sign dictated by the step, with minus for reduction (O atoms released from the solid) and plus for oxidation (O atoms absorbed by the solid). For materials with an oxygen vacancy mechanism, such as most redox-active materials studied to date ( CeO 2 and its solid solutions, perovskites, and most ferrites), the total amount of oxygen exchanged in a cycle is equal in each formulation Δ δ = Δ ϕ . More details are provided in the Supplementary Material.

The connection between the thermodynamics of a specific redox-active material and its performance in the redox cycle to an actual system is based on a previously developed model Lidor and Bulfin (2024) and is aligned with other works in the field Bulfin et al. (2021). However, one generalization presented here is the clear separation of heat and work requirements. While it is customary to convert the auxiliary work to heat equivalent in the field of solar thermochemistry Bulfin et al. (2021), there is increased interest in hybrid CST-PV systems that would be able to provide both heat and power. Moreover, the emerging field of high-temperature electric thermal energy storage (ETES) Ma et al. (2023) opens up the possibility of using renewable electricity to drive the thermochemical reaction. Hence, we have excluded the energy conversion penalties associated with the energy input source, whether via solar thermal or other sources. Instead, we have calculated the overall thermochemical conversion process efficiency based on the required heat and power per η = n ̇ fuel HHV fuel Q ̇ req + W ̇ req , (22) with η as the efficiency, n ̇ fuel as the fuel production rate, HHV as the fuel higher heating value (upper limit of a fuel thermal energy, defined as the heat released during combustion assuming condensing of products that are liquid at the initial fuel temperature and pressure), Q ̇ req as the required thermal power, and W ̇ req as the required electric power for all necessary unit operations such as pumping and separation. Equation 22 can be divided by the flow rate of the redox-active material to be turned into a general form (not tied to a specific size/flow rate): η = Δ ϕ HHV fuel q req + w req , (23) with Δ ϕ = ϕ ox − ϕ red as the specific amount of fuel produced per mole of redox-active material, q req as the specific required heat per mole of redox-active material, and w req as the specific required work per mole of redox-active material. From now on, we will use Equation 23 and the specific energy terms in this paper.

The following energy terms are considered in this work:

• Sensible heating of the redox-active material between oxidation and reduction temperatures ( q MO,h ) .

In addition, the effects of heat recovery, both in the gaseous streams (sweep gas and oxidizer gas) and in the solid stream, are evaluated by the implementation of heat recovery effectiveness values ε g and ε s , respectively. This topic has been of great importance due to the large sensible energy needed for temperature-swing operation (Lidor et al., 2023; Patankar et al., 2022; Lidor and Zimmermann, 2023). It is also important for sweep gas operated systems which require large flow rates (Hathaway et al., 2016; Lidor et al., 2021). The possibility of recovering the exothermic heat of oxidation q ox is also implemented in the same manner, using a heat recovery effectiveness ε ox . This heat can be used to meet any heating requirements at temperatures below T ox , such as oxidizer preheating or inert sweep gas heating to T ox . This heat can also drive an auxiliary power cycle that can be used to supply some of the required work. Since this model has been developed as a generic tool, capable of modeling a large variety of systems, the amount of usable exothermic heat is calculated assuming general heat recovery effectiveness, as well as heat-to-work efficiency in case this heat can cover all the thermal loads at or below T ox and still has some excess. The full details of calculating each term are presented in the Supplementary Material.

The feedstock conversion extent is given by X = 1 − n ̇ reac n ̇ reac,in (24) with n ̇ reac as the unreacted feedstock molar flow rate at the effluent and n ̇ reac,in as the feedstock molar flow rate into the reactor. Since we perform an analysis normalized by moles of redox-active material, the conversion is calculated in our model from Equation 25: X = Δ ϕ ω ox (25) which is in essence the moles of oxidizer that reacted per moles of redox-active material, over the moles of oxidizer fed into the reactor per moles of redox-active material.

The model has been developed in MATLAB Version 2023b The MathWorks Inc. (2023) and is publicly available on a GitHub repository (https://github.com/NREL/REDOTHERM), including all the main scripts and required functions. The repository includes thermodynamic functions for several redox-active materials: CeO 2 , Ce 0.8 Zr 0.2 O 2 (CeZr20), La 0.6 Ca 0.4 Mn 0.6 Al 0.4 O 3 (LCMA6464), La 0.6 Sr 0.4 MnO 3 (LSM40), and Fe 0.33 Al 0.67 O 4 (Fe33Al67). The thermodynamic data has been collected from (Bulfin et al., 2015; Bulfin et al., 2016; Carrillo and Scheffe, 2019; Warren et al., 2022), and uses either the provided functions for the reduction enthalpy Δ h red δ and reduction entropy Δ s red δ (Bulfin et al., 2015; Bulfin et al., 2016) or uses a curve fitting tool to develop polynomial correlations from the given data (Carrillo and Scheffe, 2019; Warren et al., 2022). We note that when using curve fitting methods to calculate the material thermodynamic properties, a certain degree of inaccuracy is expected, especially if the original experimental data is limited in terms of temperature and δ . Extrapolating beyond the range of the material data can be unreliable, and care should be taken in interpreting results from the sensitivity analysis and optimization capabilities described below when conditions extend beyond the limits of the measured material performance.

Thermodynamic properties of the fluids in the system are calculated using the CoolProp package (Bell et al., 2014) at the relevant process temperature and pressure, and the equilibrium composition of the oxidizer feed is calculated using the Gibbs energy minimization method implemented with Goodwin et al. (2023). The oxidizer gas can be selected as either H 2 O or CO 2 , as the necessary equilibrium constant calculation is performed for the relevant thermolysis reaction, based on the oxidizer feed of each analysis. The model is developed with an interactive input acquisition (via prompts) for the main parameters which are expected to be changed frequently, such as temperatures, and other parameters which are set within the code. There are three main scripts that are to be used when running the model with their respective code logic flowcharts in Figure 3:

• Main code: Redox_Countercurrent_Thermo_Main.m - this script is used to solve for a single data point (i.e., one combination of temperatures, pressures, molar flow rate ratios, etc.) and obtain the redox performance in terms of Δ δ , Δ ϕ , and X . It can also be used to perform parametric sweeps over combinations of T red and T ox and/or ω red and ω ox . The results can be saved in a file that can be post-processed for cycle performance and plotting via the support script Redox_Countercurrent_Thermo_Plot_Results.m. The logic flowchart of this script is presented in Figure 3a.

The reactor pressures during reduction ( p red ) and oxidation ( p ox ) are passed to the optimization problem as decision variables, but are currently fixed via constraints/bounds (depending on the optimizer used). The motivation behind this is two-fold: (1) to allow the analysis of hybrid vacuum-sweep gas systems, which requires both inert sweep gas separation as well as vacuum pumping work terms; (2) to allow future inclusion of any downstream pumping requirements, which in turn might promote the use of high pressure oxidation due to the benefits of pumping water, and obtaining high-pressure H 2 for storage or downstream processes such as Fischer-Tropsch.

Cryogenic separation was selected for inert gas separation. Pressure swing adsorption (PSA) was implemented for product separation for CO 2 , while H 2 O splitting includes options for conventional condensation and evaporation of the make-up feed and condensed water, as well as a mechanical vapor recompression (MVR) cycle separation which can recover the latent heat (based on Lidor, 2024).

All the different options for the REDOTHERM code are provided in Table 1. The full details of the calculation methods for each term are given in the Supplementary Material.

The parametric sweeps that are implemented in the main code (Redox_Countercurrent_Thermo_Main.m) are useful in investigating the effects of different operation and design parameters on all the different energy terms of the system as well as on its performance indicators. We present as an example the results of a parametric sweep over the molar flow rate ratios ω red and ω ox for CeO 2 for a case of H 2 O splitting. The temperatures are set at T red = 1550 °C and T ox = 900 °C. The inert sweep gas purity is x O 2 , in = 1 0 − 5 , and equilibrium H 2 mole fraction at inlet conditions is assumed ( x H 2 , in = 4.8 ⋅ 1 0 − 6 ) . The H 2 –  H 2 O separation method selected is MVR, and while the code analyzes the performance of both PF-PF and CF-CF configurations, we are presenting here only the CF-CF case for brevity. The trade-off between obtaining a higher Δ δ at the expense of conversion using higher gas molar flow rates (i.e., higher ω red and ω ox ) is clearly demonstrated in Figure 6. The regions which exhibit a value of X > 0.1 , which has been referred to in different studies (Bayon et al., 2022; Lidor and Bulfin, 2024) as the minimum viable conversion extent, barely overlaps with the areas of high Δ δ . Of course, this set of plots do not reveal the complete picture: having a Δ δ that is too low would reduce the efficiency since a small amount of fuel is produced.

We also calculate the system efficiency for these parametric sweep cases. The values used for the various heat recovery unit operations are ε s = 0.5 , ε g = 0.8 , and ε ox = 0.8 . The rest of the terms used in the calculation of the system efficiency, including specific energy and efficiency terms for some of the auxiliary units, are provided in the Supplementary Material. The results of the cycle efficiency analysis are presented in Figure 7 for the same sample case, using two types of product separation, conventional steam condensation (Figure 7a) and MVR-based separation (Figure 7b). A maximum efficiency of 11.2% and 13.93% is identified for condensation and MVR-based system, respectively. We note the low values compared to other predicted values from the literature. This is due to the following reasons: (a) we do not attempt to optimize the system in this example, except for the two parameters swept ( ω red and ω ox ); (b) some of our system parameter assumptions, such as heat recovery effectiveness, are more modest than other studies; (c) our analysis includes the correct thermodynamics limits, which are not implemented in most studies; and (d) we have accounted for the separation and steam generation energy terms, which are often neglected.

By examining the rest of the performance maps for the different energy fraction terms (sensible heating, etc.), one can obtain important insights into the specific thermochemical hydrogen (TCH) process that is evaluated, both from material and system perspective. The plots for this analysis are provided in the Supplementary Material.

The optimization capability of the REDOTHERM code includes a single case optimization, as well as optimization for multiple cases, in which the model sweeps over a range of values for the different system parameters, such as redox-active material type, heat recovery effectiveness, and so forth. A sample usage is given here, demonstrating how it can be used to compare different redox-active materials, under the same operating conditions and boundaries. The input parameters, as well as the ranges of the decision variables, are given in Table 2. The optimization problem is solved using a direct search method (patternsearch in MATLAB) with the Nonuniform Pattern Search (NUPS) algorithm.

Five values for ε s from 0 to 1 were included, examining all cases under both extremes, with no solid heat recovery up to an ideal (and nonphysical) full solid heat recovery. The analysis is performed for water splitting cycles and includes both condensation with reboiling and MVR as the product separation options. The analysis is demonstrated for atmospheric reduction (with sweep gas) and oxidation; however, the option to perform either step at different pressures is included as well. Performing the oxidation at elevated pressures can benefit from increased p O 2 in the oxidizer stream, which can increase the extent of re-oxidation (Tran et al., 2023). However, the major improvement is achieved from obtaining the produced H 2 at pressure, thus saving on the costly compression that is needed for any produced hydrogen. Lastly, we assume no thermal losses from the reactor ( f th,loss = 0 ) , since this analysis is performed for a generic moving oxide system without a specific reactor design.

The summary of the results from the optimization run is presented in Figure 8. We focus on a baseline case with ε s = 0.5 , an extreme case of no solid heat recovery ( ε s = 0 ) and an ideal full solid heat recovery ( ε s = 1 ) . In all of those cases the oxidation heat recovery was taken as ε ox = 0.8 . Note that the optimization in Figure 8 considers a wide range of reduction and oxidation temperature conditions, and extends beyond the range of available oxide material thermodynamic data in some cases. We present these results to demonstrate the capabilities of the model, but caution that results from the optimization should ideally be combined with sensitivity analyses to understand the behavior in the vicinity of the optimum and to verify that the material thermodynamic performance is reasonable within this region of the parameter space.

In Figure 8a we can see that the efficiency of CeO 2 is higher than all other materials. When considering the reduction enthalpy of the different materials, a clear trend is observed that materials with higher values of the enthalpy of oxygen vacancy formation exhibit higher efficiency values. As expected, the reduction temperature in all the optimal solutions converged to the upper bound (1700 °C), with less than 10 °C below this limit in all 80 solutions. This trend is identical to the results reported in (Li et al., 2018a). It implies the model can be modified to provide T red as an input, thus removing one decision variable from the optimization space, potentially lowering the required computational load. However, the temperature swing varied as can be seen from Figure 8b, with most materials converging to an optimal solution with Δ T = 800 °C while LCMA6464 benefitting from a large value of Δ T . When examining Figure 8f, one would expect that all materials will converge to the maximum Δ T possible since having a temperature swing does not require any sensible heating when ε s = 1 . However, we see that only CeO 2 exhibits the largest Δ T possible. While not included here for brevity, our hypothesis why materials with lower reduction enthalpies, that would benefit from a larger temperature swing, are limited in their Δ T is due to the high requirements of sensible heating of the oxidizer stream, that still has a value of ε g = 0.8 . This type of investigation, deeply examining all the intricate effects of the different system components and their effects on the performance for different redox materials, is what we envision REDOTHERM most suitable for. Another important finding is that CeO 2 and CeZr20 are the only materials that have a value of ω ox below 1; all other materials require excessive amounts of steam for re-oxidation, with up to ω ox = 14 . While known in the literature (Bayon et al., 2022; Lidor and Bulfin, 2024), this side-by-side comparison highlights the importance of examining all operating parameters and design considerations when comparing different redox-active materials.

Examining a specific subset of results, we present in Figure 9 the specific energy terms for each of the evaluated redox-active materials for the base case optimized solution (i.e., ε s = 0.5 and ε ox = 0.8 ). The energy terms are normalized by the amount of produced H 2 . When examining the case of using condensation as the H 2 - H 2 O separation method (Figure 9a), it can be seen that for CeO 2 and CeZr20 the specific reduction energy and specific sensible MO required heating are the dominant factors (with CeO 2 as the only case when the specific reduction energy is the largest term), while for the other materials the required heat for steam generation is dominant (the required product separation work is zero, since the load is only thermal). This supports the findings presented in Figure 8, with an inverse correlation between the ease of reduction (via the reduction enthalpy) and the performance. For the case of using MVR separation (Figure 9b) the value of the required sensible MO heat is the dominant factor across all the materials, except CeO 2 in which it is the second largest term, following the specific reduction energy (identical to the condensation case for CeO 2 ).

The example cases shown in this section provide a brief overview of the capabilities of REDOTHERM. Extending upon them for new redox-active materials, both real and hypothetical, as well as adding more technology options for the auxiliary operations, is relatively straightforward and could support the field in the search for efficient and scalable designs.

In this subsection, we demonstrate the use of the REDOTHERM repository in a few types of analysis that can benefit different research activities within the TCH and chemical looping redox fields.