diff --git a/dumux/discretization/staggered/properties.hh b/dumux/discretization/staggered/properties.hh
index 6a2816e2c968a6f950cac8fd80755f41f9ea1235..8dc80acd8f4f557946722bb574a2c6ecbc56f6c5 100644
--- a/dumux/discretization/staggered/properties.hh
+++ b/dumux/discretization/staggered/properties.hh
@@ -225,6 +225,7 @@ public:
     using type = BoundaryTypes<numEqCellCenter + numEqFace>;
 };
 
+//! A container class used to specify values for boundary/initial conditions
 SET_PROP(StaggeredModel, PrimaryVariables)
 {
 private:
@@ -234,6 +235,16 @@ public:
     using type = StaggeredPrimaryVariables<TypeTag, CellCenterPrimaryVariables, FacePrimaryVariables>;
 };
 
+//! A container class used to specify values for sources and Neumann BCs
+SET_PROP(StaggeredModel, NumEqVector)
+{
+private:
+    using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+    using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
+public:
+    using type = StaggeredPrimaryVariables<TypeTag, CellCenterPrimaryVariables, FacePrimaryVariables>;
+};
+
 SET_TYPE_PROP(StaggeredModel, GridVariables, StaggeredGridVariables<TypeTag>);
 
 //! Set one or different base epsilons for the calculations of the localJacobian's derivatives
diff --git a/dumux/freeflow/properties.hh b/dumux/freeflow/properties.hh
index c3992aef9af4b6b5829a7dc911799f7b374f6de5..71a0083b216dc51a587c81c15491222afc714f77 100644
--- a/dumux/freeflow/properties.hh
+++ b/dumux/freeflow/properties.hh
@@ -75,7 +75,7 @@ public:
 //! The variables living on the faces
 SET_TYPE_PROP(FreeFlow, FaceVariables, StaggeredFaceVariables<TypeTag>);
 
-//! A container class used to specify values for boundary/initial conditions and sources
+//! A container class used to specify values for boundary/initial conditions
 SET_PROP(FreeFlow, PrimaryVariables)
 {
 private:
@@ -87,6 +87,18 @@ public:
     using type = StaggeredPrimaryVariables<TypeTag, CellCenterBoundaryValues, FaceBoundaryValues>;
 };
 
+//! A container class used to specify values for sources and Neumann BCs
+SET_PROP(FreeFlow, NumEqVector)
+{
+private:
+    using CellCenterBoundaryValues = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using FaceBoundaryValues = Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Scalar),
+                                                 GridView::dimension>;
+public:
+    using type = StaggeredPrimaryVariables<TypeTag, CellCenterBoundaryValues, FaceBoundaryValues>;
+};
+
 //! Boundary types at a single degree of freedom
 SET_TYPE_PROP(FreeFlow,
               BoundaryTypes,
diff --git a/test/freeflow/staggered/channeltestproblem.hh b/test/freeflow/staggered/channeltestproblem.hh
index 548368cb5b0dfef60ae3e647fa1458097027bc21..7458858e06bfc5279b1fb46eb7f32959ce3a7df3 100644
--- a/test/freeflow/staggered/channeltestproblem.hh
+++ b/test/freeflow/staggered/channeltestproblem.hh
@@ -129,6 +129,7 @@ class ChannelTestProblem : public NavierStokesProblem<TypeTag>
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
 
@@ -163,9 +164,9 @@ public:
      *
      * \param globalPos The global position
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return PrimaryVariables(0.0);
+        return SourceValues(0.0);
     }
     // \}
     /*!
diff --git a/test/freeflow/staggered/closedsystemtestproblem.hh b/test/freeflow/staggered/closedsystemtestproblem.hh
index 7a39deeb389e1b4184aab344b5573f3d68f431a8..755eb3ac7e1236eac04aa981107b565338ad39c1 100644
--- a/test/freeflow/staggered/closedsystemtestproblem.hh
+++ b/test/freeflow/staggered/closedsystemtestproblem.hh
@@ -111,6 +111,7 @@ class ClosedSystemTestProblem : public NavierStokesProblem<TypeTag>
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
 
     using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
@@ -153,9 +154,9 @@ public:
      * \param values Stores the source values, acts as return value
      * \param globalPos The global position
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return PrimaryVariables(0.0);
+        return SourceValues(0.0);
     }
     // \}
     /*!
diff --git a/test/freeflow/staggered/doneatestproblem.hh b/test/freeflow/staggered/doneatestproblem.hh
index 74c192e2248f814e28a76b96880444ae92daffd6..cff61768e66e2a85a3ef9008098318ea2db76597 100644
--- a/test/freeflow/staggered/doneatestproblem.hh
+++ b/test/freeflow/staggered/doneatestproblem.hh
@@ -120,6 +120,7 @@ class DoneaTestProblem : public NavierStokesProblem<TypeTag>
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
 
@@ -188,9 +189,9 @@ public:
      *
      * \param globalPos The global position
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
     {
-        PrimaryVariables source(0.0);
+        SourceValues source(0.0);
         Scalar x = globalPos[0];
         Scalar y = globalPos[1];
 
diff --git a/test/freeflow/staggered/kovasznaytestproblem.hh b/test/freeflow/staggered/kovasznaytestproblem.hh
index 62ab2d84f8d3829c9b3818ba91c8e9528011a6fe..0bba7f84d98f610be57c0222793b819a27e92557 100644
--- a/test/freeflow/staggered/kovasznaytestproblem.hh
+++ b/test/freeflow/staggered/kovasznaytestproblem.hh
@@ -121,6 +121,7 @@ class KovasznayTestProblem : public NavierStokesProblem<TypeTag>
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
 
@@ -188,9 +189,9 @@ public:
      *
      * \param globalPos The global position
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return PrimaryVariables(0.0);
+        return SourceValues(0.0);
     }
 
     // \}
diff --git a/test/freeflow/staggerednc/channeltestproblem.hh b/test/freeflow/staggerednc/channeltestproblem.hh
index 978ee517fc31ba045b31a003e680c977da046d44..c0d69eeabd3d341a205be49d411e2d93ad8e2bbf 100644
--- a/test/freeflow/staggerednc/channeltestproblem.hh
+++ b/test/freeflow/staggerednc/channeltestproblem.hh
@@ -138,6 +138,7 @@ class ChannelNCTestProblem : public NavierStokesProblem<TypeTag>
     using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
     using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
@@ -176,9 +177,9 @@ public:
      *
      * \param globalPos The global position
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return PrimaryVariables(0.0);
+        return SourceValues(0.0);
     }
     // \}
     /*!
diff --git a/test/freeflow/staggerednc/densityflowproblem.hh b/test/freeflow/staggerednc/densityflowproblem.hh
index 498f15283747f2ceef30c5d83b20bdbfc62a7ed5..31d9da1ab34bf723d81702fb4ebfcd4c8cb125a4 100644
--- a/test/freeflow/staggerednc/densityflowproblem.hh
+++ b/test/freeflow/staggerednc/densityflowproblem.hh
@@ -123,6 +123,7 @@ class DensityDrivenFlowProblem : public NavierStokesProblem<TypeTag>
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using SourceValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
 
     using TimeLoopPtr = std::shared_ptr<CheckPointTimeLoop<Scalar>>;
     using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
@@ -164,9 +165,9 @@ public:
      *
      * \param globalPos The global position
      */
-    PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+    SourceValues sourceAtPos(const GlobalPosition &globalPos) const
     {
-        return PrimaryVariables(0.0);
+        return SourceValues(0.0);
     }
     // \}
     /*!