diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index 22673a693bc40003054f8640d51ea25e3cb03099..c32c8e1a4dcc42da83d8451e4948a02b00e6c6ab 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -235,7 +235,7 @@ struct StaggeredFaceSolution { using type = UndefinedProperty; };
 template<class TypeTag, class MyTypeTag>
 struct EnableGridFaceVariablesCache { using type = UndefinedProperty; };        //!< Switch on/off caching of face variables
 template<class TypeTag, class MyTypeTag>
-struct UpwindSchemeOrder { using type = UndefinedProperty; };                   //!< Specifies the order of the upwinding scheme (1 == basic upwinding, 2 == higher order)
+struct UpwindSchemeOrder { using type = UndefinedProperty; };                   //!< Specifies the order of the upwinding scheme (1 == first order, 2 == second order(tvd methods))
 
 /////////////////////////////////////////////////////////////
 // Properties used by the mpnc model
diff --git a/dumux/discretization/staggered/freeflow/connectivitymap.hh b/dumux/discretization/staggered/freeflow/connectivitymap.hh
index 3fe85d132fd8ba6577c31d61d2acad16c99285fd..1aa82d24fe1e852e2225648e240e1eb4ac79be2b 100644
--- a/dumux/discretization/staggered/freeflow/connectivitymap.hh
+++ b/dumux/discretization/staggered/freeflow/connectivitymap.hh
@@ -52,9 +52,11 @@ class StaggeredFreeFlowConnectivityMap
     using FaceToCellCenterMap = std::vector<std::vector<GridIndexType>>;
     using FaceToFaceMap = std::vector<std::vector<GridIndexType>>;
 
+    using SmallLocalIndex = typename IndexTraits<GridView>::SmallLocalIndex;
+
     using Stencil = std::vector<GridIndexType>;
 
-    static constexpr int upwindSchemeOrder = FVGridGeometry::upwindSchemeOrder;
+    static constexpr SmallLocalIndex upwindSchemeOrder = FVGridGeometry::upwindSchemeOrder;
     static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
 
 public:
@@ -195,9 +197,9 @@ private:
                 stencil.push_back(data.normalPair.second);
 
             // add parallel dofs
-            for (int i = 0; i < upwindSchemeOrder; i++)
+            for (SmallLocalIndex i = 0; i < upwindSchemeOrder; i++)
             {
-                if(!(data.parallelDofs[i] < 0))
+                if (data.hasParallelNeighbor[i])
                     stencil.push_back(data.parallelDofs[i]);
             }
         }
