From Trajectories to Timescales: Our Framework for Hydrogen-Based Direct Reduction

The Late-Stage Bottleneck In Hydrogen-Based Direct Reduction We Set Out to Resolve

Our team at the Max Planck Institute for Sustainable Materials has been working on a stubborn problem in green steelmaking. Hydrogen-based direct reduction works thermodynamically, but the kinetics collapse at high conversion. The last twenty percent of reduction—from wustite to iron—takes most of the time. That late tail drives up hydrogen consumption, extends residence times, and erodes the sustainability advantage of the entire process.

The classical kinetic literature tells you that mixed control exists. It does not tell you at which conversion the system switches from reaction to diffusion control. It does not tell you which design variable—temperature, pellet diameter, porosity, tortuosity, or composition—actually governs the late tail. And it certainly does not give you a constitutive map that works across different pellet composition families without refitting.

We decided to approach this backwards.

Our Inverse Kinetic Strategy to Tackle Hydrogen-Based Direct Reduction for Green Streel Making

Instead of fitting a forward shrinking-core model with prescribed rate constants, we start from the measured reduction trajectory t(X)t(X)—time as a function of conversion. We then deconvolute each trajectory into three additive contributions using the Sohn additive-time basis:

t ( X ) = τ M g M ( X ) + τ K g K ( X ) + τ D g D ( X ) t(X)=τM​gM​(X)+τK​gK​(X)+τD​gD​(X)

with the fixed dimensionless basis functions:

g M ( X ) = X , g K ( X ) = 1 − ( 1 − X ) 1 / 3 , g D ( X ) = 1 − 3 ( 1 − X ) 1 / 3 + 2 ( 1 − X ) gM​(X)=X,gK​(X)=1−(1−X)1/3,gD​(X)=1−3(1−X)1/3+2(1−X)

The ττ values are non-negative timescales extracted by non-negative least squares. They are not fitted rate constants. They are effective pellet-scale descriptors of how much each mechanism—external mass transfer, interfacial reaction, internal diffusion—contributes to the total time budget.

But a single triplet per experiment is too coarse. The pore network evolves. The product layer grows. Sintering closes pathways. So we split each trajectory at X=0.60X=0.60 into early and late stages, fitting separate ττ triplets while preserving continuity of t(X)t(X). This two-stage decomposition reduces reconstruction error systematically compared to a single-stage fit (median NRMSE improvement concentrated below 0.01) and outperforms the random pore model across 179 experiments.

The SCAM Framework: Imposing Physics on Machine Learning in Hydrogen-Based Direct Reduction for Green Streel Making

The inferred τKτK​ and τDτD​ timescales become our target variables. We then ask: how do these effective timescales depend on temperature, hydrogen partial pressure, pellet diameter, porosity, tortuosity, and composition?

We do not answer this with a black-box surrogate. Instead, we developed a scientifically constrained additive model (SCAM). The functional form is:

log ⁡ 10 τ = β 0 + ∑ j f j ( z j ) + ∑ ( j , k ) ∈ I f j k ( z j , z k ) + γ stage + θ log10​τ=β0​+j∑​fj​(zj​)+(j,k)∈I∑​fjk​(zj​,zk​)+γstage​+θ

where each fjfj​ is a smooth univariate spline with an explicit smoothness penalty, and fjkfjk​ are bivariate interaction surfaces admitted only when they survive cross-validation and bootstrap stability tests.

The constraints are the core innovation. We enforce monotonicity directly at the coefficient level where the metallurgy demands it. For example, τDτD​ must increase with Dpellet/dporeDpellet​/dpore​ and decrease with ϵ/τtortϵ/τtort​ over the sampled envelope. Non-monotonicity is permitted only where competing effects—such as sintering accelerating at high temperature while intrinsic diffusivity also increases—can plausibly generate a maximum or minimum. This is not post-hoc regularization. It is built into the basis representation using shape-constrained splines.

We implemented the SCAM using a penalized B-spline basis with second-difference roughness penalties. Monotonicity constraints are encoded as linear inequality constraints on the spline coefficients, solved via quadratic programming. The stage indicator γstageγstage​ is a binary term that allows early and late windows to have different baseline timescales while sharing the same shape-constrained response surfaces.

What the Constitutive Maps Reveal

Our fitted maps show a clean separation of control.

