diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 9936ce309ef89a2961fd69fb3c8221dafc7ac95d..6a3580b6de1d8ce1a7dfc890c3ddad1d2174384e 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -19,6 +19,23 @@
 
 namespace ChemistryLib
 {
+void PhreeqcIO::doWaterChemistryCalculation(
+    std::vector<GlobalVector*>& process_solutions, double const dt)
+{
+    setAqueousSolutionsOrUpdateProcessSolutions(
+        process_solutions, Status::SettingAqueousSolutions);
+    setTimeStep(dt);
+
+    writeInputsToFile();
+
+    execute();
+
+    readOutputsFromFile();
+
+    setAqueousSolutionsOrUpdateProcessSolutions(
+        process_solutions, Status::UpdatingProcessSolutions);
+}
+
 void PhreeqcIO::setAqueousSolutionsOrUpdateProcessSolutions(
     std::vector<GlobalVector*> const& process_solutions, Status const status)
 {
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index 230a0ba16b2f33f60b702e7971baca9f5bd57bb4..10e431d4bac321f953b642c459d5749a353ab01d 100644
--- a/ChemistryLib/PhreeqcIO.h
+++ b/ChemistryLib/PhreeqcIO.h
@@ -49,6 +49,9 @@ public:
     {
     }
 
+    void doWaterChemistryCalculation(
+        std::vector<GlobalVector*>& process_solutions, double const dt);
+
     void setAqueousSolutionsOrUpdateProcessSolutions(
         std::vector<GlobalVector*> const& process_solutions,
         Status const status);
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index efd338064c41049e84609ea894d55f7a48ca85da..ba4850a7fc03c263b445da3183bc325e809f5140 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -725,6 +725,16 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
             timestep_id, t);
     }
 
+    if (_chemical_system != nullptr)
+    {
+        // Sequential non-iterative approach applied here to perform water
+        // chemistry calculation followed by resolving component transport
+        // process.
+        // TODO: move into a global loop to consider both mass balance over
+        // space and localized chemical equilibrium between solutes.
+        _chemical_system->doWaterChemistryCalculation(_process_solutions, dt);
+    }
+
     int process_id = 0;
     for (auto& process_data : _per_process_data)
     {