diff --git a/dumux/discretization/staggered/gridfluxvariablescache.hh b/dumux/discretization/staggered/gridfluxvariablescache.hh
index 5214858f72c0a31f65d42a9e660d860274b7bb2e..3c87372776d2c95922574fc8730da2ef0249097d 100644
--- a/dumux/discretization/staggered/gridfluxvariablescache.hh
+++ b/dumux/discretization/staggered/gridfluxvariablescache.hh
@@ -28,6 +28,8 @@
 #include <dumux/discretization/localview.hh>
 #include <dumux/discretization/staggered/elementfluxvariablescache.hh>
 
+#include <dumux/freeflow/higherorderapproximation.hh>
+
 namespace Dumux {
 
 /*!
@@ -42,7 +44,6 @@ struct StaggeredDefaultGridFluxVariablesCacheTraits
     using Problem = P;
     using FluxVariablesCache = FVC;
     using FluxVariablesCacheFiller = FVCF;
-
     template<class GridFluxVariablesCache, bool cachingEnabled>
     using LocalView = StaggeredElementFluxVariablesCache<GridFluxVariablesCache, cachingEnabled>;
 };
@@ -72,6 +73,7 @@ class StaggeredGridFluxVariablesCache<P, FVC, FVCF, true, Traits>
 public:
     //! export the flux variable cache type
     using FluxVariablesCache = typename Traits::FluxVariablesCache;
+    using Scalar = typename FluxVariablesCache::Scalar;
 
     //! export the flux variable cache filler type
     using FluxVariablesCacheFiller = typename Traits::FluxVariablesCacheFiller;
@@ -82,7 +84,10 @@ public:
     //! export the type of the local view
     using LocalView = typename Traits::template LocalView<ThisType, cachingEnabled>;
 
-    StaggeredGridFluxVariablesCache(const Problem& problem) : problemPtr_(&problem) {}
+    StaggeredGridFluxVariablesCache(const Problem& problem, const std::string& paramGroup = "")
+    : problemPtr_(&problem)
+    , higherOrderApproximation_(paramGroup)
+      {}
 
     // When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
     template<class FVGridGeometry, class GridVolumeVariables, class SolutionVector>
@@ -118,6 +123,12 @@ public:
         }
     }
 
+    //! Return the HigherOrderApproximation
+    const HigherOrderApproximation<Scalar>& higherOrderApproximation() const
+    {
+        return higherOrderApproximation_;
+    }
+
     const Problem& problem() const
     { return *problemPtr_; }
 
@@ -130,6 +141,7 @@ public:
 
 private:
     const Problem* problemPtr_;
+    HigherOrderApproximation<Scalar> higherOrderApproximation_;
 
     std::vector<FluxVariablesCache> fluxVarsCache_;
     std::vector<std::size_t> globalScvfIndices_;
@@ -149,6 +161,7 @@ class StaggeredGridFluxVariablesCache<P, FVC, FVCF, false, Traits>
 public:
     //! export the flux variable cache type
     using FluxVariablesCache = typename Traits::FluxVariablesCache;
+    using Scalar = typename FluxVariablesCache::Scalar;
 
     //! export the flux variable cache filler type
     using FluxVariablesCacheFiller = typename Traits::FluxVariablesCacheFiller;
@@ -159,7 +172,10 @@ public:
     //! export the type of the local view
     using LocalView = typename Traits::template LocalView<ThisType, cachingEnabled>;
 
-    StaggeredGridFluxVariablesCache(const Problem& problem) : problemPtr_(&problem) {}
+    StaggeredGridFluxVariablesCache(const Problem& problem, const std::string& paramGroup = "")
+    : problemPtr_(&problem)
+    , higherOrderApproximation_(paramGroup)
+      {}
 
     // When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
     template<class FVGridGeometry, class GridVolumeVariables, class SolutionVector>
@@ -171,10 +187,16 @@ public:
     const Problem& problem() const
     { return *problemPtr_; }
 
-private:
-
+    //! Return the HigherOrderApproximation
+    const HigherOrderApproximation<Scalar>& higherOrderApproximation() const
+    {
+        return higherOrderApproximation_;
+    }
 
+private:
     const Problem* problemPtr_;
+    HigherOrderApproximation<Scalar> higherOrderApproximation_;
+
 };
 
 } // end namespace Dumux
diff --git a/dumux/freeflow/navierstokes/fluxvariablescache.hh b/dumux/freeflow/navierstokes/fluxvariablescache.hh
index 9f1c6ff9144b43bd215198fbaf2441212393eda0..a7651c0496c53499ea4843b6fe1b8b439956ee0f 100644
--- a/dumux/freeflow/navierstokes/fluxvariablescache.hh
+++ b/dumux/freeflow/navierstokes/fluxvariablescache.hh
@@ -59,6 +59,8 @@ class FreeFlowFluxVariablesCacheImplementation<TypeTag, DiscretizationMethod::st
     using Element = typename GridView::template Codim<0>::Entity;
 
 public:
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+
     //! Do nothing for the staggered grid specialization.
     void update(const Problem& problem,
                 const Element& element,
diff --git a/dumux/freeflow/navierstokes/problem.hh b/dumux/freeflow/navierstokes/problem.hh
index d988a03955a7e0471131bbcd694974c0bf31d05a..d16dfbfe62318301ff529aab7a11a0ca2c9f30a4 100644
--- a/dumux/freeflow/navierstokes/problem.hh
+++ b/dumux/freeflow/navierstokes/problem.hh
@@ -29,8 +29,6 @@
 #include <dumux/common/staggeredfvproblem.hh>
 #include <dumux/discretization/method.hh>
 
-#include <dumux/freeflow/higherorderapproximation.hh>
-
 #include "model.hh"
 
 namespace Dumux {
@@ -95,9 +93,8 @@ public:
      * \param paramGroup The parameter group in which to look for runtime parameters first (default is "")
      */
     NavierStokesProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, const std::string& paramGroup = "")