For τKτK​, temperature dominates. The partial dependence on 1/T1/T is approximately linear over 700–1000°C, consistent with an Arrhenius-like effective activation energy. Hydrogen partial pressure provides a secondary modulation. Pellet composition—compressed into a principal component PC₁ that increases with Fe content and decreases with oxide additions—shows a non-monotone dependence. That non-linearity is real, not a fitting artifact. It reflects competing pathways: certain gangue components can either enhance reducibility by creating microcracking or hinder it by forming stable mixed oxides such as iron silicates.

For τDτD​, the story is different. Pellet architecture governs. Porosity ϵϵ has the strongest single effect: increasing ϵϵ from 0.2 to 0.4 reduces log⁡10τDlog10​τD​ by approximately 0.6 decades. Tortuosity τtortτtort​ and pellet diameter DpelletDpellet​ provide additional, mechanistically distinct leverage. Temperature shows only a weak residual effect on τDτD​ once porosity and tortuosity are accounted for—because the dominant temperature influence on internal transport is indirect, mediated by sintering and pore coarsening, not by intrinsic diffusivity alone.

The interaction surfaces in Figure 4 tell a more nuanced story. For τKτK​, temperature and pH2pH2​​ couple strongly: the benefit of increasing hydrogen potential is largest at lower temperatures, where the interfacial reaction step remains rate-limiting over a wider conversion span. At higher temperatures, the incremental gain from more hydrogen diminishes because the bottleneck shifts toward product-layer development. For τDτD​, the strongest coupling is architectural: the increase of τDτD​ with pellet diameter is much steeper at low porosity. That is physically consistent—a large pellet with a dense, poorly connected pore network forces oxygen to diffuse through long, tortuous paths through the product layer.

Uncertainty Quantification and Validation Protocol

We do not present these maps as deterministic. Every step of our pipeline—trajectory smoothing, non-negative timescale extraction, SCAM fitting, regime boundary calculation—is propagated through bootstrap resampling at the experiment level. For each of 500 bootstrap replicates, we refit the entire pipeline and recompute all partial dependencies, interaction surfaces, and dominance probabilities.

The validation strategy goes beyond simple train-test splits. We perform composition-family holdouts. The pellet composition space is partitioned into six domains defined by basicity (CaO/SiO₂) and total Fe content. We train on five families, predict on the sixth, and repeat. This tests whether our ττ laws generalize to unseen pellet chemistries without refitting.

The results are reported in Table 1. For log⁡10τKlog10​τK​, the family-out median RMSE is 0.315, median MAE 0.282. Empirical coverage of nominal 90% prediction intervals is 0.918—slightly conservative, which is acceptable. For log⁡10τDlog10​τD​, the corresponding RMSE is 0.410, coverage 0.898. These errors are small relative to the envelope-wide spread in the targets (approximately 1.5 decades for τKτK​, 2 decades for τDτD​). The physics-only version of the model—excluding the composition correction term—performs worse on unseen families, confirming that the residual composition term captures genuine chemistry effects not mediated by the measured pore descriptors.

We also benchmarked SCAM against a lattice-based surrogate: a cubic grid interpolant in the same predictor space. The out-of-fold RMSE for log⁡10τKlog10​τK​ was 0.296 (SCAM) versus 0.317 (lattice). Regime agreement (diffusion-dominant vs reaction-dominant at X=0.8) was 0.864 with Cohen's κ=0.835κ=0.835. SCAM is not vastly more accurate—the lattice works reasonably well—but it is interpretable in a way the lattice is not. You can read the monotonicity, the interaction signs, the shape of the partial dependencies. That matters for industrial uptake.

Pathway Decomposition: Direct Versus Mediated Composition Effects

A question that often comes up in discussions with pellet producers: does composition matter directly, or only through the pore network it creates? We answered this with a causal mediation-style decomposition, adapted for our observational data.

We define the natural direct effect (NDE) as the change in log⁡10τlog10​τ when composition changes but the measured architecture descriptors—ϵϵ, dporedpore​, τtortτtort​—are held fixed at their values under the reference composition. This captures composition effects not transmitted through the pore network. The natural indirect effect (NIE) is the change transmitted through shifts in those architecture descriptors.

Figure 5 shows the result. For τDτD​, the NIE dominates, especially in the late conversion window. Pellet composition affects late-stage transport primarily by reshaping the pore network. For τKτK​, a larger residual direct effect persists. That means matching porosity and tortuosity alone will not erase composition differences in interfacial reducibility. Something else—phase assembly, local oxygen potential, impurity segregation at the reaction front—is at play.

Symbolic Regression for Reduced-Order Laws

