diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index af86f98646835dd42f826c79497c11e849bce660..2425bfefe4182adc1dfd7513ceed32606da5f894 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -916,6 +916,21 @@ public:
         return flux;
     }
 
+    std::vector<double> interpolateNodalValuesToIntegrationPoints(
+        std::vector<double> const& local_x) override
+    {
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
+
+        std::vector<double> interpolated_values(n_integration_points);
+        for (unsigned ip(0); ip < n_integration_points; ++ip)
+        {
+            NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N,
+                                             interpolated_values[ip]);
+        }
+        return interpolated_values;
+    }
+
     std::vector<double> const& getInterpolatedLocalSolution(
         const double /*t*/,
         std::vector<GlobalVector*> const& int_pt_x,
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
index e7c0af359707e9bdf39696ccd31c9796c8d78d3f..d824b58ebc6362b1f37040d0048bd7bbd7b580dc 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.cpp
@@ -174,6 +174,77 @@ void ComponentTransportProcess::
         _local_assemblers, pv.getActiveElementIDs(), _coupled_solutions);
 }
 
+std::vector<GlobalVector>
+ComponentTransportProcess::interpolateNodalValuesToIntegrationPoints(
+    std::vector<GlobalVector*> const& nodal_values_vectors) const
+{
+    // Result is for each process a vector of integration point values for each
+    // element stored consecutively.
+    auto interpolateNodalValuesToIntegrationPoints =
+        [](std::size_t mesh_item_id,
+           LocalAssemblerInterface& local_assembler,
+           std::vector<
+               std::reference_wrapper<NumLib::LocalToGlobalIndexMap>> const&
+               dof_tables,
+           std::vector<GlobalVector*> const& x,
+           std::vector<std::vector<double>>& int_pt_x) {
+            for (unsigned process_id = 0; process_id < x.size(); ++process_id)
+            {
+                auto const& dof_table = dof_tables[process_id].get();
+                auto const indices =
+                    NumLib::getIndices(mesh_item_id, dof_table);
+                auto const local_x = x[process_id]->get(indices);
+
+                std::vector<double> const interpolated_values =
+                    local_assembler.interpolateNodalValuesToIntegrationPoints(
+                        local_x);
+
+                // For each element (mesh_item_id) concatenate the integration
+                // point values.
+                int_pt_x[process_id].insert(int_pt_x[process_id].end(),
+                                            interpolated_values.begin(),
+                                            interpolated_values.end());
+            }
+        };
+
+    // Same dof table for each primary variable.
+    std::vector<std::reference_wrapper<NumLib::LocalToGlobalIndexMap>>
+        dof_tables;
+    dof_tables.reserve(nodal_values_vectors.size());
+    std::generate_n(std::back_inserter(dof_tables), nodal_values_vectors.size(),
+                    [&]() { return std::ref(*_local_to_global_index_map); });
+
+    std::vector<std::vector<double>> integration_point_values_vecs(
+        nodal_values_vectors.size());
+
+    GlobalExecutor::executeDereferenced(
+        interpolateNodalValuesToIntegrationPoints, _local_assemblers,
+        dof_tables, nodal_values_vectors, integration_point_values_vecs);
+
+    // Convert from std::vector<std::vector<double>> to
+    // std::vector<GlobalVector>
+    std::vector<GlobalVector> integration_point_values_vectors;
+    integration_point_values_vectors.reserve(nodal_values_vectors.size());
+    for (auto const& integration_point_values_vec :
+         integration_point_values_vecs)
+    {
+        GlobalIndexType const size = integration_point_values_vec.size();
+        // New vector of size (elements x integration_points_per_element)
+        GlobalVector integration_point_values_vector(size);
+
+        // Copy one by one.
+        for (GlobalIndexType i = 0; i < size; ++i)
+        {
+            integration_point_values_vector.set(
+                i, integration_point_values_vec[i]);
+        }
+        integration_point_values_vectors.push_back(
+            std::move(integration_point_values_vector));
+    }
+
+    return integration_point_values_vectors;
+}
+
 void ComponentTransportProcess::extrapolateIntegrationPointValuesToNodes(
     const double t,
     std::vector<GlobalVector*> const& integration_point_values_vectors,
diff --git a/ProcessLib/ComponentTransport/ComponentTransportProcess.h b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
index 374f771f5e707e6c4081df10e56fc146cdb3979b..d86b67547faa831a24b14322ee4da282e4636c5c 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportProcess.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportProcess.h
@@ -117,6 +117,9 @@ public:
     void setCoupledTermForTheStaggeredSchemeToLocalAssemblers(
         int const process_id) override;
 
+    std::vector<GlobalVector> interpolateNodalValuesToIntegrationPoints(
+        std::vector<GlobalVector*> const& nodal_values_vectors) const override;
+
     void extrapolateIntegrationPointValuesToNodes(
         const double t,
         std::vector<GlobalVector*> const& integration_point_values_vectors,
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index 2b09c451092607e7990315e1b440c42d584a7d71..253fa45fe5295524951debbe9d4488093fd59d8a 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -98,6 +98,12 @@ public:
                              GlobalVector const& x, double const t,
                              double const dt, bool const use_monolithic_scheme);
 
+    virtual std::vector<double> interpolateNodalValuesToIntegrationPoints(
+        std::vector<double> const& /*local_x*/)
+    {
+        return {};
+    }
+
     /// Computes the flux in the point \c p_local_coords that is given in local
     /// coordinates using the values from \c local_x.
     /// Fits to monolithic scheme.
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 5838c56e19ccc3c7575ca91667f8f33ce84b3045..43f04a22d21f63c2eca5c80a530c99e578ac1bd3 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -101,6 +101,12 @@ public:
     {
     }
 
+    virtual std::vector<GlobalVector> interpolateNodalValuesToIntegrationPoints(
+        std::vector<GlobalVector*> const& /*nodal_values_vectors*/) const
+    {
+        return {};
+    }
+
     virtual void extrapolateIntegrationPointValuesToNodes(
         const double /*t*/,
         std::vector<GlobalVector*> const& /*integration_point_values_vectors*/,
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index 8afb0961cb452e9f3f8bea43f2033680d3941375..35d7e748baba0f3d120e2ac50cd7be5aa9950d1b 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -453,6 +453,11 @@ void TimeLoop::initialize()
     {
         BaseLib::RunTime time_phreeqc;
         time_phreeqc.start();
+
+        auto& pcs = _per_process_data[0]->process;
+        auto const interpolated_process_solutions =
+            pcs.interpolateNodalValuesToIntegrationPoints(_process_solutions);
+
         _chemical_solver_interface->executeInitialCalculation(
             _process_solutions);
         INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());
@@ -815,6 +820,11 @@ TimeLoop::solveCoupledEquationSystemsByStaggeredScheme(
         // space and localized chemical equilibrium between solutes.
         BaseLib::RunTime time_phreeqc;
         time_phreeqc.start();
+
+        auto& pcs = _per_process_data[0]->process;
+        auto const interpolated_process_solutions =
+            pcs.interpolateNodalValuesToIntegrationPoints(_process_solutions);
+
         _chemical_solver_interface->doWaterChemistryCalculation(
             _process_solutions, dt);
         INFO("[time] Phreeqc took {:g} s.", time_phreeqc.elapsed());