-    : ParentType(fvGridGeometry, paramGroup),
-      gravity_(0.0),
-      higherOrderApproximation_(paramGroup)
+    : ParentType(fvGridGeometry, paramGroup)
+    , gravity_(0.0)
     {
         if (getParamFromGroup<bool>(paramGroup, "Problem.EnableGravity"))
             gravity_[dim-1]  = -9.81;
@@ -222,12 +219,6 @@ public:
         return velocitySelf / (alpha / sqrt(K) * scvf.cellCenteredParallelDistance(localSubFaceIdx,0) + 1.0);
     }
 
-    //! Return the HigherOrderApproximation
-    const HigherOrderApproximation<Scalar>& higherOrderApproximation() const
-    {
-        return higherOrderApproximation_;
-    }
-
 private:
 
     //! Returns the implementation of the problem (i.e. static polymorphism)
@@ -239,7 +230,6 @@ private:
     { return *static_cast<const Implementation *>(this); }
 
     GravityVector gravity_;
-    HigherOrderApproximation<Scalar> higherOrderApproximation_;
     bool enableInertiaTerms_;
 };
 
diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
index 448b046e25576b86f8a3b98daaebee8636e2b385..a39cd4824f2866c0e120b8c90e98ab13e3781b03 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -163,10 +163,11 @@ public:
                                              const SubControlVolumeFace& scvf,
                                              const FVElementGeometry& fvGeometry,
                                              const ElementVolumeVariables& elemVolVars,
