diff --git a/ChemistryLib/ChemicalSolverInterface.h b/ChemistryLib/ChemicalSolverInterface.h
index 21d53e9fd140591667e8aaec3cb7517ccecc76a0..f1d691ca3f73dc55a9cba4b39f6079920a493a4c 100644
--- a/ChemistryLib/ChemicalSolverInterface.h
+++ b/ChemistryLib/ChemicalSolverInterface.h
@@ -33,6 +33,12 @@ public:
         return {};
     }
 
+    virtual void computeSecondaryVariable(
+        std::size_t const /*ele_id*/,
+        std::vector<GlobalIndexType> const& /*chemical_system_indices*/)
+    {
+    }
+
     virtual ~ChemicalSolverInterface() = default;
 
 public:
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 45af3e26944053b7e5f13436ec4f3d2e0dfd90c8..8149ec0a32e31890e2071ecc6a422731cf61ce64 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -530,5 +530,35 @@ std::vector<std::string> const PhreeqcIO::getComponentList() const
 
     return component_names;
 }
+
+void PhreeqcIO::computeSecondaryVariable(
+    std::size_t const ele_id,
+    std::vector<GlobalIndexType> const& chemical_system_indices)
+{
+    auto const num_chemical_systems = chemical_system_indices.size();
+    for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
+    {
+        double amount_avg = std::accumulate(
+            chemical_system_indices.begin(), chemical_system_indices.end(), 0.0,
+            [&](double amount, auto const& id) {
+                return amount + (*kinetic_reactant.amount)[id];
+            });
+
+        amount_avg /= num_chemical_systems;
+        (*kinetic_reactant.amount_avg)[ele_id] = amount_avg;
+    }
+
+    for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
+    {
+        double amount_avg = std::accumulate(
+            chemical_system_indices.begin(), chemical_system_indices.end(), 0.0,
+            [&](double amount, auto const& id) {
+                return amount + (*equilibrium_reactant.amount)[id];
+            });
+
+        amount_avg /= num_chemical_systems;
+        (*equilibrium_reactant.amount_avg)[ele_id] = amount_avg;
+    }
+}
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index 6df10d8efbe8ba66468bc6fca197b0af6629fd12..fb7b20ef63dccc11fe599dde2160ea169b613266 100644
--- a/ChemistryLib/PhreeqcIO.h
+++ b/ChemistryLib/PhreeqcIO.h
@@ -69,6 +69,10 @@ public:
 
     friend std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io);
 
+    void computeSecondaryVariable(
+        std::size_t const ele_id,
+        std::vector<GlobalIndexType> const& chemical_system_indices) override;
+
     std::vector<std::string> const getComponentList() const override;
 
     std::string const _phreeqc_input_file;
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 957e5c8bbc8108ab7a0638f2b054e377fc635801..27bb0d6c9744d81cd1d35e19219ec2df2bc944bd 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -980,11 +980,24 @@ public:
         auto const ele_velocity_mat =
             MathLib::toMatrix(ele_velocity, GlobalDim, n_integration_points);
 
-        auto ele_id = _element.getID();
+        auto const ele_id = _element.getID();
         Eigen::Map<LocalVectorType>(
             &(*_process_data.mesh_prop_velocity)[ele_id * GlobalDim],
             GlobalDim) =
             ele_velocity_mat.rowwise().sum() / n_integration_points;
+
+        if (_process_data.chemical_solver_interface)
+        {
+            std::vector<GlobalIndexType> chemical_system_indices;
+            chemical_system_indices.reserve(n_integration_points);
+            std::transform(
+                _ip_data.begin(), _ip_data.end(),
+                std::back_inserter(chemical_system_indices),
+                [](auto const& ip_data) { return ip_data.chemical_system_id; });
+
+            _process_data.chemical_solver_interface->computeSecondaryVariable(
+                ele_id, chemical_system_indices);
+        }
     }
 
     void postTimestepConcrete(Eigen::VectorXd const& /*local_x*/,