Interpretable maps are good. Compact equations are better for engineering use. We performed symbolic regression on the SCAM-predicted timescales, restricted to physically dimensionless groups: Tref/(Tgas+273.15)Tref​/(Tgas​+273.15), pH2/prefpH2​​/pref​, Dpellet/dporeDpellet​/dpore​, and ϵ/τtortϵ/τtort​. The operator set was limited to addition, multiplication, division, and integer powers. Dimensionally inconsistent forms were gated out.

The Pareto front (complexity versus loss) gave a knee solution at complexity 11 for τKτK​ and 14 for τDτD​. The final forms are shown in Table 2. They consist of a physics-only symbolic law plus a bounded linear composition correction in basicity, Fe content, and two gauge coordinates. The composition correction is applied after the physics term—it is a residual adjustment, not a separate mechanism.

For example, the τDτD​ law reads:

log ⁡ 10 τ D = 1.588 ( T ref T + 273.15 ) 3 + ( p H 2 p ref ) 2 + 0.715 ( 0.932 − 0.215 B − 0.012 Fe − 0.02 G a + 0.034 G b ) log10​τD​=1.588(T+273.15Tref​​)3+(pref​pH2​​​)2+0.715(0.932−0.215B−0.012Fe−0.02Ga​+0.034Gb​)

with Tref=1148.15 KTref​=1148.15K, pref=1 barpref​=1bar. The cubic temperature term and quadratic hydrogen term emerged from symbolic search, not from prior assumption. They happen to fit the data better than Arrhenius or power-law forms over our limited temperature range. We do not overinterpret them as fundamental laws.

Conversion-Resolved Regime Maps

The ultimate synthesis is the regime map in Figure 7. We define the dominance ratio Λ=(ΔtD+ΔtM)/ΔtKΛ=(ΔtD​+ΔtM​)/ΔtK​ over a specified conversion window—here X=0.8X=0.8 for the late-stage comparison. Transport-dominant conditions are those with Λ>1Λ>1. The color scale shows P(Λ>1)P(Λ>1) from the bootstrap ensemble.

Two insights emerge. First, increasing the pellet-to-pore ratio Dpellet/dporeDpellet​/dpore​ pushes the system toward diffusion dominance regardless of temperature. That is a robust geometric scaling. Second, the effective diffusivity metric ϵ/τtortϵ/τtort​ acts as a master lever. At ϵ/τtort≈0.2ϵ/τtort​≈0.2, most of the envelope is diffusion-dominant. At ϵ/τtort≈0.5ϵ/τtort​≈0.5, reaction control takes over for small pellets. The boundary shifts systematically with conversion—at X=0.3X=0.3 it is mostly reaction-controlled, at X=0.5X=0.5 diffusion has already expanded substantially.

Where We Draw the Line

Our framework is validated within the sampled experimental envelope: T=700T=700–1000\,^\circ\text{C}, pH2=0.5pH2​​=0.5–1.0 bar1.0bar, Dpellet=4Dpellet​=4–18 mm18mm, porosity 0.2–0.5, tortuosity 2–8, pore size 10–50 μmμm. We do not claim extrapolation beyond these bounds. The composition family holdouts test interpolation to unseen chemistries within the same descriptor space, not extrapolation to entirely different ore types.

We also do not resolve the full time-dependent evolution of pore topology, interfacial area, or crack networks. Our τKτK​ and τDτD​ are effective pellet-scale summaries. They incorporate the integrated effect of microstructural evolution but do not provide evolution equations. That would require in situ structural measurements coupled to the reduction trajectory—a next step we are actively pursuing.

What This Means for Green Steelmaking

For the plant operator: temperature and hydrogen are early-stage levers. Once you are past 60% conversion, they matter less. Pellet architecture—high porosity, low tortuosity, small pellet-to-pore ratio—is what kills the late tail.

For the pellet designer: composition matters primarily through the pore network it creates. If you want to change late-stage kinetics, you do not need to change chemistry arbitrarily. You need to design the pore network. The regime map tells you where the diffusion-controlled region begins. Stay in the reaction-controlled region by keeping Dpellet/dporeDpellet​/dpore​ low and ϵ/τtortϵ/τtort​ high.

For the reactor modeler: the symbolic laws in Table 2 and the regime maps in Figure 7 provide experimentally anchored, conversion-resolved constitutive relations. They are not ab initio rate constants. But they are measurable, interpretable, and transferable within the stated bounds. You can plug them into shaft furnace models as pellet-scale closures.

Our code and data are open (https://github.com/bajpaianurag/Hydrogen_Reduction_Reaction_Transition_ML). We invite the community to test these constitutive maps on their own pellet datasets, within the validated envelope, and to push the boundaries where our current framework stops.