diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 33778212c4c020fac1cc4c4d7b1c5754d405594d..2f68bdeee5eba51a87d0f22aff393f150a314ff7 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -109,6 +109,8 @@ void PhreeqcIO::initialize()
         [](auto result, auto const& chemical_system_index) {
             return result + chemical_system_index.size();
         });
+
+    _chemical_system->initialize(_num_chemical_systems);
 }
 
 void PhreeqcIO::executeInitialCalculation(
diff --git a/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp b/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b4d3932bd0ea3c04ae5f68c5cc01688fcfdfbe54
--- /dev/null
+++ b/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
@@ -0,0 +1,62 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "ChemicalSystem.h"
+#include "AqueousSolution.h"
+
+#include "MathLib/LinAlg/MatrixVectorTraits.h"
+#include "MathLib/LinAlg/UnifiedMatrixSetters.h"
+
+namespace ChemistryLib
+{
+namespace PhreeqcIOData
+{
+void ChemicalSystem::initialize(std::size_t const num_chemical_systems)
+{
+    aqueous_solution->pH =
+        MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+            num_chemical_systems);
+
+    aqueous_solution->pe->resize(num_chemical_systems);
+    std::fill(aqueous_solution->pe->begin(), aqueous_solution->pe->end(),
+              aqueous_solution->pe0);
+
+    auto& components = aqueous_solution->components;
+    for (auto& component : components)
+    {
+        component.amount =
+            MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
+                num_chemical_systems);
+    }
+
+    if (!kinetic_reactants.empty())
+    {
+        for (auto& kinetic_reactant : kinetic_reactants)
+        {
+            kinetic_reactant.amount->resize(num_chemical_systems);
+            std::fill(kinetic_reactant.amount->begin(),
+                      kinetic_reactant.amount->end(),
+                      kinetic_reactant.initial_amount);
+        }
+    }
+
+    if (!equilibrium_reactants.empty())
+    {
+        for (auto& equilibrium_reactant : equilibrium_reactants)
+        {
+            equilibrium_reactant.amount->resize(num_chemical_systems);
+            std::fill(equilibrium_reactant.amount->begin(),
+                      equilibrium_reactant.amount->end(),
+                      equilibrium_reactant.initial_amount);
+        }
+    }
+}
+}  // namespace PhreeqcIOData
+}  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/ChemicalSystem.h b/ChemistryLib/PhreeqcIOData/ChemicalSystem.h
index e06f41a65eedbe33ea8eab01b3923fa930bf0390..2650f7b3d584691ce3b21731263bab1bf316d485 100644
--- a/ChemistryLib/PhreeqcIOData/ChemicalSystem.h
+++ b/ChemistryLib/PhreeqcIOData/ChemicalSystem.h
@@ -32,6 +32,8 @@ struct ChemicalSystem
     {
     }
 
+    void initialize(std::size_t const num_chemical_systems);
+
     std::unique_ptr<AqueousSolution> aqueous_solution;
     std::vector<KineticReactant> kinetic_reactants;
     std::vector<EquilibriumReactant> equilibrium_reactants;
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index 9adf4f53a87883a6dad714d1b2316bff29cebabf..8afb0961cb452e9f3f8bea43f2033680d3941375 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -417,6 +417,11 @@ double TimeLoop::computeTimeStepping(const double prev_dt, double& t,
 /// initialize output, convergence criterion, etc.
 void TimeLoop::initialize()
 {
+    if (_chemical_solver_interface != nullptr)
+    {
+        _chemical_solver_interface->initialize();
+    }
+
     for (auto& process_data : _per_process_data)
     {
         auto& pcs = process_data->process;