-                                             const ElementFaceVariables& elemFaceVars)
+                                             const ElementFaceVariables& elemFaceVars,
+                                             const GridFluxVariablesCache& gridFluxVarsCache)
     {
-        return computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars) +
-               computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars);
+        return computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache) +
+               computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache);
     }
 
     /*!
@@ -191,7 +192,8 @@ public:
                                                     const SubControlVolumeFace& scvf,
                                                     const FVElementGeometry& fvGeometry,
                                                     const ElementVolumeVariables& elemVolVars,
-                                                    const ElementFaceVariables& elemFaceVars)
+                                                    const ElementFaceVariables& elemFaceVars,
+                                                    const GridFluxVariablesCache& gridFluxVarsCache)
     {
         FacePrimaryVariables frontalFlux(0.0);
 
@@ -227,7 +229,7 @@ public:
             // distances[2]: downstream staggered cell size
             std::array<Scalar, 3> distances{0.0, 0.0, 0.0};
 
-            const auto& highOrder = problem.higherOrderApproximation();
+            const auto& highOrder = gridFluxVarsCache.higherOrderApproximation();
 
             // If a Tvd approach has been specified and I am not too near to the boundary I can use a second order
             // approximation for the velocity. In this frontal flux I use for the density always the value that I have on the scvf.
@@ -304,7 +306,8 @@ public:
                                                     const SubControlVolumeFace& scvf,
                                                     const FVElementGeometry& fvGeometry,
                                                     const ElementVolumeVariables& elemVolVars,
-                                                    const ElementFaceVariables& elemFaceVars)
+                                                    const ElementFaceVariables& elemFaceVars,
+                                                    const GridFluxVariablesCache& gridFluxVarsCache)
     {
         FacePrimaryVariables normalFlux(0.0);
         auto& faceVars = elemFaceVars[scvf];
@@ -364,7 +367,7 @@ public:
 
             // If there is no symmetry or Neumann boundary condition for the given sub face, proceed to calculate the tangential momentum flux.
             if (problem.enableInertiaTerms())
-                normalFlux += computeAdvectivePartOfLateralMomentumFlux_(problem, fvGeometry, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx, lateralFaceHasDirichletPressure, lateralFaceHasBJS);
+                normalFlux += computeAdvectivePartOfLateralMomentumFlux_(problem, fvGeometry, element, scvf, normalFace, elemVolVars, faceVars, gridFluxVarsCache, localSubFaceIdx, lateralFaceHasDirichletPressure, lateralFaceHasBJS);
 
             normalFlux += computeDiffusivePartOfLateralMomentumFlux_(problem, fvGeometry, element, scvf, normalFace, elemVolVars, faceVars, localSubFaceIdx, lateralFaceHasDirichletPressure, lateralFaceHasBJS);
         }
@@ -401,6 +404,7 @@ private:
                                                                     const SubControlVolumeFace& normalFace,
                                                                     const ElementVolumeVariables& elemVolVars,
                                                                     const FaceVariables& faceVars,
+                                                                    const GridFluxVariablesCache& gridFluxVarsCache,
                                                                     const int localSubFaceIdx,
                                                                     const bool lateralFaceHasDirichletPressure,
                                                                     const bool lateralFaceHasBJS)
@@ -432,7 +436,7 @@ private:
         // distances[2]: downstream staggered cell size
         std::array<Scalar, 3> distances{0.0, 0.0, 0.0};
 
-        const auto& highOrder = problem.higherOrderApproximation();
+        const auto& highOrder = gridFluxVarsCache.higherOrderApproximation();
 
         // If a Tvd approach has been specified and I am not too near to the boundary I can use a second order approximation.
         if (highOrder.tvdApproach() != TvdApproach::none)
diff --git a/dumux/freeflow/navierstokes/staggered/localresidual.hh b/dumux/freeflow/navierstokes/staggered/localresidual.hh
index 3ddde3953667bc5e8706645799d5d8afb4948095..fe2039b902a09824a89356f1b146f26d8ed180ea 100644
--- a/dumux/freeflow/navierstokes/staggered/localresidual.hh
+++ b/dumux/freeflow/navierstokes/staggered/localresidual.hh
@@ -179,7 +179,7 @@ public:
                                             const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         FluxVariables fluxVars;
-        return fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars);
+        return fluxVars.computeMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache.gridFluxVarsCache());
     }
 
     /*!
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh
index 77d63bdd6af87dd33bd639e8f0acb9a250af653c..d831d5917bd6b2f47a28b67994ed7a63ebcd3a73 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh
@@ -184,12 +184,13 @@ public:
                                              const SubControlVolumeFace& scvf,
                                              const FVElementGeometry& fvGeometry,
                                              const ElementVolumeVariables& elemVolVars,
-                                             const ElementFaceVariables& elemFaceVars)
+                                             const ElementFaceVariables& elemFaceVars,
+                                             const GridFluxVariablesCache& gridFluxVarsCache)
     {
         const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
 
-        return ParentType::computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars)
-               + ParentType::computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars)
+        return ParentType::computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache)
+               + ParentType::computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache)
                + 2.0 / ModelTraits::dim() * insideVolVars.density() * insideVolVars.turbulentKineticEnergy()
                  * scvf.area() * scvf.directionSign() * insideVolVars.extrusionFactor();
     }
diff --git a/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh
index bc76b71421488e78856afd28241521bf6190b271..22b659dd06e034be428cf9b0def870718e03d050 100644
--- a/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh
+++ b/dumux/freeflow/rans/twoeq/komega/staggered/fluxvariables.hh
@@ -175,12 +175,13 @@ public:
                                              const SubControlVolumeFace& scvf,
                                              const FVElementGeometry& fvGeometry,
                                              const ElementVolumeVariables& elemVolVars,
-                                             const ElementFaceVariables& elemFaceVars)
+                                             const ElementFaceVariables& elemFaceVars,
+                                             const GridFluxVariablesCache& gridFluxVarsCache)
     {
         const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
 
-        return ParentType::computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars)
-               + ParentType::computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars)
+        return ParentType::computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache)
+               + ParentType::computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache)
                + 2.0 / ModelTraits::dim() * insideVolVars.density() * insideVolVars.turbulentKineticEnergy()
                  * scvf.area() * scvf.directionSign() * insideVolVars.extrusionFactor();
     }
diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh
index 0add605a222922346538bf123f386ef98ec7a11a..a89108b327113bf87d8f0733e6fe58416be3d29e 100644
--- a/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh
+++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/staggered/fluxvariables.hh
@@ -176,12 +176,13 @@ public:
                                              const SubControlVolumeFace& scvf,
                                              const FVElementGeometry& fvGeometry,
                                              const ElementVolumeVariables& elemVolVars,
-                                             const ElementFaceVariables& elemFaceVars)
+                                             const ElementFaceVariables& elemFaceVars,
+                                             const GridFluxVariablesCache& gridFluxVarsCache)
     {
         const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
 
-        return ParentType::computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars)
-               + ParentType::computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars)
+        return ParentType::computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache)
+               + ParentType::computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache)
                + 2.0 / ModelTraits::dim() * insideVolVars.density() * insideVolVars.turbulentKineticEnergy()
                  * scvf.area() * scvf.directionSign() * insideVolVars.extrusionFactor();
     }