diff --git a/docs/REFERENCE_MANUAL_INDEX.md b/docs/REFERENCE_MANUAL_INDEX.md
index 217ce11d6e..3a27cdad4a 100644
--- a/docs/REFERENCE_MANUAL_INDEX.md
+++ b/docs/REFERENCE_MANUAL_INDEX.md
@@ -159,6 +159,7 @@ NeqSim is distributed under the Apache-2.0 license and can be used via:
 | Mercury Thermodynamics    | [docs/thermo/mercury_thermodynamics.md](thermo/mercury_thermodynamics) | Mercury-focused SRK-TwuCoon-Statoil-EOS usage, TPflash setup, and thesis-linked BIP/correlation guidance |
 | **Søreide-Whitson Model** | [docs/thermo/SoreideWhitsonModel.md](thermo/SoreideWhitsonModel)   | **Gas solubility in brine** - Modified PR EoS with salinity effects, used in NeqSimLive for produced water emissions |
 | GERG-2008                 | [docs/thermo/gerg2008_eoscg.md](thermo/gerg2008_eoscg)             | GERG-2008, GERG-2008-H2, GERG-2008-NH3, and EOS-CG equations of state                                               |
+| Thermodynamic Benchmarks  | [docs/thermo/thermodynamic_benchmarks.md](thermo/thermodynamic_benchmarks) | Experimental-data benchmarks with provenance, uncertainty-aware metrics, and H2-CO2 phase-equilibrium data |
 | GERG-2008-NH3 Notebook    | [docs/examples/GERG2008_NH3_Ammonia_Properties.ipynb](https://github.com/equinor/neqsim/blob/master/docs/examples/GERG2008_NH3_Ammonia_Properties.ipynb) | Ammonia properties with Gao EOS — density validation, mixture properties, isotherms |
 | Mixing Rules              | [docs/thermo/mixing_rules_guide.md](thermo/mixing_rules_guide)     | Mixing rules and BIPs                                                                                                |
 | Mixing Rule Package       | [docs/thermo/mixingrule/README.md](thermo/mixingrule/)             | Mixing rule implementations                                                                                          |
diff --git a/docs/thermo/thermodynamic_benchmarks.md b/docs/thermo/thermodynamic_benchmarks.md
new file mode 100644
index 0000000000..86d2115cbd
--- /dev/null
+++ b/docs/thermo/thermodynamic_benchmarks.md
@@ -0,0 +1,99 @@
+---
+title: "Thermodynamic Model Benchmarks"
+description: "Auditable comparison of NeqSim property models with published experimental data."
+---
+
+# Thermodynamic model benchmarks
+
+NeqSim provides a reusable benchmark framework in
+`neqsim.thermo.util.benchmark.ThermodynamicBenchmark`. It keeps state,
+composition, property units, experimental uncertainty, citation, DOI, and data
+reuse information together with every comparison.
+
+The framework reports:
+
+- average absolute relative deviation (AARD);
+- signed bias;
+- root-mean-square relative error;
+- maximum absolute relative error;
+- uncertainty-normalized residuals when experimental uncertainty is available;
+- all point-level experimental and predicted values.
+
+## Hydrogen-containing CO2 benchmark
+
+`H2CO2PhaseEquilibriumData.load()` reads 24 bubble- and dew-point values from
+Tables IV and V of:
+
+> Zhang et al. (2026), *Measurement of phase equilibrium characteristics and
+> equation-of-state applicability for hydrogen-containing CO2 systems*,
+> International Journal of Fluid Engineering 3, 013903.
+> [doi:10.1063/5.0288386](https://doi.org/10.1063/5.0288386)
+
+The systems cover:
+
+- 96 mol% CO2 + 4 mol% H2;
+- 96 mol% CO2 + 2 mol% H2 + 2 mol% N2;
+- temperatures from -30 to 20 °C;
+- bubble pressures from 13.8 to 64.3 bara;
+- dew pressures from 65.3 to 112.1 bara.
+
+The article does not report pointwise pressure uncertainties. NeqSim records
+these as unavailable instead of assigning an unsupported uncertainty.
+
+## Java example
+
+```java
+ThermodynamicBenchmark.Dataset dataset = H2CO2PhaseEquilibriumData.load();
+NeqSimPhaseEquilibriumPrediction prediction =
+    new NeqSimPhaseEquilibriumPrediction(
+        NeqSimPhaseEquilibriumPrediction.Model.GERG_2008_H2);
+
+ThermodynamicBenchmark.Report report =
+    ThermodynamicBenchmark.run("GERG-2008-H2", dataset, prediction);
+
+double aardPercent = report.getAverageAbsoluteRelativeDeviationPercent();
+double biasPercent = report.getBiasPercent();
+```
+
+Supported configurations are SRK, PR, standard GERG-2008, and
+GERG-2008-H2. Cubic models use NeqSim database interaction parameters and
+mixing rule 2. The GERG-2008-H2 configuration enables the hydrogen-enhanced
+binary parameters and departure functions.
+
+The experimental pressure is used only as a numerical starting point. It is
+not returned by the prediction adapter and does not alter model parameters.
+
+## Constant H2-CO2 kij regression
+
+`BinaryInteractionParameterFitter` performs bounded golden-section regression
+of one constant cubic-EOS binary interaction parameter. Its objective is the
+unweighted mean squared relative pressure error. The calculation below used
+the 12 binary CO2-H2 points for calibration and retained the 12 ternary
+CO2-H2-N2 points as a holdout. CO2-N2 and H2-N2 interaction parameters remained
+at their NeqSim database values.
+
+| Model | Fitted H2-CO2 kij | Binary calibration RMSRE | Binary calibration AARD | Ternary holdout AARD | Database-kij binary AARD |
+| --- | ---: | ---: | ---: | ---: | ---: |
+| SRK | -0.47728 | 140.68% | 115.64% | 89.06% | 167.69% |
+| PR | -0.37523 | 120.26% | 92.87% | 88.03% | 142.00% |
+
+Regression bounds were -0.5 to 0.2, the parameter tolerance was `1e-4`, and
+the maximum was 20 objective evaluations. These are in-sample fitted values,
+not literature parameters.
+
+The fit reduces the selected objective, but the remaining errors and large
+negative parameters are unacceptable for a general NeqSim default. The fitted
+values are therefore documented for reproducibility but are not installed in
+the interaction-parameter database. This result indicates that one constant
+H2-CO2 kij in the current cubic-model setup cannot represent these bubble and
+dew boundaries adequately. Before considering a database change, investigate
+the point-level residuals, flash formulation, temperature dependence, and
+additional independent binary data.
+
+## Interpretation
+
+A low aggregate deviation is not sufficient evidence of universal model
+validity. Inspect bubble and dew results separately, examine signed residuals
+against temperature and composition, and retain the source validity range.
+Do not tune binary interaction parameters without recording their provenance
+and validating them on data not used for regression.
diff --git a/src/main/java/neqsim/thermo/util/benchmark/BinaryInteractionParameterFitter.java b/src/main/java/neqsim/thermo/util/benchmark/BinaryInteractionParameterFitter.java
new file mode 100644
index 0000000000..edad4f752a
--- /dev/null
+++ b/src/main/java/neqsim/thermo/util/benchmark/BinaryInteractionParameterFitter.java
@@ -0,0 +1,163 @@
+package neqsim.thermo.util.benchmark;
+
+import java.util.ArrayList;
+import java.util.List;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Dataset;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Point;
+
+/**
+ * Bounded one-dimensional regression of a constant cubic-EOS binary interaction parameter.
+ *
+ * <p>
+ * The objective is the mean squared relative error. This gives bubble- and dew-point pressures comparable influence
+ * even when their absolute pressure levels differ.
+ * </p>
+ */
+public final class BinaryInteractionParameterFitter {
+  /** Creates model predictions for a trial binary interaction parameter. */
+  public interface PredictionFactory {
+    /**
+     * Creates a prediction function for one trial parameter.
+     *
+     * @param binaryInteractionParameter trial dimensionless binary interaction parameter
+     * @return model prediction function
+     */
+    ThermodynamicBenchmark.Prediction create(double binaryInteractionParameter);
+  }
+
+  /** Immutable scalar-regression result. */
+  public static final class Result {
+    private final double binaryInteractionParameter;
+    private final double rootMeanSquareRelativeErrorPercent;
+    private final int objectiveEvaluations;
+
+    private Result(double binaryInteractionParameter, double rootMeanSquareRelativeErrorPercent,
+        int objectiveEvaluations) {
+      this.binaryInteractionParameter = binaryInteractionParameter;
+      this.rootMeanSquareRelativeErrorPercent = rootMeanSquareRelativeErrorPercent;
+      this.objectiveEvaluations = objectiveEvaluations;
+    }
+
+    /** @return fitted dimensionless binary interaction parameter */
+    public double getBinaryInteractionParameter() {
+      return binaryInteractionParameter;
+    }
+
+    /** @return in-sample root-mean-square relative error in percent */
+    public double getRootMeanSquareRelativeErrorPercent() {
+      return rootMeanSquareRelativeErrorPercent;
+    }
+
+    /** @return number of objective evaluations */
+    public int getObjectiveEvaluations() {
+      return objectiveEvaluations;
+    }
+  }
+
+  private final Dataset dataset;
+  private final PredictionFactory predictionFactory;
+
+  /**
+   * Creates a scalar binary-interaction-parameter regression.
+   *
+   * @param dataset experimental calibration dataset
+   * @param predictionFactory model adapter factory
+   */
+  public BinaryInteractionParameterFitter(Dataset dataset, PredictionFactory predictionFactory) {
+    if (dataset == null || predictionFactory == null) {
+      throw new IllegalArgumentException("Dataset and prediction factory are required");
+    }
+    this.dataset = dataset;
+    this.predictionFactory = predictionFactory;
+  }
+
+  /**
+   * Fits within a closed interval using golden-section minimization.
+   *
+   * @param lowerBound inclusive lower parameter bound
+   * @param upperBound inclusive upper parameter bound
+   * @param tolerance parameter-space convergence tolerance
+   * @param maximumEvaluations maximum objective evaluations
+   * @return fitted parameter and in-sample objective
+   * @throws Exception when a model prediction fails
+   */
+  public Result fit(double lowerBound, double upperBound, double tolerance, int maximumEvaluations) throws Exception {
+    if (!Double.isFinite(lowerBound) || !Double.isFinite(upperBound) || lowerBound >= upperBound
+        || !Double.isFinite(tolerance) || tolerance <= 0.0 || maximumEvaluations < 3) {
+      throw new IllegalArgumentException("Valid bounds, tolerance, and evaluation limit are required");
+    }
+
+    double left = lowerBound;
+    double right = upperBound;
+    double inverseGoldenRatio = (Math.sqrt(5.0) - 1.0) / 2.0;
+    double innerLeft = right - inverseGoldenRatio * (right - left);
+    double innerRight = left + inverseGoldenRatio * (right - left);
+    double objectiveLeft = objective(innerLeft);
+    double objectiveRight = objective(innerRight);
+    int evaluations = 2;
+
+    while (right - left > tolerance && evaluations < maximumEvaluations) {
+      if (objectiveLeft <= objectiveRight) {
+        right = innerRight;
+        innerRight = innerLeft;
+        objectiveRight = objectiveLeft;
+        innerLeft = right - inverseGoldenRatio * (right - left);
+        objectiveLeft = objective(innerLeft);
+      } else {
+        left = innerLeft;
+        innerLeft = innerRight;
+        objectiveLeft = objectiveRight;
+        innerRight = left + inverseGoldenRatio * (right - left);
+        objectiveRight = objective(innerRight);
+      }
+      evaluations++;
+    }
+
+    double fittedParameter;
+    double fittedObjective;
+    if (objectiveLeft <= objectiveRight) {
+      fittedParameter = innerLeft;
+      fittedObjective = objectiveLeft;
+    } else {
+      fittedParameter = innerRight;
+      fittedObjective = objectiveRight;
+    }
+    return new Result(fittedParameter, 100.0 * Math.sqrt(fittedObjective), evaluations);
+  }
+
+  /**
+   * Returns a dataset containing only selected point indexes.
+   *
+   * @param source source dataset
+   * @param indexes zero-based point indexes
+   * @param name descriptive subset name
+   * @return immutable dataset subset retaining source provenance
+   */
+  public static Dataset subset(Dataset source, List<Integer> indexes, String name) {
+    if (source == null || indexes == null || indexes.isEmpty() || name == null || name.trim().isEmpty()) {
+      throw new IllegalArgumentException("Source, indexes, and subset name are required");
+    }
+    List<Point> points = new ArrayList<Point>();
+    for (Integer index : indexes) {
+      if (index == null || index < 0 || index >= source.getPoints().size()) {
+        throw new IllegalArgumentException("Subset index is outside the source dataset");
+      }
+      points.add(source.getPoints().get(index));
+    }
+    return new Dataset(name, source.getCitation(), source.getDoi(), source.getLicense(), points);
+  }
+
+  private double objective(double binaryInteractionParameter) throws Exception {
+    ThermodynamicBenchmark.Prediction prediction = predictionFactory.create(binaryInteractionParameter);
+    double squaredRelativeErrorSum = 0.0;
+    for (Point point : dataset.getPoints()) {
+      double predictedValue = prediction.predict(point);
+      if (!Double.isFinite(predictedValue)) {
+        throw new IllegalStateException("Non-finite prediction during parameter regression");
+      }
+      double relativeError = (predictedValue - point.getExperimentalValue()) / point.getExperimentalValue();
+      squaredRelativeErrorSum += relativeError * relativeError;
+    }
+    return squaredRelativeErrorSum / dataset.getPoints().size();
+  }
+}
diff --git a/src/main/java/neqsim/thermo/util/benchmark/H2CO2PhaseEquilibriumData.java b/src/main/java/neqsim/thermo/util/benchmark/H2CO2PhaseEquilibriumData.java
new file mode 100644
index 0000000000..abd3458052
--- /dev/null
+++ b/src/main/java/neqsim/thermo/util/benchmark/H2CO2PhaseEquilibriumData.java
@@ -0,0 +1,90 @@
+package neqsim.thermo.util.benchmark;
+
+import java.io.BufferedReader;
+import java.io.IOException;
+import java.io.InputStream;
+import java.io.InputStreamReader;
+import java.nio.charset.StandardCharsets;
+import java.util.ArrayList;
+import java.util.LinkedHashMap;
+import java.util.List;
+import java.util.Map;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Dataset;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Point;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Property;
+
+/** Published H2-CO2 and H2-N2-CO2 phase-equilibrium benchmark datasets. */
+public final class H2CO2PhaseEquilibriumData {
+  private static final String RESOURCE = "/data/thermo/benchmark/zhang2026_h2_co2_phase_equilibrium.csv";
+  private static final String DOI = "10.1063/5.0288386";
+  private static final String CITATION = "Zhang et al. (2026), Measurement of phase equilibrium characteristics and "
+      + "equation-of-state applicability for hydrogen-containing CO2 systems, "
+      + "International Journal of Fluid Engineering 3, 013903.";
+  private static final String LICENSE = "Values transcribed from Tables IV-V of the cited open article; cite the original source.";
+
+  private H2CO2PhaseEquilibriumData() {
+  }
+
+  /**
+   * Loads all published binary and ternary bubble- and dew-point data.
+   *
+   * <p>
+   * The article does not report pointwise pressure uncertainty. The framework therefore records uncertainty as
+   * unavailable rather than inventing a value.
+   * </p>
+   *
+   * @return phase-equilibrium dataset containing 24 points
+   * @throws IOException if the packaged resource cannot be read
+   */
+  public static Dataset load() throws IOException {
+    InputStream input = H2CO2PhaseEquilibriumData.class.getResourceAsStream(RESOURCE);
+    if (input == null) {
+      throw new IOException("Missing benchmark resource " + RESOURCE);
+    }
+    List<Point> points = new ArrayList<Point>();
+    try (BufferedReader reader = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
+      String line = reader.readLine();
+      if (line == null) {
+        throw new IOException("Empty benchmark resource " + RESOURCE);
+      }
+      int lineNumber = 1;
+      while ((line = reader.readLine()) != null) {
+        lineNumber++;
+        if (line.trim().isEmpty() || line.trim().startsWith("#")) {
+          continue;
+        }
+        points.add(parseCsvRow(line, lineNumber));
+      }
+    }
+    return new Dataset("Zhang 2026 hydrogen-containing CO2 phase equilibrium", CITATION, DOI, LICENSE, points);
+  }
+
+  static Point parseCsvRow(String line, int lineNumber) throws IOException {
+    String[] values = line.split(",", -1);
+    if (values.length != 8) {
+      throw new IOException("Expected 8 CSV columns at line " + lineNumber + " but found " + values.length);
+    }
+    try {
+      double temperatureK = Double.parseDouble(values[1]) + 273.15;
+      double pressureBara = Double.parseDouble(values[3]) * 10.0;
+      Map<String, Double> composition = new LinkedHashMap<String, Double>();
+      composition.put("CO2", Double.parseDouble(values[4]));
+      composition.put("hydrogen", Double.parseDouble(values[5]));
+      double nitrogenFraction = Double.parseDouble(values[6]);
+      if (nitrogenFraction > 0.0) {
+        composition.put("nitrogen", nitrogenFraction);
+      }
+      Property property;
+      if ("bubble".equals(values[2])) {
+        property = Property.BUBBLE_POINT_PRESSURE;
+      } else if ("dew".equals(values[2])) {
+        property = Property.DEW_POINT_PRESSURE;
+      } else {
+        throw new IOException("Unknown phase-equilibrium property " + values[2] + " at line " + lineNumber);
+      }
+      return new Point(property, temperatureK, pressureBara, pressureBara, Double.NaN, "bara", composition);
+    } catch (NumberFormatException exception) {
+      throw new IOException("Malformed numeric value at line " + lineNumber + ": " + line, exception);
+    }
+  }
+}
diff --git a/src/main/java/neqsim/thermo/util/benchmark/NeqSimPhaseEquilibriumPrediction.java b/src/main/java/neqsim/thermo/util/benchmark/NeqSimPhaseEquilibriumPrediction.java
new file mode 100644
index 0000000000..61af5d56a8
--- /dev/null
+++ b/src/main/java/neqsim/thermo/util/benchmark/NeqSimPhaseEquilibriumPrediction.java
@@ -0,0 +1,110 @@
+package neqsim.thermo.util.benchmark;
+
+import java.util.Map;
+import neqsim.thermo.system.SystemGERG2008Eos;
+import neqsim.thermo.system.SystemInterface;
+import neqsim.thermo.system.SystemPrEos;
+import neqsim.thermo.system.SystemSrkEos;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Point;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Prediction;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Property;
+import neqsim.thermodynamicoperations.ThermodynamicOperations;
+
+/** NeqSim bubble- and dew-point prediction adapter for experimental benchmark points. */
+public final class NeqSimPhaseEquilibriumPrediction implements Prediction {
+  /** Supported equation-of-state configurations. */
+  public enum Model {
+    /** Soave-Redlich-Kwong with database interaction parameters. */
+    SRK,
+    /** Peng-Robinson with database interaction parameters. */
+    PR,
+    /** Standard GERG-2008. */
+    GERG_2008,
+    /** GERG-2008-H2 with hydrogen-enhanced binary parameters and departure functions. */
+    GERG_2008_H2
+  }
+
+  private final Model model;
+  private final double hydrogenCarbonDioxideKij;
+
+  /**
+   * Creates a phase-equilibrium prediction adapter.
+   *
+   * @param model equation-of-state configuration
+   */
+  public NeqSimPhaseEquilibriumPrediction(Model model) {
+    this(model, Double.NaN);
+  }
+
+  /**
+   * Creates a cubic-EOS phase-equilibrium prediction adapter with a custom constant H2-CO2 interaction parameter.
+   *
+   * @param model cubic equation-of-state configuration
+   * @param hydrogenCarbonDioxideKij dimensionless constant H2-CO2 binary interaction parameter
+   */
+  public NeqSimPhaseEquilibriumPrediction(Model model, double hydrogenCarbonDioxideKij) {
+    if (model == null) {
+      throw new IllegalArgumentException("Model is required");
+    }
+    if (Double.isFinite(hydrogenCarbonDioxideKij) && model != Model.SRK && model != Model.PR) {
+      throw new IllegalArgumentException("A custom kij is supported only for SRK and PR");
+    }
+    this.model = model;
+    this.hydrogenCarbonDioxideKij = hydrogenCarbonDioxideKij;
+  }
+
+  /** @return configured equation-of-state model */
+  public Model getModel() {
+    return model;
+  }
+
+  /** @return configured H2-CO2 interaction parameter, or NaN when database values are used */
+  public double getHydrogenCarbonDioxideKij() {
+    return hydrogenCarbonDioxideKij;
+  }
+
+  /** {@inheritDoc} */
+  @Override
+  public double predict(Point point) throws Exception {
+    if (point.getProperty() != Property.BUBBLE_POINT_PRESSURE && point.getProperty() != Property.DEW_POINT_PRESSURE) {
+      throw new IllegalArgumentException("Only bubble- and dew-point pressures are supported");
+    }
+    SystemInterface system = createSystem(point.getTemperatureK(), point.getPressureBara());
+    for (Map.Entry<String, Double> component : point.getComposition().entrySet()) {
+      system.addComponent(component.getKey(), component.getValue());
+    }
+    if (model == Model.SRK || model == Model.PR) {
+      system.createDatabase(true);
+      system.setMixingRule(2);
+      if (Double.isFinite(hydrogenCarbonDioxideKij)) {
+        system.setBinaryInteractionParameter("hydrogen", "CO2", hydrogenCarbonDioxideKij);
+      }
+    }
+    ThermodynamicOperations operations = new ThermodynamicOperations(system);
+    if (point.getProperty() == Property.BUBBLE_POINT_PRESSURE) {
+      operations.bubblePointPressureFlash(false);
+    } else {
+      operations.dewPointPressureFlash();
+    }
+    double pressure = system.getPressure(point.getUnit());
+    if (!Double.isFinite(pressure) || pressure <= 0.0) {
+      throw new IllegalStateException(
+          model + " returned invalid " + point.getProperty() + " at " + point.getTemperatureK() + " K");
+    }
+    return pressure;
+  }
+
+  private SystemInterface createSystem(double temperatureK, double pressureBara) {
+    if (model == Model.SRK) {
+      return new SystemSrkEos(temperatureK, pressureBara);
+    }
+    if (model == Model.PR) {
+      return new SystemPrEos(temperatureK, pressureBara);
+    }
+    SystemGERG2008Eos system = new SystemGERG2008Eos(temperatureK, pressureBara);
+    if (model == Model.GERG_2008_H2) {
+      system.useHydrogenEnhancedModel();
+    }
+    return system;
+  }
+}
diff --git a/src/main/java/neqsim/thermo/util/benchmark/ThermodynamicBenchmark.java b/src/main/java/neqsim/thermo/util/benchmark/ThermodynamicBenchmark.java
new file mode 100644
index 0000000000..a6744df564
--- /dev/null
+++ b/src/main/java/neqsim/thermo/util/benchmark/ThermodynamicBenchmark.java
@@ -0,0 +1,328 @@
+package neqsim.thermo.util.benchmark;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.LinkedHashMap;
+import java.util.List;
+import java.util.Map;
+
+/**
+ * Reusable, uncertainty-aware comparison of thermodynamic model predictions with experimental data.
+ *
+ * <p>
+ * Experimental values, units, composition, uncertainty, and provenance remain attached to every point so benchmark
+ * results are auditable and suitable for regression testing.
+ * </p>
+ */
+public final class ThermodynamicBenchmark {
+  private ThermodynamicBenchmark() {
+  }
+
+  /** Thermodynamic properties supported by the benchmark framework. */
+  public enum Property {
+    /** Bubble-point pressure. */
+    BUBBLE_POINT_PRESSURE,
+    /** Dew-point pressure. */
+    DEW_POINT_PRESSURE,
+    /** Mass density. */
+    DENSITY,
+    /** Speed of sound. */
+    SPEED_OF_SOUND,
+    /** Dynamic viscosity. */
+    VISCOSITY,
+    /** Thermal conductivity. */
+    THERMAL_CONDUCTIVITY,
+    /** Isobaric heat capacity. */
+    HEAT_CAPACITY_CP,
+    /** Water content. */
+    WATER_CONTENT
+  }
+
+  /** Immutable experimental benchmark point. */
+  public static final class Point {
+    private final Property property;
+    private final double temperatureK;
+    private final double pressureBara;
+    private final double experimentalValue;
+    private final double standardUncertainty;
+    private final String unit;
+    private final Map<String, Double> composition;
+
+    /**
+     * Creates an experimental point.
+     *
+     * @param property measured property
+     * @param temperatureK absolute temperature in K
+     * @param pressureBara absolute pressure in bara
+     * @param experimentalValue measured value in {@code unit}
+     * @param standardUncertainty standard uncertainty in {@code unit}, or {@link Double#NaN}
+     * @param unit unit of the measured value
+     * @param composition mole-fraction composition
+     */
+    public Point(Property property, double temperatureK, double pressureBara, double experimentalValue,
+        double standardUncertainty, String unit, Map<String, Double> composition) {
+      if (property == null || unit == null || unit.trim().isEmpty() || composition == null || composition.isEmpty()) {
+        throw new IllegalArgumentException("Property, unit, and composition are required");
+      }
+      if (!Double.isFinite(temperatureK) || temperatureK <= 0.0 || !Double.isFinite(pressureBara) || pressureBara <= 0.0
+          || !Double.isFinite(experimentalValue) || experimentalValue == 0.0) {
+        throw new IllegalArgumentException("Temperature, pressure, and experimental value must be physical and finite");
+      }
+      if (!Double.isNaN(standardUncertainty) && (!Double.isFinite(standardUncertainty) || standardUncertainty < 0.0)) {
+        throw new IllegalArgumentException("Standard uncertainty must be non-negative or NaN");
+      }
+      double compositionSum = 0.0;
+      Map<String, Double> compositionCopy = new LinkedHashMap<String, Double>();
+      for (Map.Entry<String, Double> entry : composition.entrySet()) {
+        if (entry.getKey() == null || entry.getKey().trim().isEmpty() || entry.getValue() == null
+            || !Double.isFinite(entry.getValue()) || entry.getValue() < 0.0) {
+          throw new IllegalArgumentException("Composition entries must be named and non-negative");
+        }
+        compositionCopy.put(entry.getKey(), entry.getValue());
+        compositionSum += entry.getValue();
+      }
+      if (Math.abs(compositionSum - 1.0) > 1.0e-8) {
+        throw new IllegalArgumentException("Mole fractions must sum to one");
+      }
+      this.property = property;
+      this.temperatureK = temperatureK;
+      this.pressureBara = pressureBara;
+      this.experimentalValue = experimentalValue;
+      this.standardUncertainty = standardUncertainty;
+      this.unit = unit;
+      this.composition = Collections.unmodifiableMap(compositionCopy);
+    }
+
+    /** @return measured property */
+    public Property getProperty() {
+      return property;
+    }
+
+    /** @return absolute temperature in K */
+    public double getTemperatureK() {
+      return temperatureK;
+    }
+
+    /** @return absolute pressure in bara */
+    public double getPressureBara() {
+      return pressureBara;
+    }
+
+    /** @return experimental value */
+    public double getExperimentalValue() {
+      return experimentalValue;
+    }
+
+    /** @return reported standard uncertainty, or NaN when unavailable */
+    public double getStandardUncertainty() {
+      return standardUncertainty;
+    }
+
+    /** @return property unit */
+    public String getUnit() {
+      return unit;
+    }
+
+    /** @return immutable mole-fraction composition */
+    public Map<String, Double> getComposition() {
+      return composition;
+    }
+  }
+
+  /** Immutable experimental dataset with source provenance. */
+  public static final class Dataset {
+    private final String name;
+    private final String citation;
+    private final String doi;
+    private final String license;
+    private final List<Point> points;
+
+    /**
+     * Creates a dataset.
+     *
+     * @param name dataset name
+     * @param citation full source citation
+     * @param doi source DOI
+     * @param license reuse statement for the encoded data
+     * @param points experimental points
+     */
+    public Dataset(String name, String citation, String doi, String license, List<Point> points) {
+      if (name == null || name.trim().isEmpty() || citation == null || doi == null || license == null || points == null
+          || points.isEmpty()) {
+        throw new IllegalArgumentException("Dataset metadata and points are required");
+      }
+      this.name = name;
+      this.citation = citation;
+      this.doi = doi;
+      this.license = license;
+      this.points = Collections.unmodifiableList(new ArrayList<Point>(points));
+    }
+
+    /** @return dataset name */
+    public String getName() {
+      return name;
+    }
+
+    /** @return source citation */
+    public String getCitation() {
+      return citation;
+    }
+
+    /** @return source DOI */
+    public String getDoi() {
+      return doi;
+    }
+
+    /** @return data reuse statement */
+    public String getLicense() {
+      return license;
+    }
+
+    /** @return immutable experimental points */
+    public List<Point> getPoints() {
+      return points;
+    }
+  }
+
+  /** Supplies a model prediction for one experimental point. */
+  public interface Prediction {
+    /**
+     * Predicts the property represented by a point.
+     *
+     * @param point experimental point defining state and composition
+     * @return predicted value in the point's unit
+     * @throws Exception when the model cannot calculate the point
+     */
+    double predict(Point point) throws Exception;
+  }
+
+  /** Result for one experimental point. */
+  public static final class Row {
+    private final Point point;
+    private final double predictedValue;
+    private final double signedRelativeErrorPercent;
+    private final double uncertaintyNormalizedResidual;
+
+    private Row(Point point, double predictedValue) {
+      this.point = point;
+      this.predictedValue = predictedValue;
+      this.signedRelativeErrorPercent = 100.0 * (predictedValue - point.getExperimentalValue())
+          / point.getExperimentalValue();
+      double uncertainty = point.getStandardUncertainty();
+      this.uncertaintyNormalizedResidual = Double.isFinite(uncertainty) && uncertainty > 0.0
+          ? (predictedValue - point.getExperimentalValue()) / uncertainty
+          : Double.NaN;
+    }
+
+    /** @return experimental point */
+    public Point getPoint() {
+      return point;
+    }
+
+    /** @return predicted value */
+    public double getPredictedValue() {
+      return predictedValue;
+    }
+
+    /** @return signed relative error in percent */
+    public double getSignedRelativeErrorPercent() {
+      return signedRelativeErrorPercent;
+    }
+
+    /** @return residual divided by standard uncertainty, or NaN when unavailable */
+    public double getUncertaintyNormalizedResidual() {
+      return uncertaintyNormalizedResidual;
+    }
+  }
+
+  /** Immutable aggregate benchmark report. */
+  public static final class Report {
+    private final String modelName;
+    private final Dataset dataset;
+    private final List<Row> rows;
+    private final double averageAbsoluteRelativeDeviationPercent;
+    private final double biasPercent;
+    private final double rootMeanSquareRelativeErrorPercent;
+    private final double maximumAbsoluteRelativeErrorPercent;
+
+    private Report(String modelName, Dataset dataset, List<Row> rows) {
+      this.modelName = modelName;
+      this.dataset = dataset;
+      this.rows = Collections.unmodifiableList(new ArrayList<Row>(rows));
+      double absoluteErrorSum = 0.0;
+      double signedErrorSum = 0.0;
+      double squaredErrorSum = 0.0;
+      double maximumError = 0.0;
+      for (Row row : rows) {
+        double error = row.getSignedRelativeErrorPercent();
+        absoluteErrorSum += Math.abs(error);
+        signedErrorSum += error;
+        squaredErrorSum += error * error;
+        maximumError = Math.max(maximumError, Math.abs(error));
+      }
+      this.averageAbsoluteRelativeDeviationPercent = absoluteErrorSum / rows.size();
+      this.biasPercent = signedErrorSum / rows.size();
+      this.rootMeanSquareRelativeErrorPercent = Math.sqrt(squaredErrorSum / rows.size());
+      this.maximumAbsoluteRelativeErrorPercent = maximumError;
+    }
+
+    /** @return model name */
+    public String getModelName() {
+      return modelName;
+    }
+
+    /** @return source dataset */
+    public Dataset getDataset() {
+      return dataset;
+    }
+
+    /** @return immutable point results */
+    public List<Row> getRows() {
+      return rows;
+    }
+
+    /** @return average absolute relative deviation in percent */
+    public double getAverageAbsoluteRelativeDeviationPercent() {
+      return averageAbsoluteRelativeDeviationPercent;
+    }
+
+    /** @return mean signed relative error in percent */
+    public double getBiasPercent() {
+      return biasPercent;
+    }
+
+    /** @return root mean square relative error in percent */
+    public double getRootMeanSquareRelativeErrorPercent() {
+      return rootMeanSquareRelativeErrorPercent;
+    }
+
+    /** @return maximum absolute relative error in percent */
+    public double getMaximumAbsoluteRelativeErrorPercent() {
+      return maximumAbsoluteRelativeErrorPercent;
+    }
+  }
+
+  /**
+   * Runs a model over all points in a dataset.
+   *
+   * @param modelName auditable model name
+   * @param dataset experimental dataset
+   * @param prediction prediction implementation
+   * @return aggregate report
+   * @throws Exception when any prediction fails or is non-finite
+   */
+  public static Report run(String modelName, Dataset dataset, Prediction prediction) throws Exception {
+    if (modelName == null || modelName.trim().isEmpty() || dataset == null || prediction == null) {
+      throw new IllegalArgumentException("Model name, dataset, and prediction are required");
+    }
+    List<Row> rows = new ArrayList<Row>();
+    for (Point point : dataset.getPoints()) {
+      double predictedValue = prediction.predict(point);
+      if (!Double.isFinite(predictedValue)) {
+        throw new IllegalStateException("Non-finite prediction for " + point.getProperty());
+      }
+      rows.add(new Row(point, predictedValue));
+    }
+    return new Report(modelName, dataset, rows);
+  }
+}
diff --git a/src/main/resources/data/thermo/benchmark/zhang2026_h2_co2_phase_equilibrium.csv b/src/main/resources/data/thermo/benchmark/zhang2026_h2_co2_phase_equilibrium.csv
new file mode 100644
index 0000000000..1cae18319c
--- /dev/null
+++ b/src/main/resources/data/thermo/benchmark/zhang2026_h2_co2_phase_equilibrium.csv
@@ -0,0 +1,25 @@
+system,temperature_C,property,pressure_MPa,x_CO2,x_H2,x_N2,source_table
+CO2_H2_96_4,-30,bubble,1.45,0.96,0.04,0.00,IV
+CO2_H2_96_4,-30,dew,11.21,0.96,0.04,0.00,IV
+CO2_H2_96_4,-20,bubble,2.14,0.96,0.04,0.00,IV
+CO2_H2_96_4,-20,dew,9.72,0.96,0.04,0.00,IV
+CO2_H2_96_4,-10,bubble,2.76,0.96,0.04,0.00,IV
+CO2_H2_96_4,-10,dew,9.11,0.96,0.04,0.00,IV
+CO2_H2_96_4,0,bubble,3.65,0.96,0.04,0.00,IV
+CO2_H2_96_4,0,dew,8.87,0.96,0.04,0.00,IV
+CO2_H2_96_4,10,bubble,4.85,0.96,0.04,0.00,IV
+CO2_H2_96_4,10,dew,8.61,0.96,0.04,0.00,IV
+CO2_H2_96_4,20,bubble,6.43,0.96,0.04,0.00,IV
+CO2_H2_96_4,20,dew,8.51,0.96,0.04,0.00,IV
+CO2_H2_N2_96_2_2,-30,bubble,1.38,0.96,0.02,0.02,V
+CO2_H2_N2_96_2_2,-30,dew,6.67,0.96,0.02,0.02,V
+CO2_H2_N2_96_2_2,-20,bubble,1.95,0.96,0.02,0.02,V
+CO2_H2_N2_96_2_2,-20,dew,6.53,0.96,0.02,0.02,V
+CO2_H2_N2_96_2_2,-10,bubble,2.86,0.96,0.02,0.02,V
+CO2_H2_N2_96_2_2,-10,dew,6.56,0.96,0.02,0.02,V
+CO2_H2_N2_96_2_2,0,bubble,3.97,0.96,0.02,0.02,V
+CO2_H2_N2_96_2_2,0,dew,6.85,0.96,0.02,0.02,V
+CO2_H2_N2_96_2_2,10,bubble,4.92,0.96,0.02,0.02,V
+CO2_H2_N2_96_2_2,10,dew,7.21,0.96,0.02,0.02,V
+CO2_H2_N2_96_2_2,20,bubble,6.28,0.96,0.02,0.02,V
+CO2_H2_N2_96_2_2,20,dew,8.09,0.96,0.02,0.02,V
diff --git a/src/test/java/neqsim/thermo/util/benchmark/ThermodynamicBenchmarkTest.java b/src/test/java/neqsim/thermo/util/benchmark/ThermodynamicBenchmarkTest.java
new file mode 100644
index 0000000000..9a56ed33fd
--- /dev/null
+++ b/src/test/java/neqsim/thermo/util/benchmark/ThermodynamicBenchmarkTest.java
@@ -0,0 +1,147 @@
+package neqsim.thermo.util.benchmark;
+
+import static org.junit.jupiter.api.Assertions.assertEquals;
+import static org.junit.jupiter.api.Assertions.assertThrows;
+import static org.junit.jupiter.api.Assertions.assertTrue;
+
+import java.io.IOException;
+import java.util.LinkedHashMap;
+import java.util.Map;
+import org.junit.jupiter.api.Test;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Dataset;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Point;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Property;
+import neqsim.thermo.util.benchmark.ThermodynamicBenchmark.Report;
+
+/** Tests for reusable thermodynamic benchmark infrastructure and published H2-CO2 data. */
+class ThermodynamicBenchmarkTest {
+  @Test
+  void loadsPublishedH2CO2DatasetWithProvenance() throws Exception {
+    Dataset dataset = H2CO2PhaseEquilibriumData.load();
+
+    assertEquals(24, dataset.getPoints().size());
+    assertEquals("10.1063/5.0288386", dataset.getDoi());
+    assertTrue(dataset.getCitation().contains("Zhang"));
+
+    Point firstPoint = dataset.getPoints().get(0);
+    assertEquals(Property.BUBBLE_POINT_PRESSURE, firstPoint.getProperty());
+    assertEquals(243.15, firstPoint.getTemperatureK(), 1.0e-12);
+    assertEquals(14.5, firstPoint.getExperimentalValue(), 1.0e-12);
+    assertEquals(0.96, firstPoint.getComposition().get("CO2"), 1.0e-12);
+    assertEquals(0.04, firstPoint.getComposition().get("hydrogen"), 1.0e-12);
+
+    Point finalPoint = dataset.getPoints().get(23);
+    assertEquals(Property.DEW_POINT_PRESSURE, finalPoint.getProperty());
+    assertEquals(80.9, finalPoint.getExperimentalValue(), 1.0e-12);
+    assertEquals(0.02, finalPoint.getComposition().get("nitrogen"), 1.0e-12);
+  }
+
+  @Test
+  void reportsMalformedNumericCsvRowsAsIoFailures() {
+    IOException exception = assertThrows(IOException.class,
+        () -> H2CO2PhaseEquilibriumData.parseCsvRow("binary,bad,bubble,1.45,0.96,0.04,0.0,bara", 12));
+
+    assertTrue(exception.getMessage().contains("line 12"));
+    assertTrue(exception.getCause() instanceof NumberFormatException);
+  }
+
+  @Test
+  void calculatesAardBiasRmsAndMaximumError() throws Exception {
+    Dataset dataset = H2CO2PhaseEquilibriumData.load();
+    Report report = ThermodynamicBenchmark.run("synthetic +2 percent", dataset,
+        point -> point.getExperimentalValue() * 1.02);
+
+    assertEquals(2.0, report.getAverageAbsoluteRelativeDeviationPercent(), 1.0e-12);
+    assertEquals(2.0, report.getBiasPercent(), 1.0e-12);
+    assertEquals(2.0, report.getRootMeanSquareRelativeErrorPercent(), 1.0e-12);
+    assertEquals(2.0, report.getMaximumAbsoluteRelativeErrorPercent(), 1.0e-12);
+    assertEquals(24, report.getRows().size());
+  }
+
+  @Test
+  void calculatesUncertaintyNormalizedResidualWhenAvailable() throws Exception {
+    Map<String, Double> composition = new LinkedHashMap<String, Double>();
+    composition.put("CO2", 0.96);
+    composition.put("hydrogen", 0.04);
+    Point point = new Point(Property.BUBBLE_POINT_PRESSURE, 273.15, 36.5, 36.5, 0.5, "bara", composition);
+    Dataset dataset = new Dataset("uncertainty test", "test citation", "10.0000/test", "test data",
+        java.util.Collections.singletonList(point));
+
+    Report report = ThermodynamicBenchmark.run("test model", dataset, pointValue -> 37.5);
+
+    assertEquals(2.0, report.getRows().get(0).getUncertaintyNormalizedResidual(), 1.0e-12);
+  }
+
+  @Test
+  void rejectsCompositionThatDoesNotSumToOne() {
+    Map<String, Double> composition = new LinkedHashMap<String, Double>();
+    composition.put("CO2", 0.90);
+    composition.put("hydrogen", 0.04);
+
+    assertThrows(IllegalArgumentException.class,
+        () -> new Point(Property.BUBBLE_POINT_PRESSURE, 273.15, 36.5, 36.5, Double.NaN, "bara", composition));
+  }
+
+  @Test
+  void rejectsInvalidStateValueAndUncertainty() {
+    Map<String, Double> composition = new LinkedHashMap<String, Double>();
+    composition.put("CO2", 0.96);
+    composition.put("hydrogen", 0.04);
+
+    assertThrows(IllegalArgumentException.class,
+        () -> new Point(Property.BUBBLE_POINT_PRESSURE, 273.15, 0.0, 36.5, Double.NaN, "bara", composition));
+    assertThrows(IllegalArgumentException.class,
+        () -> new Point(Property.BUBBLE_POINT_PRESSURE, 273.15, 36.5, 0.0, Double.NaN, "bara", composition));
+    assertThrows(IllegalArgumentException.class,
+        () -> new Point(Property.BUBBLE_POINT_PRESSURE, 273.15, 36.5, 36.5, -0.1, "bara", composition));
+  }
+
+  @Test
+  void exposesConfiguredNeqSimModel() {
+    NeqSimPhaseEquilibriumPrediction prediction = new NeqSimPhaseEquilibriumPrediction(
+        NeqSimPhaseEquilibriumPrediction.Model.GERG_2008_H2);
+
+    assertEquals(NeqSimPhaseEquilibriumPrediction.Model.GERG_2008_H2, prediction.getModel());
+  }
+
+  @Test
+  void returnsPredictionInPointPressureUnit() throws Exception {
+    Point source = H2CO2PhaseEquilibriumData.load().getPoints().get(0);
+    Point pressureInMegapascal = new Point(source.getProperty(), source.getTemperatureK(), source.getPressureBara(),
+        source.getExperimentalValue() / 10.0, source.getStandardUncertainty(), "MPa", source.getComposition());
+    NeqSimPhaseEquilibriumPrediction prediction = new NeqSimPhaseEquilibriumPrediction(
+        NeqSimPhaseEquilibriumPrediction.Model.SRK);
+
+    double predictionBara = prediction.predict(source);
+    double predictionMegapascal = prediction.predict(pressureInMegapascal);
+
+    assertTrue(Double.isFinite(predictionMegapascal));
+    assertEquals(predictionBara, predictionMegapascal * 10.0, 1.0e-8);
+  }
+
+  @Test
+  void fitsScalarObjectiveWithinBounds() throws Exception {
+    Dataset dataset = H2CO2PhaseEquilibriumData.load();
+    BinaryInteractionParameterFitter fitter = new BinaryInteractionParameterFitter(dataset,
+        parameter -> point -> point.getExperimentalValue() * (1.0 + parameter - 0.125));
+
+    BinaryInteractionParameterFitter.Result result = fitter.fit(-0.2, 0.3, 1.0e-8, 50);
+
+    assertEquals(0.125, result.getBinaryInteractionParameter(), 1.0e-6);
+    assertTrue(result.getRootMeanSquareRelativeErrorPercent() < 1.0e-5);
+  }
+
+  @Test
+  void rejectsCustomKijForNonCubicModel() {
+    assertThrows(IllegalArgumentException.class,
+        () -> new NeqSimPhaseEquilibriumPrediction(NeqSimPhaseEquilibriumPrediction.Model.GERG_2008_H2, 0.1));
+  }
+
+  @Test
+  void exposesCustomKijForCubicModel() {
+    NeqSimPhaseEquilibriumPrediction prediction = new NeqSimPhaseEquilibriumPrediction(
+        NeqSimPhaseEquilibriumPrediction.Model.PR, -0.1);
+
+    assertEquals(-0.1, prediction.getHydrogenCarbonDioxideKij(), 1.0e-12);
+  }
+}