diff --git a/ChemistryLib/ChemicalSolverInterface.h b/ChemistryLib/ChemicalSolverInterface.h
index b212a5e8f80037712a4c02eff58a8a449a475ecf..2a50dc04b593a7a27938bf61591a1cf626c3a069 100644
--- a/ChemistryLib/ChemicalSolverInterface.h
+++ b/ChemistryLib/ChemicalSolverInterface.h
@@ -12,6 +12,16 @@
 
 #include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
 
+namespace MaterialPropertyLib
+{
+class Medium;
+}
+
+namespace ParameterLib
+{
+class SpatialPosition;
+}
+
 namespace ChemistryLib
 {
 class ChemicalSolverInterface
@@ -19,6 +29,14 @@ class ChemicalSolverInterface
 public:
     virtual void initialize() {}
 
+    virtual void initializeChemicalSystemConcrete(
+        GlobalIndexType const& /*chemical_system_id*/,
+        MaterialPropertyLib::Medium const* /*medium*/,
+        ParameterLib::SpatialPosition const& /*pos*/,
+        double const /*t*/)
+    {
+    }
+
     virtual void executeInitialCalculation(
         std::vector<GlobalVector> const& interpolated_process_solutions) = 0;
 
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 69de25b11c7c2b75909262087f2f07e4ba39af25..038ff4a778af2d1b220570812fe8d2bd8b735ef5 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -20,6 +20,7 @@
 
 #include "BaseLib/Algorithm.h"
 #include "BaseLib/ConfigTreeUtil.h"
+#include "MaterialLib/MPL/Medium.h"
 #include "MathLib/LinAlg/Eigen/EigenVector.h"
 #include "MathLib/LinAlg/LinAlg.h"
 #include "MeshLib/Mesh.h"
@@ -48,6 +49,21 @@ std::ostream& operator<<(std::ostream& os,
               std::ostream_iterator<DataBlock>(os));
     return os;
 }
+
+template <typename Reactant>
+void setReactantAmount(Reactant& reactant,
+                       GlobalIndexType const& chemical_system_id,
+                       MaterialPropertyLib::Phase const& solid_phase,
+                       ParameterLib::SpatialPosition const& pos, double const t)
+{
+    auto const& solid_constituent = solid_phase.component(reactant.name);
+
+    auto const amount =
+        solid_constituent.property(MaterialPropertyLib::PropertyType::amount)
+            .template initialValue<double>(pos, t);
+
+    (*reactant.amount)[chemical_system_id] = amount;
+}
 }  // namespace
 
 PhreeqcIO::PhreeqcIO(std::string const project_file_name,
@@ -114,6 +130,25 @@ void PhreeqcIO::initialize()
     }
 }
 
+void PhreeqcIO::initializeChemicalSystemConcrete(
+    GlobalIndexType const& chemical_system_id,
+    MaterialPropertyLib::Medium const* medium,
+    ParameterLib::SpatialPosition const& pos, double const t)
+{
+    auto const& solid_phase = medium->phase("Solid");
+    for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
+    {
+        setReactantAmount(kinetic_reactant, chemical_system_id, solid_phase,
+                          pos, t);
+    }
+
+    for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
+    {
+        setReactantAmount(equilibrium_reactant, chemical_system_id, solid_phase,
+                          pos, t);
+    }
+}
+
 void PhreeqcIO::executeInitialCalculation(
     std::vector<GlobalVector> const& interpolated_process_solutions)
 {
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index dd403a7a7380804265768ba0b5bcb99fe607533a..f0eacf1b7300cd10f89f76f212fe4a9fe766fe15 100644
--- a/ChemistryLib/PhreeqcIO.h
+++ b/ChemistryLib/PhreeqcIO.h
@@ -46,6 +46,12 @@ public:
 
     void initialize() override;
 
+    void initializeChemicalSystemConcrete(
+        GlobalIndexType const& chemical_system_id,
+        MaterialPropertyLib::Medium const* medium,
+        ParameterLib::SpatialPosition const& pos,
+        double const t) override;
+
     void executeInitialCalculation(std::vector<GlobalVector> const&
                                        interpolated_process_solutions) override;
 
diff --git a/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp b/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
index 6cc79850a58b29768e647540c7ce174e8179acd0..669c57d3b88816b9e901e3fbde076554adc0dc88 100644
--- a/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
+++ b/ChemistryLib/PhreeqcIOData/ChemicalSystem.cpp
@@ -39,17 +39,11 @@ void ChemicalSystem::initialize(std::size_t const num_chemical_systems)
     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);
     }
 
     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