@@ -207,7 +209,7 @@ private:
 
     void addHigherOrderInAxisDofs_(const SubControlVolumeFace& scvf, Stencil& stencil, std::true_type)
     {
-        for (int i = 0; i < upwindSchemeOrder - 1; i++)
+        for (SmallLocalIndex i = 0; i < upwindSchemeOrder - 1; i++)
         {
             if (scvf.hasBackwardNeighbor(i))
                 stencil.push_back(scvf.axisData().inAxisBackwardDofs[i]);
diff --git a/dumux/discretization/staggered/freeflow/facevariables.hh b/dumux/discretization/staggered/freeflow/facevariables.hh
index 159018d08c95c92185f58858e8d6dce942458346..4645f66c52013ef06d216e19ec2e33b9b1fdd660 100644
--- a/dumux/discretization/staggered/freeflow/facevariables.hh
+++ b/dumux/discretization/staggered/freeflow/facevariables.hh
@@ -53,8 +53,11 @@ struct InAxisVelocities<Scalar, 1>
  * \ingroup StaggeredDiscretization
  * \brief The face variables class for free flow staggered grid models.
  *        Contains all relevant velocities for the assembly of the momentum balance.
+ *        When the upwindSchemeOrder is set to 2, additional velocities located at Dofs
+ *        further from the central stencil will be added and used when calculating the
+ *        advective term. When the order remains at 1, these velocities will not be provided.
  */
-template<class FacePrimaryVariables, int dim, int upwindSchemeOrder> // TODO doc
+template<class FacePrimaryVariables, int dim, int upwindSchemeOrder>
 class StaggeredFaceVariables
 {
     static constexpr int numPairs = (dim == 2) ? 2 : 4;
@@ -203,7 +206,7 @@ private:
         // forward wrt the self face by degree i
         for (int i = 0; i < scvf.axisData().inAxisForwardDofs.size(); i++)
         {
-             if (!(scvf.axisData().inAxisForwardDofs[i] < 0))
+             if (scvf.hasForwardNeighbor(i))
                  inAxisVelocities_.forward[i]= faceSol[scvf.axisData().inAxisForwardDofs[i]];
         }
 
@@ -211,7 +214,7 @@ private:
         // behind the opposite face by degree i
         for (int i = 0; i < scvf.axisData().inAxisBackwardDofs.size(); i++)
         {
-             if (!(scvf.axisData().inAxisBackwardDofs[i] < 0))
+             if (scvf.hasBackwardNeighbor(i))
                  inAxisVelocities_.backward[i] = faceSol[scvf.axisData().inAxisBackwardDofs[i]];
         }
     }
diff --git a/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh b/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
index 0032f2cb009c30087e78874b74ca181d16bacc85..c9e88535066c1de09024b420d6427001c08d74ea 100644
--- a/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
+++ b/dumux/discretization/staggered/freeflow/fvgridgeometrytraits.hh
@@ -39,15 +39,15 @@ namespace Dumux {
  * \ingroup StaggeredDiscretization
  * \brief Default traits for the finite volume grid geometry.
  */
-template<class GridView, int USO, class MapperTraits = DefaultMapperTraits<GridView>>
+template<class GridView, int upwOrder, class MapperTraits = DefaultMapperTraits<GridView>>
 struct StaggeredFreeFlowDefaultFVGridGeometryTraits
 : public MapperTraits
 {
     using SubControlVolume = CCSubControlVolume<GridView>;
-    using SubControlVolumeFace = FreeFlowStaggeredSubControlVolumeFace<GridView, USO>;
+    using SubControlVolumeFace = FreeFlowStaggeredSubControlVolumeFace<GridView, upwOrder>;
     using IntersectionMapper = ConformingGridIntersectionMapper<GridView>;
-    using GeometryHelper = FreeFlowStaggeredGeometryHelper<GridView, USO>;
-    static constexpr int upwindSchemeOrder = USO;
+    using GeometryHelper = FreeFlowStaggeredGeometryHelper<GridView, upwOrder>;
+    static constexpr int upwindSchemeOrder = upwOrder;
 
     struct DofTypeIndices
     {
diff --git a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
index 09c568e7a904c736034d3f6442fa90d9d5f3a3a0..071888dfddb0381fe66e3ae4b540fe40dadb856d 100644
--- a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
+++ b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh
@@ -25,27 +25,37 @@
 #define DUMUX_DISCRETIZATION_STAGGERED_GEOMETRY_HELPER_HH
 
 #include <dune/geometry/multilineargeometry.hh>
-#include <dune/geometry/referenceelements.hh>
 
 #include <dumux/common/math.hh>
+#include <dumux/common/indextraits.hh>
 #include <type_traits>
 #include <algorithm>
 #include <array>
+#include <bitset>
 
 namespace Dumux
 {
 
+namespace Detail {
+
 /*!
  * \ingroup StaggeredDiscretization
  * \brief Parallel Data stored per sub face
  */
-template<class Scalar, class GlobalPosition, int upwindSchemeOrder>
+template<class GridView, int upwindSchemeOrder>
 struct PairData
 {
-    std::array<int, upwindSchemeOrder> parallelDofs;
+    using Scalar = typename GridView::ctype;
+    using GlobalPosition = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
+    using GridIndexType = typename IndexTraits<GridView>::GridIndex;
+    using SmallLocalIndexType = typename IndexTraits<GridView>::SmallLocalIndex;
+
+    std::bitset<upwindSchemeOrder> hasParallelNeighbor;
+    std::array<GridIndexType, upwindSchemeOrder> parallelDofs;
     std::array<Scalar, upwindSchemeOrder+1> parallelDistances; // TODO store only two distances, not three.
-    std::pair<signed int, signed int> normalPair;
-    int localNormalFaceIdx;
+    bool hasNormalNeighbor = false;
+    std::pair<GridIndexType, GridIndexType> normalPair;
+    SmallLocalIndexType localNormalFaceIdx;
     Scalar normalDistance;
     GlobalPosition virtualFirstParallelFaceDofPos;
 };
@@ -55,27 +65,24 @@ struct PairData
  * \ingroup StaggeredDiscretization
  * \brief In Axis Data stored per sub face
  */
-template<class Scalar, int upwindSchemeOrder>
+template<class GridView, int upwindSchemeOrder>
 struct AxisData;
 /*!
  * \ingroup StaggeredDiscretization
  * \brief In Axis Data stored per sub face
  */
-template<class Scalar, int upwindSchemeOrder>
+template<class GridView, int upwindSchemeOrder>
 struct AxisData
 {
-    AxisData()
-    {
-        inAxisForwardDofs.fill(-1);
-        inAxisBackwardDofs.fill(-1);
-        inAxisForwardDistances.fill(0);
-        inAxisBackwardDistances.fill(0);
-    }
-
-    int selfDof;
-    int oppositeDof;
-    std::array<int, upwindSchemeOrder-1> inAxisForwardDofs;
-    std::array<int, upwindSchemeOrder-1> inAxisBackwardDofs;
+    using Scalar = typename GridView::ctype;
+    using GridIndexType = typename IndexTraits<GridView>::GridIndex;
+
+    GridIndexType selfDof;
+    GridIndexType oppositeDof;
+    std::bitset<upwindSchemeOrder-1> hasForwardNeighbor;
+    std::bitset<upwindSchemeOrder-1> hasBackwardNeighbor;
+    std::array<GridIndexType, upwindSchemeOrder-1> inAxisForwardDofs;
+    std::array<GridIndexType, upwindSchemeOrder-1> inAxisBackwardDofs;
     Scalar selfToOppositeDistance;
     std::array<Scalar, upwindSchemeOrder-1> inAxisForwardDistances;
     std::array<Scalar, upwindSchemeOrder-1> inAxisBackwardDistances;
@@ -85,16 +92,18 @@ struct AxisData
  * \ingroup StaggeredDiscretization
  * \brief In Axis Data stored per sub face
  */
-template<class Scalar>
-struct AxisData<Scalar, 1>
+template<class GridView>
+struct AxisData<GridView, 1>
 {
-    int selfDof;
-    int oppositeDof;
+    using Scalar = typename GridView::ctype;
+    using GridIndexType = typename IndexTraits<GridView>::GridIndex;
+
+    GridIndexType selfDof;
+    GridIndexType oppositeDof;
     Scalar selfToOppositeDistance;
 };
 
-// template<class Scalar, int upwindSchemeOrder>
-// using AxisData
+} // namespace Detail
 
 /*!
  * \ingroup StaggeredDiscretization
@@ -104,8 +113,7 @@ template<class Vector>
 inline static unsigned int directionIndex(Vector&& vector)
 {
     const auto eps = 1e-8;
-    const int idx = std::find_if(vector.begin(), vector.end(), [eps](const auto& x) { return std::abs(x) > eps; } ) - vector.begin();
-    return idx;
+    return std::find_if(vector.begin(), vector.end(), [eps](const auto& x) { return std::abs(x) > eps; } ) - vector.begin();
 }
 
 /*!
@@ -117,30 +125,28 @@ template<class GridView, int upwindSchemeOrder>
 class FreeFlowStaggeredGeometryHelper
 {
     using Scalar = typename GridView::ctype;
-    static constexpr int dim = GridView::dimension;
-    static constexpr int dimWorld = GridView::dimensionworld;
-
-    using ScvGeometry = Dune::CachedMultiLinearGeometry<Scalar, dim, dimWorld>;
-    using ScvfGeometry = Dune::CachedMultiLinearGeometry<Scalar, dim-1, dimWorld>;
-
-    using GlobalPosition = typename ScvGeometry::GlobalCoordinate;
+    static constexpr auto dim = GridView::dimension;
+    static constexpr auto dimWorld = GridView::dimensionworld;
 
     using Element = typename GridView::template Codim<0>::Entity;
     using Intersection = typename GridView::Intersection;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
-    using ReferenceElements = typename Dune::ReferenceElements<Scalar, dim>;
+    using GridIndexType = typename IndexTraits<GridView>::GridIndex;
+    using SmallLocalIndexType = typename IndexTraits<GridView>::SmallLocalIndex;
 
     //TODO include assert that checks for quad geometry
-    static constexpr int codimIntersection =  1;
-    static constexpr int codimCommonEntity = 2;
-    static constexpr int numFacetSubEntities = (dim == 2) ? 2 : 4;
-    static constexpr int numfacets = dimWorld * 2;
-    static constexpr int numPairs = 2 * (dimWorld - 1);
+    static constexpr auto codimIntersection =  1;
+    static constexpr auto codimCommonEntity = 2;
+    static constexpr auto numFacetSubEntities = (dim == 2) ? 2 : 4;
+    static constexpr auto numfacets = dimWorld * 2;
+    static constexpr auto numPairs = 2 * (dimWorld - 1);
 
     static constexpr bool useHigherOrder = upwindSchemeOrder > 1;
 
-
 public:
+    using PairData = Detail::PairData<GridView, upwindSchemeOrder>;
+    using AxisData = Detail::AxisData<GridView, upwindSchemeOrder>;
 
     FreeFlowStaggeredGeometryHelper(const Element& element, const GridView& gridView) : element_(element), elementGeometry_(element.geometry()), gridView_(gridView)
     { }
@@ -157,7 +163,7 @@ public:
     /*!
      * \brief Returns the local index of the face (i.e. the intersection)
      */
-    int localFaceIndex() const
+    SmallLocalIndexType localFaceIndex() const
     {
         return intersection_.indexInInside();
     }
@@ -165,7 +171,7 @@ public:
     /*!
      * \brief Returns a copy of the axis data
      */
-    auto axisData() const
+    AxisData axisData() const
     {
         return axisData_;
     }
@@ -173,7 +179,7 @@ public:
     /*!
      * \brief Returns a copy of the pair data
      */
-    auto pairData() const
+    std::array<PairData, numPairs> pairData() const
     {
         return pairData_;
     }
@@ -181,15 +187,15 @@ public:
     /*!
      * \brief Returns the dirction index of the primary facet (0 = x, 1 = y, 2 = z)
      */
-    int directionIndex() const
+    unsigned int directionIndex() const
     {
-        return Dumux::directionIndex(std::move(intersection_.centerUnitOuterNormal()));
+        return Dumux::directionIndex(intersection_.centerUnitOuterNormal());
     }
 
     /*!
      * \brief Returns the dirction index of the facet passed as an argument (0 = x, 1 = y, 2 = z)
      */
-    int directionIndex(const Intersection& intersection) const
+    unsigned int directionIndex(const Intersection& intersection) const
     {
         return Dumux::directionIndex(std::move(intersection.centerUnitOuterNormal()));
     }
@@ -202,6 +208,8 @@ private:
      */
     void fillAxisData_()
     {
+        fillOuterAxisData_(std::integral_constant<bool, useHigherOrder>{});
+
         const auto inIdx = intersection_.indexInInside();
         const auto oppIdx = localOppositeIdx_(inIdx);
 
@@ -216,7 +224,6 @@ private:
         const auto opposite = getFacet_(oppIdx, element_);
         axisData_.selfToOppositeDistance = (self.geometry().center() - opposite.geometry().center()).two_norm();
 
-        fillOuterAxisData_(std::integral_constant<bool, useHigherOrder>{});
     }
 
     /*!
@@ -231,6 +238,9 @@ private:
      */
     void fillOuterAxisData_(std::true_type)
     {
+        // reset the axis data struct
+        axisData_ = {};
+
         const auto numForwardBackwardAxisDofs = axisData_.inAxisForwardDofs.size();
 
         // Set the Forward Dofs
@@ -266,6 +276,7 @@ private:
         while(!inAxisForwardElementStack.empty())
         {
             int forwardIdx = inAxisForwardElementStack.size()-1;
+            axisData_.hasForwardNeighbor.set(forwardIdx, true);
             axisData_.inAxisForwardDofs[forwardIdx] = gridView_.indexSet().subIndex(inAxisForwardElementStack.top(), inIdx, codimIntersection);
             selfFace = getFacet_(inIdx, inAxisForwardElementStack.top());
             forwardFaceCoordinates[forwardIdx] = selfFace.geometry().center();
@@ -273,16 +284,12 @@ private:
         }
 
         const auto self = getFacet_(inIdx, element_);
-        for(int i = 0; i< forwardFaceCoordinates.size(); i++)
+        for (int i = 0; i< forwardFaceCoordinates.size(); i++)
         {
-            if(i == 0)
-            {
+            if (i == 0)
                 axisData_.inAxisForwardDistances[i] = (forwardFaceCoordinates[i] - self.geometry().center()).two_norm();
-            }
             else
-            {
                 axisData_.inAxisForwardDistances[i] = (forwardFaceCoordinates[i]- forwardFaceCoordinates[i-1]).two_norm();
-            }
         }
 
         // Set the Backward Dofs
@@ -325,23 +332,20 @@ private:
         while(!inAxisBackwardElementStack.empty())
         {
             int backwardIdx = inAxisBackwardElementStack.size()-1;
-            this->axisData_.inAxisBackwardDofs[backwardIdx] = gridView_.indexSet().subIndex(inAxisBackwardElementStack.top(), oppIdx, codimIntersection);
+            axisData_.hasBackwardNeighbor.set(backwardIdx, true);
+            axisData_.inAxisBackwardDofs[backwardIdx] = gridView_.indexSet().subIndex(inAxisBackwardElementStack.top(), oppIdx, codimIntersection);
             oppFace = getFacet_(oppIdx, inAxisBackwardElementStack.top());
             backwardFaceCoordinates[backwardIdx] = oppFace.geometry().center();
             inAxisBackwardElementStack.pop();
         }
 
         const auto opposite = getFacet_(oppIdx, element_);
-        for(int i = 0; i< backwardFaceCoordinates.size(); i++)
+        for (int i = 0; i< backwardFaceCoordinates.size(); i++)
         {
-            if(i == 0)
-            {
-                this->axisData_.inAxisBackwardDistances[i] = (backwardFaceCoordinates[i] - opposite.geometry().center()).two_norm();
-            }
+            if (i == 0)
+                axisData_.inAxisBackwardDistances[i] = (backwardFaceCoordinates[i] - opposite.geometry().center()).two_norm();
             else
-            {
-                this->axisData_.inAxisBackwardDistances[i] = (backwardFaceCoordinates[i] - backwardFaceCoordinates[i-1]).two_norm();
-            }
+                axisData_.inAxisBackwardDistances[i] = (backwardFaceCoordinates[i] - backwardFaceCoordinates[i-1]).two_norm();
         }
     }
 
@@ -351,42 +355,38 @@ private:
      */
     void fillPairData_()
     {
-        // initialize values that could remain unitialized if the intersection lies on a boundary
-        for(auto& data : pairData_)
-        {
-            // outer normals
-            data.normalPair.second = -1;
-        }
+        // reset the pair data structs
+        pairData_ = {};
 
         // set basic global positions
-        const auto& elementCenter = this->element_.geometry().center();
+        const auto& elementCenter = element_.geometry().center();
         const auto& FacetCenter = intersection_.geometry().center();
 
         // get the inner normal Dof Index
-        int numPairInnerNormalIdx = 0;
-        for(const auto& innerElementIntersection : intersections(gridView_, element_))
+        SmallLocalIndexType numPairInnerNormalIdx = 0;
+        for (const auto& innerElementIntersection : intersections(gridView_, element_))
         {
-            if(facetIsNormal_(innerElementIntersection.indexInInside(), intersection_.indexInInside()))
+            if (facetIsNormal_(innerElementIntersection.indexInInside(), intersection_.indexInInside()))
             {
-                const int innerElementIntersectionIdx = innerElementIntersection.indexInInside();
+                const auto innerElementIntersectionIdx = innerElementIntersection.indexInInside();
                 setNormalPairFirstInfo_(innerElementIntersectionIdx, element_, numPairInnerNormalIdx);
                 numPairInnerNormalIdx++;
             }
         }
 
         // get the outer normal Dof Index
-        int numPairOuterNormalIdx = 0;
-        if(intersection_.neighbor())
+        SmallLocalIndexType numPairOuterNormalIdx = 0;
+        if (intersection_.neighbor())
         {
             // the direct neighbor element and the respective intersection index
             const auto& directNeighborElement = intersection_.outside();
 
-            for(const auto& directNeighborElementIntersection : intersections(gridView_, directNeighborElement))
+            for (const auto& directNeighborElementIntersection : intersections(gridView_, directNeighborElement))
             {
                 // skip the directly neighboring face itself and its opposing one
-                if(facetIsNormal_(directNeighborElementIntersection.indexInInside(), intersection_.indexInOutside()))
+                if (facetIsNormal_(directNeighborElementIntersection.indexInInside(), intersection_.indexInOutside()))
                 {
-                    const int directNeighborElemIsIdx = directNeighborElementIntersection.indexInInside();
+                    const auto directNeighborElemIsIdx = directNeighborElementIntersection.indexInInside();
                     setNormalPairSecondInfo_(directNeighborElemIsIdx, directNeighborElement, numPairOuterNormalIdx);
                     numPairOuterNormalIdx++;
                 }
@@ -395,11 +395,11 @@ private:
         else // intersection is on boundary
         {
             // fill the normal pair entries
-            for(int pairIdx = 0; pairIdx < numPairs; ++pairIdx)
+            for(SmallLocalIndexType pairIdx = 0; pairIdx < numPairs; ++pairIdx)
             {
-                assert(pairData_[pairIdx].normalPair.second == -1);
+                assert(!pairData_[pairIdx].hasNormalNeighbor);
                 const auto distance = FacetCenter - elementCenter;
-                this->pairData_[pairIdx].normalDistance = std::move(distance.two_norm());
+                pairData_[pairIdx].normalDistance = std::move(distance.two_norm());
                 numPairOuterNormalIdx++;
             }
         }
@@ -413,32 +413,26 @@ private:
      */
     void fillParallelPairData_(std::false_type)
     {
-        // initialize values that could remain unitialized if the intersection lies on a boundary
-        for(auto& data : pairData_)
-        {
-            // parallel Dofs and Distances
-            data.parallelDofs.fill(-1);
-            data.parallelDistances.fill(0.0);
-        }
-
         // set basic global positions and stencil size definitions
         const auto& elementCenter = element_.geometry().center();
 
         // get the parallel Dofs
-        const int parallelLocalIdx = intersection_.indexInInside();
-        int numPairParallelIdx = 0;
-        for(const auto& intersection : intersections(gridView_, element_))
+        const auto parallelLocalIdx = intersection_.indexInInside();
+        SmallLocalIndexType numPairParallelIdx = 0;
+
+        for (const auto& intersection : intersections(gridView_, element_))
         {
-            if( facetIsNormal_(intersection.indexInInside(), parallelLocalIdx) )
+            if (facetIsNormal_(intersection.indexInInside(), parallelLocalIdx))
             {
                 // Store the parallel dimension of self cell in the direction of the axis
                 auto parallelAxisIdx = directionIndex(intersection);
                 pairData_[numPairParallelIdx].parallelDistances[0] = setParallelPairDistances_(element_, parallelAxisIdx);
 
-                if( intersection.neighbor() )
+                if (intersection.neighbor())
                 {
                     // If the normal intersection has a neighboring cell, go in and store the parallel information.
-                    auto outerElement = intersection.outside();
+                    const auto& outerElement = intersection.outside();
+                    pairData_[numPairParallelIdx].hasParallelNeighbor.set(0, true);
                     pairData_[numPairParallelIdx].parallelDofs[0] = gridView_.indexSet().subIndex(outerElement, parallelLocalIdx, codimIntersection);
                     pairData_[numPairParallelIdx].parallelDistances[1] = setParallelPairDistances_(outerElement, parallelAxisIdx);
                 }
@@ -447,7 +441,7 @@ private:
                     // If the intersection has no neighbor we have to deal with the virtual outer parallel dof
                     const auto& boundaryFacetCenter = intersection.geometry().center();
                     const auto distance = boundaryFacetCenter - elementCenter;
-                    const auto virtualFirstParallelFaceDofPos = this->intersection_.geometry().center() + distance;
+                    const auto virtualFirstParallelFaceDofPos = intersection_.geometry().center() + distance;
 
                     pairData_[numPairParallelIdx].virtualFirstParallelFaceDofPos = std::move(virtualFirstParallelFaceDofPos);
                 }
@@ -462,21 +456,13 @@ private:
      */
     void fillParallelPairData_(std::true_type)
     {
-        // initialize values that could remain unitialized if the intersection lies on a boundary
-        for(auto& data : pairData_)
-        {
-            // parallel Dofs and Distances
-            data.parallelDofs.fill(-1);
-            data.parallelDistances.fill(0.0);
-        }
-
         // set basic global positions and stencil size definitions
-        unsigned int numParallelFaces = pairData_[0].parallelDistances.size();
-        const auto& elementCenter = this->element_.geometry().center();
+        const auto numParallelFaces = pairData_[0].parallelDistances.size();
+        const auto& elementCenter = element_.geometry().center();
 
         // get the parallel Dofs
-        const int parallelLocalIdx = intersection_.indexInInside();
-        int numPairParallelIdx = 0;
+        const auto parallelLocalIdx = intersection_.indexInInside();
+        SmallLocalIndexType numPairParallelIdx = 0;
         std::stack<Element> parallelElementStack;
         for(const auto& intersection : intersections(gridView_, element_))
         {
@@ -513,9 +499,10 @@ private:
                     {
                         if(parallelElementStack.size() > 1)
                         {
-                            this->pairData_[numPairParallelIdx].parallelDofs[parallelElementStack.size()-2] = gridView_.indexSet().subIndex(parallelElementStack.top(), parallelLocalIdx, codimIntersection);
+                            pairData_[numPairParallelIdx].hasParallelNeighbor.set(parallelElementStack.size()-2, true);
+                            pairData_[numPairParallelIdx].parallelDofs[parallelElementStack.size()-2] = gridView_.indexSet().subIndex(parallelElementStack.top(), parallelLocalIdx, codimIntersection);
                         }
-                        this->pairData_[numPairParallelIdx].parallelDistances[parallelElementStack.size()-1] = setParallelPairDistances_(parallelElementStack.top(), parallelAxisIdx);
+                        pairData_[numPairParallelIdx].parallelDistances[parallelElementStack.size()-1] = setParallelPairDistances_(parallelElementStack.top(), parallelAxisIdx);
                         parallelElementStack.pop();
                     }
 
@@ -526,13 +513,13 @@ private:
                     const auto& boundaryFacetCenter = intersection.geometry().center();
 
                     const auto distance = boundaryFacetCenter - elementCenter;
-                    const auto virtualFirstParallelFaceDofPos = this->intersection_.geometry().center() + distance;
+                    const auto virtualFirstParallelFaceDofPos = intersection_.geometry().center() + distance;
 
-                    this->pairData_[numPairParallelIdx].virtualFirstParallelFaceDofPos = std::move(virtualFirstParallelFaceDofPos);
+                    pairData_[numPairParallelIdx].virtualFirstParallelFaceDofPos = std::move(virtualFirstParallelFaceDofPos);
 
                     // The distance is saved doubled because with scvf.cellCenteredSelfToFirstParallelDistance
                     // an average between parallelDistances[0] and parallelDistances[1] will be computed
-                    this->pairData_[numPairParallelIdx].parallelDistances[0] = std::move(distance.two_norm() * 2);
+                    pairData_[numPairParallelIdx].parallelDistances[0] = std::move(distance.two_norm() * 2);
                 }
                 numPairParallelIdx++;
             }
@@ -569,26 +556,25 @@ private:
     void setNormalPairFirstInfo_(const int isIdx, const Element& element, const int numPairsIdx)
     {
         // store the inner normal dofIdx
-        auto& dofIdx = pairData_[numPairsIdx].normalPair.first;
-        dofIdx = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
+        pairData_[numPairsIdx].normalPair.first = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
 
         // store the local normal facet index
-        this->pairData_[numPairsIdx].localNormalFaceIdx = isIdx;
+        pairData_[numPairsIdx].localNormalFaceIdx = isIdx;
     }
 
     //! Sets the information about the normal faces (in the neighbor element)
     void setNormalPairSecondInfo_(const int isIdx, const Element& element, const int numPairsIdx)
     {
         // store the dofIdx
-        auto& dofIdx = pairData_[numPairsIdx].normalPair.second;
-        dofIdx = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
+        pairData_[numPairsIdx].hasNormalNeighbor = true;
+        pairData_[numPairsIdx].normalPair.second = gridView_.indexSet().subIndex(element, isIdx, codimIntersection);
 
         // store the element distance
         const auto& outerNormalFacet = getFacet_(isIdx, element);
         const auto outerNormalFacetPos = outerNormalFacet.geometry().center();
         const auto& innerNormalFacet = getFacet_(isIdx, element_);
         const auto innerNormalFacetPos = innerNormalFacet.geometry().center();
-        this->pairData_[numPairsIdx].normalDistance = (innerNormalFacetPos - outerNormalFacetPos).two_norm();
+        pairData_[numPairsIdx].normalDistance = (innerNormalFacetPos - outerNormalFacetPos).two_norm();
     }
 
     //! Sets the information about the parallel distances
@@ -631,8 +617,8 @@ private:
     const Element element_; //!< The respective element
     const typename Element::Geometry elementGeometry_; //!< Reference to the element geometry
     const GridView gridView_; //!< The grid view
-    AxisData<Scalar, upwindSchemeOrder> axisData_; //!< Data related to forward and backward faces
-    std::array<PairData<Scalar, GlobalPosition, upwindSchemeOrder>, numPairs> pairData_; //!< Collection of pair information related to normal and parallel faces
+    AxisData axisData_; //!< Data related to forward and backward faces
+    std::array<PairData, numPairs> pairData_; //!< Collection of pair information related to normal and parallel faces
 };
 
 } // end namespace Dumux
diff --git a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
index fe7800a2caf8c1f402294b94b6e41a2641ba6469..9f163725c141730cb9a1906fa099adbddbb668d5 100644
--- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
@@ -72,6 +72,9 @@ class FreeFlowStaggeredSubControlVolumeFace
     using GridIndexType = typename IndexTraits<GV>::GridIndex;
     using LocalIndexType = typename IndexTraits<GV>::LocalIndex;
 
+    using PairData = typename T::PairData;
+    using AxisData = typename T::AxisData;
+
     using Scalar = typename T::Scalar;
     static const int dim = Geometry::mydimension;
     static const int dimworld = Geometry::coorddimension;
@@ -106,7 +109,7 @@ public:
       boundary_(is.boundary()),
 
       axisData_(geometryHelper.axisData()),
-      pairData_(geometryHelper.pairData()),
+      pairData_(std::move(geometryHelper.pairData())),
       localFaceIdx_(geometryHelper.localFaceIndex()),
       dirIdx_(geometryHelper.directionIndex()),
       outerNormalSign_(sign(unitOuterNormal_[directionIndex()])),
@@ -226,19 +229,19 @@ public:
     }
 
     //! Returns the data for one sub face
-    const PairData<Scalar, GlobalPosition, upwindSchemeOrder>& pairData(const int idx) const
+    const PairData& pairData(const int idx) const
     {
         return pairData_[idx];
     }
 
     //! Return an array of all pair data
-    const auto& pairData() const
+    const std::array<PairData, numPairs>& pairData() const
     {
         return pairData_;
     }
 
     //! Return an array of all pair data
-    const auto& axisData() const
+    const AxisData& axisData() const
     {
         return axisData_;
     }
@@ -257,7 +260,7 @@ public:
     */
     bool hasParallelNeighbor(const int localSubFaceIdx, const int parallelDegreeIdx) const
     {
-        return !(pairData(localSubFaceIdx).parallelDofs[parallelDegreeIdx] < 0);
+        return pairData(localSubFaceIdx).hasParallelNeighbor[parallelDegreeIdx];
     }
 
    /*!
@@ -267,7 +270,7 @@ public:
     */
     bool hasOuterNormal(const int localSubFaceIdx) const
     {
-        return !(pairData_[localSubFaceIdx].normalPair.second < 0);
+        return pairData(localSubFaceIdx).hasNormalNeighbor;
     }
 
    /*!
@@ -278,7 +281,7 @@ public:
     template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
     bool hasBackwardNeighbor(const int backwardIdx) const
     {
-        return !(axisData().inAxisBackwardDofs[backwardIdx] < 0);
+        return axisData().hasBackwardNeighbor[backwardIdx];
     }
 
    /*!
@@ -289,7 +292,7 @@ public:
     template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
     bool hasForwardNeighbor(const int forwardIdx) const
     {
-        return !(axisData().inAxisForwardDofs[forwardIdx] < 0);
+        return axisData().hasForwardNeighbor[forwardIdx];
     }
 
     //! Returns the dof of the face
@@ -347,8 +350,8 @@ private:
 
     int dofIdx_;
     Scalar selfToOppositeDistance_;
-    AxisData<Scalar, upwindSchemeOrder> axisData_;
-    std::array<PairData<Scalar, GlobalPosition, upwindSchemeOrder>, numPairs> pairData_;
+    AxisData axisData_;
+    std::array<PairData, numPairs> pairData_;
 
     int localFaceIdx_;
     unsigned int dirIdx_;
diff --git a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
index 62a49aaa3a6972e3b5c3ee1b722b7a3650425a1f..39ae0b74f5c856bf8352fcd5d316e51b0499a879 100644
--- a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
@@ -84,11 +84,11 @@ public:
      *        Then the corresponding set of momenta are collected and the prescribed
      *        upwinding method is used to calculate the momentum.
      */
-    FacePrimaryVariables computeUpwindedFrontalMomentum(const SubControlVolumeFace& scvf,
-                                                        const ElementFaceVariables& elemFaceVars,
-                                                        const ElementVolumeVariables& elemVolVars,
-                                                        const GridFluxVariablesCache& gridFluxVarsCache,
-                                                        const Scalar transportingVelocity)
+    static FacePrimaryVariables computeUpwindedFrontalMomentum(const SubControlVolumeFace& scvf,
+                                                               const ElementFaceVariables& elemFaceVars,
+                                                               const ElementVolumeVariables& elemVolVars,
+                                                               const GridFluxVariablesCache& gridFluxVarsCache,
+                                                               const Scalar transportingVelocity)
     {
         Scalar momentum(0.0);
         const bool canHigherOrder = canFrontalSecondOrder_(scvf, transportingVelocity, std::integral_constant<bool, useHigherOrder>{});
@@ -118,17 +118,17 @@ public:
      *        Then the corresponding set of momenta are collected and the prescribed
      *        upwinding method is used to calculate the momentum.
      */
-    FacePrimaryVariables computeUpwindedLateralMomentum(const Problem& problem,
-                                                        const FVElementGeometry& fvGeometry,
-                                                        const Element& element,
-                                                        const SubControlVolumeFace& scvf,
-                                                        const SubControlVolumeFace& normalFace,
-                                                        const ElementVolumeVariables& elemVolVars,
-                                                        const FaceVariables& faceVars,
-                                                        const GridFluxVariablesCache& gridFluxVarsCache,
-                                                        const int localSubFaceIdx,
-                                                        const bool lateralFaceHasDirichletPressure,
-                                                        const bool lateralFaceHasBJS)
+    static FacePrimaryVariables computeUpwindedLateralMomentum(const Problem& problem,
+                                                               const FVElementGeometry& fvGeometry,
+                                                               const Element& element,
+                                                               const SubControlVolumeFace& scvf,
+                                                               const SubControlVolumeFace& normalFace,
+                                                               const ElementVolumeVariables& elemVolVars,
+                                                               const FaceVariables& faceVars,
+                                                               const GridFluxVariablesCache& gridFluxVarsCache,
+                                                               const int localSubFaceIdx,
+                                                               const bool lateralFaceHasDirichletPressure,
+                                                               const bool lateralFaceHasBJS)
     {
         // Get the transporting velocity, located at the scvf perpendicular to the current scvf where the dof
         // of interest is located.
@@ -171,9 +171,9 @@ private:
      * \param transportingVelocity The average of the self and opposite velocities.
      * \param false_type False if higher order methods are not enabled.
      */
-    bool canFrontalSecondOrder_(const SubControlVolumeFace& ownScvf,
-                                const Scalar transportingVelocity,
-                                std::false_type) const
+    static bool canFrontalSecondOrder_(const SubControlVolumeFace& ownScvf,
+                                       const Scalar transportingVelocity,
+                                       std::false_type)
     {
         return false;
     }
@@ -191,9 +191,9 @@ private:
      * \param transportingVelocity The average of the self and opposite velocities.
      * \param true_type True if higher order methods are enabled.
      */
-    bool canFrontalSecondOrder_(const SubControlVolumeFace& ownScvf,
-                                const Scalar transportingVelocity,
-                                std::true_type) const
+    static bool canFrontalSecondOrder_(const SubControlVolumeFace& ownScvf,
+                                       const Scalar transportingVelocity,
+                                       std::true_type)
     {
         // Depending on selfIsUpstream I have to check if I have a forward or a backward neighbor to retrieve
         const bool selfIsUpstream = ownScvf.directionSign() != sign(transportingVelocity);
@@ -223,11 +223,11 @@ private:
      * \param transportingVelocity The average of the self and opposite velocities.
      * \param false_type False if higher order methods are not enabled.
      */
-    std::array<Scalar, 2> getFrontalUpwindingMomenta_(const SubControlVolumeFace& scvf,
-                                                      const ElementFaceVariables& elemFaceVars,
-                                                      const Scalar density,
-                                                      const Scalar transportingVelocity,
-                                                      std::false_type) const
+    static std::array<Scalar, 2> getFrontalUpwindingMomenta_(const SubControlVolumeFace& scvf,
+                                                             const ElementFaceVariables& elemFaceVars,
+                                                             const Scalar density,
+                                                             const Scalar transportingVelocity,
+                                                             std::false_type)
     {
         const Scalar momentumSelf = elemFaceVars[scvf].velocitySelf() * density;
         const Scalar momentumOpposite = elemFaceVars[scvf].velocityOpposite() * density;
@@ -246,11 +246,11 @@ private:
      * \param transportingVelocity The average of the self and opposite velocities.
      * \param true_type True if higher order methods are enabled.
      */
-    std::array<Scalar, 3> getFrontalUpwindingMomenta_(const SubControlVolumeFace& scvf,
-                                                      const ElementFaceVariables& elemFaceVars,
-                                                      const Scalar density,
-                                                      const Scalar transportingVelocity,
-                                                      std::true_type) const
+    static std::array<Scalar, 3> getFrontalUpwindingMomenta_(const SubControlVolumeFace& scvf,
+                                                             const ElementFaceVariables& elemFaceVars,
+                                                             const Scalar density,
+                                                             const Scalar transportingVelocity,
+                                                             std::true_type)
     {
         const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
         std::array<Scalar, 3> momenta{0.0, 0.0, 0.0};
@@ -279,10 +279,10 @@ private:
      * \param transportingVelocity The average of the self and opposite velocities.
      * \param gridFluxVarsCache The grid flux variables cache
      */
-    Scalar doFrontalMomentumUpwinding_(const SubControlVolumeFace& scvf,
-                                const std::array<Scalar, 2>& momenta,
-                                const Scalar transportingVelocity,
-                                const GridFluxVariablesCache& gridFluxVarsCache) const
+    static Scalar doFrontalMomentumUpwinding_(const SubControlVolumeFace& scvf,
+                                              const std::array<Scalar, 2>& momenta,
+                                              const Scalar transportingVelocity,
+                                              const GridFluxVariablesCache& gridFluxVarsCache)
     {
         const auto& upwindScheme = gridFluxVarsCache.staggeredUpwindMethods();
         return upwindScheme.upwind(momenta[0], momenta[1]);
@@ -297,10 +297,10 @@ private:
      * \param gridFluxVarsCache The grid flux variables cache
      */
     template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
-    Scalar doFrontalMomentumUpwinding_(const SubControlVolumeFace& scvf,
-                                       const std::array<Scalar, 3>& momenta,
-                                       const Scalar transportingVelocity,
-                                       const GridFluxVariablesCache& gridFluxVarsCache) const
+    static Scalar doFrontalMomentumUpwinding_(const SubControlVolumeFace& scvf,
+                                              const std::array<Scalar, 3>& momenta,
+                                              const Scalar transportingVelocity,
+                                              const GridFluxVariablesCache& gridFluxVarsCache)
     {
         const auto& upwindScheme = gridFluxVarsCache.staggeredUpwindMethods();
         const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity);
@@ -315,8 +315,8 @@ private:
      * \param selfIsUpstream bool describing upstream face.
      */
     template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
-    std::array<Scalar, 3> getFrontalDistances_(const SubControlVolumeFace& ownScvf,
-                                               const bool selfIsUpstream) const
+    static std::array<Scalar, 3> getFrontalDistances_(const SubControlVolumeFace& ownScvf,
+                                                      const bool selfIsUpstream)
     {
         // Depending on selfIsUpstream the downstream and the (up)upstream distances are saved.
         // distances {upstream to downstream distance, up-upstream to upstream distance, downstream staggered cell size}
@@ -348,10 +348,10 @@ private:
      * \param localSubFaceIdx The local subface index
      * \param true_type True when second order is enabled.
      */
-    bool canLateralSecondOrder_(const SubControlVolumeFace& ownScvf,
-                                const bool selfIsUpstream,
-                                const int localSubFaceIdx,
-                                std::true_type) const
+    static bool canLateralSecondOrder_(const SubControlVolumeFace& ownScvf,
+                                       const bool selfIsUpstream,
+                                       const int localSubFaceIdx,
+                                       std::true_type)
     {
         if (selfIsUpstream)
         {
@@ -373,10 +373,10 @@ private:
      * \brief Check if a second order approximation for the lateral part of the advective term can be used
      *        If higher order methods are not prescribed, this is called and false is returned.
      */
-    bool canLateralSecondOrder_(const SubControlVolumeFace& ownScvf,
-                                const bool selfIsUpstream,
-                                const int localSubFaceIdx,
-                                std::false_type) const
+    static bool canLateralSecondOrder_(const SubControlVolumeFace& ownScvf,
+                                       const bool selfIsUpstream,
+                                       const int localSubFaceIdx,
+                                       std::false_type)
     {   return false; }
 
 
@@ -385,17 +385,17 @@ private:
      *
      *  Only called if higher order methods are enabled and the scvf can use higher order methods.
      */
-    std::array<Scalar, 3> getLateralUpwindingMomenta_(const Problem& problem,
-                                                      const FVElementGeometry& fvGeometry,
-                                                      const Element& element,
-                                                      const SubControlVolumeFace& ownScvf,
-                                                      const ElementVolumeVariables& elemVolVars,
-                                                      const FaceVariables& faceVars,
-                                                      const Scalar transportingVelocity,
-                                                      const int localSubFaceIdx,
-                                                      bool lateralFaceHasDirichletPressure,
-                                                      bool lateralFaceHasBJS,
-                                                      std::true_type) const
+    static std::array<Scalar, 3> getLateralUpwindingMomenta_(const Problem& problem,
+                                                             const FVElementGeometry& fvGeometry,
+                                                             const Element& element,
+                                                             const SubControlVolumeFace& ownScvf,
+                                                             const ElementVolumeVariables& elemVolVars,
+                                                             const FaceVariables& faceVars,
+                                                             const Scalar transportingVelocity,
+                                                             const int localSubFaceIdx,
+                                                             bool lateralFaceHasDirichletPressure,
+                                                             bool lateralFaceHasBJS,
+                                                             std::true_type)
     {
         std::array<Scalar, 3> momenta{0.0, 0.0, 0.0};
         const SubControlVolumeFace& lateralFace = fvGeometry.scvf(ownScvf.insideScvIdx(), ownScvf.pairData(localSubFaceIdx).localNormalFaceIdx);
@@ -456,17 +456,17 @@ private:
      *
      *  Called if higher order methods are not enabled of if the scvf can not use higher order methods.
      */
-    std::array<Scalar, 2> getLateralUpwindingMomenta_(const Problem& problem,
-                                                       const FVElementGeometry& fvGeometry,
-                                                       const Element& element,
-                                                       const SubControlVolumeFace& ownScvf,
-                                                       const ElementVolumeVariables& elemVolVars,
-                                                       const FaceVariables& faceVars,
-                                                       const Scalar transportingVelocity,
-                                                       const int localSubFaceIdx,
-                                                       bool lateralFaceHasDirichletPressure,
-                                                       bool lateralFaceHasBJS,
-                                                       std::false_type) const
+    static std::array<Scalar, 2> getLateralUpwindingMomenta_(const Problem& problem,
+                                                             const FVElementGeometry& fvGeometry,
+                                                             const Element& element,
+                                                             const SubControlVolumeFace& ownScvf,
+                                                             const ElementVolumeVariables& elemVolVars,
+                                                             const FaceVariables& faceVars,
+                                                             const Scalar transportingVelocity,
+                                                             const int localSubFaceIdx,
+                                                             bool lateralFaceHasDirichletPressure,
+                                                             bool lateralFaceHasBJS,
+                                                             std::false_type)
     {
          // Check whether the own or the neighboring element is upstream.
         const SubControlVolumeFace& lateralFace = fvGeometry.scvf(ownScvf.insideScvIdx(), ownScvf.pairData(localSubFaceIdx).localNormalFaceIdx);
@@ -494,12 +494,12 @@ private:
      *
      *        Fowards to Frontal Momentum Upwinding method
      */
-    Scalar doLateralMomentumUpwinding_(const SubControlVolumeFace& scvf,
-                                       const SubControlVolumeFace& normalScvf,
-                                       const std::array<Scalar, 2>& momenta,
-                                       const Scalar transportingVelocity,
-                                       const int localSubFaceIdx,
-                                       const GridFluxVariablesCache& gridFluxVarsCache) const
+    static Scalar doLateralMomentumUpwinding_(const SubControlVolumeFace& scvf,
+                                              const SubControlVolumeFace& normalScvf,
+                                              const std::array<Scalar, 2>& momenta,
+                                              const Scalar transportingVelocity,
+                                              const int localSubFaceIdx,
+                                              const GridFluxVariablesCache& gridFluxVarsCache)
     {
         return doFrontalMomentumUpwinding_(scvf, momenta, transportingVelocity, gridFluxVarsCache);
     }
@@ -514,12 +514,12 @@ private:
      * \param localSubFaceIdx  The local index of the subface
      * \param gridFluxVarsCache The grid flux variables cache
      */
-    Scalar doLateralMomentumUpwinding_(const SubControlVolumeFace& scvf,
-                                       const SubControlVolumeFace& normalScvf,
-                                       const std::array<Scalar, 3>& momenta,
-                                       const Scalar transportingVelocity,
-                                       const int localSubFaceIdx,
-                                       const GridFluxVariablesCache& gridFluxVarsCache) const
+    static Scalar doLateralMomentumUpwinding_(const SubControlVolumeFace& scvf,
+                                              const SubControlVolumeFace& normalScvf,
+                                              const std::array<Scalar, 3>& momenta,
+                                              const Scalar transportingVelocity,
+                                              const int localSubFaceIdx,
+                                              const GridFluxVariablesCache& gridFluxVarsCache)
     {
         const bool selfIsUpstream = ( normalScvf.directionSign() == sign(transportingVelocity) );
         const auto& upwindScheme = gridFluxVarsCache.staggeredUpwindMethods();
@@ -535,9 +535,9 @@ private:
      * \param selfIsUpstream bool describing upstream face.
      */
     template<bool enable = useHigherOrder, std::enable_if_t<enable, int> = 0>
-    std::array<Scalar, 3> getLateralDistances_(const SubControlVolumeFace& ownScvf,
-                                               const int localSubFaceIdx,
-                                               const bool selfIsUpstream) const
+    static std::array<Scalar, 3> getLateralDistances_(const SubControlVolumeFace& ownScvf,
+                                                      const int localSubFaceIdx,
+                                                      const bool selfIsUpstream)
     {
         // distances {upstream to downstream distance, up-upstream to upstream distance, downstream staggered cell size}
         std::array<Scalar, 3> distances{0.0, 0.0, 0.0};
@@ -578,14 +578,14 @@ private:
      * \param lateralFaceHasDirichletPressure @c true if there is a dirichlet condition for the pressure on the boundary
      * \param lateralFaceHasBJS @c true if there is a BJS condition fot the velocity on the boundary
      */
-    Scalar getParallelVelocityFromBoundary_(const Problem& problem,
-                                            const SubControlVolumeFace& scvf,
-                                            const SubControlVolumeFace& normalFace,
-                                            const Scalar velocitySelf,
-                                            const int localSubFaceIdx,
-                                            const Element& element,
-                                            const bool lateralFaceHasDirichletPressure,
-                                            const bool lateralFaceHasBJS) const
+    static Scalar getParallelVelocityFromBoundary_(const Problem& problem,
+                                                   const SubControlVolumeFace& scvf,
+                                                   const SubControlVolumeFace& normalFace,
+                                                   const Scalar velocitySelf,
+                                                   const int localSubFaceIdx,
+                                                   const Element& element,
+                                                   const bool lateralFaceHasDirichletPressure,
+                                                   const bool lateralFaceHasBJS)
     {
         // If there is a Dirichlet condition for the pressure we assume zero gradient for the velocity,
         // so the velocity at the boundary equal to that on the scvf.
@@ -607,12 +607,12 @@ private:
      * \param boundaryElement The element that is on the boundary
      * \param parallelVelocity The velocity over scvf
      */
-    Scalar getParallelVelocityFromOtherBoundary_(const Problem& problem,
-                                                 const FVElementGeometry& fvGeometry,
-                                                 const SubControlVolumeFace& scvf,
-                                                 const int localIdx,
-                                                 const Element& boundaryElement,
-                                                 const Scalar parallelVelocity) const
+    static Scalar getParallelVelocityFromOtherBoundary_(const Problem& problem,
+                                                        const FVElementGeometry& fvGeometry,
+                                                        const SubControlVolumeFace& scvf,
+                                                        const int localIdx,
+                                                        const Element& boundaryElement,
+                                                        const Scalar parallelVelocity)
     {
         // A ghost subface at the boundary is created, featuring the location of the sub face's center
         const SubControlVolumeFace& boundaryNormalFace = fvGeometry.scvf(scvf.insideScvIdx(), scvf.pairData(localIdx).localNormalFaceIdx);
@@ -646,14 +646,14 @@ private:
 
 
     //! helper function to conveniently create a ghost face used to retrieve boundary values from the problem
-    SubControlVolumeFace makeGhostFace_(const SubControlVolumeFace& ownScvf, const GlobalPosition& pos) const
+    static SubControlVolumeFace makeGhostFace_(const SubControlVolumeFace& ownScvf, const GlobalPosition& pos)
     {
         return SubControlVolumeFace(pos, std::vector<unsigned int>{ownScvf.insideScvIdx(), ownScvf.outsideScvIdx()},
                                     ownScvf.directionIndex(), ownScvf.axisData().selfDof, ownScvf.index());
     };
 
     //! helper function to conveniently create a ghost face which is outside the domain, parallel to the scvf of interest
-    SubControlVolumeFace makeParallelGhostFace_(const SubControlVolumeFace& ownScvf, const int localSubFaceIdx) const
+    static SubControlVolumeFace makeParallelGhostFace_(const SubControlVolumeFace& ownScvf, const int localSubFaceIdx)
     {
         return makeGhostFace_(ownScvf, ownScvf.pairData(localSubFaceIdx).virtualFirstParallelFaceDofPos);
     };