diff --git a/ChemistryLib/ChemicalSolverInterface.h b/ChemistryLib/ChemicalSolverInterface.h
index 79443449e6540fafb40ef956058c3a44cc1c47a4..492e88e79c033055c18e3ea91455ec3fac2b7d3a 100644
--- a/ChemistryLib/ChemicalSolverInterface.h
+++ b/ChemistryLib/ChemicalSolverInterface.h
@@ -10,6 +10,7 @@
 
 #pragma once
 
+#include "MaterialLib/MPL/VariableType.h"
 #include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
 
 namespace MaterialPropertyLib
@@ -39,7 +40,11 @@ public:
 
     virtual void setChemicalSystemConcrete(
         std::vector<double> const& /*concentrations*/,
-        GlobalIndexType const& /*chemical_system_id*/)
+        GlobalIndexType const& /*chemical_system_id*/,
+        MaterialPropertyLib::Medium const* /*medium*/,
+        MaterialPropertyLib::VariableArray const& /*vars*/,
+        ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+        double const /*dt*/)
     {
     }
 
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index be287e4f6632b3ee83b944124f2f217d02d0ba4d..9ecc4807c5de52e7d9ab5a14a1270f178c088e39 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -116,6 +116,41 @@ void initializeReactantMolality(Reactant& reactant,
     }
 }
 
+template <typename Reactant>
+void setReactantMolality(Reactant& reactant,
+                         GlobalIndexType const& chemical_system_id,
+                         MaterialPropertyLib::Medium const* medium,
+                         MaterialPropertyLib::VariableArray const& vars,
+                         ParameterLib::SpatialPosition const& pos,
+                         double const t, double const dt)
+{
+    auto const& solid_phase = medium->phase("Solid");
+    auto const& solid_constituent = solid_phase.component(reactant.name);
+
+    if (!solid_constituent.hasProperty(
+            MaterialPropertyLib::PropertyType::molality))
+    {
+        auto const volume_fraction =
+            (*reactant.volume_fraction)[chemical_system_id];
+
+        auto const& liquid_phase = medium->phase("AqueousLiquid");
+        auto const fluid_density =
+            liquid_phase.property(MaterialPropertyLib::PropertyType::density)
+                .template value<double>(vars, pos, t, dt);
+
+        auto const porosity = std::get<double>(
+            vars[static_cast<int>(MaterialPropertyLib::Variable::porosity)]);
+
+        auto const molar_volume =
+            solid_constituent
+                .property(MaterialPropertyLib::PropertyType::molar_volume)
+                .template value<double>(vars, pos, t, dt);
+
+        (*reactant.molality)[chemical_system_id] =
+            volume_fraction / fluid_density / porosity / molar_volume;
+    }
+}
+
 template <typename Reactant>
 static double averageReactantMolality(
     Reactant const& reactant,
@@ -219,10 +254,25 @@ void PhreeqcIO::initializeChemicalSystemConcrete(
 
 void PhreeqcIO::setChemicalSystemConcrete(
     std::vector<double> const& concentrations,
-    GlobalIndexType const& chemical_system_id)
+    GlobalIndexType const& chemical_system_id,
+    MaterialPropertyLib::Medium const* medium,
+    MaterialPropertyLib::VariableArray const& vars,
+    ParameterLib::SpatialPosition const& pos, double const t, double const dt)
 {
     setAqueousSolution(concentrations, chemical_system_id,
                        *_chemical_system->aqueous_solution);
+
+    for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
+    {
+        setReactantMolality(kinetic_reactant, chemical_system_id, medium, vars,
+                            pos, t, dt);
+    }
+
+    for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
+    {
+        setReactantMolality(equilibrium_reactant, chemical_system_id, medium,
+                            vars, pos, t, dt);
+    }
 }
 
 void PhreeqcIO::executeInitialCalculation()
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index be513a73635ee3b4b666cc57a64dfa6de462b1d4..376f8aaa712bb07d0b2b75c58b85578a5719033a 100644
--- a/ChemistryLib/PhreeqcIO.h
+++ b/ChemistryLib/PhreeqcIO.h
@@ -55,7 +55,11 @@ public:
 
     void setChemicalSystemConcrete(
         std::vector<double> const& concentrations,
-        GlobalIndexType const& chemical_system_id) override;
+        GlobalIndexType const& chemical_system_id,
+        MaterialPropertyLib::Medium const* medium,
+        MaterialPropertyLib::VariableArray const& vars,
+        ParameterLib::SpatialPosition const& pos, double const t,
+        double const dt) override;
 
     void executeInitialCalculation() override;
 
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 9f3b8d049a75c986c55f8b5d4eeaa6d30f89747e..2d345fa4de207e5661954bd9ee162bf330e4a97a 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -286,16 +286,25 @@ public:
     }
 
     void setChemicalSystemConcrete(Eigen::VectorXd const& local_x,
-                                   double const /*t*/, double /*dt*/) override
+                                   double const t, double dt) override
     {
         assert(_process_data.chemical_solver_interface);
 
+        auto const& medium =
+            _process_data.media_map->getMedium(_element.getID());
+
+        MaterialPropertyLib::VariableArray vars;
+
+        ParameterLib::SpatialPosition pos;
+        pos.setElementID(_element.getID());
+
         unsigned const n_integration_points =
             _integration_method.getNumberOfPoints();
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             auto& ip_data = _ip_data[ip];
             auto const& N = ip_data.N;
+            auto& porosity = ip_data.porosity;
             auto const& chemical_system_id = ip_data.chemical_system_id;
 
             auto const n_component = _transport_process_variables.size();
@@ -314,8 +323,11 @@ public:
                                                  C_int_pt[component_id]);
             }
 
+            vars[static_cast<int>(MaterialPropertyLib::Variable::porosity)] =
+                porosity;
+
             _process_data.chemical_solver_interface->setChemicalSystemConcrete(
-                C_int_pt, chemical_system_id);
+                C_int_pt, chemical_system_id, medium, vars, pos, t, dt);
         }
     }