diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
index 6f309631377062dae369d2b72ca48164944cbf04..425c20b8b9cdabc1812964aed3aea406902b0410 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
@@ -26,18 +26,22 @@ namespace ProcessLib
 namespace GroundwaterFlow
 {
 
-enum class IntegrationPointValue {
-    DarcyVelocityX,
-    DarcyVelocityY,
-    DarcyVelocityZ
-};
-
 const unsigned NUM_NODAL_DOF = 1;
 
 class GroundwaterFlowLocalAssemblerInterface
         : public ProcessLib::LocalAssemblerInterface
-        , public NumLib::Extrapolatable<IntegrationPointValue>
-{};
+        , public NumLib::ExtrapolatableElement
+{
+public:
+    virtual std::vector<double> const& getIntPtDarcyVelocityX(
+        std::vector<double>& /*cache*/) const = 0;
+
+    virtual std::vector<double> const& getIntPtDarcyVelocityY(
+        std::vector<double>& /*cache*/) const = 0;
+
+    virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+        std::vector<double>& /*cache*/) const = 0;
+};
 
 template <typename ShapeFunction,
          typename IntegrationMethod,
@@ -118,22 +122,21 @@ public:
     }
 
     std::vector<double> const&
-    getIntegrationPointValues(IntegrationPointValue const property,
-                              std::vector<double>& /*cache*/) const override
+    getIntPtDarcyVelocityX(std::vector<double>& /*cache*/) const override
     {
-        switch (property)
-        {
-        case IntegrationPointValue::DarcyVelocityX:
             return _darcy_velocities[0];
-        case IntegrationPointValue::DarcyVelocityY:
-            assert(GlobalDim > 1);
-            return _darcy_velocities[1];
-        case IntegrationPointValue::DarcyVelocityZ:
-            assert(GlobalDim > 2);
-            return _darcy_velocities[2];
-        }
+    }
 
-        OGS_FATAL("");
+    std::vector<double> const&
+    getIntPtDarcyVelocityY(std::vector<double>& /*cache*/) const override
+    {
+            return _darcy_velocities[0];
+    }
+
+    std::vector<double> const&
+    getIntPtDarcyVelocityZ(std::vector<double>& /*cache*/) const override
+    {
+            return _darcy_velocities[0];
     }
 
 private:
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
index e25c277069c31cd3c1e787efa2944df7be25888f..5caf1007c06cf6a4e7af6e4e9c0ea1930dde5377 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.cpp
@@ -57,24 +57,25 @@ void GroundwaterFlowProcess::initializeConcreteProcess(
     _extrapolator.reset(new ExtrapolatorImplementation(
         Base::getMatrixSpecifications(), *Base::_local_to_global_index_map));
 
-    Base::_secondary_variables.addSecondaryVariable(
+    _secondary_variables.addSecondaryVariable(
         "darcy_velocity_x", 1,
-        makeExtrapolator(IntegrationPointValue::DarcyVelocityX, *_extrapolator,
-                         _local_assemblers));
+        makeExtrapolator(
+            *_extrapolator, _local_assemblers,
+            &GroundwaterFlowLocalAssemblerInterface::getIntPtDarcyVelocityX));
 
-    if (mesh.getDimension() > 1)
-    {
-        Base::_secondary_variables.addSecondaryVariable(
+    if (mesh.getDimension() > 1) {
+        _secondary_variables.addSecondaryVariable(
             "darcy_velocity_y", 1,
-            makeExtrapolator(IntegrationPointValue::DarcyVelocityY,
-                             *_extrapolator, _local_assemblers));
+            makeExtrapolator(*_extrapolator, _local_assemblers,
+                             &GroundwaterFlowLocalAssemblerInterface::
+                                 getIntPtDarcyVelocityY));
     }
-    if (mesh.getDimension() > 2)
-    {
-        Base::_secondary_variables.addSecondaryVariable(
+    if (mesh.getDimension() > 2) {
+        _secondary_variables.addSecondaryVariable(
             "darcy_velocity_z", 1,
-            makeExtrapolator(IntegrationPointValue::DarcyVelocityZ,
-                             *_extrapolator, _local_assemblers));
+            makeExtrapolator(*_extrapolator, _local_assemblers,
+                             &GroundwaterFlowLocalAssemblerInterface::
+                                 getIntPtDarcyVelocityZ));
     }
 }
 
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index adaaf66d597dcaac9aee5336b5b892d4e445ee45..8ea98af4e6c86744520189811dcc6492f694f9ba 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -43,11 +43,9 @@ public:
 
 private:
     using ExtrapolatorInterface =
-        NumLib::Extrapolator<IntegrationPointValue,
-                             GroundwaterFlowLocalAssemblerInterface>;
+        NumLib::Extrapolator;
     using ExtrapolatorImplementation =
-        NumLib::LocalLinearLeastSquaresExtrapolator<
-            IntegrationPointValue, GroundwaterFlowLocalAssemblerInterface>;
+        NumLib::LocalLinearLeastSquaresExtrapolator;
 
     void initializeConcreteProcess(
         NumLib::LocalToGlobalIndexMap const& dof_table,