diff --git a/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp b/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
index 6a7286fb1c9012b5e6fd995076e727dd8aa12b7a..1382ba202359b5fd0fc9316c5a2da996de6b35f3 100644
--- a/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
+++ b/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
@@ -39,11 +39,13 @@ void ChemicalSystem::initialize(std::size_t const num_chemical_systems)
     for (auto& kinetic_reactant : kinetic_reactants)
     {
         kinetic_reactant.molality->resize(num_chemical_systems);
+        kinetic_reactant.volume_fraction->resize(num_chemical_systems);
     }
 
     for (auto& equilibrium_reactant : equilibrium_reactants)
     {
         equilibrium_reactant.molality->resize(num_chemical_systems);
+        equilibrium_reactant.volume_fraction->resize(num_chemical_systems);
     }
 }
 }  // namespace PhreeqcIOData
diff --git a/ChemistryLib/PhreeqcIOData/CreateEquilibriumReactants.cpp b/ChemistryLib/PhreeqcIOData/CreateEquilibriumReactants.cpp
index 2694bdd872d4272f898444ab6364197b88fda5f7..366c3e0a75015f31fbe4275b7fce8416c97cb46b 100644
--- a/ChemistryLib/PhreeqcIOData/CreateEquilibriumReactants.cpp
+++ b/ChemistryLib/PhreeqcIOData/CreateEquilibriumReactants.cpp
@@ -46,12 +46,15 @@ std::vector<EquilibriumReactant> createEquilibriumReactants(
         auto molality = MeshLib::getOrCreateMeshProperty<double>(
             mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
 
+        auto volume_fraction = MeshLib::getOrCreateMeshProperty<double>(
+            mesh, "phi_" + name, MeshLib::MeshItemType::IntegrationPoint, 1);
+
         auto mesh_prop_molality = MeshLib::getOrCreateMeshProperty<double>(
             mesh, name + "_avg", MeshLib::MeshItemType::Cell, 1);
         mesh_prop_molality->resize(mesh.getNumberOfElements());
 
         equilibrium_reactants.emplace_back(
-            std::move(name), molality, mesh_prop_molality, saturation_index);
+            std::move(name), molality, volume_fraction, mesh_prop_molality, saturation_index);
     }
 
     return equilibrium_reactants;
diff --git a/ChemistryLib/PhreeqcIOData/CreateKineticReactant.cpp b/ChemistryLib/PhreeqcIOData/CreateKineticReactant.cpp
index 042bb5d256ac511288a949b01f9685188dd61da9..75e23b633c537e95cd7d1d9a847dc53b68ad8576 100644
--- a/ChemistryLib/PhreeqcIOData/CreateKineticReactant.cpp
+++ b/ChemistryLib/PhreeqcIOData/CreateKineticReactant.cpp
@@ -54,6 +54,9 @@ std::vector<KineticReactant> createKineticReactants(
         auto molality = MeshLib::getOrCreateMeshProperty<double>(
             mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
 
+        auto volume_fraction = MeshLib::getOrCreateMeshProperty<double>(
+            mesh, "phi_" + name, MeshLib::MeshItemType::IntegrationPoint, 1);
+
         auto mesh_prop_molality = MeshLib::getOrCreateMeshProperty<double>(
             mesh, name + "_avg", MeshLib::MeshItemType::Cell, 1);
         mesh_prop_molality->resize(mesh.getNumberOfElements());
@@ -68,6 +71,7 @@ std::vector<KineticReactant> createKineticReactants(
         kinetic_reactants.emplace_back(std::move(name),
                                        std::move(chemical_formula),
                                        molality,
+                                       volume_fraction,
                                        mesh_prop_molality,
                                        std::move(parameters),
                                        fix_amount);
diff --git a/ChemistryLib/PhreeqcIOData/EquilibriumReactant.h b/ChemistryLib/PhreeqcIOData/EquilibriumReactant.h
index 03a1b7dd14b7b33b2fc7a9d393eacd5c8aa5df0a..198f5e672bc1ba315c43b683890b2fc70a089e44 100644
--- a/ChemistryLib/PhreeqcIOData/EquilibriumReactant.h
+++ b/ChemistryLib/PhreeqcIOData/EquilibriumReactant.h
@@ -30,10 +30,12 @@ struct EquilibriumReactant
 {
     EquilibriumReactant(std::string name_,
                         MeshLib::PropertyVector<double>* molality_,
+                        MeshLib::PropertyVector<double>* volume_fraction_,
                         MeshLib::PropertyVector<double>* mesh_prop_molality_,
                         double saturation_index_)
         : name(std::move(name_)),
           molality(molality_),
+          volume_fraction(volume_fraction_),
           mesh_prop_molality(mesh_prop_molality_),
           saturation_index(saturation_index_)
     {
@@ -43,6 +45,7 @@ struct EquilibriumReactant
 
     std::string const name;
     MeshLib::PropertyVector<double>* molality;
+    MeshLib::PropertyVector<double>* volume_fraction;
     MeshLib::PropertyVector<double>* mesh_prop_molality;
     double const saturation_index;
     static const ItemType item_type = ItemType::EquilibriumReactant;
diff --git a/ChemistryLib/PhreeqcIOData/KineticReactant.h b/ChemistryLib/PhreeqcIOData/KineticReactant.h
index 93dd3ef9fac7e15f38ee34e3af9dec9891e16123..ddfbd460524ab75f726dfc0e2edb40b9b65772fd 100644
--- a/ChemistryLib/PhreeqcIOData/KineticReactant.h
+++ b/ChemistryLib/PhreeqcIOData/KineticReactant.h
@@ -27,12 +27,14 @@ struct KineticReactant
     KineticReactant(std::string name_,
                     std::string chemical_formula_,
                     MeshLib::PropertyVector<double>* molality_,
+                    MeshLib::PropertyVector<double>* volume_fraction_,
                     MeshLib::PropertyVector<double>* mesh_prop_molality_,
                     std::vector<double>&& parameters_,
                     bool const fix_amount_)
         : name(std::move(name_)),
           chemical_formula(std::move(chemical_formula_)),
           molality(molality_),
+          volume_fraction(volume_fraction_),
           mesh_prop_molality(mesh_prop_molality_),
           parameters(std::move(parameters_)),
           fix_amount(fix_amount_)
@@ -44,6 +46,7 @@ struct KineticReactant
     std::string const name;
     std::string const chemical_formula;
     MeshLib::PropertyVector<double>* molality;
+    MeshLib::PropertyVector<double>* volume_fraction;
     MeshLib::PropertyVector<double>* mesh_prop_molality;
     std::vector<double> const parameters;
     bool const fix_amount;