diff --git a/dumux/discretization/staggered/freeflow/velocityoutput.hh b/dumux/discretization/staggered/freeflow/velocityoutput.hh
index 1d8683a4f6a569f1edc20168d76e64ca0a7ee4da..fadcbd7115836ea63febccab7ee5a1611df33cdf 100644
--- a/dumux/discretization/staggered/freeflow/velocityoutput.hh
+++ b/dumux/discretization/staggered/freeflow/velocityoutput.hh
@@ -82,6 +82,10 @@ public:
     static std::string phaseName(int phaseIdx)
     { return GET_PROP_TYPE(TypeTag, FluidSystem)::phaseName(phaseIdx); }
 
+    // returns the number of phase velocities computed by this class
+    static constexpr int numPhaseVelocities()
+    { return GET_PROP_TYPE(TypeTag, ModelTraits)::numPhases(); }
+
     //! Return the problem boundary types
     auto problemBoundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const
     { return problem_.boundaryTypes(element, scvf); }
diff --git a/dumux/io/vtkoutputmodule.hh b/dumux/io/vtkoutputmodule.hh
index 63522f9d8970d18a6406abf6d4a311e42dc062d4..da5cb25c25e4050995c44f9c03f81e6c6bf683a2 100644
--- a/dumux/io/vtkoutputmodule.hh
+++ b/dumux/io/vtkoutputmodule.hh
@@ -68,8 +68,7 @@ class VtkOutputModule
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using VelocityOutput = typename GET_PROP_TYPE(TypeTag, VelocityOutput);
-    using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
-    static constexpr int numPhases = ModelTraits::numPhases();
+    static constexpr int numPhaseVelocities = VelocityOutput::numPhaseVelocities();
 
     using VV = typename GridVariables::VolumeVariables;
     using Scalar = typename GridVariables::Scalar;
@@ -263,7 +262,7 @@ private:
 
         // instatiate the velocity output
         VelocityOutput velocityOutput(problem_, gridGeom_, gridVariables_, sol_);
-        std::array<std::vector<VelocityVector>, numPhases> velocity;
+        std::array<std::vector<VelocityVector>, numPhaseVelocities> velocity;
 
         // process rank
         static bool addProcessRank = getParamFromGroup<bool>(paramGroup_, "Vtk.AddProcessRank");
@@ -291,7 +290,7 @@ private:
 
             if (velocityOutput.enableOutput())
             {
-                for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
                 {
                     if(isBox && dim == 1)
                         velocity[phaseIdx].resize(numCells);
@@ -343,7 +342,7 @@ private:
 
                 // velocity output
                 if (velocityOutput.enableOutput())
-                    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                    for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
                         velocityOutput.calculateVelocity(velocity[phaseIdx], elemVolVars, fvGeometry, element, phaseIdx);
 
                 //! the rank
@@ -378,7 +377,7 @@ private:
             {
                 if (isBox && dim > 1)
                 {
-                    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                    for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
                         sequenceWriter_.addVertexData( Field(gridGeom_.gridView(), gridGeom_.vertexMapper(), velocity[phaseIdx],
                                                              "velocity_" + velocityOutput.phaseName(phaseIdx+phaseIdxOffset) + " (m/s)",
                                                              /*numComp*/dimWorld, /*codim*/dim).get() );
@@ -386,7 +385,7 @@ private:
                 // cell-centered models
                 else
                 {
-                    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                    for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
                         sequenceWriter_.addCellData( Field(gridGeom_.gridView(), gridGeom_.elementMapper(), velocity[phaseIdx],
                                                            "velocity_" + velocityOutput.phaseName(phaseIdx+phaseIdxOffset) + " (m/s)",
                                                            /*numComp*/dimWorld, /*codim*/0).get() );
@@ -432,7 +431,7 @@ private:
 
         // instatiate the velocity output
         VelocityOutput velocityOutput(problem_, gridGeom_, gridVariables_, sol_);
-        std::array<std::vector<VelocityVector>, numPhases> velocity;
+        std::array<std::vector<VelocityVector>, numPhaseVelocities> velocity;
 
         // process rank
         static bool addProcessRank = getParamFromGroup<bool>(paramGroup_, "Vtk.AddProcessRank");
@@ -462,7 +461,7 @@ private:
 
             if (velocityOutput.enableOutput())
             {
-                for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
                 {
                     if(isBox && dim == 1)
                         velocity[phaseIdx].resize(numCells);
@@ -520,7 +519,7 @@ private:
 
                 // velocity output
                 if (velocityOutput.enableOutput())
-                    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                    for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
                         velocityOutput.calculateVelocity(velocity[phaseIdx], elemVolVars, fvGeometry, element, phaseIdx);
 
                 //! the rank
@@ -546,14 +545,14 @@ private:
             {
                 // node-wise velocities
                 if (dim > 1)
-                    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                    for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
                         sequenceWriter_.addVertexData( Field(gridGeom_.gridView(), gridGeom_.vertexMapper(), velocity[phaseIdx],
                                                              "velocity_" + velocityOutput.phaseName(phaseIdx+phaseIdxOffset) + " (m/s)",
                                                              /*numComp*/dimWorld, /*codim*/dim).get() );
 
                 // cell-wise velocities
                 else
-                    for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                    for (int phaseIdx = 0; phaseIdx < numPhaseVelocities; ++phaseIdx)
                         sequenceWriter_.addCellData( Field(gridGeom_.gridView(), gridGeom_.elementMapper(), velocity[phaseIdx],
                                                            "velocity_" + velocityOutput.phaseName(phaseIdx+phaseIdxOffset) + " (m/s)",
                                                            /*numComp*/dimWorld, /*codim*/0).get() );
diff --git a/dumux/porousmediumflow/velocityoutput.hh b/dumux/porousmediumflow/velocityoutput.hh
index b7877c5e2a07e59ef446f7b3fb1bde92ee238ed4..a80fbcb5a5cabf2e9f063aeb512ff7d3a9f9b71e 100644
--- a/dumux/porousmediumflow/velocityoutput.hh
+++ b/dumux/porousmediumflow/velocityoutput.hh
@@ -111,6 +111,10 @@ public:
     static std::string phaseName(int phaseIdx)
     { return FluidSystem::phaseName(phaseIdx); }
 
+    // returns the number of phase velocities computed by this class
+    static constexpr int numPhaseVelocities()
+    { return GET_PROP_TYPE(TypeTag, ModelTraits)::numPhases(); }
+
     // The following SFINAE enable_if usage allows compilation, even if only a
     //
     // boundaryTypes(const Element&, const scv&)