diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 3cda89519f8f9f99a53e67db287beaaee2c926f3..e1dfd16138dfaa739a170eca054a376fd3e68efc 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -373,13 +373,21 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
     auto const& surface = phreeqc_io._surface;
     int const num_skipped_lines = surface.empty() ? 1 : 2;
 
+    auto& kinetic_reactants = phreeqc_io._chemical_system->kinetic_reactants;
+    for (auto& kinetic_reactant : kinetic_reactants)
+    {
+        std::fill(kinetic_reactant.amount_avg->begin(),
+                  kinetic_reactant.amount_avg->end(), 0.);
+    }
+
     for (std::size_t mesh_item_id = 0;
          mesh_item_id < phreeqc_io.chemical_system_index_map.size();
          ++mesh_item_id)
     {
         auto const& local_chemical_system =
             phreeqc_io.chemical_system_index_map[mesh_item_id];
-        for (std::size_t ip = 0; ip < local_chemical_system.size(); ++ip)
+        auto const num_local_chemical_system = local_chemical_system.size();
+        for (std::size_t ip = 0; ip < num_local_chemical_system; ++ip)
         {
             auto const& chemical_system_id = local_chemical_system[ip];
 
@@ -443,8 +451,6 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
             auto& components = aqueous_solution->components;
             auto& equilibrium_reactants =
                 phreeqc_io._chemical_system->equilibrium_reactants;
-            auto& kinetic_reactants =
-                phreeqc_io._chemical_system->kinetic_reactants;
             auto& user_punch = phreeqc_io._user_punch;
             for (int item_id = 0;
                  item_id < static_cast<int>(accepted_items.size());
@@ -508,6 +514,8 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
                                 "'.");
                         (*kinetic_reactant.amount)[chemical_system_id] =
                             accepted_items[item_id];
+                        (*kinetic_reactant.amount_avg)[mesh_item_id] +=
+                            accepted_items[item_id];
                         break;
                     }
                     case ItemType::SecondaryVariable:
@@ -528,6 +536,12 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
                 }
             }
         }
+
+        for (auto& kinetic_reactant : kinetic_reactants)
+        {
+            (*kinetic_reactant.amount_avg)[mesh_item_id] /=
+                num_local_chemical_system;
+        }
     }
 
     return in;
diff --git a/ChemistryLib/PhreeqcIOData/CreateKineticReactant.cpp b/ChemistryLib/PhreeqcIOData/CreateKineticReactant.cpp
index 1c0b122f4331ba82df1668297951a6019eb47b66..bd39ebe6db2dfe8391fb1dd42e799fa1591a760b 100644
--- a/ChemistryLib/PhreeqcIOData/CreateKineticReactant.cpp
+++ b/ChemistryLib/PhreeqcIOData/CreateKineticReactant.cpp
@@ -57,6 +57,10 @@ std::vector<KineticReactant> createKineticReactants(
         auto amount = MeshLib::getOrCreateMeshProperty<double>(
             mesh, name, MeshLib::MeshItemType::IntegrationPoint, 1);
 
+        auto mesh_prop_amount = MeshLib::getOrCreateMeshProperty<double>(
+            mesh, name + "_avg", MeshLib::MeshItemType::Cell, 1);
+        mesh_prop_amount->resize(mesh.getNumberOfElements());
+
         if (chemical_formula.empty() && fix_amount)
         {
             OGS_FATAL(
@@ -67,6 +71,7 @@ std::vector<KineticReactant> createKineticReactants(
         kinetic_reactants.emplace_back(std::move(name),
                                        std::move(chemical_formula),
                                        amount,
+                                       mesh_prop_amount,
                                        initial_amount,
                                        std::move(parameters),
                                        fix_amount);
diff --git a/ChemistryLib/PhreeqcIOData/KineticReactant.h b/ChemistryLib/PhreeqcIOData/KineticReactant.h
index 1141683e84df79c70e414b6bbd91045b1e2cae64..14cd01dcc3d5ed7e92cb6c940a8e5fb9c167e256 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>* amount_,
+                    MeshLib::PropertyVector<double>* amount_avg_,
                     double const initial_amount_,
                     std::vector<double>&& parameters_,
                     bool const fix_amount_)
         : name(std::move(name_)),
           chemical_formula(std::move(chemical_formula_)),
           amount(amount_),
+          amount_avg(amount_avg_),
           initial_amount(initial_amount_),
           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>* amount;
+    MeshLib::PropertyVector<double>* amount_avg;
     double const initial_amount;
     std::vector<double> const parameters;
     bool const fix_amount;